Bio.bgzf模块¶

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

SAM/BAM文件格式(序列比对/映射)以纯文本格式(SAM)和压缩二进制格式(BAM)提供。后者使用修改后的gzip压缩形式，称为BGZF(Block GNU Zip Format)，它可以应用于任何文件格式，以提供高效随机访问的压缩。在http://samtools.sourceforge.net/SAM1.pdf上描述了bgzf和sam/bam文件格式

在使用BGZF文件进行随机访问之前，请阅读下面有关“虚拟偏移”的文本。

本模块的目的¶

Python gzip库可用于读取BGZF文件，因为对于解压缩而言，它们只是(专门的)gzip文件。本模块旨在促进对BGZF文件的随机访问(使用“虚拟偏移”的概念)，以及写入BGZF文件(这意味着使用大小合适的gzip块，并在gzip头文件中写入额外的“bc”字段)。与gzip库一样，zlib库在内部使用。

除了随机访问和写入BAM文件所需之外，BGZF格式还可以用于其他顺序数据(从一个记录接着另一个记录的意义上讲)，例如Bio.SeqIO中支持的大多数序列数据格式(如FASTA、FASTQ、GenBank等)或大型MAF比对。

Bio.SeqIO索引函数使用此模块支持BGZF文件。

BGZF技术简介¶

gzip文件格式允许多个压缩块，每个压缩块可以是一个独立的gzip文件。有趣的是，这意味着您可以使用Unix cat 通过连接两个或多个gzip文件将它们合并为一个文件。此外，每个挡路可以有几个压缩级别中的一个(包括未压缩，由于gzip头文件的缘故，这实际上会占用更多一点的空间)。

bam的设计者意识到，虽然随机访问存储在传统gzip文件中的数据很慢，但将文件分成gzip块将允许快速随机访问每个挡路。要访问特定的解压缩数据，您只需要知道它从哪个挡路开始(gzip挡路开始的偏移量)，以及需要阅读挡路(解压缩后的)内容的程度。

这样做的一个问题是高效地找到gzip挡路大小。您可以使用标准的gzip文件进行解压，但这需要对每个挡路进行解压缩--这会相当慢。此外，典型的gzip文件可能使用非常大的块。

bgzf的不同之处在于，每个gzip挡路的压缩大小被限制为2^16字节，并且gzip头部中的一个额外的‘BC’字段记录了这个大小。传统的解压缩工具可以忽略这一点，并像解压缩任何其他gzip文件一样解压缩该文件。

关键是您可以查看第一个bgzf挡路，从这个‘bc’头找出它有多大，从而立即搜索到第二个挡路，依此类推。

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

这将您限制为bam文件，其中最后一个挡路以2^48字节(或256PB)开始，每个挡路的解压缩大小最多为2^16字节(64KB)。注意，这与bgzf‘bc’字段大小匹配，后者将每个挡路的压缩大小限制为2^16字节，允许bam文件使用不带gzip压缩的bgzf(对于内存中的中间文件很有用，以减少cpu负载)。

有关命名空间的警告¶

使用“from XXX import”被认为是个坏主意 * “在Python中，因为它污染了命名空间。这是Bio.bgzf(和标准Python库gzip)的一个真正问题，因为它们包含一个名为open的函数，即假设您这样做：

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

或,

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

请注意，OPEN函数已被替换。如果需要，您可以通过导入内置的OPEN函数来“修复”此问题：

>>> from builtins import open

但是，我们建议改用显式名称空间，例如

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

示例¶

这是一个用BGZF压缩的普通GenBank文件，所以可以用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()

我们也可以使用BGZF阅读器访问该文件-但请注意文件偏移量，下面将对此进行说明：

>>> 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()

请注意，句柄的偏移量看起来与BGZF文件不同。这就把我们带到了BGZF的关键一点，那就是挡路的架构：

>>> 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个字节，用作‘文件结束’(EOF)标记。如果缺少此文件，则可能是您的BGZF文件不完整。

通过提前读取70,000字节，我们进入第二个BGZF挡路，此时，BGZF虚拟偏移量开始看起来与gzip库公开的解压缩数据中的简单偏移量不同。

例如，考虑寻找解压缩位置196734。由于196734=65536+65536+65536+126=65536*3+126，这相当于跳过前三个块(在此特定示例中，解压缩后的大小均为65536-这并不总是成立)并从第四个挡路的第126字节开始(解压缩后)。对于BGZF，我们需要知道第四个挡路的偏移量55074和挡路内部的偏移量126后才能得到BGZF的虚拟偏移量。

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

因此，对于该bgzf文件，解压缩位置196734对应于虚拟偏移3609329790。但是，另一个具有不同内容的BGZF文件将或多或少地进行压缩，因此压缩块的大小将不同。这意味着未压缩偏移和压缩虚拟偏移之间的映射取决于您使用的BGZF文件。

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

BGZF虚拟偏移的问题是，虽然可以比较它们(哪个偏移在文件中最先出现)，但您不能安全地减去它们以获得它们之间的数据大小，也不能添加/减去相对偏移。

