转自：https://gitee.com/wangshun1121/TaxonomyPickUp
基因注释时，常常需要用blast比对一些大型数据库，例如NT、NR、TrEMBL等。这些数据库中的序列数目非常多，且包罗万象。于是，考虑根据先验信息，blast注释的时候，从这群数据库当中只选择特定物种的蛋白或核酸序列，即提高效率，又增加注释精度。
前期数据调研
无论NCBI的数据库和TrEMBL数据库，都提供了序列物种信息，而且物种表示方式统一为NCBI的TaxID。
NCBI提供的物种信息位于ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/，其中既有蛋白库物种信息和核酸库物种信息。蛋白库的物种信息prot.accession2taxid.gz。格式如下：
accession accession.version taxid gi
P26567 P26567.2 4577 1168978
P12208 P12208.1 3197 116525
P12210 P12210.1 4097 116527
P24064 P24064.2 4565 17374148
UniProt的数据库包括TrEMBL和SwissProt两部分，提供的注释信息文件idmapping_selected.tab.gz，表头信息有说明，摘要如下：
idmapping_selected.tab
We also provide this tab-delimited table which includes
the following mappings delimited by tab:
1\. UniProtKB-AC
2\. UniProtKB-ID
3\. GeneID (EntrezGene)
4\. RefSeq
5\. GI
6\. PDB
7\. GO
8\. UniRef100
9\. UniRef90
10\. UniRef50
11\. UniParc
12\. PIR
13\. NCBI-taxon
14\. MIM
15\. UniGene
16\. PubMed
17\. EMBL
18\. EMBL-CDS
19\. Ensembl
20\. Ensembl_TRS
21\. Ensembl_PRO
22\. Additional PubMed
或者idmapping.dat.gz，更是简单粗暴，直接3列内容：
1) idmapping.dat
This file has three columns, delimited by tab:
1\. UniProtKB-AC
2\. ID_type
3\. ID
where ID_type is the database name as appearing in UniProtKB cross-references,
and as supported by the ID mapping tool on the UniProt web site,
http://www.uniprot.org/mapping and where ID is the identifier in
that cross-referenced database.
制作物种信息总表
首先，从NCBI的服务器上下载taxdmp.zip。下载后解压之：
wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdmp.zip
unzip taxdmp.zip
该压缩包当中有好几个文件，但只需要其中两个names.dmp和nodes.dmp。用自己撰写的perl脚本处理之：
perl NodeExtract.pl |gzip -c >TaxInfo.txt.gz
得到的文件为TaxInfo.txt.gz，格式如下：
1 root 1 no rank
2 Bacteria 2 superkingdom-131567:no rank
6 Azorhizobium 7 genus-335928:family-356:order-28211:class-1224:phylum-2:superkingdom-131567:no rank
7 Azorhizobium caulinodans 8 species-6:genus-335928:family-356:order-28211:class-1224:phylum-2:superkingdom-131567:no rank
9 Buchnera aphidicola 8 species-32199:genus-1903409:family-91347:order-1236:class-1224:phylum-2:superkingdom-131567:no rank
第一列为TaxID，第二列为该TaxID的拉丁学名。第三列表示该TaxID所处的层级：即经过追溯多少次Parents到达分类级别的最高层级。第四列则是TaxID向上追溯的详细信息。不同分类层级之间以“-”间隔，最左边为该TaxID代表的层级，例如7 “Azorhizobium caulinodans”表示species，其上一级的TaxID为6，层级为genus，再向上为335928的family，依次类推。
这张基本表存储了NCBI的Taxonomy物种分类的详细层级信息，有了它，当指定了高级分类单位的时候，即可方便提取出下属所有低级分类单位，进而实现序列的筛选。
指定特定物种中的蛋白序列
UniProt库为2018年4月的版本，SwissProt和TrEMBL的蛋白数量与下面截图中完全一致： [图片上传失败...(image-cd764e-1692440765273)]
依据下面的命令，准备NR库和UniProt库中蛋白到TaxID的对应关系
gunzip -c prot.accession2taxid.gz|cut -f 1,3|gzip -c >NR2Tax.txt.gz
zgrep "NCBI_TaxID" idmapping.dat.gz|cut -f 1,3|gzip -c >UniProt2Tax.txt.gz
最终，我们需要的表格格式，只需要两列即可：accessionID和taxid。
我撰写了一份脚本 GetSeqFromTaxonomy.pl，用它来提取双翅目昆虫的蛋白序列。首先在NCBI的Taxonomy数据库输入双翅目的拉丁名Diptera，拿到该科的taxid为7147，然后输入下面的命令：
perl GetSeqFromTaxonomy.pl -taxid 7147 -outdir Diptera
需要等待一会儿，就会在./Diptera文件夹当中拿到双翅目昆虫的物种列表，以及NR库和UniProt库中的序列ID了。每个文件格式如下：
- Diptera.taxid.xls。记录的是双翅目所述的全部低级分类单位的拉丁名，taxID以及分类层级信息。第一列为TaxID，第二列为拉丁名，第三列为分类层级信息。居然有87833个分类单位之多！
- Diptera.NR.txt。NR库中双翅目的蛋白ID与其对应的taxID表，一共2487593条蛋白。
- Diptera.NR.stat.xls。双翅目物种在NR库中蛋白数目汇总统计。第一列为TaxID，第二列为蛋白数目，第三列拉丁名，第四列为分类层级。69675个双翅目物种在NR库中能够找得到蛋白序列。
- Diptera.UniProt.txt。TrEMBL+SwissProt库中双翅目的蛋白ID与其对应的taxID表，一共1297510条蛋白。
- Diptera.UniProt.stat.xls。双翅目物种在TrEMBL+SwissProt库中蛋白数目汇总统计。第一列为TaxID，第二列为蛋白数目，第三列拉丁名，第四列为分类层级。51446个双翅目物种在TrEMBL+SwissProt库中能够找得到蛋白序列。
将来，若提取更详细的序列信息，该脚本仍然可以改进。但是，若提高效率，需要学习数据库的知识了。
接下来，就可以将这批ID列表送入blast中做蛋白比对了。Diptera.UniProt.txt为TrEMBL和SwissProt两库公用的蛋白列表，注意：若该表取的物种太偏，会导致SwissProt的注释率严重下降，毕竟SwissProt为人工注释核查的蛋白库，以模式物种为主。
