相信做生物的小伙伴都接觸過序列比對(duì)，不論是想知道測(cè)序得到的序列還是你手中的由ATCG組成的一條序列想是什么，都需要進(jìn)行進(jìn)行序列的比對(duì)，才能確定它的身份。 序列比對(duì)是一種用于比較和分析序列之間的相似性和差異性的方法。它在生物學(xué)和計(jì)算機(jī)科學(xué)領(lǐng)域被廣泛應(yīng)用于DNA、RNA和蛋白質(zhì)序列的比較分析。

序列比對(duì)簡(jiǎn)介

序列比對(duì)的原理基于兩個(gè)基本假設(shè)：保守性和演化性。保守性假設(shè)認(rèn)為在不同物種或個(gè)體之間的相似序列可能在進(jìn)化過程中被保留下來，而相似的序列可能具有相似的功能。演化性假設(shè)認(rèn)為序列之間會(huì)發(fā)生變異和突變，通過比較序列之間的差異可以推斷它們之間的共同祖先和演化關(guān)系。根據(jù)比對(duì)的目的和序列的長(zhǎng)度，選擇適當(dāng)?shù)谋葘?duì)算法。常用的比對(duì)算法包括全局比對(duì)算法（如Needleman-Wunsch算法）和局部比對(duì)算法（如Smith-Waterman算法）。 需要注意的是，序列比對(duì)是一個(gè)計(jì)算復(fù)雜度較高的任務(wù)，特別是當(dāng)處理較長(zhǎng)的序列或大規(guī)模的多序列比對(duì)時(shí)。因此，研究人員通常利用高性能計(jì)算和優(yōu)化算法來加快比對(duì)的速度和提高比對(duì)的準(zhǔn)確性。

序列比對(duì)定義

簡(jiǎn)單的進(jìn)行匹配，可以看到僅有兩個(gè)堿基匹配

如果引入一個(gè)空格或空格（gap），匹配上的狀況就可以明顯改善

因此，小果根據(jù)查閱文獻(xiàn)認(rèn)為以下說法是序列比對(duì)的最好解釋。 根據(jù)特定的計(jì)分規(guī)則，通過一定的算法對(duì)兩條或者多條DNA或者蛋白序列進(jìn)行比較，找出他們之間的最優(yōu)匹配或者最大相似度匹配。序列匹配的意義在于，通過序列比對(duì)可以獲得一個(gè)序列相似性度比值，通過這個(gè)值可以獲得兩條序列相似度或者同源性的統(tǒng)計(jì)學(xué)判定結(jié)果。

動(dòng)態(tài)規(guī)劃

想要了解blastn或者序列比對(duì)的底層邏輯，我們需要把動(dòng)態(tài)規(guī)劃弄明白，講清楚。

打分矩陣

首先，需要了解什么是打分矩陣。 對(duì)堿基匹配情況進(jìn)行一個(gè)標(biāo)準(zhǔn)化的矩陣，可以根據(jù) 1、替代矩陣： DNA替代矩陣：用于評(píng)估DNA序列中不同堿基之間的替代得分。常見的DNA替代矩陣包括簡(jiǎn)單的等分矩陣（例如，A和T之間得分為1，A和G之間得分為-1）以及更復(fù)雜的模型，如PAM（Point Accepted Mutation）和BLOSUM（Blocks Substitution Matrix）矩陣。 蛋白質(zhì)替代矩陣：用于評(píng)估蛋白質(zhì)序列中不同氨基酸之間的替代得分。常見的蛋白質(zhì)替代矩陣包括BLOSUM系列和PAM系列矩陣。這些矩陣根據(jù)蛋白質(zhì)家族的多序列比對(duì)結(jié)果建立，考慮了氨基酸的進(jìn)化關(guān)系和替代頻率。 2、缺失和插入罰分： Gap罰分矩陣：用于懲罰序列比對(duì)中的缺失（Gap）和插入。引入Gap會(huì)導(dǎo)致序列不匹配和不對(duì)齊，為了懲罰Gap的引入，通常在打分矩陣中引入了開放(Gap Opening Penalty)和擴(kuò)展(Gap Extension Penalty)處罰得分。這些得分會(huì)隨著Gap的長(zhǎng)度逐漸增加。 3、直接得分矩陣： 匹配-不匹配得分：在某些情況下，打分矩陣中的得分可以根據(jù)堿基或氨基酸之間的匹配或不匹配來定義。例如在DNA序列比對(duì)中，匹配得分可以設(shè)為正值1。 ?這里對(duì)各種矩陣不理解可以直接跳過，直接了解下邊例如中的兩個(gè)矩陣。 例如：下邊兩個(gè)矩陣，第一個(gè)只考慮匹配上的得分為1，第二個(gè)只考慮了考慮顛換和置換的情況（序列中嘌呤被嘧啶替換或反之，稱為顛換（transversion），如A→C、G→T；嘌呤或嘧啶自己互換稱為置換（transition），如A→G、C→T?？蓪?duì)核苷酸的顛換與置換給予不同的評(píng)分。）在序列比對(duì)中打分矩陣可以根據(jù)需求進(jìn)行自定義～

