用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

首先为什么不验证=错误改善运行时间

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