使用 scikit-bio 读取 fastq 的最快方法

Posted 2023-02-23

【中文标题】使用 scikit-bio 读取 fastq 的最快方法

【英文标题】：Fastest way to read a fastq with scikit-bio

【发布时间】：2016-08-25 17:02:05

【问题描述】：

我正在尝试使用 scikit-bio 读取 fastq 格式的文本文件。

鉴于它是一个相当大的文件，执行操作很慢。

最终，我试图将 fastq 文件复制到字典中：

f = 'Undetermined_S0_L001_I1_001.fastq'
seqs = skbio.io.read(f, format='fastq')

seq_dic = 
for seq in seqs:
    seq = str(seq)
    if seq in seq_dic.keys():
        seq_dic[seq] +=1
    else:
        seq_dic[seq] = 1

这里大部分时间是在读取文件的时候使用的：

%%time
f = 'Undetermined_S0_L001_I1_001.fastq'
seqs = skbio.io.read(f, format='fastq')

for seq in itertools.islice(seqs, 100000):
    seq

CPU times: user 46.2 s, sys: 334 ms, total: 46.5 s
Wall time: 47.8 s

我的理解是不验证序列会提高运行时间，但情况似乎并非如此：

%%time
f = 'Undetermined_S0_L001_I1_001.fastq'
seqs = skbio.io.read(f, format='fastq', verify=False, variant='illumina1.8')

for seq in itertools.islice(seqs, 100000):
    seq

CPU times: user 47 s, sys: 369 ms, total: 47.4 s
Wall time: 48.9 s

所以我的问题是，首先为什么 verify=False 不能提高运行时间，其次是否有更快的方法使用 scikit-bio 读取序列？

【参考方案1】：

首先为什么 verify=False 不能提高运行时间

verify=False 是 scikit-bio 的 I/O API 一般接受的参数。它并不特定于特定的文件格式。 verify=False 告诉 scikit-bio 不要调用文件格式的嗅探器来仔细检查文件是否为用户指定的格式。来自文档 [1]：

verify : bool, optional
    When True, will double check the format if provided.

所以verify=False 不会关闭序列数据验证；它关闭文件格式嗅探器验证。使用verify=False，您将获得最小的性能提升。

seqs = skbio.io.read(f, format='fastq', verify=False, variant='illumina1.8') 将生成skbio.Sequence 对象的生成器。不执行序列字母验证，因为skbio.Sequence 没有字母，所以这不是您的性能瓶颈所在。请注意，如果您想将 FASTQ 文件读入特定类型的生物序列（DNA、RNA 或蛋白质），您可以传递 constructor=skbio.DNA（例如）。这将对相关序列类型执行字母验证，目前在阅读时无法关闭。由于您遇到性能问题，我不建议您传递constructor，因为这只会进一步减慢速度。

其次，有没有更快的方法使用 scikit-bio 读取序列？

没有更快的方法来使用 scikit-bio 读取 FASTQ 文件。存在一个问题 [2] 探索可以加快这一进程的想法，但这些想法尚未实施。

scikit-bio 读取 FASTQ 文件的速度很慢，因为它支持读取可能跨越多行的序列数据和质量分数。这使读取逻辑复杂化并影响性能。然而，具有多行数据的 FASTQ 文件不再常见； Illumina 曾经生成这些文件，但他们现在更喜欢/建议编写正好四行的 FASTQ 记录（序列标题、序列数据、质量标题、质量分数）。事实上，这就是 scikit-bio 写入 FASTQ 数据的方式。使用这种更简单的记录格式，读取 FASTQ 文件变得更快、更容易。 scikit-bio 读取 FASTQ 文件的速度也很慢，因为它会解码和验证质量分数。它还将序列数据和质量分数存储在skbio.Sequence 对象中，这会产生性能开销。

在您的情况下，您不需要解码质量分数，并且您可能有一个包含简单四行记录的 FASTQ 文件。这是一个兼容 Python 3 的生成器，它读取 FASTQ 文件并将序列数据生成为 Python 字符串：

import itertools

def read_fastq_seqs(filepath):
    with open(filepath, 'r') as fh:
        for seq_header, seq, qual_header, qual in itertools.zip_longest(*[fh] * 4):
            if any(line is None for line in (seq_header, seq, qual_header, qual)):
                raise Exception(
                    "Number of lines in FASTQ file must be multiple of four "
                    "(i.e., each record must be exactly four lines long).")
            if not seq_header.startswith('@'):
                raise Exception("Invalid FASTQ sequence header: %r" % seq_header)
            if qual_header != '+\n':
                raise Exception("Invalid FASTQ quality header: %r" % qual_header)
            if qual == '\n':
                raise Exception("FASTQ record is missing quality scores.")

            yield seq.rstrip('\n')

如果您确定您的文件是有效的 FASTQ 文件，其中包含正好四行长的记录，则可以在此处删除验证检查。

这与您的问题无关，但我想指出您的计数器逻辑中可能存在错误。当你第一次看到一个序列时，你的计数器设置为零而不是 1。我认为逻辑应该是：

if seq in seq_dic:  # calling .keys() is necessary
    seq_dic[seq] +=1
else:
    seq_dic[seq] = 1

[1]http://scikit-bio.org/docs/latest/generated/skbio.io.registry.read.html

[2]https://github.com/biocore/scikit-bio/issues/907

【讨论】：

感谢@jairideout，这是一个很好的回应/解决方案！

【参考方案2】：

如果你想计算fastq文件中每个唯一序列的出现次数，我建议尝试pyGATB的Bank解析器，以及collections中方便的Counter对象标准库的模块。

#!/usr/bin/env python3

from collections import Counter
from gatb import Bank

# (gunzipping is done transparently)
seq_counter = Counter(seq.sequence for seq in Bank("SRR077487_2.filt.fastq.gz"))

这对于 python 模块来说非常有效（根据我在 Bioinformatics SE 中的基准：https://bioinformatics.stackexchange.com/a/380/292）。

这个计数器的行为应该像你的seq_dic、plus some convenient methods。

例如，如果我想打印 10 个最频繁的序列及其计数：

print(*seq_counter.most_common(10), sep="\n")