如何利用NR库快速进行物种鉴定
作者：charon
审稿：童蒙
编辑：angelica
引言
在做基因或者蛋白质功能注释的时候，NR数据库包含非常重要的信息。但是有没有人跟我一样遇到相同的苦恼，就是在运行的时候特别耗时，经过一番搜索后，发现可以根据实际数据建立自己的NR子库，这样运行效率将大幅提高。那么怎么创建NR子库呢，跟着我一起来做吧。
01 NR 数 据 库 介 绍
目前我们可以从很多数据库中获取蛋白质的序列信息，比如GenPept, Swissprot, PIR, PDF, PDB 和 NCBI RefSeq等。但是在这些数据库之间，蛋白质序列存在冗余性，为了解决这个问题，NCBI构建了一个非冗余的蛋白质序列数据库，即NR（non-redundant proteins）。
我们可以看到最新版的NR库是刚刚更新的，压缩后的文件存储已经达到了73G，这是相当恐怖的。跟着来建立自己的NR子库吧~
02 准 备 工 作
首先我们下载NR数据库及需要用到的软件。
2.1 下载NR数据库
https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz
下载后需要解压，fasta格式文件，存放路径 mypath/nr 。
既然是非冗余数据库，那么NR数据库中序列合并必须同时满足以下两个原则：
- 两条序列有相同的长度；
- 两条序列上的每个残基都是相同的。
2.2 下载分类数据库
-
物种及分类信息的文件：
https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz -
解压文件 taxdump.tar.gz：
1 mkdir mypath/taxdump
2 tar -zxvf mypath/taxdump.tar.gz -C mypath/taxdump
- taxdump 目录中有两个重要文件：
names.dmp：记录物种名及其分类编号
nodes.dmp：记录分类编号的分类节点信息
2.3 下载accession与taxid的对应关系
下载链接为：
https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz
1gzip -dc prot.accession2taxid.gz >prot.accession2taxid
2.4 其他工具
- taxonkit:
https://github.com/shenwei356/taxonkit/releases/download/v0.6.0/taxonkit_linux_amd64.tar.gz - csvtk:
https://github.com/shenwei356/csvtk/releases/download/v0.20.0/csvtk_linux_amd64.tar.gz - ncbi-blast:
ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.10.0+-x64-linux.tar.gz - 这几个工具都是下载后解压可以直接用的，解压目录路分别为 ：
mypath/taxonkit/taxonkit;mypath/csvtk/csvtk;mypath/ncbi-blast/ 。
03 构 建 命 令
我们以构建病毒子库为例，首先要知道病毒的tax id为10239，这个需要要先验知识，或者在taxdump/names.dmp文件中查找。
1.提取病毒包含的所有物种编号
1 mypath/taxonkit/taxonkit list -j 2 --ids 10239 --indent "" --data-dir mypath/taxdump/ >mypath/virus.list
这行命令运行很快，大概十秒钟，看下结果，总共提取出来205975个taxid。
-j ，线程，给2足以；
--ids，需要提取的分类的taxid；
--data-dir，该目录下必须包含文件names.dmp和nodes.dmp；
--indent，提取的物种编号缩进位置，这个参数很重要，记得一定要设置为空 ""
2.根据taxid获取对应accession.version
1 zcat mypath/prot.accession2taxid.gz | mypath/csvtk/csvtk -t grep -f taxid -P mypath/virus.list | mypath/csvtk/csvtk -t cut -f accession.version > mypath/virus.taxid.acc.txt
这一步运行也比较快，我这不到10分钟就运行结束了，这一步提取出来的accession就能对应到NR数据库中的标识符，总共提取出来了6065494行数据。
很显然这个数目比上一步taxid的数目多了很多，这是因为一个accession可能会对应多个taxid，这个我们在mypath/prot.accession2taxid.gz 中就可以看到。
3.NR数据库建索引
1 mypath/ncbi-blast/bin/makeblastdb -in nr -dbtype prot -parse_seqids
建索引的步骤非常耗时，可能会花费好几个小时，需要耐心等待。
4.提取NR子库
1 mypath/ncbi-blast/bin/blastdb_aliastool -seqidlist mypath/virus.taxid.acc.txt -db nr -out nr_virus -title nr_virus
最后这一步也会比较耗时，大概运行了2个多小时。运行结束后只生成了文件 mypath/nr_virus.pal，不要怀疑自己，确实只生成一个文件。
04 功 能 注 释
对于NR数据库，一般都是blast直接比对，获取基因或者蛋白质的功能注释信息，更近一步，如果有Taxonomy 数据库的数据，我们还可以做物种注释。
接下来我们试试用提取的病毒子库做功能注释吧，测试数据我选取了100条蛋白质序列，命名为test.fa。
运行命令：
1 mypath/ncbi-blast/bin/blastp -query mypath/test.fa -db mypath/nr_virus -evalue 1e-5 -num_alignments 5 -outfmt 6 -out mypath/test.blastp.out
该文件包含12列，每列信息简单介绍如下：
- qseqid：查询序列ID标识
- sseqid：比对上的目标序列ID标识
- pident：序列比对的一致性百分比
- length：符合比对的比对区域的长度
- mismatch：比对区域的错配数
- gapopen：比对区域的gap数目
- qstart：比对区域在查询序列(qseqid)上的起始位点
- qend：比对区域在查询序列(qseqid)上的终止位点
- sstart：比对区域在目标序列(sseqid)上的起始位点
10.send：比对区域在目标序列(sseqid)上的终止位点
11.evalue：比对结果的期望值，evalue越小，越有可能是真实的相似序列
12.bitscore：比对结果的bit score值
注 意 ! ! !
- 因为我们用的ncbi-blast+中的blastp，-outfmt 要设置为6，如果你用的是blast中的blastp工具，要设置 -m 8；
- 一般我们比较关注第3、11、12列信息；
- 一条序列比对可能有多条结果，我们可以按照bitscore筛选得分最高的一条即可。
子库运行时间不到50分钟，而选择整个NR数据库比对时，将近4个小时才结束，很显然，选择子库比对大大提升了比对效率。因此，如果我们对研究的物种比较确定，提取子库再进行比对，将会提升运行效率降低集群压力。
05 物 种 注 释
5.1 合并物种分类文件
names.dmp，nodes.dmp 这两个文件的内容我们上文已经介绍过，现在我们先来处理一下。将这两个文件的信息整合在一起，便于后续应用。
这里介绍网上一个脚本，首先运行
1 git clone https://gitee.com/wangshun1121/TaxonomyPickUp.git
命令后会在当前路径下生成目录mypath/TaxonomyPickUp/，然后将names.dmp，nodes.dmp 两个文件拷贝到目录mypath/TaxonomyPickUp/。
最后运行
1 perl NodeExtract.pl >TaxInfo.txt
大概一两分钟就得到结果，它存储了NCBI的Taxonomy物种分类的详细层级信息，格式如下：
- 第一列为TaxID，
- 第二列为该TaxID的拉丁学名。
- 第三列表示该TaxID所处的层级：即经过追溯多少次Parents到达分类级别的最高层级。
- 第四列则是TaxID向上追溯的详细信息。不同分类层级之间以“-”间隔，最左边为该TaxID代表的层级，例如7 “Azorhizobium caulinodans”表示species，其上一级的TaxID为6，层级为genus，再向上为335928的family，依次类推。
5.2 进行物种注释
比对之后我们可以继续做物种分类注释，需要用到上一步blastp的结果，prot.accession2taxid.gz，names.dmp，nodes.dmp 这几个文件了。接下来，我们只需要写个简单的脚本，就能从blastp比对结果mypath/test.blastp.out获取每条reads对应的物种分类信息了。这里简单说下逻辑：
- 从mypath/test.blastp.out文件前两列获取reads id 和accession对应关系；
- prot.accession2taxid.gz文件包含accession和taxid的对应关系，结合1可以得到每条reads id和taxid的对应关系；
- 上述TaxInfo.txt包含了taxid和物种的对应关系，结合2就能获取reads id和物种的关系了，还可以方便筛选所有层级的分类单位。
通过这个详细的描述，相信大家都已经掌握了整个过程，快动手试试吧。
参 考 资 料
https://gitee.com/wangshun1121/TaxonomyPickUp/tree/master
https://ftp.ncbi.nlm.nih.gov
该文来源于“生信阿拉丁”，关注公众号，第一时间查收“新款”生信学习干货。
