按物种拆分blast本地库

按物种拆分blast本地库
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"


三、不用分库

简单的比对,省事的话,blast支持指定物种 taxids 来做比对,这样效果也不错:


## blastn 查询需要比对的污染物种库  taxids
# get_species_taxids.sh -n Bacteria #2 
# get_species_taxids.sh -n Archaea #2157
# get_species_taxids.sh -n Viruses #10239
# 得到物种taxids 列表
# get_species_taxids.sh -t 2 >Bacteria.taxids
# get_species_taxids.sh -t 2157 >Archaea.taxids
# get_species_taxids.sh -t 10239 >Viruses.taxids
# 合并 
# cat Bacteria.taxids Archaea.taxids  Viruses.taxids >all.taxids
#blast 比对 
blastn -query  novel.fasta \
  -db  /share/work/database/blastdb/nt/nt \
  -taxidlist all.taxids \
  -out novel.nt.blast \
  -evalue 1e-5  \
  -perc_identity 0.8 \
  -task megablast \
-outfmt '6 std qcovs stitle staxid' \
-max_target_seqs 5 -num_threads 10


  • 发表于 2023-08-15 15:23
  • 阅读 ( 1835 )
  • 分类:软件工具

2 条评论

请先 登录 后评论
omicsgene
omicsgene

生物信息

693 篇文章

作家榜 »

  1. omicsgene 693 文章
  2. 安生水 341 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. 红橙子 78 文章
  6. xun 76 文章
  7. rzx 74 文章
  8. CORNERSTONE 72 文章