Bio.bgzf模块

读写BGZR压缩文件(BAM中使用的GZip变体)。

Sam/BAM文件格式(序列对齐/映射)有纯文本格式(Sam)和压缩二进制格式(BAM)。后者使用一种名为BGZR(Blocked GNU Zip Form)的修改后的gZip压缩形式,它可以应用于任何文件格式,以提供高效随机访问的压缩。BGZR与Sam/BAM文件格式一起在https://samtools.sourceforge.net/SAM1.pdf上描述

在使用BGZR文件进行随机访问之前,请阅读下面有关“虚拟补偿”的文本。

本模块的目标

Python gZip库可用于读取BGZR文件,因为对于解压来说,它们只是(专门的)gZip文件。该模块的目的是促进随机访问BGZR文件(使用“虚拟偏置”想法)并编写BGZR文件(这意味着使用适当大小的gZip块并在gZip头中写入额外的“BC”字段)。与gZip库一样,zlib库是在内部使用的。

除了随机访问和写入BAM文件所需外,BGZR格式还可以用于其他顺序数据(从一条又一条记录的意义上说),例如Bio.SeqIO中支持的大多数序列数据格式(如FASTA、FASTQ、SEN等)或大型APM比对。

Bio.SeqIO索引功能使用此模块来支持BGZR文件。

BGZR技术介绍

gZip文件格式允许多个压缩块,每个块都可以是独立的gZip文件。作为一个有趣的奖励,这意味着您可以使用Unix cat 通过将两个或多个gZip文件连接起来将它们组合成一个。此外,每个块都可以具有几个压缩级别之一(包括未压缩的,由于gZip标头,这实际上会占用更多的空间)。

BAM设计者意识到,虽然随机访问存储在传统gZip文件中的数据很慢,但将文件分解为gZip块将允许对每个块进行快速随机访问。要访问特定的解压缩数据,您只需要知道它从哪个块开始(gZip块开始的偏差),以及您需要阅读的块(解压缩)内容有多远。

其中一个问题是有效地找到gZip块大小。您可以使用标准的gZip文件来实现这一点,但它需要对每个块进行解压--而且这会相当慢。此外,典型的gZip文件可能使用非常大的块。

BGZR中的唯一不同之处在于每个gZip块的压缩大小限制为2 ' 16字节,gZip标头中额外的“BC”字段记录了此大小。传统的解压工具可以忽略这一点,并像任何其他gZip文件一样解压该文件。

这样做的目的是,你可以查看第一个BGZF块,从这个“BC”头中找出它有多大,从而立即查找第二个块,等等。

BAM索引方案使用64位“虚拟偏置”记录读取位置,包括 coffset << 16 | uoffset ,在哪里 coffset 是包含读取开始的BGZR块的文件偏置量(使用最多64-16 = 48位的无符号整数),以及 uoffset 是(解压缩)块内的偏移量(无符号16位整数)。

这将您限制在最后一个块以2#48字节(即256 PB)开头的ADAM文件中,并且每个块的解压缩大小最多为2#16字节(即64 KB)。请注意,这与BGZR“BC”字段大小相匹配,该字段将每个块的压缩大小限制为2,16字节,允许BAM文件使用BGZR,而无需gZip压缩(对于内存中的中间文件有用,以减少处理器负载)。

有关名称空间的警告

使用“从XXX导入”被认为是一个坏主意 * “在Python中,因为它会污染命名空间。这是Bio.bgzf(以及标准Python库gZip)的一个真正问题,因为它们包含一个名为Open的函数。假设您这样做:

>>> from Bio.bgzf import *
>>> print(open.__module__)
Bio.bgzf

或者,

>>> from gzip import *
>>> print(open.__module__)
gzip

请注意,打开函数已被替换。如果需要,您可以通过导入内置的打开函数来“修复”这个问题:

>>> from builtins import open

然而,我们建议使用显式命名空间,例如

>>> from Bio import bgzf
>>> print(bgzf.open.__module__)
Bio.bgzf

示例

这是一个使用BGZR压缩的普通基因库文件,因此可以使用gZip进行解压,