当然，您可以使用Bio.SeqIO使用BgzfReader解析此文件，尽管使用gzip.open(.)没有任何好处，除非您希望索引BGZF压缩序列文件：

>>> 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(.)一样，BGZF代码默认以二进制模式打开文件。

您可以请求以文本模式打开文件，但要注意这是硬编码为简单的“latin1”(也称为“iso-8859-1”)编码(包括所有ASCII字符)，该编码可以很好地与大多数西欧语言配合使用。但是，它与更广泛使用的UTF-8编码并不完全兼容。

在像UTF-8这样的可变宽度编码中，Unicode文本输出中的一些单个字符由原始二进制形式的多个字节表示。这在bgzf中是有问题的，因为我们不能总是孤立地解码每个挡路-一个unicode字符可能会被分成两个块。这甚至可能发生在固定宽度的Unicode编码中，因为BGZF挡路的大小不是固定的。

因此，此模块目前仅限于支持单字节Unicode编码，如ASCII、“latin1”(它是ASCII的超集)或其他可能的字符映射(未实现)。

此外，与Python3上的默认文本模式不同，我们不会尝试实现通用的换行模式。这会将各种操作系统新行约定，如Windows(CRLF或“rn”)、Unix(只有LF，“n”)或旧Mac(只有CR，“r”)转换为只有LF(“n”)。这里我们有同样的问题--挡路末尾的“r”是不完整的视窗风格的新行吗？

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

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

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

打开BGZF文件进行读取、写入或追加。

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

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

Bio.bgzf.make_virtual_offset(block_start_offset, within_block_offset)¶

从挡路起点和挡路偏移内计算bgzf虚拟偏移。

BAM索引方案使用64位“虚拟偏移量”记录读取位置，该虚拟偏移量在C术语中包括：

挡路_START_OFFSET<<16|在_挡路_OFFSET内

这里挡路_START_OFFSET是BGZF挡路START(无符号整数，最大使用64-16=48位)的文件偏移量，在(解压缩)挡路内的_挡路_OFFSET(无符号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位BGZF虚拟偏移划分为挡路开始&在挡路偏移内。

>>> (100000, 0) == split_virtual_offset(6553600000)
True
>>> (100000, 10) == split_virtual_offset(6553600010)
True

Bio.bgzf.BgzfBlocks(handle)¶

用于检查BGZF块的低级调试功能。

需要使用内置OPEN函数以二进制读取模式打开的BGZF压缩文件。不要使用此bgzf模块或gzip模块的open函数中的句柄，因为这将对文件进行解压缩。

作为迭代器返回挡路开始偏移量(请参见虚拟偏移量)、挡路长度(添加这些以作为下一个挡路的开始)和块内容的解压缩长度(在bgzf中限制为65536)，作为迭代器-每个bgzf挡路返回一个元组。

>>> 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，因为所有的块(除了最后一个块和虚拟的空EOF标记挡路)都是65536字节。以后的版本避免在两个块之间拆分读取，并为报头提供自己的挡路(对于加快替换报头很有用)。您可以在使用Samtools 0.1.18创建的ex1_fresh h.bam中看到：

Samtools view-b ex1.bam>ex1_fresh h.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

BGZF读取器，作用类似于只读句柄，但查找/告知不同。

让我们使用BgzfBlock函数在示例BAM文件中的BGZF块处达到峰值，

>>> 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，是吗？这是一个BGZF 64位虚拟偏移量，这意味着：

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

你应该会点18239作为第二届博鳌亚洲论坛挡路的起点，而4则是对这一挡路的补偿。另请参见make_virtual_Offset，

>>> make_virtual_offset(18239, 4)
1195311108

让我们回到几乎是文件开头的地方，

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

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

__init__(filename=None, mode='r', fileobj=None, max_cache=100)¶: 初始化类。

tell()¶: 返回64位无符号BGZF虚拟偏移量。

seek(virtual_offset)¶: 查找64位无符号BGZF虚拟偏移量。

read(size=-1)¶: BGZF模块的读取方法。

readline()¶: 读取BGZF文件的一行。

__next__()¶: 返回下一行。

__iter__()¶: 迭代BGZF文件中的各行。

close()¶: 关闭BGZF文件。

seekable()¶: 返回True，表示BGZF支持随机访问。

isatty()¶: 如果连接到TTY设备，则返回True。

fileno()¶: 返回整数文件描述符。

__enter__()¶: 打开可使用WITH语句操作的文件。

__exit__(type, value, traceback)¶: 使用WITH语句关闭文件。

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字节的BGZF EOF标记，关闭BGZF文件。

Samtools将查找一个神奇的EOF标记，它只是一个28字节的空BGZF挡路，如果它丢失了，就会警告bam文件可能被截断。除了Samtools编写这个挡路之外，bgzip也是如此-这个实现也是如此。

tell()¶: 返回BGZF 64位虚拟偏移量。

seekable()¶: 返回True，表示BGZF支持随机访问。

isatty()¶: 如果连接到TTY设备，则返回True。

fileno()¶: 返回整数文件描述符。

__enter__()¶: 打开可使用WITH语句操作的文件。

__exit__(type, value, traceback)¶: 使用WITH语句关闭文件。