比對(duì)算法

全局比對(duì)：給定兩個(gè)序列A＝a1a2…an 和B＝b1b2…bm，S（i，j）表示兩個(gè)序列所有比對(duì)的最好分?jǐn)?shù)。在設(shè)定初值后可以用以下遞歸關(guān)系計(jì)算：

局部比對(duì)的具體算法與全局比對(duì)幾乎一樣，給定兩個(gè)序列A＝a1a2…an 和B＝b1b2…bm，S（i，j）表示兩個(gè)序列任何比對(duì)的最好分?jǐn)?shù)。在設(shè)定初值后可以用以下遞歸關(guān)系計(jì)算：

上述的算法中，不論是全局比對(duì)還是局部比對(duì)都是來尋找最大值的過程，我們需要得到Mij的分?jǐn)?shù)，是取了三個(gè)方向中的最大值。

例如：現(xiàn)在還是用開頭時(shí)候我們用的序列，我們現(xiàn)在用動(dòng)態(tài)規(guī)劃的方式來得到下圖中的匹配方式：

規(guī)定的打分矩陣為：匹配得分5分；空位罰分-1分；錯(cuò)配 -3分 小果建議可以在excel嘗試自己算出以下小果填入的數(shù)字。 這里注意，是mismatch還是gap，如果你是從對(duì)角線移動(dòng)的為錯(cuò)配，從上邊或者右邊的格子移動(dòng)的即為gap，這兩種情況的罰分是不一樣的。

得到上述格子后，查閱到得分最高的格子，追溯到是怎么來的序列即可。

blast比對(duì)原理

1、編譯一個(gè)由查詢序列生成的長(zhǎng)度固定的字段編譯列表“words” 2、在數(shù)據(jù)庫中掃描獲得與編譯列表中的字段匹配的序列記錄； 3、以編譯列表中的字段對(duì)為中心向兩端延伸以尋找超過閾值分?jǐn)?shù)S的高分值片段HSP。 4、對(duì)于每部分HSP使用Smith-Watermans算法進(jìn)行局部性比對(duì)，為每部分打分，作為最終結(jié)果。 下圖，是blast的文章中的原理 E值衡量了在隨機(jī)情況下，數(shù)據(jù)庫存在的比當(dāng)前匹配分?jǐn)?shù)更好的比對(duì)的數(shù)目，因此可以用該值作為指標(biāo)評(píng)價(jià)HSP比對(duì)序列的可信度，E值計(jì)算公式：E＝Kmne－λS

了解了blast的原理后，下面介紹一些blast的應(yīng)用和參數(shù)。

blast比對(duì)實(shí)際操作

比對(duì)：

blastn -query input.fasta -db database -out output.txt -evalue 0.001 -outfmt 6 blastn -query input.fasta?-db database -out output.txt -outfmt " 7 qaccver saccver pident length mismatch qlen qstart qend sstart send qseq evalue bitscore" -num_threads 10 -word_size#也可以嘗試7的格式。

在這個(gè)示例中，blastn是BLAST程序的名稱，input.fasta是輸入文件，database是BLAST數(shù)據(jù)庫，output.txt是輸出結(jié)果文件，-evalue參數(shù)指定過濾結(jié)果的期望值閾值，-outfmt指定輸出格式為表格格式（6） ?

構(gòu)建數(shù)據(jù)庫：

makeblastdb -dbtype nucl??-in ${name}.fa?-out ${name} #完成后會(huì)顯示以下提示

#會(huì)生成以下的幾個(gè)文件

參數(shù)解讀：

