阅读fastq的最快方法

Question

我正在尝试使用科基特-生物读取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

首先，为什么verify=False不改进运行时？

verify=False是scikit的I/O普遍接受的参数。它不特定于特定的文件格式。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，因为这样只会使事情更慢。

第二，是否有一种更快的方法使用科学生物来读取序列？

没有更快的方法来读取FASTQ文件与scikit-bio。有一个问题2探索的想法，可以加快这一点，但这些想法尚未落实。

scikit bio在读取FASTQ文件方面速度缓慢，因为它支持读取序列数据和可能跨越多行的质量分数。这会使阅读逻辑复杂化，并对性能产生影响。然而，包含多行数据的FASTQ文件已不再常见；Illumina过去用于生成这些文件，但他们现在更喜欢/建议编写FASTQ记录，这些记录正好是四行(序列头、序列数据、质量标头、质量分数)。事实上，这就是科学生物如何编写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

2

Answer 2

如果您想要计算fastq文件中每个唯一序列的出现次数，我建议尝试pyGATB的pyGATB解析器，以及标准库的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模块非常有效(根据我在生物信息学SE：https://bioinformatics.stackexchange.com/a/380/292中的基准测试)。

这个计数器应该像你的seq_dic，加上一些方便的方法。

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

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