序列比對工具 | BLAST、 BLAT、 diamond

基本概念

相似性(similarity)

  • 一種很直接的數量關系,比如部分相同或相似的百分比或其他一些合適的度量
  • 如:A序列和B序列的相似性是80%

同源性(homology)

  • 從一些數據中推斷出的兩個基因或者蛋白序列具有共同祖先的結論,屬于質的判斷
  • 可以說A序列和B序列是同源序列,但不能說同源性80%

常用工具

  • BLAST
  • BLAT
    序列比對的常用工具:BLAST,但是其運行速度慢的令人捉急。

一、BLAST(Basic Local Alignment Search Tool,局部相似性基本查詢工具)

BLAST(Basic Local Alignment Search Tool)是一套在蛋白質數據庫或DNA數據庫中進行相似性比較的分析工具。BLAST程序能迅速與公開數據庫進行相似性序列比較。BLAST結果中的得分是對一種對相似性的統(tǒng)計說明。

https://www.biomart.cn/experiment/599/608/19912_0.htm
-F Filter query sequence (DUST with blastn, SEG with others) [String] default = T
查詢序列過濾,將那些 給出影響比對結果的低復雜度區(qū)域過濾掉。用blastn進行查詢的序列用DUST程序過濾,其他的用SEG過濾 。對DUST和SEG的詳細情況,用戶可以自己查詢資料。
使用此參數 -F F 即可獲得完整的比對

在對核苷酸或蛋白質序列數據庫進行Blast搜索之前,必須要對所使用的序列數據庫進行formatdb, 即對序列數 據庫進行格式化,這是所有使用BLAST所必須的一步。

1.格式化序列數據庫— —formatdb

formatdb 簡單介紹:formatdb處理的都是格式為 ASN.1和 FASTA,而且不論是核苷酸序列數據庫,還是蛋白質序列數據庫;不論是使用Blastall ,還是Blastpgp,Mega Blast應用程序,這一步都是不可少的。

主要參數的說明
-i 輸入需要格式化的源數據庫名稱 Optional
-p 文件類型,是核苷酸序列數據庫,還是蛋白質序列數據庫 T:protein , F:nucleotide
-n 重命名數據庫文件的名稱 ;
-t 數據庫的標題【可選】;
-l 日志文件名,formatdb.log
-o 解析選項. T - True: 解析序列標識并且建立目錄,F - False: 與上相反,[T/F] Optional default = F

示例:

formatdb -i uniref100.fasta -n uniref100 -t uniref100 -l uniref100.log -p T
formatdb -i uniref90.fasta -n uniref90 -t uniref90 -l uniref90.log -p T
formatdb -i uniref50.fasta -n uniref50 -t uniref50 -l uniref50.log -p T

2.運行blastall

blastall -i query.fasta -d uniref50 -o blast.out -p blastn -F F -e 1e-5 -m 8

1.-e參數
用來過濾比對較差的結果的,用"-e"參數指定一個實數,blast 會過濾掉期望值大于這個數的比對結果。這樣不但簡化了結果,還縮短了運行時間和結果占用的空間。

2. -p參數
-p blastp:蛋白序列與蛋白庫做比對。
-p blastx:核酸序列對蛋白庫的比對。
-p blastn:核酸序列對核酸庫的比對。
-p tblastn:蛋白序列對核酸庫的比對。

3.-F 參數
-F(T/F)參數是用來屏蔽簡單重復和低復雜度序列的。如果選“T”,程序在比對過程中會屏蔽query中的簡單重復和低復雜度序列;選"F"則不會屏蔽。缺省值為"T"。

4. -m

m8格式如下圖,12列:

1、Query id:查詢序列ID標識
2、Subject id:比對上的目標序列ID標識
3、% identity:序列比對的一致性百分比
4、alignment length:符合比對的比對區(qū)域的長度
5、mismatches:比對區(qū)域的錯配數
6、gap openings:比對區(qū)域的gap數目
7、q. start:比對區(qū)域在查詢序列(Query id)上的起始位點
8、q. end:比對區(qū)域在查詢序列(Query id)上的終止位點
9、s. start:比對區(qū)域在目標序列(Subject id)上的起始位點
10、s. end:比對區(qū)域在目標序列(Subject id)上的終止位點
11、e-value:比對結果的期望值,解釋是大概多少次隨機比對才能出現一次這個score,Evalue越小,表明這種情況從概率上越不可能發(fā)生,那么這個比對的可靠性越高。
12、bit score:比對結果的bit score值
一般情況我們看第3、11、12兩列,e值越小越可靠。
blastall 老版本對應的參數是 -m 8
blast+對應的參數是-outfmt 6