1、E-value（期望值）： E-value是一個(gè)衡量比對(duì)結(jié)果的期望隨機(jī)發(fā)生次數(shù)的指標(biāo)。它表示在進(jìn)行大量隨機(jī)比對(duì)的情況下，期望得到與查詢序列的比對(duì)相似性至少一樣好（或更好）的結(jié)果的次數(shù)。 較小的E-value表示比對(duì)結(jié)果的相似性更顯著和顯著。通常，E-value小于某個(gè)預(yù)設(shè)的閾值（例如0.001或0.05）被認(rèn)為是顯著的比對(duì)結(jié)果。在BLAST輸出中，E-value一般以科學(xué)計(jì)數(shù)法的形式給出，例如1.0e-10，表示10的負(fù)10次方（即0.0000000001）。 2、Score是基于比對(duì)中的匹配和不匹配進(jìn)行計(jì)算的一個(gè)值，用于衡量比對(duì)的質(zhì)量和相似性。具體的得分計(jì)算方式取決于所使用的打分矩陣。一般來說，較高的Score表示比對(duì)結(jié)果具有更好的相似性和匹配程度。通常情況下，Score越高，比對(duì)的質(zhì)量越好。在BLAST輸出中，Score是以具體數(shù)值的形式給出。 3、“Identity”（相同性）是輸出結(jié)果中的一個(gè)重要指標(biāo)，用于描述比對(duì)序列的相似性程度。它表示在比對(duì)過程中，兩個(gè)序列中存在相同的堿基或氨基酸的百分比。相同性（Identity）：在比對(duì)結(jié)果中，相同性通常以百分比的形式表示，反映了兩個(gè)序列中相同位置的堿基或氨基酸的百分比。較高的相同性表示比對(duì)結(jié)果具有更高的相似性。相同性的計(jì)算：相同性的計(jì)算是通過比對(duì)過程中的匹配（Matches）和不匹配（Mismatches）數(shù)量來得出的。計(jì)算公式：Identity（相同性）= (Matches / Alignment Length) × 100% 其中，Matches表示比對(duì)中完全匹配的位置數(shù)，Alignment Length表示比對(duì)的序列長(zhǎng)度。 例如，假設(shè)有一個(gè)比對(duì)結(jié)果，序列的長(zhǎng)度為100，其中有80個(gè)位置是完全匹配的，那么相同性為： Identity = (80 / 100) × 100% = 80% 請(qǐng)注意，相同性只考慮了完全匹配的位置，不考慮缺失（Gaps）和插入（Insertions）。缺失和插入會(huì)影響相似性和比對(duì)質(zhì)量的其他方面，例如Gap Coverage和Gap Length。 Identity、E-value、Score、Coverage等。這些指標(biāo)的綜合分析可以提供更全面的了解和解釋比對(duì)結(jié)果的相似性和可靠性。 小果附上blasn命令的詳細(xì)解釋： -query ：指定輸入文件，其中包含待比對(duì)的查詢序列。 -db ：指定比對(duì)所使用的數(shù)據(jù)庫。 -out ：指定輸出文件的名稱或路徑。 -evalue ：設(shè)置過濾結(jié)果的期望值（E-value）閾值。默認(rèn)值為10。 -outfmt ：指定輸出結(jié)果的格式。常用的格式參數(shù)包括： 0：Pairwise（默認(rèn)格式） 5：XML格式 6：Tabular格式 7：Tabular格式（注釋） 更多格式參數(shù)可以通過blastn --help查看。 除了上述常見的參數(shù)外，還有一些其他參數(shù)可以進(jìn)一步調(diào)整BLAST比對(duì)的行為： -num_threads ：設(shè)置使用的線程數(shù)（并行處理）。默認(rèn)值為1。 -task ：指定BLAST任務(wù)類型，例如blastn、blastn-short、blastn-fast等。不同的任務(wù)類型針對(duì)不同的查詢和目標(biāo)序列類型進(jìn)行了優(yōu)化。 -word_size ：設(shè)置匹配單詞的長(zhǎng)度。默認(rèn)值為11。 -gapopen ：設(shè)置開啟一個(gè)gap（間隙）的懲罰分?jǐn)?shù)。默認(rèn)值為5。 -gapextend ：設(shè)置擴(kuò)展一個(gè)gap（間隙）的懲罰分?jǐn)?shù)。默認(rèn)值為2。 -max_target_seqs ：設(shè)置從數(shù)據(jù)庫中返回的最大比對(duì)序列數(shù)。默認(rèn)值為500。 -perc_identity ：設(shè)置比對(duì)結(jié)果中的最低相似度百分比。默認(rèn)值為0。 好了今天的blast就講到這里，歡迎大家有問題與小果一起討論哦～