如何在blast输出结果中添加物种名称

最近做一个项目需要利用blastn结果来画出进化树，这样就需要有物种名称。一种方法是利用blastn输出的gi去NCBI查询获取到物种名称，虽然也是可行的，但是有没有简单一点的方法呢？笔者经过各种Google终于找到了一种方法。

1. 所需要的基础知识

首先有几个基础知识是需要掌握的：

第一，用于构建blast数据库的fasta序列文件里面必须包含NCBI的gi信息，这个一般在NCBI下载的序列文件里面都会包含，例如我下载的nt的序列文件一般都是这样的：

>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.gz

第三，构建blast索引时可以把映射文件包含进去，以便输出时能输出物种的ID。

$ ./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.gz

第五，blastn等程序可以自定义输出格式

blastn(2.6.0+)的输出格式有19种，通过outfmt来指定：

$ ./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的输出条目是可以修改的：

$ ./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加入到输出条目中去了。

2. 实战

那么，根据先前讲到的基础知识很容易，就可以在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

可以看到物种名称已经在输入结果里面啦！