method by id是什么_根据id快速提取fastq序列

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

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

$ wc -l id.name426648 id.name$ wc -l raw.fastq 41867248 raw.fastq

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

extract_fastq_reads_by_bam_id.py：

import pandas as pdimport sysfrom Bio import SeqIOdf_id=pd.read_csv(sys.argv[1],names=["name"])#input id file id.name name=sys.argv[1].split(".")[0]#prefix of output filename_list=set(df_id["name"].values)#all idsout_handle = open(name+".ext.fq", "w+")#open a file handle to write in readsfor seq_record in SeqIO.parse(sys.argv[2],"fastq"):#input raw.fastq if seq_record.id in name_list:  SeqIO.write(seq_record,out_handle,"fastq")#write in readsout_handle.close()

运行时间：

$ time python3 extract_fastq_reads_by_bam_id.py id.name raw.fastqpython3 extract_fastq_reads_by_bam_id.py id.name 156.89s user 4.10s system 102% cpu 2:37.37 total

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

extract_fastq_reads_by_bam_id_v1.py：

import pandas as pdimport sysfrom Bio import SeqIOfrom Bio.SeqIO.QualityIO import FastqGeneralIteratordf_id=pd.read_csv(sys.argv[1],names=["name"])name=sys.argv[1].split(".")[0]name_list=set(df_id["name"].values)out_handle = open(name+".MT.fq", "w+")for title,seq,qual in FastqGeneralIterator(open(sys.argv[2])):#different title1=title.split(" ")[0] if title1 in name_list:  seq_record="\n".join(["@"+title, seq, "+", qual])  out_handle.write(seq_record)  out_handle.write("\n")out_handle.close()

运行时间：

$ time python3 extract_fastq_reads_by_bam_id_v1.py id.name raw.fastqpython3 extract_fastq_reads_by_bam_id_v1.py id.name 24.79s user 2.85s system 108% cpu 25.532 total

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

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

(膜拜大神Brian Bushnell三秒)：

$ time ~/tool/bbmap/filterbyname.sh in=raw.fastq out=raw.ext.fq names=id.name include=t

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

运行时间：

Time:  11.618 seconds.Reads Processed:      10466k   900.95k reads/secBases Processed:       1046m   90.10m bases/secReads Out:          426648Bases Out:          42664800

11秒！

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

今天的"科普"就到这里，希望大家多多支持关注我们，如果有什么问题的话就给我留言吧。ps：(这期实在是太忙了，下期有时间的话，给大家补上Numba的向量运算以及bbmap的安装和小用法，谢谢大家！) 公众号ID：生物信息学扫码关注最新动态