Bio.bgzf 模組

讀取和寫入 BGZF 壓縮檔案 (BAM 中使用的 GZIP 變體)。

SAM/BAM 檔案格式 (序列比對/映射) 有純文字格式 (SAM) 和壓縮二進位格式 (BAM)。後者使用一種修改過的 gzip 壓縮形式，稱為 BGZF (區塊 GNU Zip 格式)，它可以應用於任何檔案格式，以提供高效隨機存取的壓縮。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 位整數）。

這限制了您只能使用最後一個區塊從 2^48 位元組開始的 BAM 檔案，或 256 PB，並且每個區塊的解壓縮大小最大為 2^16 位元組，或 64 KB。請注意，這與 BGZF「BC」欄位大小相符，該欄位將每個區塊的壓縮大小限制為 2^16 位元組，允許 BAM 檔案使用沒有 gzip 壓縮的 BGZF（這對於記憶體中的中間檔案很有用，以減少 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 虛擬偏移量的問題在於，雖然它們可以比較（哪個偏移量在檔案中先出現），但您不能安全地減去它們以取得它們之間資料的大小，也不能加上/減去相對偏移量。

當然，您可以使用 BgzfReader 使用 Bio.SeqIO 解析此檔案，儘管除非您想索引 BGZF 壓縮的序列檔案，否則使用 gzip.open(…) 沒有任何好處

>>> 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 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 陳述式關閉檔案。