创建NCBI的Nr数据库的子集数据库

1. 为什么要做Nr子集数据库

NCBI官网仅提供Nr全数据库。该数据库太大,将物种的蛋白序列使用Blastp比对到Nr数据库非常消耗计算和时间。对1个蛋白序列可能需要1个CPU计算半个小时。若对全基因组2万个基因分析,普通台式机8个CPU要计算2000/(2*8*24)=52天。这主要是由于Nr数据库太大导致的。为了能尽快得到Nr注释结果,可以按物种分类将Nr数据库分割成子集数据库,能得到更快的比对速度。以下是创建Nr子集数据库的步骤。

2. 创建Nr子集数据库的步骤

从NCBI下载Nr数据库和分类数据库文件

cd /opt/biosoft/Nr_database
# 下载Nr数据库(FASTA文件)
ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --host=ftp.ncbi.nih.gov --user=anonftp --mode=recv /blast/db/FASTA/nr.gz ./

# 下载NCBI的分类数据库文件
ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --host=ftp.ncbi.nih.gov --user=anonftp --mode=recv /pub/taxonomy/taxdump.tar.gz ./
ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --host=ftp.ncbi.nih.gov --user=anonftp --mode=recv /pub/taxonomy/accession2taxid/prot.accession2taxid.gz ./

# 解压缩两个NCBI的分类数据库文件
gzip -dc prot.accession2taxid.gz > prot.accession2taxid
mkdir ~/.taxonkit
tar zxf taxdump.tar.gz -C ~/.taxonkit
# 其主要有效文件有两个:
# names.dmp 记录物种名及其分类编号
# nodes.dmp 记录分类编号的节点信息
# 查看~/.taxonkit/names.dmp文件,使用关键词检索得到目标类的分类编号,例如:
# fungi 4751             # grep -P "\|\s+[fF]ungi\w*\s*\|" ~/.taxonkit/names.dmp
# plants 3193            # grep -P "\|\s+[pP]lant\w*\s*\|" ~/.taxonkit/names.dmp
# animals 33208          # grep -P "\|\s+[aA]nimal\w*\s*\|" ~/.taxonkit/names.dmp

下载并安装NCBI分类数据库解析软件TaxonKitTaxonKit,并解析nodes.dmp文件的物种节点信息,得到指定类的所有物种列表信息。

# 下载并安装NCBI分类数据库解析软件TaxonKit
wget https://github.com/shenwei356/taxonkit/releases/download/v0.2.4/taxonkit_linux_amd64.tar.gz
tar zxvf taxonkit_linux_amd64.tar.gz

# 提取古菌(2157)、细菌(2)和病毒(10239)这几个大类下的所有物种编号。
./taxonkit list -j 8 --ids 2,2157,10239 > sub.meta.list

# 再编写程序extract_sub_data_from_Nr.pl获得列表中物种在Nr数据库中的序列信息。
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.meta.list --acc2taxid prot.accession2taxid - > nr_meta.fasta

提取fungi/plants/animals子集

./taxonkit list -j 8 --ids 4751 > sub.fungi.list
./taxonkit list -j 8 --ids 3193 > sub.plants.list
./taxonkit list -j 8 --ids 33208 > sub.animals.list
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.fungi.list --acc2taxid prot.accession2taxid - > nr_fungi.fasta
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.plants.list --acc2taxid prot.accession2taxid - > nr_plants.fasta
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.animals.list --acc2taxid prot.accession2taxid - > nr_animals.fasta

使用makeblastdb创建blast本地数据库

makeblastdb -in nr_fungi.fasta -dbtype prot -title nr_fungi -parse_seqids -out nr_fungi_`date +%Y%m%d` -logfile nr_fungi_`date +%Y%m%d`.log
makeblastdb -in nr_fungi.fasta -dbtype prot -title nr_plants -parse_seqids -out nr_plants_`date +%Y%m%d` -logfile nr_plants_`date +%Y%m%d`.log
makeblastdb -in nr_fungi.fasta -dbtype prot -title nr_animals -parse_seqids -out nr_animals_`date +%Y%m%d` -logfile nr_animals_`date +%Y%m%d`.log