最近作一个项目须要利用blastn
结果来画出进化树,这样就须要有物种名称。一种方法是利用blastn
输出的gi
去NCBI
查询获取到物种名称,虽然也是可行的,可是有没有简单一点的方法呢?笔者通过各类Google
终于找到了一种方法。node
首先有几个基础知识是须要掌握的:数据库
第一,用于构建blast
数据库的fasta
序列文件里面必须包含NCBI
的gi
信息,这个通常在NCBI
下载的序列文件里面都会包含,例如我下载的nt
的序列文件通常都是这样的:bash
>gi|1092|emb|X60725.1| Feline immunodeficiency virus ENV gene for envelope precursor glycoprotein ATGGCAGAAGGGTTTGTAGCCAATGGACAATGGATAGGACCAGAAGAAGCTGAAGAGTTAGTAGATTTTGAAATAGCAAC ACAAATGAATGAAGAAGGGCCACTAAATCCAGGAATAAACCCATTTAGGGTACCTGGAATAACAAAACAAGAAAAGCAGG AATATTGTAGCACAATGCAACCCAAATTACAAGCTCTAAGGAATGAAATTCAAGAGGTAAAACTGGAAGAAGGAAATGCA GGTAAGTTTAGAAGAGCAAGATTTTTAAGATACTCTGATGAAACTATATTGTCTCTGATTTACTTGTTCATAGGATATTT TAGATATTTAGTAGATAGAAAAAGGTTTGGGTCCTTAAGACATGACATAGATATAGAAGCACCTCAAGAAGAGTGTTATA
第二,NCBI
能够下载到gi
到物种ID
的映射文件,核酸序列的映射文件地址是:ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/gi_taxid_nucl.dmp.gzapp
第三,构建blast
索引时能够把映射文件包含进去,以便输出时能输出物种的ID
。less
$ ./makeblastdb -help -taxid_map <File_In> Text file mapping sequence IDs to taxonomy IDs. Format:<SequenceId> <TaxonomyId><newline> * Requires: parse_seqids * Incompatible with: taxid
第四,NCBI
能够下载到物种ID
到物种名称的映射文件,文件地址为:ftp://ftp.ncbi.nlm.nih.gov/blast/db/taxdb.tar.gzide
第五,blastn
等程序能够自定义输出格式ui
blastn(2.6.0+)
的输出格式有19
种,经过outfmt
来指定:code
$ ./blastn -help -outfmt <String> alignment view options: 0 = Pairwise, 1 = Query-anchored showing identities, 2 = Query-anchored no identities, 3 = Flat query-anchored showing identities, 4 = Flat query-anchored no identities, 5 = BLAST XML, 6 = Tabular, 7 = Tabular with comment lines, 8 = Seqalign (Text ASN.1), 9 = Seqalign (Binary ASN.1), 10 = Comma-separated values, 11 = BLAST archive (ASN.1), 12 = Seqalign (JSON), 13 = Multiple-file BLAST JSON, 14 = Multiple-file BLAST XML2, 15 = Single-file BLAST JSON, 16 = Single-file BLAST XML2, 17 = Sequence Alignment/Map (SAM), 18 = Organism Report
其中,格式6
、格式7
、格式10
、格式17
的输出条目是能够修改的:orm
$ ./blastn -help The supported format specifiers for options 6, 7 and 10 are: qseqid means Query Seq-id qgi means Query GI qacc means Query accesion qaccver means Query accesion.version qlen means Query sequence length sseqid means Subject Seq-id sallseqid means All subject Seq-id(s), separated by a ';' sgi means Subject GI sallgi means All subject GIs sacc means Subject accession saccver means Subject accession.version sallacc means All subject accessions slen means Subject sequence length qstart means Start of alignment in query qend means End of alignment in query sstart means Start of alignment in subject send means End of alignment in subject qseq means Aligned part of query sequence sseq means Aligned part of subject sequence evalue means Expect value bitscore means Bit score ...
对于这些格式,若是不提供输出条目的话,默认是这样的:索引
qaccver saccver pident length mismatch gapopen qstart qend sstart send evalue bitscore
也就是说,
$ ./blastn -outfmt 6
和
$ ./blastn -outfmt '6 qaccver saccver pident length mismatch gapopen qstart qend sstart send evalue bitscore'
是同样的,有了这个基础知识就能够显而易见地将ssciname
加入到输出条目中去了。
那么,根据先前讲到的基础知识很容易,就能够在blast
时,输出物种名称,如下以Linux
平台nt
的细菌、病毒基因序列为例,下载安装blast
没必要多说:
第一步:构建blast
数据库:
注意,这一步是很是吃内存的,据我估计直接使用
gi_taxid_nucl_dmp
这个文件构建blast
数据库的话至少要有100G
左右的可用内存。若是机器内存不够用的话,就要想办法把你的输入序列文件的gi
找出来,而后再根据gi
抽取出gi_taxid_nucl_dmp
相应的部分。这一步也是至关吃内存的,反正要具体对待,最终的目的是找到一个文本文件,使用空格或者tab键
分割,有两列,第一列是gi
,第二列是taxid
。
$ ./makeblastdb -in nt_BCT_VRL -dbtype nucl -parse_seqids -taxid_map gi_taxid_nucl.dmp # -in 指定输入序列 # -dbtype nucl表明这是核酸序列 # -parse_seqids表示须要从输入文件nt_BCT_VRL的每一条记录中获取gi # -taxid_map 指定gi到物种ID的映射文件
把生成的以nt_BCT_VRL.
开头的文件所有挪到$BLASTDB
里面去,不知道什么是$BLASTDB
的,请自行查看blast
文档。
第二步:把下载下来的taxdb.tar.gz
放入$BLASTDB
中。
taxdb.tar.gz
解压出来共有两个文件taxdb.btd
、taxdb.bti
两个文件,都放入$BLASTDB
便可。
第三步:运行blastn
$ ./blastn -db nt_BCT_VRL -query /tmp/a.fa -out /tmp/a.txt -outfmt '6 qseqid qlen sseqid sgi slen pident length mismatch gapopen qstart qend sstart send evalue bitscore staxid ssciname'
$ less /tmp/a.txt NODE_1_length_317_cov_19.776026 333 gi|384474605|emb|HE793683.1| 384474605 3360 96.396 333 12 0 1 333 389 721 1.98e-153 549 12305 Cucumber mosaic virus NODE_1_length_317_cov_19.776026 333 gi|25173579|emb|AJ511988.1| 25173579 3359 96.396 333 12 0 1 333 388 720 1.98e-153 549 12305 Cucumber mosaic virus NODE_1_length_317_cov_19.776026 333 gi|222028|dbj|D00356.1|MCVFRNA1 222028 3357 96.396 333 12 0 1 333 388 720 1.98e-153 549 12305 Cucumber mosaic virus NODE_1_length_317_cov_19.776026 333 gi|85677263|emb|AM183117.1| 85677263 3360 96.396 333 12 0 1 333 389 721 1.98e-153 549 367798 Cucumber mosaic virus (strain Ri-8) NODE_1_length_317_cov_19.776026 333 gi|60649875|emb|AJ888943.1| 60649875 732 96.096 333 13 0 1 333 216 548 9.19e-152 544 12305 Cucumber mosaic virus
能够看到物种名称已经在输入结果里面啦!