m6 格式

查找了一下,列名分別為:

  1. qseqid query (e.g., unknown gene) sequence id
  2. sseqid subject (e.g., reference genome) sequence id
  3. pident percentage of identical matches
  4. length alignment length (sequence overlap)
  5. mismatch number of mismatches
  6. gapopen number of gap openings
  7. qstart start of alignment in query
  8. qend end of alignment in query
  9. sstart start of alignment in subject
  10. send end of alignment in subject
  11. evalue expect value
  12. bitscore bit score

二、BLAT

https://blog.csdn.net/alnx37271/article/details/101358955?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2aggregatepagefirst_rank_ecpm_v1~rank_v31_ecpm-2-101358955.pc_agg_new_rank&utm_term=blastall+%E4%BD%BF%E7%94%A8&spm=1000.2123.3001.4430

1. 簡介

Blat,全稱 The BLAST-Like Alignment Tool,可以稱為"類BLAST 比對工具"。

  • 對于DNA序列,BLAT是用來設計尋找95%及以上相似至少40個堿基的序列。
  • 對于蛋白序列,BLAT是用來設計尋找80%及以上相似至少20個氨基酸的序列。

特點:

  • 速度快(直接把數據庫索引讀入內存,無需訪問硬盤);
  • 共線性輸出結果簡單易讀;
  • 對于比較小的序列和大基因組的比對,BLAT是首選;

對于比較小的序列(如 cDNA 等)對大基因組的blat與blast比較比對,blat 無疑是首選。Blat 把相關的呈共線性的比對結果連接成為更大的 比對結果,從中也可以很容易的找到 exons 和 introns。因此,在相近物種的基因同源性分析和EST 分析中,blat 得到了廣泛的應用。

Blat 的比對速度之所以能比Blast快幾百倍,是因為此兩者之間的比對機制有著本質的差別。

  • Blast是將查詢序列索引化,然后線性搜索龐大的目標數據庫,期間頻繁地訪問硬盤數據,時間和空間上的數據相關性較??;
  • Blast 則將龐大的目標數據庫索引化,然后線性搜索查詢序列,這種搜索方式在時間和空間上的數據相關性比較大。

Blat將數據庫索引一次性讀入內存,可以反復地高速調用,無需訪問硬盤,占用的系統(tǒng)資源很少。只要索引建立,查詢序列的量越大,Blat的優(yōu)勢就越明顯。

2. 原理

首先 blat 將參考序列拆分成tiles/kmers,其拆分的方式取決于兩個參數:

  • -tileSize 決定tiles/kmers的大小,一般設定范圍是:8-12,預設DNA為11,蛋白質為5;
  • -stepSize 決定tiles/kmers移動的步長;
image

3. 軟件下載與安裝

簡單方便,只需無腦輸入以下命令:

wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/blat/blat
chmod u+x blat

4. 軟件運行(比對)

命令如下:

用法:
blat database query [-ooc=11.ooc] output.psl

blat nt test.fa test.out -out=blast8

說明:

  • blat有很多參數可以選擇,但大部分時候我們按照默認的即可。
  • blat的輸入文件和待查詢數據庫的格式,可以是fa/nib/2bit三種格式中的任意一種。運行時十分簡單,不需要進行建庫就可以直接比對。
  • 軟件運行時,依次輸入軟件名(軟件執(zhí)行路徑),待比對的數據庫,待搜索的序列,輸出結果。順序前后不能顛倒。這樣就可以開始比對了
  • 程序開始運行后,會在讀完database 中的所有subject 序列時在屏幕輸出database的統(tǒng)計結果。

以下是一些常用的參數:

  • -noHead: 不輸出表頭信息,有助于結果文件的后續(xù)處理
  • -out: 指定輸出的文件格式,此處使用的是blast的m8格式
  • -t: 指定數據庫的類型,dna/prot/dnax
  • -q: 指定序列類型,dna/rna/prot/dnax/rna

blat詳細參數

