生信之旅

NCBI的nr、nt库包含了基本上所有现有已知物种的序列信息，在进行blast时，由于数据库较大，搜索时间会比较长，同时可能会有其他物种信息干扰。因此我们可以根据NCBI提供的物种分类号对其进行拆分，下面以nr库进行举例说明。

一、准备

1.1、数据：

1. nr库：https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/
2. accession2taxid: https://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz
3. taxdump: https://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz

以上文件下载完后最好进行md5校验，检查文件完整性。

1.2、软件：

1. blast: https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/
2. taxonkit: https://bioinf.shenwei.me/taxonkit/download/

二、构建本地库

2.1、解压

gzip -d nr.gz
mv nr nr.fa
gzip -d prot.accession2taxid.gz
gzip -d taxdump.tar.gz

2.2、nr本地库构建

makeblastdb -in nr.fa -dbtype prot -title nr -parse_seqids -hash_index -out nr -logfile nr.log

2.3、提取特定物种分类信息

这里我们以病毒为例，病毒对应的Taxonomy id是10239

taxonkit list --ids 10239 --indent "" > Viruses.taxid   # 提取病毒分类号
cat prot.accession2taxid |csvtk -t grep -f taxid -P Viruses.taxid |csvtk -t cut -f accession.version > Viruses.acc  # 获取对应的accession号

# 也可以使用blast自带的脚本提取
get_species_taxids.sh -n Viruses  # 获取Viruses的taxid，即10239
get_species_taxids.sh -t 10239 > Viruses.taxid

2.4、构建子库

目前构建子库有两种方法：

1. 提取子库对应的fa文件，然后进行建库；
2. 使用blastdb_aliastool创建数据库别名；

方法一：提取子库fa文件
对于该方法，我们也有两种方法进行提取。

• 1. 从构建好的nr库中提取
• 2.从nr.fa文件中提取

# 方法一
blastdbcmd -db nr -entry_batch Viruses_nr.acc  -out  Viruses_nr.fa

# 方法二
seqtk subseq nr.fa Viruses_nr.acc > Viruses_nr.fa

提取完后使用makeblastdb构建对应的子库即可。

但是这两种方法提取寄生虫（Parasitus，Taxonomy id：708403）的分类数据时，发现第一种方法提取的序列不全，不知道什么情况。

方法二：使用blastdb_aliastool创建数据库别名

blastdb_aliastool -db nr  -seqidlist Viruses_nr.acc -out Viruses -title "Viruses_nr"

三、错误回顾

在提取的过程中我发现了几个问题：

1. 使用blastdbcmd提取寄生虫（Parasitus，Taxonomy id：708403）的fa文件时，发现有很多id提取不出，报错  Error: [blastdbcmd] Skipped xxxxx，这些ID在ncbi中能搜索到对应的序列，还有就是有个ID（AHM18782.1）不在种属分类号中，但是提取的序列中存在该条序列，该条序列的命名信息中包含QCF46254.1这个信息，该ID在种属分类号中；
2. 使用seqtk直接从nr.fa文件中提取对应的序列，发现能提取的数据比问题1的多，其中提取不出的序列ID，在nr.fa中找不到，但是NCBI中存在；
3. 使用blastdb_aliastool创建寄生虫库的别名时，报错 BLAST Database error: BLASTDB alias file creation failed.  Some referenced files may be missing；

下载的nr、nt库等均使用了md5校验，文件是完整的。使用的blast版本为 blast 2.11.0, build Oct  6 2020 03:24:05

针对于以上问题：

1. 对于第一个问题，暂时没能找到具体的原因(发现使用ncbi建好的nr库，提取得到的序列多点，感觉自己用fa文件构建的库有点问题)，对于后面的小问题，我发现在ncbi上AHM18782.1与QCF46254.1的序列是一致的，在nr库中应该为同个序列，其中AHM18782.1在nr库中能找到该id，QCF46254.1则找不到，其仅存在于AHM18782.1 id对应的命名信息中。这个问题将会影响到第二个问题的序列提取；
2. 对于第二个问题，由于nr库是非冗余蛋白库，所以ncbi中存在的序列不一定在nr库中存在；
3. 对于第三个问题，也暂时没找到具体原因，猜测和第一个问题有关；

如果大家知道具体错误原因的，欢迎交流学习！

四、总结

目前建议直接下载ncbi建好的库直接进行提取fa文件会好点，自己建的库可能会有点问题（尤其是数据量大的情况下，不太好发现问题所在）；或者直接根据ID在nr.fa文件中进行提取也问题不大；使用blastdb_aliastool命令目前暂时没能成功。