食谱

分割大型文件

问题

随着现代测序技术的出现，生成非常大的数据集变得相对便宜且容易。事实上，有时一个文件中的数据量可能太大，一些在线资源，如 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) 来实现类似的功能，但如果总计数是批大小的精确倍数，则此函数会在末尾生成一个空文件。