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

约翰切斯

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

空中飞人

首先,为什么不验证=假,从而提高了运行时间

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

本文收集自互联网,转载请注明来源。

如有侵权,请联系 [email protected] 删除。

编辑于
0

我来说两句

0 条评论
登录 后参与评论

相关文章

性能:使用Python读取文件的最快方法

使用 Numpy / Opencv / scikit-image 更改 RGB 圖像像素值的最快方法是什麼

使用 OpenSSL 内存 BIO 写入和读取以空字符结尾的字符串的正确方法

使用C ++从文件读取时,切换字节序的最快方法是什么?

使用Python读取大型二进制文件的最快方法

使用Spark从大型hdfs目录中读取几行的最快方法是什么?

使用 Bio 库在 django 中读取上传的 fasta 文件

scikit-learn 的 .fit() 方法中使用了什么优化算法?

如何在scikit中使用内核密度估计作为一维聚类方法?

使用keras模型解决多标签问题的scikit学习链分类器的拟合方法错误

如何使用numpy.compress方法减少图像部分?(numpy + scikit-图像)

使用scikit-learn计算AUC的正确方法是哪一种?

我应该如何使用它的 `.components` 编写代码 scikit-learn PCA `.transform()` 方法?

读取网络数据的最快方法?

从文件读取long []的最快方法?

逐行读取STDIN的最快方法?

读取word文件的最快方法

使用点屏蔽栅格的最快方法

使用常见OpenOption组合的最快方法

在for循环中使用条件的最快方法

使用JDBC插入数据的最快方法

使用Python提取tar文件的最快方法

使用Java扫描端口的最快方法

使用JDBC遍历大表的最快方法

使用ffmpeg提取帧的最快方法?

使用mysql什么是最快的联接表的方法

使用Javascript ArrayBuffer执行XOR的最快方法

使用常量的最佳和最快的方法

有没有一种方法可以使用scikit-learn“重新创建”数据?