如何利用NR库快速进行物种鉴定

引言

在做基因或者蛋白质功能注释的时候，NR数据库包含非常重要的信息。但是有没有人跟我一样遇到相同的苦恼，就是在运行的时候特别耗时，经过一番搜索后，发现可以根据实际数据建立自己的NR子库，这样运行效率将大幅提高。那么怎么创建NR子库呢，跟着我一起来做吧。

01

NR 数 据 库 介 绍

目前我们可以从很多数据库中获取蛋白质的序列信息，比如GenPept, Swissprot, PIR, PDF, PDB 和 NCBI RefSeq等。但是在这些数据库之间，蛋白质序列存在冗余性，为了解决这个问题，NCBI构建了一个非冗余的蛋白质序列数据库，即NR（non-redundant proteins）。

我们可以看到最新版的NR库是刚刚更新的，压缩后的文件存储已经达到了73G，这是相当恐怖的。跟着来建立自己的NR子库吧~

02

准 备 工 作

首先我们下载NR数据库及需要用到的软件。

下载NR数据库

• https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz
• 下载后需要解压，fasta格式文件，存放路径 mypath/nr 。
• 既然是非冗余数据库，那么NR数据库中序列合并必须同时满足以下两个原则：

      两条序列有相同的长度；

      两条序列上的每个残基都是相同的。

下载分类数据库

• 物种及分类信息的文件：
  https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz
• 解压文件 taxdump.tar.gz：
  1mkdir  mypath/taxdump
  2tar  -zxvf  mypath/taxdump.tar.gz   -C mypath/taxdump
  

• taxdump 目录中有两个重要文件：
  names.dmp：记录物种名及其分类编号
  nodes.dmp：记录分类编号的分类节点信息

下载accession与taxid的对应关系

下载链接为：https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz

1gzip  -dc  prot.accession2taxid.gz >prot.accession2taxid

其他工具

taxonkit:https://github.com/shenwei356/taxonkit/releases/download/v0.6.0/taxonkit_linux_amd64.tar.gzcsvtk:https://github.com/shenwei356/csvtk/releases/download/v0.20.0/csvtk_linux_amd64.tar.gzncbi-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.提取病毒包含的所有物种编号

1mypath/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

1zcat 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数据库建索引

1mypath/ncbi-blast/bin/makeblastdb -in nr -dbtype prot -parse_seqids

建索引的步骤非常耗时，可能会花费好几个小时，需要耐心等待。

4.提取NR子库

1mypath/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。
运行命令：

1mypath/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

运行结果是m8格式的文件如下图：

该文件包含12列，每列信息简单介绍如下：1. qseqid：查询序列ID标识
2. sseqid：比对上的目标序列ID标识
3. pident：序列比对的一致性百分比
4. length：符合比对的比对区域的长度
5. mismatch：比对区域的错配数
6. gapopen：比对区域的gap数目
7. qstart：比对区域在查询序列(qseqid)上的起始位点
8. qend：比对区域在查询序列(qseqid)上的终止位点
9. 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

物 种 注 释

合并物种分类文件

names.dmp，nodes.dmp 这两个文件的内容我们上文已经介绍过，现在我们先来处理一下。将这两个文件的信息整合在一起，便于后续应用。
这里介绍网上一个脚本，首先运行

1git clone  https://gitee.com/wangshun1121/TaxonomyPickUp.git

命令后会在当前路径下生成目录mypath/TaxonomyPickUp/，然后将names.dmp，nodes.dmp 两个文件拷贝到目录mypath/TaxonomyPickUp/。最后运行

1perl NodeExtract.pl >TaxInfo.txt

大概一两分钟就得到结果，它存储了NCBI的Taxonomy物种分类的详细层级信息，格式如下：

• 第一列为TaxID，
• 第二列为该TaxID的拉丁学名。
• 第三列表示该TaxID所处的层级：即经过追溯多少次Parents到达分类级别的最高层级。
• 第四列则是TaxID向上追溯的详细信息。不同分类层级之间以“-”间隔，最左边为该TaxID代表的层级，例如7 “Azorhizobium caulinodans”表示species，其上一级的TaxID为6，层级为genus，再向上为335928的family，依次类推。

进行物种注释

比对之后我们可以继续做物种分类注释，需要用到上一步blastp的结果，prot.accession2taxid.gz，names.dmp，nodes.dmp 这几个文件了。接下来，我们只需要写个简单的脚本，就能从blastp比对结果mypath/test.blastp.out获取每条reads对应的物种分类信息了。
这里简单说下逻辑：

1. 从mypath/test.blastp.out文件前两列获取reads id 和accession对应关系；
2. prot.accession2taxid.gz文件包含accession和taxid的对应关系，结合1可以得到每条reads id和taxid的对应关系；
3. 上述TaxInfo.txt包含了taxid和物种的对应关系，结合2就能获取reads id和物种的关系了，还可以方便筛选所有层级的分类单位。

通过这个详细的描述，相信大家都已经掌握了整个过程，快动手试试吧。

参 考 资 料

https://gitee.com/wangshun1121/TaxonomyPickUp/tree/master

https://ftp.ncbi.nlm.nih.gov

作者：charon审稿：童蒙编辑：angelica

往期精选

udp -葡萄糖6-脱氢酶功能缺失突变可引起隐性发展性癫痫性脑病
如何拆分pacbio的数据
如何用matplotlib绘制图片
超快组装软件的使用--hifiasm软件
ARGO-工作流部署与管理工具
全长转录本的鉴定
DNA-RNA互作检测技术
ORF的前世今生

来都来了 点个赞再走吧~~~