Bio.bgzf 模块

读取和写入 BGZF 压缩文件（BAM 中使用的 GZIP 变体）。

SAM/BAM 文件格式（序列比对/映射）采用纯文本格式（SAM）和压缩二进制格式（BAM）。后者使用 gzip 压缩的修改形式，称为 BGZF（分块 GNU 压缩格式），可以应用于任何文件格式以提供高效随机访问的压缩。BGZF 与 SAM/BAM 文件格式一起在 https://samtools.sourceforge.net/SAM1.pdf 中有描述。

在使用 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 字节或 256 PB，每个块的解压缩大小最多为 2^16 字节或 64KB。请注意，这与 BGZF“BC”字段大小相匹配，该大小将每个块的压缩大小限制为 2^16 字节，允许 BAM 文件使用 BGZF 而不进行 gzip 压缩（这对内存中的中间文件很有用，可以减少 CPU 负载）。

关于命名空间的警告

在 Python 中使用“from XXX import *”被认为是一个坏习惯，因为它会污染命名空间。这对于 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 的超集），或可能的其他字符映射（未实现）。

此外，与 Python 3 上的默认文本模式不同，我们不尝试实现通用换行模式。这将 Windows（CR LF 或“rn”）、Unix（仅 LF，“n”）或旧 Mac（仅 CR，“r”）等各种操作系统换行约定转换为仅 LF（”n”）。这里我们有同样的问题 - 块末尾的“r”是否是不完整的 Windows 风格换行符？

相反，您将获得 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 术语中包含

block_start_offset << 16 | within_block_offset

这里 block_start_offset 是 BGZF 块起始位置的文件偏移量（使用最多 64-16 = 48 位的无符号整数），而 within_block_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 避免将一个读分割到两个块中，并为头部分配一个独立的块（有助于加速头部替换）。您可以在使用 samtools 0.1.18 创建的 ex1_refresh.bam 中看到这一点

samtools view -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

BGZF 读取器，像只读句柄一样，但 seek/tell 不同。

让我们使用 BgzfBlocks 函数来窥视一个示例 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 是第二个 BGZF 块的起始位置，而 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 块的数量。默认值为 100，由于每个块最多可以为 64kb，因此默认缓存最多可能占用 6MB 的 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 控制内存中缓存的 BGZF 块的最大数量。每个块最多可以为 64kb，因此默认的 100 个块可能最多占用 6MB 的 RAM。这对于高效的随机访问很重要，对于一次性读取整个文件，较小的值就足够了。

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 语句关闭文件。