>>> import gzip
>>> handle = gzip.open("GenBank/NC_000932.gb.bgz", "r")
>>> assert 0 == handle.tell()
>>> line = handle.readline()
>>> assert 80 == handle.tell()
>>> line = handle.readline()
>>> assert 143 == handle.tell()
>>> data = handle.read(70000)
>>> assert 70143 == handle.tell()
>>> handle.close()

我们还可以使用BGZR阅读器访问文件-但请注意文件偏差,下文将对此进行解释:

>>> handle = BgzfReader("GenBank/NC_000932.gb.bgz", "r")
>>> assert 0 == handle.tell()
>>> print(handle.readline().rstrip())
LOCUS       NC_000932             154478 bp    DNA     circular PLN 15-APR-2009
>>> assert 80 == handle.tell()
>>> print(handle.readline().rstrip())
DEFINITION  Arabidopsis thaliana chloroplast, complete genome.
>>> assert 143 == handle.tell()
>>> data = handle.read(70000)
>>> assert 987828735 == handle.tell()
>>> print(handle.readline().rstrip())
f="GeneID:844718"
>>> print(handle.readline().rstrip())
     CDS             complement(join(84337..84771,85454..85843))
>>> offset = handle.seek(make_virtual_offset(55074, 126))
>>> print(handle.readline().rstrip())
    68521 tatgtcattc gaaattgtat aaagacaact cctatttaat agagctattt gtgcaagtat
>>> handle.close()

请注意,手柄的偏差看起来与BGZR文件不同。这让我们来到了BGZR的关键点,即区块结构:

>>> handle = open("GenBank/NC_000932.gb.bgz", "rb")
>>> for values in BgzfBlocks(handle):
...     print("Raw start %i, raw length %i; data start %i, data length %i" % values)
Raw start 0, raw length 15073; data start 0, data length 65536
Raw start 15073, raw length 17857; data start 65536, data length 65536
Raw start 32930, raw length 22144; data start 131072, data length 65536
Raw start 55074, raw length 22230; data start 196608, data length 65536
Raw start 77304, raw length 14939; data start 262144, data length 43478
Raw start 92243, raw length 28; data start 305622, data length 0
>>> handle.close()

在本例中,前三个块是“完整的”,并保存65536字节的未压缩数据。第四个块未满,包含43478个字节。最后,有一个特殊的空第五个块,它占用磁盘上的28个字节,并充当“文件结束”(RST)标记。如果缺少该内容,则您的BGZR文件可能不完整。

通过提前读取70,000字节,我们移动到第二个BGZF块中,此时BGZF虚拟偏移量开始与gzip库公开的解压缩数据中的简单偏移量不同。

例如,考虑寻找解压位置196734。由于196734 = 65536 + 65536 + 65536 + 126 = 65536*3 + 126,这相当于跳过前三个块(在这个特定示例中,解压后大小均为65536-这并不总是保持)并从第四个块的字节126开始(解压后)。对于BGZR,我们需要知道第四个块的偏差55074和块内的偏差126才能获得BGZR虚拟偏差。

>>> print(55074 << 16 | 126)
3609329790
>>> print(bgzf.make_virtual_offset(55074, 126))
3609329790

因此,对于该BGZR文件,解压缩位置196734对应于虚拟偏差3609329790。然而,具有不同内容的另一个BGZR文件的压缩效率会或多或少,因此压缩的块的大小会不同。这意味着未压缩的偏差和压缩的虚拟偏差之间的映射取决于您正在使用的BGZR文件。

如果您要通过此模块访问BGZR文件,只需使用handle.tell()方法来记录您稍后可能想要使用handle.seek()返回的位置的虚拟偏差。

BGZR虚拟偏差的问题是,虽然可以比较它们(哪个偏差在文件中是第一个),但您无法安全地减去它们来获得它们之间的数据大小,也无法添加/减去相对偏差。

当然,您可以使用BgzfReader使用Bio.SeqIO解析此文件,尽管与使用gzip.open(.)相比没有任何好处,除非您想对BGZR压缩序列文件进行索引:

>>> from Bio import SeqIO
>>> handle = BgzfReader("GenBank/NC_000932.gb.bgz")
>>> record = SeqIO.read(handle, "genbank")
>>> handle.close()
>>> print(record.id)
NC_000932.1

