我第一次聽說BLAST是在大二的結束的那個暑假。我不知道BLAST是幹嘛的,但是當自己在NCBI上輸入一些序列之後,點擊確定,過了幾分鐘跳出一頁常常的結果時,感覺自己突然就很厲害了。
這就是傳說中的好用的工具,你不需要知道他是什麼原理,你只要隨便敲敲,隨便摸索一下,就能學會使用了,而且結果還符合自己的需要。
關於BLAST的教程,其實在連載系列Biostar: 課程11、12 中就介紹過了,內容都其實差不多,其實還比我更全面,但是沒辦法,我的排版好看呀。
Search may take place in nucleotide and/or protein space or translated spaces where nucleotides are translated into proteins.
Searches may implement search 「strategies」: optimizations to a certain task. Different search strategies will return different alignments.
Searches use alignments that rely on scoring matrices
Searches may be customized with many additional parameters. BLAST has many subtle functions that most users never need.
用makeblastdb為BLAST提供資料庫
選擇blast工具,blastn,blastp
運行工具,有必要的還可以對輸出結果進行修飾
第一步:建立檢索所需資料庫BLAST資料庫分為兩類,核酸資料庫和胺基酸資料庫,可以用makeblastbd創建。可以用help參數簡單看下說明。
$ makeblastdb -helpUSAGE makeblastdb [-h] [-help] [-in input_file] [-input_type type] -dbtype molecule_type [-title database_title] [-parse_seqids] [-hash_index] [-mask_data mask_data_files] [-mask_id mask_algo_ids] [-mask_desc mask_algo_descriptions] [-gi_mask] [-gi_mask_name gi_based_mask_names] [-out database_name] [-max_file_sz number_of_bytes] [-logfile File_Name] [-taxid TaxID] [-taxid_map TaxIDMapFile] [-version]-dbtype <String, `nucl', `prot'>
具體以擬南芥基因組作為案例,介紹使用方法:
注: 擬南芥的基因組可以在TAIR上下在,也可在ensemblgenomes下載。後者還可以下載其他植物的基因組
# 下載基因組wget ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gzgzip -d Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz# 構建核酸BLAST資料庫makeblastdb -in Arabidopsis_thaliana.TAIR10.dna.toplevel.fa -dbtype nucl -out TAIR10 -parse_seqids# 下載擬南芥protein數據wget ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz# 構建蛋白BLAST資料庫gzip -dArabidopsis_thaliana.TAIR10.pep.all.fa.gzmakeblastdb -in Arabidopsis_thaliana.TAIR10.pep.all.fa -dbtype prot -out TAIR10 -parse_seqids
如果你從NCBI或者其他渠道下載了格式化過的資料庫,那麼可以用blastdbcmd去檢索blast資料庫,參數很多,常用就如下幾個
db string : string表示資料庫所在路徑
dbtype string,: string在(guess, nucl, prot)中選擇一個
檢索相關參數
輸出格式 -outfmt %f %s %a %g …默認是%f
out 輸出文件
show_blastdb_search_path: blast檢索資料庫路徑
使用案例
# 查看信息blastdbcmd -db TAIR10 -dbtype nucl -info# 所有數據blastdbcmd -db TAIR10 -dbtype nucl -entry all | head# 具體關鍵字,如GI號blastdbcmd -db TAIR10 -dbtype nucl -entry 3 | head
還有其他有意思的參數,可以看幫助文件了解
第二步:選擇blast工具根據不同的需求,比如說你用的序列是胺基酸還是核苷酸,你要查找的數據是核甘酸還是幹計算,選擇合適的blast工具。不同需求的對應關係可以見下圖(來自biostars handbook)
不同工具的應用範圍雖然不同,但是基本參數都是一致的,學會blastn,基本上其他諸如blastp,blastx也都會了。
blastn的使用參數很多blastn [-h],但是比較常用是如下幾個
-db : 資料庫在本地的位置,或者是NCBI上資料庫的類型,如 -db nr
-query: 檢索文件
-query_loc : 指定檢索的位置
-strand: 搜索正義鏈還是反義鏈,還是都要
out : 輸出文件
-remote: 可以用NCBI的遠程資料庫, 一般與 -db nr
-evalue 科學計數法,比如說1e3,定義期望值閾值。E值表明在隨機的情況下,其它序列與目標序列相似度要大於這條顯示的序列的可能性。 與S值有關,S值表示兩序列的同源性,分值越高表明它們之間相似的程度越大
E值總結: 1.E值適合於有一定長度,而且複雜度不能太低的序列。2. 當E值小於10-5時,表明兩序列有較高的同源性,而不是因為計算錯誤。3. 當E值小於10-6時,表時兩序列的同源性非常高,幾乎沒有必要再做確認。
與聯配計分相關參數: -gapopen,打開gap的代價;-gapextend, gap延伸的代價;-penalty:核算錯配的懲罰; -reward, 核酸正確匹配的獎勵;
結果過濾:-perc_identity, 根據相似度
注 BLAST還提供一個task參數,感覺很有用的樣子,好像會針對任務進行優化速度。
第三步:運行blast,調整輸出格式。我隨機找了一段序列進行檢索
echo '>test' > query.fa echo TGAAAGCAAGAAGAGCGTTTGGTGGTTTCTTAACAAATCATTGCAACTCCACAAGGCGCCTGTAATAGACAGCTTGTGCATGGAACTTGGTCCACAGTGCCCTACCACTGATGATGTTGATATCGGAAAGTGGGTTGCAAAAGCTGTTGATTGTTTGGTGATGACGCTAACAATCAAGCTCCTCTGGT >> query.fa
用的是blastn 檢索核酸資料庫。最簡單的用法就是提供資料庫所在位置和需要檢索的序列文件
blastn -db BLAST/TAIR10 -query query.fa# 還可以指定檢索序列的位置blastn -db BLAST/TAIR10 -query query.fa -query_loc 20-100# 或者使用遠程資料庫blastn -db nr -remote -query query.fa blastn -db nt -remote -query query.fa
以上是默認輸出,blast的-outfmt選項提供個性化的選擇。一共有18個選擇,默認是0。
0 = Pairwise,
1 = Query-anchored showing identities,
2 = Query-anchored no identities,
3 = Flat query-anchored showing identities,
4 = Flat query-anchored no identities,
5 = BLAST XML,
6 = Tabular,
7 = Tabular with comment lines,
8 = Seqalign (Text ASN.1),
9 = Seqalign (Binary ASN.1),
10 = Comma-separated values,
11 = BLAST archive (ASN.1),
12 = Seqalign (JSON),
13 = Multiple-file BLAST JSON,
14 = Multiple-file BLAST XML2,
15 = Single-file BLAST JSON,
16 = Single-file BLAST XML2,
17 = Sequence Alignment/Map (SAM),
18 = Organism Report
其中6,7,10,17可以自定輸出格式。默認是
qaccver saccver pident length mismatch gapopen qstart qend sstart send evalue bitscore
簡寫含義qaccver查詢的AC版本(與此類似的還有qseqid,qgi,qacc,與序列命名有關)saccver目標的AC版本(於此類似的還有sseqid,sallseqid,sgi,sacc,sallacc,也是序列命名相關)pident完全匹配百分比 (響應的nident則是匹配數)length聯配長度(另外slen表示查詢序列總長度,qlen表示目標序列總長度)mismatch錯配數目gapopengap的數目qstart查詢序列起始qstart查詢序列結束sstart目標序列起始send目標序列結束evalue期望值bitscoreBit得分score原始得分AC: accession
以格式7為實例進行輸出,並且對在線資料庫進行檢索
blastn -task blastn -remote -db nr -query query.fa -outfmt 7 -out query.txthead -n 15 query.txt
PS: 這是一篇欠了我師兄好久的教程,現在終於補上了
更多更全的內容點擊原文連結,看看biostars handbook系列。