Fastq序列提取

喜欢就点关注吧!

根据fastq序列的id,从原始fastq中提取序列这个操作,应该是大家在处理序列文件的过程中经常遇到的。如果大家用过Biopython,应该知道Bio模块在做fastq这些文件的处理时非常方便。但是有时序列达到几百万几千万条的时候,Bio的速度可能就无法满足要求了。

还是举个例子比较好,我从比对筛选过滤之后的bam文件中提取了第一列序列名,保存为id.name文件,想根据这个id文件从原始的fastq文件(单端)raw.fastq中把序列提出来。这里id.name中id数目42万左右,raw.fastq序列数1000万左右:

我首先写了一个脚本:(这里要用到biopython模块以及pandas模块,如果没安装的话可以装一下anaconda,它已经集成了这些常用包了,安装教程及使用见这里Anaconda:解决你装包的烦恼)

extract_fastq_reads_by_bam_id.py:

运行时间:

两分钟,感觉有点久,然后我查了下Bio中其实有针对fastq快速处理的FastqGeneralIterator,于是我使用FastqGeneralIterator又写了一个脚本:

extract_fastq_reads_by_bam_id_v1.py:

运行时间:

25秒!是不是感觉到人生苦短,我用Python这句话的真谛了,不过还没完。

其实还有更快的,只是不是Python了,而是Brian Bushnell用java写的工具包BBmap

(膜拜大神Brian Bushnell三秒):

这里很多参数的意义都很明了,include=t是提取id.name中的序列,include=f是提取非id.name中的序列,这里我们应该用t。filterbyname.sh还可以操作双端fastq,具体的可以用-h查看参数,作者写的很清楚呢。

运行时间:

11秒!

所以如果大家觉得麻烦也可以装一下bbmap,但其实Biopython已经很优秀了!

今天的"科普"就到这里,希望大家多多支持关注我们,如果有什么问题的话就给我留言吧。ps:(这期实在是太忙了,下期有时间的话,给大家补上Numba的向量运算以及bbmap的安装和小用法,谢谢大家!)

感谢您的阅读,欢迎点赞、评论、支持和转发!!

觉得本文好看,请点这里

  • 发表于:
  • 原文链接https://kuaibao.qq.com/s/20191113A03QE400?refer=cp_1026
  • 腾讯「云+社区」是腾讯内容开放平台帐号(企鹅号)传播渠道之一,根据《腾讯内容开放平台服务协议》转载发布内容。
  • 如有侵权,请联系 yunjia_community@tencent.com 删除。

扫码关注云+社区

领取腾讯云代金券