文本模式

就像标准库gzip.open(.)一样,BGZR代码默认以二进制模式打开文件。

您可以请求以文本模式打开该文件,但请注意,这是硬编码为简单的“latin 1”(又名“iso-8859-1”)编码(其中包括所有ASC字符),该编码与大多数西欧语言都很好用。然而,它与更广泛使用的UTF-8编码并不完全兼容。

在像UTF-8这样的可变宽度编码中,unicode文本输出中的一些单个字符由原始二进制形式的多个字节表示。这对于BGZR来说是有问题的,因为我们不能总是孤立地解码每个块-一个unicode字符可能会分裂到两个块上。这种情况甚至可能发生在固定宽度的unicode编码中,因为BGZR块大小不固定。

因此,该模块目前仅限于仅支持单字节Unicode编码,例如ASCI、“latin 1”(这是ASCI的超集)或潜在的其他字符映射(未实现)。

Furthermore, unlike the default text mode on Python 3, we do not attempt to implement universal new line mode. This transforms the various operating system new line conventions like Windows (CR LF or "rn"), Unix (just LF, "n"), or old Macs (just CR, "r"), into just LF ("n"). Here we have the same problem - is "r" at the end of a block an incomplete Windows style new line?

相反,您将按原样获得CR(“r”)和LF(“n”)字符。

如果你的数据是UTF-8或任何其他不兼容的编码,你必须使用二进制模式,并自己解码适当的片段。

Bio.bgzf.open(filename, mode='rb')

打开BGZR文件进行读取、写入或附加。

如果请求文本模式,为了避免多字节字符,将硬编码为使用“latin1”编码,并按原样传递“r”和“n”(不实现通用新行模式)。

如果你的数据是UTF-8或任何其他不兼容的编码,你必须使用二进制模式,并自己解码适当的片段。

Bio.bgzf.make_virtual_offset(block_start_offset, within_block_offset)

计算从块开始和块偏移内的BGZF虚拟偏移。

BAM索引方案使用64位“虚拟偏置”记录读取位置,用C术语来说包括:

块_start_off << 16| in_bank_off

这里,bank_start_oft是BGZR块开始(使用最多64-16 = 48位的无符号整数)和(解压缩)块内的_bank_oft(无符号16位整数)的文件偏置。

>>> make_virtual_offset(0, 0)
0
>>> make_virtual_offset(0, 1)
1
>>> make_virtual_offset(0, 2**16 - 1)
65535
>>> make_virtual_offset(0, 2**16)
Traceback (most recent call last):
...
ValueError: Require 0 <= within_block_offset < 2**16, got 65536
>>> 65536 == make_virtual_offset(1, 0)
True
>>> 65537 == make_virtual_offset(1, 1)
True
>>> 131071 == make_virtual_offset(1, 2**16 - 1)
True
>>> 6553600000 == make_virtual_offset(100000, 0)
True
>>> 6553600001 == make_virtual_offset(100000, 1)
True
>>> 6553600010 == make_virtual_offset(100000, 10)
True
>>> make_virtual_offset(2**48, 0)
Traceback (most recent call last):
...
ValueError: Require 0 <= block_start_offset < 2**48, got 281474976710656
Bio.bgzf.split_virtual_offset(virtual_offset)

将64位BGZR虚拟偏差分为块开始和块内偏差。

>>> (100000, 0) == split_virtual_offset(6553600000)
True
>>> (100000, 10) == split_virtual_offset(6553600010)
True
Bio.bgzf.BgzfBlocks(handle)

检查BGZR块的低级调试功能。

预计使用内置打开函数以二进制读取模式打开BGZR压缩文件。请勿使用此bgzf模块或gZip模块的打开函数中的handle,这些函数将初始化该文件。

作为迭代器返回块开始偏置(请参阅虚拟偏置)、块长度(将这些添加到下一个块的开始)和块内容的解压缩长度(在BGZR中限制为65536)-每个BGZR块一个字节组。

