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

Question

我正在尝试使用 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 读取序列的更快方法？

Answer 1

first why isn't verify=False improving run time

verify=False一般是scikit-bio的I/OAPI接受的参数。它不特定于特定的文件格式。 verify=False 告诉 scikit-bio 不要调用文件格式的嗅探器来 double-check 该文件是用户指定的格式。来自文档 [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 object 的生成器。未执行序列字母表验证，因为 skbio.Sequence 没有字母表，因此这不是您的性能瓶颈所在。请注意，如果您想将 FASTQ 文件读入特定类型的生物序列（DNA、RNA 或蛋白质），您可以传递 constructor=skbio.DNA（例如）。此将对相关序列类型执行字母表验证，目前无法在阅读时将其关闭。由于您遇到性能问题，我不建议传递 constructor，因为这只会进一步降低速度。

and second is there a faster way using scikit-bio to read sequences?

没有比 scikit-bio 更快的读取 FASTQ 文件的方法了。有一个问题 [2] 探索可以加快速度的想法，但这些想法尚未实施。

scikit-bio 在读取 FASTQ 文件时速度较慢，因为它支持读取可能跨越多行的序列数据和质量分数。这会使读取逻辑复杂化并且会影响性能。然而，具有 multi-line 数据的 FASTQ 文件不再常见； Illumina 过去常常生成这些文件，但他们现在 prefer/recommend 写入恰好四行的 FASTQ 记录（序列 header、序列数据、质量 header、质量分数）。事实上，scikit-bio就是这样写FASTQ数据的。使用这种更简单的记录格式，可以更快更轻松地读取 FASTQ 文件。 scikit-bio 读取 FASTQ 文件的速度也很慢，因为它会解码和验证质量分数。它还将序列数据和质量分数存储在 skbio.Sequence object 中，这会产生性能开销。

在您的情况下，您不需要解码质量分数，并且您可能有一个包含简单 four-line 记录的 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

Answer 2

如果你想计算 fastq 文件中每个唯一序列的出现次数，我建议尝试使用 Bank parser of pyGATB，以及 [=] 中方便的 Counter 对象13=] 标准库的模块。

#!/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")

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

Fastest way to read a fastq with scikit-bio

bioinformatics

python-3.x

skbio