python-3.x - 使用 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/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，因为这只会进一步降低速度。

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

没有比使用 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