>>> from builtins import open
>>> handle = open("SamBam/ex1.bam", "rb")
>>> for values in BgzfBlocks(handle):
...     print("Raw start %i, raw length %i; data start %i, data length %i" % values)
Raw start 0, raw length 18239; data start 0, data length 65536
Raw start 18239, raw length 18223; data start 65536, data length 65536
Raw start 36462, raw length 18017; data start 131072, data length 65536
Raw start 54479, raw length 17342; data start 196608, data length 65536
Raw start 71821, raw length 17715; data start 262144, data length 65536
Raw start 89536, raw length 17728; data start 327680, data length 65536
Raw start 107264, raw length 17292; data start 393216, data length 63398
Raw start 124556, raw length 28; data start 456614, data length 0
>>> handle.close()

我们可以间接地看出这个文件来自旧版本的samtools,因为所有块(除了最后一个和虚拟空的收件箱标记块)都是65536字节。 后来的版本避免在两个块之间拆分读取,并为标头赋予自己的块(有助于加快替换标头)。您可以在使用samtools 0.1.18创建的ex1_refresh.bam中看到这一点:

samtools视图-b ex1.bam > ex1_refresh.bam

>>> handle = open("SamBam/ex1_refresh.bam", "rb")
>>> for values in BgzfBlocks(handle):
...     print("Raw start %i, raw length %i; data start %i, data length %i" % values)
Raw start 0, raw length 53; data start 0, data length 38
Raw start 53, raw length 18195; data start 38, data length 65434
Raw start 18248, raw length 18190; data start 65472, data length 65409
Raw start 36438, raw length 18004; data start 130881, data length 65483
Raw start 54442, raw length 17353; data start 196364, data length 65519
Raw start 71795, raw length 17708; data start 261883, data length 65411
Raw start 89503, raw length 17709; data start 327294, data length 65466
Raw start 107212, raw length 17390; data start 392760, data length 63854
Raw start 124602, raw length 28; data start 456614, data length 0
>>> handle.close()

上面的示例没有嵌入式的Sam标头(因此第一个块非常小,只有38个字节的解压缩数据),而下一个示例则有(更大的103个字节的块)。请注意,其余块显示相同的大小(它们包含相同的读取数据):

>>> handle = open("SamBam/ex1_header.bam", "rb")
>>> for values in BgzfBlocks(handle):
...     print("Raw start %i, raw length %i; data start %i, data length %i" % values)
Raw start 0, raw length 104; data start 0, data length 103
Raw start 104, raw length 18195; data start 103, data length 65434
Raw start 18299, raw length 18190; data start 65537, data length 65409
Raw start 36489, raw length 18004; data start 130946, data length 65483
Raw start 54493, raw length 17353; data start 196429, data length 65519
Raw start 71846, raw length 17708; data start 261948, data length 65411
Raw start 89554, raw length 17709; data start 327359, data length 65466
Raw start 107263, raw length 17390; data start 392825, data length 63854
Raw start 124653, raw length 28; data start 456679, data length 0
>>> handle.close()
class Bio.bgzf.BgzfReader(filename=None, mode='r', fileobj=None, max_cache=100)

基类:object

BGZR阅读器,就像只读手柄,但寻求/告诉不同。

让我们使用BgzfBlocks函数来浏览一个示例BAM文件中的BGZR块,

>>> from builtins import open
>>> handle = open("SamBam/ex1.bam", "rb")
>>> for values in BgzfBlocks(handle):
...     print("Raw start %i, raw length %i; data start %i, data length %i" % values)
Raw start 0, raw length 18239; data start 0, data length 65536
Raw start 18239, raw length 18223; data start 65536, data length 65536
Raw start 36462, raw length 18017; data start 131072, data length 65536
Raw start 54479, raw length 17342; data start 196608, data length 65536
Raw start 71821, raw length 17715; data start 262144, data length 65536
Raw start 89536, raw length 17728; data start 327680, data length 65536
Raw start 107264, raw length 17292; data start 393216, data length 63398
Raw start 124556, raw length 28; data start 456614, data length 0
>>> handle.close()

现在让我们看看如何使用此块信息跳转到解压缩的BAM文件的特定部分:

>>> handle = BgzfReader("SamBam/ex1.bam", "rb")
>>> assert 0 == handle.tell()
>>> magic = handle.read(4)
>>> assert 4 == handle.tell()

