blast本地?cái)?shù)據(jù)庫構(gòu)建已經(jīng)是一個(gè)很繁瑣的事情了，比對NT或NR數(shù)據(jù)庫還是會(huì)遇到各種問題：
廢話少說，先記錄一下本地構(gòu)建BLAST數(shù)據(jù)庫：

下載安裝BLAST+的過程，并添加環(huán)境變量

參考鏈接：徐洲更大神的
還有這位大神的
還有提取各類子庫的方法，有空搞一下試試：
http://www.itdecent.cn/p/5a72f42e0412
http://www.itdecent.cn/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

第一步：下載相關(guān)數(shù)據(jù)包，解壓并建立數(shù)據(jù)庫索引

#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數(shù)據(jù)庫并解壓的意思
#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可能會(huì)報(bào)錯(cuò)，如果缺少模塊直接conda安裝，然后再運(yùn)行perl 
#另外一種方法：aspera高速下載nt,nr數(shù)據(jù)庫 參考：https://zhuanlan.zhihu.com/p/245450890；https://cloud.tencent.com/developer/article/1654480
數(shù)據(jù)庫網(wǎng)址：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目錄下，添加環(huán)境變量
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解壓出來共有兩個(gè)文件taxdb.btd、taxdb.bti兩個(gè)文件，都放入$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 #輸入下載的數(shù)據(jù)序列；-dbtype #數(shù)據(jù)庫序列類型； -out 輸出自定義數(shù)據(jù)庫名用于-db參數(shù)；-parse_seqids輸入對應(yīng)的序列ID或者說gi
#檢查是否創(chuàng)建成功
blastdbcmd -db nt -info

第二步：可以本地配置數(shù)據(jù)庫索引，隨時(shí)可直接調(diào)用本地?cái)?shù)據(jù)庫

vi ~/.ncbirc #家目錄下創(chuàng)建一個(gè)配置文件，內(nèi)容如下：

; 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的結(jié)果](http://www.itdecent.cn/p/c75dca4607d9)
#diamond的簡單用法：
diamond makedb --in nr.fa -d nr
diamond blastx -d nr -q test.merged.transcript.fasta -o test.matches.m8
# 但是diamond使用有限制，只能用于比對蛋白數(shù)據(jù)庫

#如果要提取數(shù)據(jù)庫里面的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

但是，當(dāng)我有一次用了megablast模塊且需要taxonomyID信息時(shí)，blast結(jié)果輸出有問題，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沒有物種注釋信息，但是不可能所有的都沒有注釋信息，因此肯定不是數(shù)據(jù)庫的問題，后來檢查各種以前的腳本，發(fā)現(xiàn)用nt.00數(shù)據(jù)庫（以前只下載了這一個(gè)子集）是有注釋信息staxids的，然后突然發(fā)現(xiàn)以前創(chuàng)建了一個(gè)megablast的nt.00數(shù)據(jù)庫的索引

##Indexed megaBLAST search## 此為用于megablast的索引，如果使用megablast的話一定要運(yùn)行這一步?。?！
makembindex -input nt.00 -iformat blastdb  #已建立名為nt.00的數(shù)據(jù)庫，用于-db后面的參數(shù)

因此原來是我用了megablast卻沒有創(chuàng)建相應(yīng)的索引?。。?/h1>

makembindex -input nt -iformat blastdb

數(shù)據(jù)庫很大，創(chuàng)建起來很麻煩！

image.png

最后有一個(gè)問題是如果blast過程中忘了輸出物種信息，該怎么加上去呢？
見鏈接：https://segmentfault.com/a/1190000012055972
先下載一個(gè)gi_taxid_nucl.dmp，然后makeblastdb的時(shí)候加上參數(shù)-taxid_map gi_taxid_nucl.dmp，然后再進(jìn)行以下操作。但下載網(wǎng)址好像有問題，可能改名了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'

今天又是無所事事、浪費(fèi)感情的一天?。?！