# TaxonomyPickUp

**Repository Path**: wangshun1121/TaxonomyPickUp

## Basic Information

- **Project Name**: TaxonomyPickUp
- **Description**: 从NR库或UniProt库中寻找指定物种中的蛋白序列ID，用于后续blast比对的加速
- **Primary Language**: Unknown
- **License**: Not specified
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 6
- **Forks**: 0
- **Created**: 2018-04-20
- **Last Updated**: 2024-11-23

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

基因注释时，常常需要用blast比对一些大型数据库，例如NT、NR、TrEMBL等。这些数据库中的序列数目非常多，且包罗万象。于是，考虑根据先验信息，blast注释的时候，从这群数据库当中只选择特定物种的蛋白或核酸序列，即提高效率，又增加注释精度。

## 前期数据调研

无论NCBI的数据库和TrEMBL数据库，都提供了序列物种信息，而且物种表示方式统一为NCBI的TaxID。

NCBI提供的物种信息位于ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/，其中既有蛋白库物种信息和核酸库物种信息。蛋白库的物种信息[**prot.accession2taxid.gz**](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**](ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/idmapping/idmapping_selected.tab.gz)，表头信息有[说明](ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/idmapping/README)，摘要如下：

```
2) 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**](ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/idmapping/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**](ftp://ftp.ncbi.nih.gov/pub/taxonomy/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的蛋白数量与下面截图中完全一致：
![20180423UniProt](http://ot0uvr7n3.bkt.clouddn.com/18-4-23/382369.jpg)

依据下面的命令，准备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为人工注释核查的蛋白库，以模式物种为主。