美文网首页生信
BLAST数据库构建过程中sseqid无法关联到staxids的

BLAST数据库构建过程中sseqid无法关联到staxids的

作者: 徐诗芬 | 来源:发表于2021-07-02 17:21 被阅读0次

    blast本地数据库构建已经是一个很繁琐的事情了,比对NT或NR数据库还是会遇到各种问题:
    废话少说,先记录一下本地构建BLAST数据库:

    下载安装BLAST+的过程,并添加环境变量

    参考链接:徐洲更大神的
    还有这位大神的
    还有提取各类子库的方法,有空搞一下试试:
    https://www.jianshu.com/p/5a72f42e0412
    https://www.jianshu.com/p/ec7773da27b9

    wget https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.10.1/ncbi-blast-2.10.1+-x64-linux.tar.gz
    ###BLAST:
    tar zxvpf ncbi-blast-2.10.1+-x64-linux.tar.gz
    HOME = /Storage/data002/qiaogx/Software/
    export PATH=$PATH:$HOME/ncbi-blast-2.10.1+/bin
    mkdir /Storage/data002/qiaogx/Software/blastdb
    export BLASTDB=/Storage/data002/qiaogx/Software/blastdb
    

    第一步:下载相关数据包,解压并建立数据库索引

    #Database Download#
    cd /Storage/data002/qiaogx/Software/blastdb
    perl /Storage/data002/qiaogx/Software/ncbi-blast-2.10.1+/bin/update_blastdb.pl --passive --decompress nt ## --decompress下载nt数据库并解压的意思
    #perl /Storage/data002/qiaogx/Software/ncbi-blast-2.10.1+/bin/update_blastdb.pl --passive --decompress 16S_ribosomal_RNA
    #perl /Storage/data002/qiaogx/Software/ncbi-blast-2.10.1+/bin/update_blastdb.pl --passive --decompress ref_prok_rep_genomes
    #perl /Storage/data002/qiaogx/Software/ncbi-blast-2.10.1+/bin/update_blastdb.pl --passive --decompress nr
    #以上方法需要更新,而且update_blastdb.pl可能会报错,如果缺少模块直接conda安装,然后再运行perl 
    #另外一种方法:aspera高速下载nt,nr数据库 参考:https://zhuanlan.zhihu.com/p/245450890;https://cloud.tencent.com/developer/article/1654480
    数据库网址:https://ftp.ncbi.nlm.nih.gov/blast/db/  #注:在FASTA文件夹里面是最全的NT和NR库
    wget https://ak-delivery04-mul.dhe.ibm.com/sar/CMA/OSA/092u0/0/ibm-aspera-connect-3.10.0.180973-linux-g2.12-64.tar.gz
    tar xvf ibm-aspera-connect-3.10.0.180973-linux-g2.12-64.tar.gz
    sh ibm-aspera-connect-3.10.0.180973-linux-g2.12-64.sh
    /Storage/data002/qiaogx/.aspera/connect
    # 所有安装文件都在~/.aspera/connect目录下,添加环境变量
    echo 'PATH=$PATH:~/.aspera/connect/bin' >> ~/.bashrc
    which ascp # ~/.aspera/connect/bin/ascp; ~/.aspera/connect/etc/asperaweb_id_dsa.openssh
    ascp -v -k 1 -T -l 500m -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/blast/db/FASTA/nt.gz ./
    ascp -v -k 1 -T -l 500m -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/blast/db/FASTA/nr.gz ./
    #生成md5
    nohup md5sum nt.gz &
    #解压nt.gz文件
    nohup gunzip nt.gz &
    #需要下载物种注释信息,把下载的taxdb.tar.gz放入$BLASTDB中
    wget https://ftp.ncbi.nlm.nih.gov/blast/db/taxdb.tar.gz
    tar zxvf taxdb.tar.gz
    taxdb.tar.gz解压出来共有两个文件taxdb.btd、taxdb.bti两个文件,都放入$BLASTDB即可
    #Create a custom database from a multi-FASTA file,  build a blast database by a reference genome sequence
    makeblastdb -in nt -dbtype nucl -out nt -parse_seqids
    #-in #输入下载的数据序列;-dbtype #数据库序列类型; -out 输出自定义数据库名用于-db参数;-parse_seqids输入对应的序列ID或者说gi
    #检查是否创建成功
    blastdbcmd -db nt -info
    

    第二步:可以本地配置数据库索引,随时可直接调用本地数据库

    vi ~/.ncbirc #家目录下创建一个配置文件,内容如下:

    ; Start the section for BLAST configuration
    
    [BLAST]
    
    ; Specifies the path where BLAST databases are installed
    BLASTDB=/Storage/data002/qiaogx/Software/blastdb
    
    ; Specifies the data sources to use for automatic resolution
    
    ; for sequence identifiers
    DATA_LOADERS=blastdb
    
    ; Specifies the BLAST database to use resolve protein sequences
    BLASTDB_PROT_DATA_LOADER=/Storage/data002/qiaogx/Software/blastdb/nr
    
    ; Specifies the BLAST database to use resolve nucleic acid sequences
    BLASTDB_NUCL_DATA_LOADER=/Storage/data002/qiaogx/Software/blastdb/nt
    
    BATCH_SIZE=10G
    
    ; Windowmasker settings
    
    [WINDOW_MASKER]
    WINDOW_MASKER_PATH=/Storage/data002/qiaogx/Software/blastdb/windowmasker
    
    ; end of file
    

    BLAST的一些简单用法

    blastn -db nt -query test.fasta -use_index true -out megablast.out -outfmt '6 std' -evalue 1e-5 -num_threads 8
    #-best_hit_score_edge 0.01 -best_hit_overhang 0.25 # best BLAST match method
    #best hit的筛选参考:[选择blast比对第一条best hit的结果](https://www.jianshu.com/p/c75dca4607d9)
    #diamond的简单用法:
    diamond makedb --in nr.fa -d nr
    diamond blastx -d nr -q test.merged.transcript.fasta -o test.matches.m8
    # 但是diamond使用有限制,只能用于比对蛋白数据库
    
    #如果要提取数据库里面的fasta文件
    blastdbcmd -entry all -db nr -out nr.fsa
    blastdbcmd -entry u00001 -db nr -out u00001.fsa # 序列名(accession) u00001
    
    #如果要获取taxonomy信息,输入staxids后 获取taxonomy信息
    get_species_taxids.sh -t 91347 > 91347.txids
    

    但是,当我有一次用了megablast模块且需要taxonomyID信息时,blast结果输出有问题,staxids这一列都是0

    blastn -task megablast -word_size 30 -db nt -query changed.fa \
    -outfmt '6 qseqid staxids bitscore std' -perc_identity 70 -num_threads 16 -evalue 1e-25 \
    -out blast.70.1e25.out 
    
    1625216558(1).png

    一开始以为下载的NT库里面匹配上的hit没有物种注释信息,但是不可能所有的都没有注释信息,因此肯定不是数据库的问题,后来检查各种以前的脚本,发现用nt.00数据库(以前只下载了这一个子集)是有注释信息staxids的,然后突然发现以前创建了一个megablast的nt.00数据库的索引

    ##Indexed megaBLAST search## 此为用于megablast的索引,如果使用megablast的话一定要运行这一步!!!
    makembindex -input nt.00 -iformat blastdb  #已建立名为nt.00的数据库,用于-db后面的参数
    

    因此原来是我用了megablast却没有创建相应的索引!!!

    makembindex -input nt -iformat blastdb
    

    数据库很大,创建起来很麻烦!


    image.png

    最后有一个问题是如果blast过程中忘了输出物种信息,该怎么加上去呢?
    见链接:https://segmentfault.com/a/1190000012055972
    先下载一个gi_taxid_nucl.dmp,然后makeblastdb的时候加上参数-taxid_map gi_taxid_nucl.dmp,然后再进行以下操作。但下载网址好像有问题,可能改名了ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/gi_taxid_nucl.dmp.gz

    blastn -db nt -query /tmp/a.fa -out /tmp/a.txt -outfmt '6 qseqid qlen sseqid sgi slen pident length mismatch gapopen qstart qend sstart send evalue bitscore staxid ssciname'
    

    今天又是无所事事、浪费感情的一天!!!

    相关文章

      网友评论

        本文标题:BLAST数据库构建过程中sseqid无法关联到staxids的

        本文链接:https://www.haomeiwen.com/subject/lmzeultx.html