按物种拆分blast本地库
NCBI的nr、nt库包含了基本上所有现有已知物种的序列信息,在进行blast时,由于数据库较大,搜索时间会比较长,同时可能会有其他物种信息干扰。因此我们可以根据NCBI提供的物种分类号对其进行拆分,下面以nr库进行举例说明。
一、准备
1.1、数据:
- nr库:https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/
- accession2taxid: https://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz
- taxdump: https://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz
以上文件下载完后最好进行md5校验,检查文件完整性。
1.2、软件:
- blast: https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/
- 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、构建子库
目前构建子库有两种方法:
- 提取子库对应的fa文件,然后进行建库;
- 使用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"
三、错误回顾
在提取的过程中我发现了几个问题:
- 使用blastdbcmd提取寄生虫(Parasitus,Taxonomy id:708403)的fa文件时,发现有很多id提取不出,报错 Error: [blastdbcmd] Skipped xxxxx,这些ID在ncbi中能搜索到对应的序列,还有就是有个ID(AHM18782.1)不在种属分类号中,但是提取的序列中存在该条序列,该条序列的命名信息中包含QCF46254.1这个信息,该ID在种属分类号中;
- 使用seqtk直接从nr.fa文件中提取对应的序列,发现能提取的数据比问题1的多,其中提取不出的序列ID,在nr.fa中找不到,但是NCBI中存在;
- 使用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
针对于以上问题:
对于第一个问题,暂时没能找到具体的原因(发现使用ncbi建好的nr库,提取得到的序列多点,感觉自己用fa文件构建的库有点问题),对于后面的小问题,我发现在ncbi上AHM18782.1与QCF46254.1的序列是一致的,在nr库中应该为同个序列,其中AHM18782.1在nr库中能找到该id,QCF46254.1则找不到,其仅存在于AHM18782.1 id对应的命名信息中。这个问题将会影响到第二个问题的序列提取;- 对于第二个问题,由于nr库是非冗余蛋白库,所以ncbi中存在的序列不一定在nr库中存在;
- 对于第三个问题,也暂时没找到具体原因,猜测和第一个问题有关;
如果大家知道具体错误原因的,欢迎交流学习!
四、总结
目前建议直接下载ncbi建好的库直接进行提取fa文件会好点,自己建的库可能会有点问题(尤其是数据量大的情况下,不太好发现问题所在);或者直接根据ID在nr.fa文件中进行提取也问题不大;使用blastdb_aliastool命令目前暂时没能成功。
版权声明:本文转载请注明出处!
最新评论:
RLiu 2022年1月29日 04:25
(4) 但我也遇到了一些问题, 我的acc.txt文件中一共才52 millions 个条目,但是建出的库居然有469 millions 个条目, 居然多出了一个数量级! 难道一个accession id会对应多个条目吗? 我所以我打算按照你的步骤试一下, 有什么问题再来请教.
RLiu 2022年1月29日 04:25
(3) 后来我尝试用 blastdbcmd -db nr -entry_batch acc.txt -out seq.fa 建库, 与你分享的构建fa库的原理算是异曲同工吧,
RLiu 2022年1月29日 04:24
(2) 我用英文google了有无人遇到相似的问题, 发现有人早在16年就发现blastdb_aliastool -gilist 报错. 可我找的那些中文教程都是发布与20年以后的, 真是无语, 发教程前自己都不试一下吗...
RLiu 2022年1月29日 04:23
十分感谢你的分享. 终于找到一篇靠谱的教程了. 我跟着网上的众多教程尝试用specific organisms accession id 经由blastdb_aliastool -gilist命令运行会报错, 后来我自己提取了gilist再试这个命令还是会报错.
发表评论
电子邮件地址不会被公开。 必填项已用*标注