如何在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

能夠看到物種名稱已經在輸入結果裏面啦！