食譜

分割大型檔案

問題

隨著現代定序技術的發展，產生非常大的數據集變得相對便宜和容易。事實上，有時一個檔案中的資料太多了，像 PLAN 或 TMHMM 這樣的線上資源會限制使用者查詢的大小。在這種情況下，能夠將序列檔案分割成一組較小的檔案非常有用，每個檔案都包含原始檔案序列的子集。

解決方案

有很多可能的方法可以解決這個普遍的問題，這個食譜使用一個產生器函數，以避免一次將所有數據都放在記憶體中。

def batch_iterator(iterator, batch_size):
    """Returns lists of length batch_size.

    This can be used on any iterator, for example to batch up
    SeqRecord objects from Bio.SeqIO.parse(...), or to batch
    Alignment objects from Bio.Align.parse(...), or simply
    lines from a file handle.

    This is a generator function, and it returns lists of the
    entries from the supplied iterator.  Each list will have
    batch_size entries, although the final list may be shorter.
    """
    batch = []
    for entry in iterator:
        batch.append(entry)
        if len(batch) == batch_size:
            yield batch
            batch = []
    if batch:
        yield batch

這是一個使用此函數分割大型 FASTQ 檔案的範例，SRR014849_1.fastq（來自此壓縮檔案）

from Bio import SeqIO

record_iter = SeqIO.parse(open("SRR014849_1.fastq"), "fastq")
for i, batch in enumerate(batch_iterator(record_iter, 10000)):
    filename = "group_%i.fastq" % (i + 1)
    with open(filename, "w") as handle:
        count = SeqIO.write(batch, handle, "fastq")
    print("Wrote %i records to %s" % (count, filename))

以及輸出結果

Wrote 10000 records to group_1.fastq
Wrote 10000 records to group_2.fastq
Wrote 10000 records to group_3.fastq
Wrote 10000 records to group_4.fastq
Wrote 10000 records to group_5.fastq
Wrote 10000 records to group_6.fastq
...
Wrote 7251 records to group_27.fastq

您可以修改此食譜以使用 Bio.SeqIO 支援的任何輸入和輸出格式，例如將大型 FASTA 檔案分割為 1000 個序列的單位

from Bio import SeqIO

record_iter = SeqIO.parse(open("large.fasta"), "fasta")
for i, batch in enumerate(batch_iterator(record_iter, 1000)):
    filename = "group_%i.fasta" % (i + 1)
    with open(filename, "w") as handle:
        count = SeqIO.write(batch, handle, "fasta")
    print("Wrote %i records to %s" % (count, filename))

運作方式

可以使用 list(SeqIO.parse(...)) 將檔案的全部內容讀取到記憶體中，然後將列表的切片寫為較小的檔案。對於大型檔案（如本食譜所討論的檔案），這將佔用大量的記憶體，相反地，我們可以定義一個產生器函數 batch_iterator()，它一次載入一個記錄，然後將其附加到列表中，重複此過程直到列表包含一個檔案份量的序列。

定義好該函數後，只需給它一個迭代器（在此範例中是 SeqIO.parse(...) 實例）並寫出它產生的記錄批次。

請注意，您可以使用內建的 Python 函數 itertools.islice(record_iter, batch_size) 來實現此功能，但如果總計數是批次大小的整數倍，則最後會得到一個空檔案。