BLAST的大名對(duì)于做生物的同學(xué)可以說是如雷貫耳，哪怕非生信的同學(xué)也或多或少接觸過這個(gè)東西。我們通常會(huì)使用ncbi的blast在線工具進(jìn)行比對(duì)，但具有一個(gè)缺點(diǎn)就是很慢，因此我們常常會(huì)搭建一個(gè)本地blast系統(tǒng)進(jìn)行比對(duì)分析。

NCBI提供了一套用于運(yùn)行BLAST的命令行工具，稱為BLAST+。這允許用戶在自己的服務(wù)器上執(zhí)行BLAST搜索，而不受大小、容量和數(shù)據(jù)庫的限制。BLAST+可以使用命令行運(yùn)行，對(duì)于linux用戶可以說是相當(dāng)友好了。

blast+的軟件安裝

wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.10.1+-x64-linux.tar.gz
tar -zxvf ncbi-blast-2.10.1+-x64-linux.tar.gz

# 加入環(huán)境變量
export PATH=$PATH:$PWD/ncbi-blast-2.10.1+/bin/

nr/nt數(shù)據(jù)庫下載和使用

nr指的是protein sequence，而nt指的是nucleotide。所以nr/nt數(shù)據(jù)庫指的是兩個(gè)數(shù)據(jù)庫，前者是蛋白質(zhì)數(shù)據(jù)庫，后者是核酸序列數(shù)據(jù)庫。
計(jì)算資源足夠而帶寬不足時(shí)，建議下載fasta序列，然后使用blast的本地命令建立index。命令如下

wget -c ftp://ftp.ncbi.nih.gov/blast/db/FASTA/nr.gz # nr數(shù)據(jù)庫 
wget -c ftp://ftp.ncbi.nih.gov/blast/db/FASTA/nt.gz # nt數(shù)據(jù)庫

# 解壓縮
gunzip nr.gz
gunzip nt.gz

# 構(gòu)建本地blast nr/nt數(shù)據(jù)庫
mkdir nr_db; mkdir nt_db
makeblastdb -in nr -dbtype prot -title my_nr -parse_seqids -out ./nr_db/nr -logfile make_nr.log
makeblastdb -in nt -dbtype nucl -title my_nt -parse_seqids -out ./nt_db/nt -logfile make_nt.log

makeblastdb常用參數(shù)的解釋：

• -in: 后面接輸入文件，即我們要格式的fasta序列
• -dbtype: 后接序列類型，nucl為核酸，prot為蛋白
• -title: 給生成的blast數(shù)據(jù)庫起一個(gè)別名
• -parse_seqids: 幫助我們解析fa文件中，“>”后面的id信息
• -out: 后接數(shù)據(jù)庫名稱，起一個(gè)有意義的名稱，后續(xù)進(jìn)行blast比對(duì)時(shí)，-db參數(shù)后跟這個(gè)名稱
• -logfile: 日志文件，如果沒有則輸出到標(biāo)準(zhǔn)輸出（屏幕）上

帶寬足夠而計(jì)算資源不足，建議下載已經(jīng)建立好索引的nr/nt庫，命令如下：

mkdir nr_db; cd nr_db
wget -c ftp://ftp.ncbi.nih.gov/blast/db/nr*
for i in *gz
do
tar -zxvf $i
done

mkdir ../nt_db; cd ../nt_db
wget -c ftp://ftp.ncbi.nih.gov/blast/db/nt*
for i in *gz
do
tar -zxvf $i
done

除了可以使用wget下載已經(jīng)建立好索引的nr/nt庫還可以使用blast+自帶的update_blastdb.pl腳本下載。

update_blastdb.pl --decompress nr [*]
update_blastdb.pl --decompress nt [*]
# 下載壓縮后的nr/nt索引數(shù)據(jù)庫，并進(jìn)行解壓

更推薦wget下載方法，可以斷點(diǎn)繼續(xù)下載。

序列比對(duì)

以核酸序列為例進(jìn)行序列比對(duì)

blastn -query test.fa -out test.result -db nr -outfmt 6 -evalue 1e-5 -num_threads 4

常用參數(shù)介紹：

• -query: 需要比對(duì)的數(shù)據(jù)的文件路徑以及文件名
• -out: 輸出比對(duì)結(jié)果
• -db: 后面跟著建立好的索引數(shù)據(jù)庫
• -outfmt: 指定輸出結(jié)果的格式，此處指定為6，即常見的m8格式
• -evalue: 設(shè)置輸出結(jié)果的最小e-value值
• -num_threads: 指定比對(duì)時(shí)使用的線程數(shù)

此外，如果我們的目的是看測序數(shù)據(jù)的污染物來源，推薦加上下面兩個(gè)參數(shù)：

• -subject_besthit：數(shù)據(jù)庫中中的所有sequence，只輸出blast最優(yōu)的那個(gè)結(jié)果
• -max_target_seqs：設(shè)置每條 query reads所能比對(duì)上的最多的個(gè)數(shù)，這里推薦設(shè)置為1。

輸出結(jié)果的格式解析

輸出的比對(duì)結(jié)果，一共有12列，分別代表：

1. Query id：查詢序列ID標(biāo)識(shí)
2. Subject id：比對(duì)上的目標(biāo)序列ID標(biāo)識(shí)
3. % identity：序列比對(duì)的一致性百分比
4. alignment length：符合比對(duì)的比對(duì)區(qū)域的長度
5. mismatches：比對(duì)區(qū)域的錯(cuò)配數(shù)
6. gap openings：比對(duì)區(qū)域的gap數(shù)目
7. q. start：比對(duì)區(qū)域在查詢序列(Query id)上的起始位點(diǎn)
8. q. end：比對(duì)區(qū)域在查詢序列(Query id)上的終止位點(diǎn)
9. s. start：比對(duì)區(qū)域在目標(biāo)序列(Subject id)上的起始位點(diǎn)
10. s. end：比對(duì)區(qū)域在目標(biāo)序列(Subject id)上的終止位點(diǎn)
11. e-value：比對(duì)結(jié)果的期望值
12. bit score：比對(duì)結(jié)果的bit score值

我們一般看第1，3，11，12列。其中第一列告訴我們比對(duì)到了哪個(gè)物種，剩余三列告訴我們比對(duì)的可信程度。

導(dǎo)師愛說，垃圾進(jìn)垃圾出……但哪怕是垃圾，你也要知道它是什么垃圾。濕垃圾or干垃圾，或者像我一樣是個(gè)學(xué)術(shù)垃圾~