Blast 單機版語法

最近在用blast比對TM protein database,
但下載了win單機版的blast[2.2.29+]版本後卻苦於不清楚指令 囧,

參考這篇文章 終於找到方案了,
而更詳細的blast指令說明文件可以參考官網的資料
(網址:http://www.ncbi.nlm.nih.gov/books/NBK1763/)

Blast比對的概念是這樣的,你輸入一個要比對的資料庫然後系統會跟預先建好的參考資料庫作比對,所以進行比對前要先建立好自己的參考資料庫。

而實際執行的方法是這樣的:

0. 很重要,不要再點那個執行檔了,沒有反應的,請先開一個命令列視窗然後轉移到安裝blast的資料夾,後面以 ...\blast---\bin\ 代表,而所有的檔案都請放在這個資料夾,之後的輸入命令就在命令列視窗進行 (我居然困在這一步兩天阿 囧rz)。

1. 要進行比對前要有一個比對的資料庫,其指令為

makeblastdb  -in [輸入檔案] -dbtype [prot 為蛋白質資料庫; nucl為核酸資料庫] -out [輸出資料庫名]

例如:要從 seq_test.txt 這樣一個蛋白質序列資料庫檔生成一個叫 ooxx 的參考資料庫, 只要開命令列視窗在blast下的...\blast---\bin\ 中輸入

makeblastdb -in seq_test.txt -dbtype prot -out ooxx


執行就會根據seq_test.txt這個序列資料跑出 ooxx.phr ; ooxx.pin ooxx.psq 三個檔案

要注意的是seq_test.txt 的格式要是fasta檔, 我用的格式是由matlab讀檔生成,一行是檔名,下一行沒有分隔符號,全序列結尾也沒有空行,然後也放在...\blast---\bin\資料夾

fasta例:
======================================
>1ppj_I 
MLSVAARSGPFAPVLSATSRGVAGALRPLVQAAVPATSESPVLDLKLSVLCRESLRG
>1ppj_H 
GDPKEEEEEEEELVDPLTTVREQCEQLEKCVKARERLELCDERVSSR
>1ppj_G 
GRQFGHLTRVRHVITYSLSPFEQRAFPHYFSKGIPNVLRRTRACILRVAPPFVAFYL

======================================

然後要做blastp的時候呢.. 可以用下面的指令就可以跑出來了

blastp -query [輸入檔] -db [參考資料庫名] -task [指定] -evalue [Evalue值] -outfmt [Num] -out [輸出檔名]

要注意的是參考資料庫名不需輸入副檔名, 然後outfmt是輸出格式0~9格各有不同格式如下表格

而 task可指定一些特別的功能 如:

例如: 針對一個叫 seq_all_test.txt 的待測序列檔 使用TM_5_test這個資料庫測試, 並指定對短序列最佳化,並設 Evalue>10^10 輸出格式為含表格介紹的表格(7) 存成ooxx.txt檔, 指令可以這樣下

blastp -query seq_all_test.txt -db TM_5_test -task blastp-short -evalue 10000000000 -outfmt 7 -out ooxx.txt

輸出的ooxx.txt檔的形式為

================================================

# BLASTP 2.2.29+

# Query: 1a0t_P

# Database: TM_5_test

# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score

# 23 hits found

1a0t_P 1a0t_P 100.00 413 0 0 1 413 1 413 0.0  867

1a0t_P 1af6_A 27.05 451 257 22 4 413 2 421 2e-019 80.8

1a0t_P 1bcc_A 62.50 8 3 0 273 280 371 378   12 15.5

1a0t_P 1bcc_A 47.37 19 4 1 59 77 19 31   41 13.8

1a0t_P 1bcc_A 21.11 90 48 4 91 157 205 294   48 13.4

                                                       ...

=================================================

                                                   2014 0220 小耿



註:目前跑出來的結果好像會把一些序列切成很多段分開比對, 我還找不到怎麼把兩個序列完整比對的參數 QQ

========
% 2014 0413加註:

如果對於單機版的blast使用不熟練也可以去myblast網站,那邊能夠線上執行blast,
(myblast 網址: http://mybioweb.nhri.org.tw/myblast/ )

只要註冊帳號後就可以用視窗界面上傳參考與目標資料庫即可進行大量序列的比對(當然使用前還是要先了解blast比對的參考與目標資料庫的概念跟格式就是),

myblast的優點是沒有blast 2.2.29版的序列切割問題,數據呈現較於簡潔,但是要注意的是至少在2014 0220時測試的結果該網站使用的還是舊版的blast, 版本為2.2.19。


而當下載大量的比對資料後,要把資料讀進matlab生成E value矩陣運算時,未免一次讀進太大的檔案造成matlab當機,我針對myblast跑出來未設定格式的輸出結果寫了一個逐行讀檔的讀取e value程式,可以在此下載( http://1drv.ms/1kTVQ5D ), 使用時用matlab打開該m檔案然後根據m檔案前面的說明1.設定要運算檔案的正確路徑後(預設資料夾為C:\blast-2.2.29+\bin\),執行即可自動抓檔讀取E value跟序列順序名單,輸出結果為:

序列名單: IDname_meta
E value矩陣: E_value_matrix  

要注意的是此處為未處理的E value 矩陣,此處的E value不合於距離定義下該符合的 Eij =Eji, 為了達到此步驟我設定的修改後的E值為  E'ij=(Eij*Eji)^0.5 。