用法:blat database query [-ooc=11.ooc] output.psl
database      輸入文件必須是其中一種類型:a .fa , .nib or .2bit file
query         輸入文件必須是其中一種類型:a .fa , .nib or .2bit file
output.psl    輸出文件
-t=type       數據庫類型,可選項: dna/prot/dnax
-q=type       查詢序列的類型,可選項:dna/prot/dnax/rnax
-prot         等同于 -t=prot -q=prot
-ooc=N.ooc Use overused tile file N.ooc.  N should correspond to the tileSize
-tileSize=N   設定tiles/kmers的大小
-stepSize=N   設定tiles/kmers在比對時移動的步長,即兩個相鄰tiles/kmers之間的距離,預設值是tileSize
-oneOff=N     如果設定為 1 ,則表示在比對到tile上允許有一個錯配堿基(mismatch),預設值是0
-minMatch=N   設定至少匹配的tile的個數,一般設置值的范圍是2-4,通常核苷酸的預設值為2,蛋白質的預設值為1
-minScore=N   設定最小分值。 由于indel通常會對序列的功能產生影響,所以空位在比對過程中總是對應于一個負分,也就是所謂的空位罰分(Gap penalty)。根據打分機制,這個分值等于匹配堿基分值減去替換分值(mismatch)和空位罰分。預設值為30
-minIdentity=N 設置序列相似度(sequence identity)最小百分比。通常核苷酸(nucleotide searches)預設值為90,蛋白質和翻譯蛋白(protein or translated protein searches)預設值為25
-maxGap=N     在一定長度序列中,設定兩個tiles/kmers之間的允許最大的空位(gap)大小。通常設定范圍是0-3,預設值為2,且僅在minMatch > 1時搭配使用
-noHead       抑制.psl頭文件的輸出,內容全部均是以制表符為分隔符的文件
-makeOoc=N.ooc Make overused tile file. Target needs to be complete genome.
-repMatch=N   在一段序列被標記為overused之前,設定允許tiles/kmers重復次數。如果超過設定值,該tiles/kmers將會被標記為overused。通常當tileSize設定為12時,repMatch則設定為256;當tileSize設定為11時,repMatch則設定為1024;當tileSize設定為10時,repMatch則設定為4096。
-mask=type Mask out repeats.  Alignments won't be started in masked region but may extend through it in nucleotide searches.  Masked areas are ignored entirely in protein or translated searches. Types are
            lower - mask out lower cased sequence
            upper - mask out upper cased sequence
            out   - mask according to database.out RepeatMasker .out file
            file.out - mask database according to RepeatMasker file.out
-qMask=type Mask out repeats in query sequence. 類型選擇與參數-mask相同
-repeats=type 類型選擇與參數-mask相同。無論如何重復堿基不會被掩蓋(masked),但是在匹配重復區(qū)域時將會在psl輸出文件中會單獨展示其匹配結果,即與其他區(qū)域的匹配結果是分開的。
-minRepDivergence=NN - minimum percent divergence of repeats to allow them to be unmasked.  Default is 15.  Only relevant for masking using RepeatMasker .out files.
-dots=N     每N個序列就輸出一個點,用于展示程序運行的進度
-trimT      剪切首部的poly-T
-noTrimA    不剪切尾部的poly-A
-trimHardA  從psl輸出文件中的qSize和alignments中移除poly-A尾巴
-fastMap    快速的DNA/DNA remapping,要求查詢序列長度不超過5000、高相似度和不進行內含子的比對
-out=type  輸出文件格式,格式如下:
                   psl - Default.  Tab separated format, no sequence
                   pslx - Tab separated format with sequence
                   axt - blastz-associated axt format
                   maf - multiz-associated maf format
                   sim4 - similar to sim4 format
                   wublast - similar to wublast format
                   blast - similar to NCBI blast format 
                   blast8- NCBI blast tabular format
                   blast9 - NCBI blast tabular format with comments
-fine           對于高質量的mRNAs搜索small initial和terminal exons更為嚴苛。此選項不推薦應用于ESTs  
                For high quality mRNAs look harder for small initial and terminal exons.
-maxIntron=N    設定內含子最大的序列長度. Default is 750000
-extendThroughN 允許序列的比對可以從大段N區(qū)域延伸

使用案例

# 處理單個job
blat chr11.fa human/test.fa test.psl                # 輸出不含序列
blat chr11.fa human/test.fa -out=pslx test.pslx     # 輸出含序列
blat chr11.fa human/test.fa -out=blast test.blast   # 輸出格式同NABI的blast格式


# 并行處理多個jobs
time parallel blat chr{}.fa human/human.fa test_{}.psl ::: {1..22} X Y M

5. 問題

1. 內存溢出
值得關注的是,因為blat的算法原理,它是將整個帶查詢數據庫全部讀入內存進行比對。因此如果你的服務器內存不大,不建議使用blat進行nt/nr庫的比對。

