TrEMBL

转自：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为人工注释核查的蛋白库，以模式物种为主。

https://gitee.com/wangshun1121/TaxonomyPickUp