从NT_NR数据库中提取子库
原创
1 首先下载NCBI的taxonomy数据库
下载完解压缩，其中names.dmp和nodes.dmp两个文件很重要，是后续提取子库的基础
2
下载NCBI的TaxonKit软件，http://bioinf.shenwei.me/taxonkit/download/，linux系统直接解压
而后把names.dmp和nodes.dmp两个文件直接cp到~/.taxonkit下，其余的.dmp也可一并cp到~/.taxonkit下
cp taxdump/* ~/.taxonkit
3
下载NCBI的csvtk软件，http://bioinf.shenwei.me/csvtk/download/，linux系统也是直接解压，即可使用
4 （选择性步骤）NCBI
taxonomy数据库下还有accession2taxid库，这个库里面也有蛋白以及核酸的accession以及对应的分类id，但是经过尝试，采取这种方法提取的子库序列往往出乎意料的少，很可能是该库的accession与NT/NR库的accession不一致，前者可能冗余更多，因此该方法可忽略，见仁见智吧，下面给个例子，例如：
#从taxonomy数据库中的nucl_wgs.accession2taxid提取accession号
pigz -dc prot.accession2taxid.gz \
| csvtk grep -t -f taxid -P $id.taxid.txt \
| csvtk cut -t -f accession.version,taxid \
| sed 1d \
$id.acc2taxid.txt
cut -f 1 $id.acc2taxid.txt > $id.acc.txt5 查看，并提取动物的taxid
grep -P "|\s+aAnimal\w\s|" ~/.taxonkit/names.dmp
33208
taxonkit list --ids 33208 --indent "" > $id.taxid.txt6 从下载好的NT库提取对应的accession
blastdbcmd -db $NT -entry all -outfmt "%a %T" | csvtk grep -d ' ' -D ' ' -f
2 -P $id.taxid.txt \
| cut -d ' ' -f 1 \
$id.acc.txt7 从NT库提取完整的nt序列，并提取子库序列
blastdbcmd -db $NT -dbtype nucl -entry all -outfmt "%f" -out - | pigz -c >
nt.fa.gz
time cat <(echo) <(pigz -dc nt.fa.gz) \
| perl -e 'BEGIN{ $/ = "\n>"; <>; } while(<>){s/>$//; $i = index $, "\n";
$h = substr $, 0, $i; $s = substr $, $i+1; if ($h !~ />/) { print ">$";
next; }; $h = ">$h"; while($h =~ />(^ + .+?) ?(?=>|$)/g){ $h1 = $1; $h1 =~
s/^\W+//; print ">$h1\n$s";} } ' \
| seqkit grep -f $id.acc.txt -o nt.$id.fa.gz需要注意的是，这里又使用了seqkit软件。这种从NT库中还原的nt.fa序列里面有很多重复的头文件，例如
所以使用的话，还需要写个perl把这些序列拆开，最终形成nt.anmail.fa.gz
8 如果直接想构建子库，那么没必要搞序列，直接运行
blastdb_aliastool -gilist $id.acc.txt -db $NT -out nt_animal -title
nt_animal原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
腾讯云开发者
