快速序列比對之BLAST

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

NCBI提供了一套用于運行BLAST的命令行工具泛源，稱為BLAST+拔妥。這允許用戶在自己的服務(wù)器上執(zhí)行BLAST搜索，而不受大小达箍、容量和數(shù)據(jù)庫的限制没龙。BLAST+可以使用命令行運行，對于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ù)庫指的是兩個數(shù)據(jù)庫赃磨，前者是蛋白質(zhì)數(shù)據(jù)庫筝家，后者是核酸序列數(shù)據(jù)庫。
計算資源足夠而帶寬不足時邻辉，建議下載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ù)庫起一個別名
-parse_seqids: 幫助我們解析fa文件中吱瘩，“>”后面的id信息
-out: 后接數(shù)據(jù)庫名稱道伟，起一個有意義的名稱，后續(xù)進行blast比對時使碾，-db參數(shù)后跟這個名稱
-logfile: 日志文件蜜徽，如果沒有則輸出到標(biāo)準(zhǔn)輸出（屏幕）上

帶寬足夠而計算資源不足祝懂，建議下載已經(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ù)庫嫂易，并進行解壓

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

序列比對

以核酸序列為例進行序列比對

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

常用參數(shù)介紹：

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

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

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

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

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

Query id：查詢序列ID標(biāo)識
Subject id：比對上的目標(biāo)序列ID標(biāo)識
% identity：序列比對的一致性百分比
alignment length：符合比對的比對區(qū)域的長度
mismatches：比對區(qū)域的錯配數(shù)
gap openings：比對區(qū)域的gap數(shù)目
q. start：比對區(qū)域在查詢序列(Query id)上的起始位點
q. end：比對區(qū)域在查詢序列(Query id)上的終止位點
s. start：比對區(qū)域在目標(biāo)序列(Subject id)上的起始位點
s. end：比對區(qū)域在目標(biāo)序列(Subject id)上的終止位點
e-value：比對結(jié)果的期望值
bit score：比對結(jié)果的bit score值

我們一般看第1驾霜，3，11买置，12列粪糙。其中第一列告訴我們比對到了哪個物種，剩余三列告訴我們比對的可信程度忿项。

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