到目前为止,还没有什么奇怪的,我们在解压缩的BAM文件的开头使用了神奇标记,而且手柄位置是有意义的。然而,现在让我们通过读取65536个字节跳转到这个块的结尾,并将4个字节转入下一个块,

>>> data = handle.read(65536)
>>> len(data)
65536
>>> assert 1195311108 == handle.tell()

你期待4 + 65536 = 65540吗?这是一个BGZR 64位虚拟失调,这意味着:

>>> split_virtual_offset(1195311108)
(18239, 4)

您应该将18239指定为第二个BGZR块的开始,而4则是该块的偏差。另请参阅make_virtual_Offset,

>>> make_virtual_offset(18239, 4)
1195311108

让我们几乎跳回到文件的开头,

>>> make_virtual_offset(0, 2)
2
>>> handle.seek(2)
2
>>> handle.close()

请注意,您可以使用max_ache参数来限制内存中缓存的BGZR块的数量。默认值为100,并且由于每个块最多可达64 KB,因此默认缓存最多可占用6 MB的RAM。缓存对于一次通读文件并不重要,但对于提高随机访问的性能很重要。

__init__(filename=None, mode='r', fileobj=None, max_cache=100)

初始化类以读取BGZF文件。

您通常会使用顶级 bgzf.open(...) 函数将在内部调用这个类。不鼓励直接使用。

要么 filename (字符串)或 fileobj (二进制模式下的输入文件对象)必须提供参数,但不能同时提供参数。

论点 mode 控制数据是否以文本模式(“RT”、“TR”或默认“r”)或字节二进制模式(“rb”或“br”)的字符串形式返回。参数名称与内置的匹配 open(...) 和标准库 gzip.open(...) 功能

如果请求文本模式,为了避免多字节字符,将硬编码为使用“latin1”编码,并按原样传递“r”和“n”(不实现通用新行模式)。没有 encoding 论点

如果你的数据是UTF-8或任何其他不兼容的编码,你必须使用二进制模式,并自己解码适当的片段。

论点 max_cache 控制内存中要缓存的BGZR块的最大数量。每个块最多可达64 KB,因此默认的100个块可能占用最多6 MB的RAM。这对于高效的随机访问很重要,对于一次读取文件来说,小值是可以的。

tell()

返回64位无符号BGZR虚拟偏差。

seek(virtual_offset)

寻求64位无符号BGZR虚拟偏差。

read(size=-1)

BGZR模块的读取方法。

readline()

读取BGZR文件的一行。

__next__()

返回下一步。

__iter__()

迭代BGZR文件中的行。

close()

关闭BGZR文件。

seekable()

返回True表示BGZR支持随机访问。

isatty()

如果连接到ty设备,则返回True。

fileno()

返回integer文件描述符。

__enter__()

打开一个可使用WHITE声明操作的文件。

__exit__(type, value, traceback)

使用AND声明关闭文件。

__firstlineno__ = 493
__static_attributes__ = ('_block_raw_length', '_block_start_offset', '_buffer', '_buffers', '_handle', '_newline', '_text', '_within_block_offset', 'max_cache')
class Bio.bgzf.BgzfWriter(filename=None, mode='w', fileobj=None, compresslevel=6)

基类:object

定义BGZFWriter对象。

__init__(filename=None, mode='w', fileobj=None, compresslevel=6)

初始化课程。

write(data)

为类编写方法。

flush()

显式刷新数据。

close()

刷新数据,写入28字节BGZR收件箱标记,并关闭BGZR文件。

samtools将寻找一个神奇的收件箱标记,只是一个28字节的空BGZR块,如果缺少它,则警告BAM文件可能会被截断。除了samtools编写这个块之外,bgZip也是如此-所以这个实现也是如此。

tell()

返回BGZR 64位虚拟偏差。

seekable()

返回True表示BGZR支持随机访问。

isatty()

如果连接到ty设备,则返回True。

fileno()

返回integer文件描述符。

__enter__()

打开一个可使用WHITE声明操作的文件。

__exit__(type, value, traceback)

使用AND声明关闭文件。

__firstlineno__ = 794
__static_attributes__ = ('_buffer', '_handle', '_text', 'compresslevel')