如何在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 
   Text file mapping sequence IDs to taxonomy IDs.
   Format: 
    * 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 
   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代表这是核酸序列