我正在尝试使用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
本文收集自互联网,转载请注明来源。
如有侵权,请联系 [email protected] 删除。
我来说两句