下面給出一個簡單的計算方法,用于評估自己的服務器是否可以順溜的run blat。對于帶查詢基因組,平均每個堿基,需要2bytes的內存。wiki給出的人類基因組大小為3100Mb,因此我們大概需要6.2G的內存才可以順利的對人類基因組進行blat查詢。

年少無知的我,用128G內存的服務器跑nt數據庫,理所當然的把我們課題組的服務器跑崩了~

https://www.biostars.org/p/9479310/#9479314

參考鏈接:Using blat

參考鏈接:Blat-The BLAST-Like Alignment Tool (詳細的使用教程)

6. GNU Parallel

安裝編譯

wget ftp://ftp.gnu.org/gnu/parallel/parallel-20170822.tar.bz2
tar -jxvf parallel-20170822.tar.bz2 
cd parallel-20170822/
cat README 
./configure && make && sudo make install

使用

7. 網絡版

鏈接 :http://genome.ucsc.edu/cgi-bin/hgBlat

操作方法

  • Genome:選擇物種,比如人
  • Assembly:版本號
  • Query type:用于查詢的序列類型(DNA/蛋白質)
  • Sort output:結果排序方式
  • Output type:輸出格式
    • hyperlink:指向結果的超鏈接,便于可視化
    • psl:制表符分隔的表格,便于數據處理

查詢結果(hyperlink)

ACTIONS QUERY SCORE START END QSIZE IDENTITY CHROM STRAND START END SPAN
browser details CRP_HUMAN 671 1 224 224 100.0% chr1 +- 159713528 159714485 958
browser details CRP_HUMAN 105 119 183 224 77.0% chr1 +- 159705131 159705325 195
browser details CRP_HUMAN 54 117 188 224 62.5% chr1 ++ 159276797 159277012 216

詳情

點擊Browser可以進入詳情界面

image

查詢結果(psl)

match mismatch rep. match N's Q gap count Q gap bases T gap count T gap bases strand Q name Q size Q start Q end T name T size T start T end block count blockSizes qStarts tStarts
224 0 0 0 0 0 1 286 +- CRP_HUMAN 224 0 224 chr1 248956422 159713527 159714485 2 19,205, 0,19, 89241937,89242280,
50 15 0 0 0 0 0 0 +- CRP_HUMAN 224 118 183 chr1 248956422 159705130 159705325 1 65, 118, 89251097,
45 27 0 0 0 0 0 0 ++ CRP_HUMAN 224 116 188 chr1 248956422 159276796 159277012 1 72, 116, 159276796,

三、diamond

diamond作為一個和BLAST具有相似功能的軟件,具有以下特點:

  • 比BLAST快500到20,000倍
  • 長序列的移框聯配分析(frameshift alignment)
  • 資源消耗小,普通臺式機和筆記本都能運行
  • 輸出格式多樣

軟件的安裝很簡單,以下是具體命令:

wget http://github.com/bbuchfink/diamond/releases/download/v0.9.25/diamond-linux64.tar.gz
tar -zxzf diamond-linux64.tar.gz

diamond 的功能是將蛋白序列或者其翻譯后的核苷酸和蛋白質數據庫進行比對,與blast相比功能單一,但也讓它的使用格外的簡單。

不推薦將核苷酸序列與蛋白質數據庫進行比對,因為可能有許多序列是比對到非編碼序列上的,利用蛋白質數據庫進行比對,將導致假陰性過高。

下載nr數據庫并建庫

具體命令如下:

wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz
gunzip nr.gz

diamond makedb --in nr --db nr

  • --in: 后面跟蛋白質數據庫
  • --db: 指定生成的diamond數據庫名稱

比對搜索

相當簡單,只有兩個子命令,blastx和blastp,前者比對DNA序列,后者比對蛋白

diamond blastx --db nr -q reads.fna -o dna_matches_fmt6.txt
diamond blastp --db nr -q reads.faa -o protein_matches_fmt6.txt

參數簡單介紹:

  • -q: 輸入的待檢索序列
  • -o:指定輸出文件,默認以 --outfmt 6 格式進行輸出
  • --db: 后面跟著我們建立好的diamond 蛋白數據庫

參考鏈接
http://m.itdecent.cn/p/f6f868010949
http://m.itdecent.cn/p/7d536d8da3fb

最后編輯于
?著作權歸作者所有,轉載或內容合作請聯系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成,瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明:文章內容(如有圖片或視頻亦包括在內)由作者上傳并發(fā)布,文章內容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務。

相關閱讀更多精彩內容

友情鏈接更多精彩內容