生物計(jì)算技術(shù)多重序列比對分析

生物計(jì)算技術(shù)多重序列比對分析第1頁/共98頁8條免疫球蛋白序列片段的多重比對：半光氨酸色氨酸疏水殘基保守區(qū)域SP得分Biocomputingtechnology—Multiplesequencealignment第2頁/共98頁Biocomputingtechnology—Multiplesequencealignment

通過序列的多重比對，可以得到一個(gè)序列家族的序列特征。當(dāng)給定一個(gè)新序列時(shí)，根據(jù)序列特征，可以判斷這個(gè)序列是否屬于該家族。

對于多序列比對，現(xiàn)有的大多數(shù)算法都基于漸進(jìn)比對的思想，在序列兩兩比對的基礎(chǔ)上逐步優(yōu)化多序列比對的結(jié)果。進(jìn)行多序列比對后，可以對比對結(jié)果進(jìn)行進(jìn)一步處理，例如構(gòu)建序列的特征模式，將序列聚類、構(gòu)建分子進(jìn)化樹等。第3頁/共98頁4.2多重序列比對算法原理4.2.1SP模型4.2.2多重比對的動(dòng)態(tài)規(guī)劃算法4.2.3優(yōu)化算法4.2.4星型比對4.2.5樹形比對4.2.6CLUSTALW算法4.2.7隱馬爾可夫模型Biocomputingtechnology—Multiplesequencealignment第4頁/共98頁4.2.1SP模型

(Sum-of-Pairs)

逐對加和函數(shù)作用：評價(jià)多重序列比對的結(jié)果SP計(jì)算的兩種方法：Biocomputingtechnology—Multiplesequencealignment第5頁/共98頁方法1：先計(jì)算多重比對結(jié)果的每一列字符的得分，

然后求整體多重比對得分Biocomputingtechnology—Multiplesequencealignment假設(shè):得分函數(shù)(代價(jià)函數(shù))具有加和性，即多重比對的得分是各列得分總和。思路:如何給比對的每一列打分，然后將各列的和加起來，成為一個(gè)總得分。每一列的處理方式:

尋找一個(gè)具有k個(gè)變量的打分函數(shù)，每一個(gè)變量或者是一個(gè)來自特定字母表中的字符，或者是一個(gè)空位。

k是參與多重比對的序列的個(gè)數(shù)。第6頁/共98頁Biocomputingtechnology—Multiplesequencealignment顯式函數(shù)應(yīng)滿足如下條件：函數(shù)形式簡單，具有統(tǒng)一的形式，不隨序列的個(gè)數(shù)而發(fā)生形式的變化。2.根據(jù)得分函數(shù)的意義，函數(shù)值應(yīng)獨(dú)立于各參數(shù)的順序，即與待比較的序列先后次序無關(guān)。3.對相同的或相似字符的比對，獎(jiǎng)勵(lì)的得分值高，而對于不相關(guān)的字符比對或空白，則進(jìn)行懲罰（得分為負(fù)值）。滿足上述條件的一個(gè)函數(shù)就是常用的逐對加和函數(shù)，SP函數(shù)。第7頁/共98頁方法1：先計(jì)算多重比對結(jié)果的每一列字符的得分，

然后求整體多重比對得分其中，c1,c2,…,ck是一列中的k個(gè)字符，

p是關(guān)于一對字符相似性的打分函數(shù)。

SP_score(c1,c2,…,ck)是多重序列比對中某一列的得分.Biocomputingtechnology—Multiplesequencealignment第8頁/共98頁例：圖4.1多重比對的倒數(shù)第3列的SP得分：打分函數(shù)：

P（a,a）=0P（a,b）=-1（a≠b）

P（a,-）=P（-,b）=-1

P（-,-）=0逐對計(jì)算p(1,2)，p(1,3)，...，p(1,8)，p(2,3)，p(2,4)，...p(2,8)...，p(7,8)的所有得分：（-7-6-5-4-3-2-1）+2=-26然后將一個(gè)多重比對所有列的得分全部加起來，其和即為該多重比對的得分。將所有多重比對的得分計(jì)算出來進(jìn)行比較，得分最高的，應(yīng)該是最好的。

Biocomputingtechnology—Multiplesequencealignment第9頁/共98頁多重比對在兩條特定序列上的投影Biocomputingtechnology—Multiplesequencealignment第10頁/共98頁方法2：先計(jì)算多重序列結(jié)果的序列兩兩比對得分，

然后計(jì)算整體多重比對得分。是一個(gè)多重比對ij是由推演出來的序列si和sj的兩兩比對方法1和方法2等價(jià)的條件：P(-,-)=0Biocomputingtechnology—Multiplesequencealignment第11頁/共98頁4.2.2多重比對的動(dòng)態(tài)規(guī)劃算法

多重序列比對的最終目標(biāo)是通過處理得到一個(gè)得分最高（或代價(jià)最小）的序列對比排列，從而分析各序列之間的相似性和差異。Biocomputingtechnology—Multiplesequencealignment第12頁/共98頁4.2.2多重比對的動(dòng)態(tài)規(guī)劃算法s1:VSN-Ss2:-SNA-s3:---ASBiocomputingtechnology—Multiplesequencealignment第13頁/共98頁前趨節(jié)點(diǎn)的個(gè)數(shù)等于2k-1問題： 計(jì)算量巨大 時(shí)間復(fù)雜度為O(2ki=1,...,k

si)

↓

O(2kNk)

Biocomputingtechnology—Multiplesequencealignment第14頁/共98頁NP-完全問題的定義Biocomputingtechnology—MultiplesequencealignmentP類問題為多項(xiàng)式界的問題；

NP類問題是這樣一類問題，如果有一個(gè)復(fù)雜度為多項(xiàng)式的算法解決其中的某個(gè)問題，則所有這些問題都在P類中；而NP-完全問題是這樣一類問題，如果其中的某個(gè)問題存在多項(xiàng)式界的算法，則NP類中的每一個(gè)問題都存在一個(gè)多項(xiàng)式界算法。

NP完全問題通常被認(rèn)為是一些人們難以在有限的時(shí)間、空間內(nèi)對問題求出最佳解的問題，幾乎所有專家都認(rèn)為不可能在多項(xiàng)式時(shí)間內(nèi)準(zhǔn)確求解NP-完全問題。第15頁/共98頁NP-完全問題的近似求解方法Biocomputingtechnology—Multiplesequencealignment舍去尋找最優(yōu)解的要求，尋找對一般問題比較接近最優(yōu)解的近似解；2.利用非常規(guī)的求解技術(shù)求解，例如，利用神經(jīng)網(wǎng)絡(luò)、遺傳算法等方法進(jìn)行問題求解。第16頁/共98頁生物信息學(xué)中NP-完全問題的近似求解方法Biocomputingtechnology—Multiplesequencealignment1.只求解規(guī)模比較小的問題；2.利用動(dòng)態(tài)規(guī)劃、分支約束等技術(shù)減小搜索空間，提高求解問題的效率；3.針對具體問題的特點(diǎn)，根據(jù)實(shí)際輸入情況，設(shè)計(jì)實(shí)用求解算法，這樣的算法雖然在最壞的情況下其時(shí)間復(fù)雜度是非多項(xiàng)式的，但是算法執(zhí)行的平均效率和復(fù)雜度與多項(xiàng)式的算法相當(dāng)；4.采用近似算法或者啟發(fā)式方法，如局部搜索、模擬退火、遺傳算法等。對基于SP模型尋找最優(yōu)多重序列比對這樣一個(gè)問題，可以用近似的方法求解，其算法的時(shí)間復(fù)雜度可用多項(xiàng)式表示。第17頁/共98頁4.2.3優(yōu)化計(jì)算方法標(biāo)準(zhǔn)動(dòng)態(tài)規(guī)劃算法存在的問題：搜索空間大剪枝技術(shù)：將搜索空間限定在一個(gè)較小的區(qū)域范圍內(nèi)。若問題是搜索一條得分最高（或代價(jià)最?。┑穆窂?，則在搜索時(shí)如果當(dāng)前路徑的得分低于某個(gè)下限（或累積代價(jià)已經(jīng)超過某個(gè)上限），則對當(dāng)前路徑進(jìn)行剪枝，即不再搜索當(dāng)前路徑的后續(xù)空間。Biocomputingtechnology—Multiplesequencealignment第18頁/共98頁Biocomputingtechnology—Multiplesequencealignment

在序列兩兩比對中,Fickett和Ukkonen設(shè)計(jì)了一種稱為定界約束過程

(boundingprocedure)的方法來縮小搜索空間,減少計(jì)算量,其中距離矩陣的上界和下界可以預(yù)先確定或動(dòng)態(tài)變化。為了在多維空間上使用動(dòng)態(tài)規(guī)劃算法，Carrillo和Lipman將這種思想引入到多重序列比對，即先進(jìn)行初步的序列雙重比對，以限制進(jìn)一步做多重序列全面比對所需要的多維空間的大小和計(jì)算量，從而克服了多序列的維數(shù)、空間和運(yùn)算量之間的矛盾第19頁/共98頁Carrillo-Lipman的優(yōu)化計(jì)算方法Biocomputingtechnology—Multiplesequencealignment

設(shè)k條序列的長度分別為n1n2……nk，按照SP得分模型計(jì)算這些序列的最優(yōu)比對。依然采用動(dòng)態(tài)規(guī)劃方法，但并不計(jì)算超晶格空間中所有的節(jié)點(diǎn)，而是僅處理與最優(yōu)路徑“相關(guān)”的節(jié)點(diǎn)。但是，哪些節(jié)點(diǎn)是相關(guān)的呢？這需要觀察節(jié)點(diǎn)在兩條序列上的投影。確定相關(guān)節(jié)點(diǎn)的方法：假設(shè)α是關(guān)于k條序列s1s2……sk

的最優(yōu)多重比對。從某個(gè)節(jié)點(diǎn)向任何兩條序列所在的平面投影，如果該投影是這兩條序列兩兩最優(yōu)比對的一部分（前面一部分），則該節(jié)點(diǎn)是與最優(yōu)比對相關(guān)的節(jié)點(diǎn)。問題的提出:第20頁/共98頁一種計(jì)算兩條序列經(jīng)過特定斷點(diǎn)的最優(yōu)比對的算法Biocomputingtechnology—Multiplesequencealignment設(shè)有兩條序列s、t,|s|=m,|t|=n;位點(diǎn)i將序列s一分為二,位點(diǎn)j將序列t一分為二,則:序列s、t對于經(jīng)過特定斷點(diǎn)（i、j）的最優(yōu)比對可分為兩個(gè)部分:①對應(yīng)于兩條序列前綴0:s:i

與0:t:j

的最優(yōu)比對;②對應(yīng)于兩條序列后綴i:s:m與j:t:n

的最優(yōu)比對;第21頁/共98頁例:Biocomputingtechnology—Multiplesequencealignment比對兩條序列:s=ATTCGG,t=GATTC打分函數(shù):p(a,a)=1,p(a,b)=-1,p(a,-)=p(-,b)=-1第22頁/共98頁Biocomputingtechnology—Multiplesequencealignment

如果超晶格空間中的一個(gè)節(jié)點(diǎn)想任意兩條序列所在的平面投影,投影在這些”斷點(diǎn)”中,則超晶格空間中的這個(gè)節(jié)點(diǎn)就是與最優(yōu)路徑相關(guān)的節(jié)點(diǎn),否則不是相關(guān)節(jié)點(diǎn).小結(jié):在進(jìn)行多重序列比對時(shí),首先要進(jìn)行序列的兩兩比對,其目的就是要找到任意兩條序列通過特定斷點(diǎn)的最優(yōu)比對,找到這些斷點(diǎn),然后,將多重比對中的超晶格空間的節(jié)點(diǎn)向任意兩條序列所在的平面投影,看看投影是否在這些斷點(diǎn)上,如果節(jié)點(diǎn)向各個(gè)平面的投影均在相應(yīng)的斷點(diǎn)上,則這個(gè)節(jié)點(diǎn)是與多重序列比對的最優(yōu)路徑相關(guān)的節(jié)點(diǎn),否則,就不是相關(guān)節(jié)點(diǎn),要進(jìn)行”剪枝”處理.第23頁/共98頁Biocomputingtechnology—Multiplesequencealignment(1)設(shè)是關(guān)于s1s2……sk

的最優(yōu)比對，如果

SP_score()≥L，則

score(i,j)≥Li,j(4-6)其中，①score(i,j

)是在si和sj所在平面投影的得分,②

這里,L實(shí)際上是最優(yōu)多重序列比對的一個(gè)下限,Lij是序列si

和序列sj

比對得分的一個(gè)下限

雖然最優(yōu)多重比對的投影不一定是兩兩最優(yōu)比對，但是我們可以為投影的得分設(shè)立一個(gè)下限。判斷超晶格空間中一個(gè)節(jié)點(diǎn)是否是相關(guān)節(jié)點(diǎn)的方法：第24頁/共98頁Biocomputingtechnology—Multiplesequencealignment(2)

現(xiàn)在的問題:

需要判斷超晶格中的一個(gè)節(jié)點(diǎn)i=(i1,i2,……ik)是否在L的限制下與最優(yōu)比對相關(guān)。(3)簡單地說,如果一個(gè)節(jié)點(diǎn)的所有投影滿足(4-6)式的條件，則該節(jié)點(diǎn)是相關(guān)的。若條件不滿足，即score(i,j)小，則該節(jié)點(diǎn)不可能是相關(guān)的，因此，i肯定不會(huì)處于最優(yōu)路徑上。(4)公式(4-6)的條件只是一個(gè)必要條件，但不是充分條件。滿足條件只是說明i可能處于最優(yōu)路徑，但不一定處于最優(yōu)路徑。公式(4-6)條件的作用是限制搜索空間，提高算法的實(shí)施效率。第25頁/共98頁Biocomputingtechnology—Multiplesequencealignment(5)將判斷條件公式(4-6)進(jìn)一步具體化，則對于所有的1≤x＜y≤k，如果i滿足

Cx,y[ix,iy]≥Lx,y（4-7）則i是相關(guān)的。這里，Cx,y是序列sx

和sy

的總得分矩陣;Cx,y[ix,iy]表示在點(diǎn)[ix,iy]處的值,即經(jīng)過[ix,iy]斷點(diǎn)的sx和sy的最優(yōu)比對得分

;[ix,iy]是斷點(diǎn);Lx,y是sx

和sy

的比對得分的下限

注意:盡管不是所有的相關(guān)節(jié)點(diǎn)均參與最優(yōu)比對,

但只有與最優(yōu)路徑相關(guān)的節(jié)點(diǎn)才參與最優(yōu)比對.第26頁/共98頁Biocomputingtechnology—Multiplesequencealignment(6)判斷非相關(guān)節(jié)點(diǎn)的方法:

假設(shè)α是關(guān)于s1s2……sk

的最優(yōu)比對，其所對應(yīng)的路徑通過節(jié)點(diǎn)i=(i1,i2,…ix,…iy,…ik),則:Cx,y[ix,iy]≥Score(i,j)≥

Lx,y反之,如果Cx,y[ix,iy]＜Lx,y,則多重序列最優(yōu)比對路徑不會(huì)經(jīng)過節(jié)點(diǎn)i=(i1,i2,…ix,…iy,…ik

),因而,該超晶格節(jié)點(diǎn)是非相關(guān)節(jié)點(diǎn).

第27頁/共98頁Biocomputingtechnology—Multiplesequencealignment

為了得到一個(gè)合理的下限L，我們可以任選一個(gè)包含所有序列的多重比對，計(jì)算其得分，以此作為L。若選取的L接近于最優(yōu)值，算法速度將大大提高。注意：L的值不能超過最優(yōu)值，否則算法出錯(cuò)。在實(shí)現(xiàn)上述優(yōu)化計(jì)算方法時(shí)，必須非常仔細(xì)。不可能對超晶格中的每一個(gè)節(jié)點(diǎn)都進(jìn)行上述判斷，否則，計(jì)算時(shí)間不會(huì)有多大的減少。我們需要一種完全消除無關(guān)單元的辦法,以便不需再處理它們。下面介紹一種可能的策略:第28頁/共98頁Biocomputingtechnology—Multiplesequencealignment

在計(jì)算機(jī)中實(shí)現(xiàn)“剪枝”技術(shù)的措施-1:

從超晶格的零點(diǎn)0=(0,0,…,0)出發(fā),此節(jié)點(diǎn)總是相關(guān)的,并影響依賴于它的節(jié)點(diǎn).每個(gè)這樣的節(jié)點(diǎn)又有它自己的受影響的節(jié)點(diǎn),如此展開,直至達(dá)到在最終角落的節(jié)點(diǎn)(n1n2…nk).(2)以數(shù)組a[]保存各節(jié)點(diǎn)的計(jì)算結(jié)果.

如果在計(jì)算a[j]時(shí)用到i,稱節(jié)點(diǎn)i影響另一個(gè)節(jié)點(diǎn)j,

又稱,節(jié)點(diǎn)j依賴于節(jié)點(diǎn)i。每個(gè)節(jié)點(diǎn)依賴于至多2k-1個(gè)其它節(jié)點(diǎn),也至多影響

2k-1個(gè)其它節(jié)點(diǎn).(3)節(jié)點(diǎn)i和節(jié)點(diǎn)j關(guān)系的另一特征是：b=j-ib是一個(gè)非零的二進(jìn)向量.第29頁/共98頁Biocomputingtechnology—Multiplesequencealignment在計(jì)算機(jī)中實(shí)現(xiàn)“剪枝”技術(shù)的措施-2:(4)為了便于處理，設(shè)置一個(gè)緩沖區(qū)，該緩沖區(qū)內(nèi)僅存放

相關(guān)節(jié)點(diǎn)的后續(xù)節(jié)點(diǎn)。(5)首先將0放入其中。(6)當(dāng)節(jié)點(diǎn)i進(jìn)入緩沖區(qū)時(shí)，其對應(yīng)的值a[i]被初始化，然后a[i]的值在隨后的階段中被更新。當(dāng)節(jié)點(diǎn)i離開緩沖區(qū)時(shí)，其值即為該點(diǎn)真正的值，并用于其他節(jié)點(diǎn)(依賴于此節(jié)點(diǎn))的計(jì)算。(7)節(jié)點(diǎn)i的后續(xù)節(jié)點(diǎn)是否要計(jì)算，還取決于i是否為相關(guān)節(jié)點(diǎn)，若不是，則不再計(jì)算其后續(xù)的其他節(jié)點(diǎn)。第30頁/共98頁具體實(shí)現(xiàn)過程:Biocomputingtechnology—Multiplesequencealignment設(shè)節(jié)點(diǎn)j是一個(gè)依賴于節(jié)點(diǎn)i的相關(guān)節(jié)點(diǎn)，如果j不在緩沖區(qū)內(nèi)，則將其放入緩沖區(qū)，并計(jì)算

a[j]a[i]+SP_Score(Colum(s,i,b))(3)如果j早已在緩沖區(qū)中，則按下式更新

a[j]max(a[j],a[i]+SP_Score(Colum(s,i,b)))注意:Carrilo-Lipman算法要求待比較的多個(gè)序列具有較大的相似性，并且序列數(shù)不能太多。第31頁/共98頁4.2.4星形比對Biocomputingtechnology—Multiplesequencealignment*啟發(fā)式方法作為首選。*啟發(fā)式方法不一定保證最終能得到最優(yōu)解，但在大多數(shù)情況下，其計(jì)算結(jié)果接近于最優(yōu)結(jié)果。*啟發(fā)式這類方法能夠大大減少所需的計(jì)算時(shí)間，加快計(jì)算速度。*漸進(jìn)法是多重序列比對中常用到的啟發(fā)式方法。其基本思想是將序列多重比對轉(zhuǎn)化為序列兩兩比對，逐漸將兩兩比對組合起來，最終形成完整的多序列比對。*組合兩兩序列比對的方法有：

星形結(jié)構(gòu)或者樹形結(jié)構(gòu)。第32頁/共98頁1.星形比對的基本思想：

星形比對是一種啟發(fā)式方法，首先由Gusfield提出。

在給定的若干序列中，選擇一個(gè)核心序列，通過該序列與其它序列的兩兩比對形成所有序列的多重比對，從而使得在核心序列和任何一個(gè)其它序列方向的投影是最優(yōu)的兩兩比對。Biocomputingtechnology—Multiplesequencealignment第33頁/共98頁2.星形比對算法概述-1Biocomputingtechnology—Multiplesequencealignment*設(shè)s1,s2,…sk是k條待比對的序列。*假設(shè)已知核心序列是sc，1≤c≤k，則可以利用標(biāo)準(zhǔn)的動(dòng)態(tài)規(guī)劃算法求出所有sc

和si

的最優(yōu)兩兩比對。這個(gè)工作的時(shí)間復(fù)雜度為O(kn2),假設(shè)所有序列的長度約為n。*將這些序列的兩兩比對聚集起來，并采用“只要是空位，則永遠(yuǎn)是空位”的原則。*聚集過程從某一個(gè)兩兩比對開始，然后逐步加上其他的兩兩比對。在這個(gè)過程中，逐步增加sc中的空位字符，以適應(yīng)其他的比對，但決不刪除sc中已存在的空位字符。第34頁/共98頁2.星形比對算法概述-2Biocomputingtechnology—Multiplesequencealignment*隨著后續(xù)比對的不斷加入,一方面我們有一個(gè)由sc指導(dǎo)的、已經(jīng)建立好的部分序列的多重比對，另一方面我們有sc與一個(gè)新序列之間的比對。在兩種比對中我們插入盡可能少而必要的空位,以保持與擴(kuò)展的sc序列相匹配。然后,將新擴(kuò)展的序列加入序列群中,且它和其它擴(kuò)展的序列具有相同的長度。*這個(gè)過程一直進(jìn)行到所有的兩兩比對都加入以后結(jié)束，這一步所需的計(jì)算量為O(k2n)。*算法總的時(shí)間復(fù)雜度為O(kn2+k2n)。第35頁/共98頁scs1s2…sk(sc,s1)(sc,s2)…(sc,sk)兩兩比對多重比對Biocomputingtechnology—Multiplesequencealignment第36頁/共98頁方法1：嘗試將每一個(gè)序列分別作為核心序列，進(jìn)行星形多重序列比對，取比對結(jié)果最好的一個(gè)。即SP得分最高的為最好。方法2：是計(jì)算所有的兩兩比對，取下式值最大的一個(gè)：3.如何選擇核心序列？Biocomputingtechnology—Multiplesequencealignment第37頁/共98頁4.舉例

有5個(gè)序列：

s1=ATTGCCATTs2=ATGGCCATTs3=ATCCAATTTTs4=ATCTTCTTs5=ACTGACCsc=s1ATTGCCATTATTGCCATT--ATTGCCATTATTGCCATTATGGCCATTATC-CAATTTTATCTTC-TTACTGACC--

(s1,s2)(s1,s3)(s1,s4)(s1,s5)ATTGCCATT--ATGGCCATT--ATC-CAATTTTATCTTC-TT--ACTGACC----Biocomputingtechnology—Multiplesequencealignment第38頁/共98頁

星形比對是一種多重序列比對近似的方法，然而是一種比較好的近似方法。如果用代價(jià)來評判多重序列的比對結(jié)果，則可以證明，用該方法所得到多重序列比對的代價(jià)不會(huì)大于最優(yōu)多重序列比對代價(jià)的兩倍。Biocomputingtechnology—Multiplesequencealignment第39頁/共98頁4.2.5樹形比對k個(gè)待比對的序列→具有k個(gè)葉節(jié)點(diǎn)的樹每個(gè)葉節(jié)點(diǎn)對應(yīng)一條序列將序列賦予樹的內(nèi)部節(jié)點(diǎn)，可以計(jì)算樹中每個(gè)分支的權(quán)值。權(quán)值代表對應(yīng)分支連接的兩個(gè)序列之間的相似性。所有權(quán)值的和就是這棵樹的得分樹形比對的問題：尋找一種樹的內(nèi)部節(jié)點(diǎn)序列賦予方式，使得樹的得分最大。星型比對可以看作是樹形比對的一種特例。Biocomputingtechnology—Multiplesequencealignment第40頁/共98頁將CT、CG、CT分別賦予節(jié)點(diǎn)x、y、z，則樹的得分為8。CTCGCT多重序列比對→兩兩序列比對 →合并兩個(gè)比對（比對的比對）

AligmentofAlignment——AA算法打分函數(shù)：

P（a,a）=1P（a,b）=0（a≠b）

P（a,-）=P（-,b）=-1111122G-TCATC-GCTG比對結(jié)果：Biocomputingtechnology—Multiplesequencealignment第41頁/共98頁AA算法概述-1Biocomputingtechnology—Multiplesequencealignment假設(shè)有兩個(gè)多重比對α1和α2

α1代表序列s1,s2,…si的多重比對

α2代表序列t1,t2,…tj的多重比對并且,(s1,s2,…si)∩(t1,t2,…tj)=Φ

α代表s1和t1的兩兩比對,則計(jì)算與α相一致的α1

和α2比對的算法如下:第42頁/共98頁AA算法概述-2Biocomputingtechnology—Multiplesequencealignment標(biāo)定α1的各列,如果s1在比對中對應(yīng)位置的編輯操作不是插入或刪除,則這些列分別標(biāo)記為s1對應(yīng)位置上的字符:a1,a2,…al1︱s1︱=l1②

標(biāo)定α2的各列,如果t1在比對中對應(yīng)位置的編輯操作不是插入或刪除,則這些列分別標(biāo)記為t1對應(yīng)位置上的字符:b1,b2,…bl2︱t1︱=l2

對a1,a2,…al1和b1,b2,…bl2進(jìn)行比對,

在所得到的比對中,對于α1、α2和α中原來有插入或刪除操作的位置,恢復(fù)其原有的實(shí)際字符或空位字符”-”.第43頁/共98頁Biocomputingtechnology—Multiplesequencealignmenta1a2a3a4b1b2b3b4b5第44頁/共98頁對于n個(gè)序列的樹形比對的基本算法過程如下：（1）初始化，對于每個(gè)序列，生成一個(gè)葉節(jié)點(diǎn)（2）利用AA算法合并兩個(gè)節(jié)點(diǎn)，形成一個(gè)新節(jié)點(diǎn)，合并的結(jié)果放在新節(jié)點(diǎn)中，原來的兩個(gè)節(jié)點(diǎn)作為新節(jié)點(diǎn)的子節(jié)點(diǎn)（3）反復(fù)執(zhí)行（2），直到形成n個(gè)葉節(jié)點(diǎn)的樹根為止，根節(jié)點(diǎn)中的序列即為最終的多重比對結(jié)果。s1s2s3s4α1α2αBiocomputingtechnology—Multiplesequencealignment第45頁/共98頁4.2.6CLUSTAL算法Biocomputingtechnology—MultiplesequencealignmentCLUSTAL算法是一種漸進(jìn)的比對方法，它是由D.G.Higgins和P.M.Sharp于1988年首次提出的。。目的：對來自不同物種的功能相同或相似的蛋白進(jìn)行比對和聚類分析，可對這些物種的親緣關(guān)系進(jìn)行判斷，并且對這些蛋白在生物進(jìn)化過程中的保守性作出估計(jì)。

第46頁/共98頁Clustal算法包括三個(gè)階段：1.先將多重序列進(jìn)行兩兩比對?；谶@些比對，計(jì)算得到一個(gè)距離矩陣，該矩陣反映每對序列之間的關(guān)系。2.根據(jù)距離矩陣計(jì)算產(chǎn)生系統(tǒng)發(fā)育樹，以此來確定被比較序列間相似的程度3.有了這棵系統(tǒng)發(fā)育樹的指導(dǎo)，從關(guān)系最緊密的兩條序列開始，逐步引入鄰近的序列或序列組，并不斷重新構(gòu)建比對，直到所有序列都被加入為止。Biocomputingtechnology—Multiplesequencealignment第47頁/共98頁舉例：Biocomputingtechnology—Multiplesequencealignment已知三級結(jié)構(gòu)的七個(gè)球蛋白序列，分別為：Hbb_Human:人的β-球蛋白Hbb_Horse:馬的β-球蛋白Hba_Human:人的α-球蛋白Hba_Horse:馬的α-球蛋白Myg_phyca:抹香鯨的血紅蛋白Glb5_petma:七鰓鰻藍(lán)血紅蛋白Lgb2_Luplu:羽扇豆的豆血紅蛋白第48頁/共98頁第一階段：兩兩比對產(chǎn)生距離矩陣Biocomputingtechnology—Multiplesequencealignment

用完全動(dòng)態(tài)規(guī)劃法計(jì)算的7個(gè)球蛋白序列之間的7×7的距離矩陣第49頁/共98頁第二階段：產(chǎn)生指導(dǎo)樹Biocomputingtechnology—Multiplesequencealignment

根據(jù)第一階段結(jié)果中的距離矩陣，用鄰近歸并法來計(jì)算產(chǎn)生一棵有分枝長度和序列權(quán)重的有根樹。第50頁/共98頁第三階段：漸進(jìn)的比對Biocomputingtechnology—Multiplesequencealignment

這一階段基本的步驟是通過一系列兩兩比對來構(gòu)建更大的比對序列組。按照指導(dǎo)樹中的分支順序，從有根樹的末梢到樹根逐漸進(jìn)行。根據(jù)圖4.22的有根樹，通過下面的順序?qū)π蛄羞M(jìn)行比對：(1)對(2)→(3)對(4)→(8)對(9)

→(5)對(10)→(6)對(11)→(7)對(12)。*每一階段使用了有殘基權(quán)重矩陣和空位開放及空位擴(kuò)展罰分的完全動(dòng)態(tài)規(guī)劃算法。*每一階段都由對兩個(gè)已經(jīng)存在的排列或序列進(jìn)行比對組成。*在先前比對中出現(xiàn)的空位仍然是固定的第51頁/共98頁圖4.23Biocomputingtechnology—Multiplesequencealignment第52頁/共98頁序列權(quán)重的作用:計(jì)算多重序列比對得分Biocomputingtechnology—Multiplesequencealignmentpeeksavtalgeekaavlalpadktnvkaaaadktnvkaaegewqlvlhvaaektkirsa第53頁/共98頁多重序列比對的統(tǒng)計(jì)特征分析

對于所得到的多重序列比對，我們往往需要進(jìn)行歸納分析，總結(jié)這些序列的特征，或者給出這些序列共性的表示。表示方式1：保守序列表示方式

表示出序列每個(gè)位置上最可能出現(xiàn)的字符或所有可能出現(xiàn)的字符表示方式2；特征統(tǒng)計(jì)圖譜（profile）表示方式（或特征統(tǒng)計(jì)矩陣表示方法）Biocomputingtechnology—Multiplesequencealignment第54頁/共98頁令P=(P1,P2,…Pj,…,PL)；P表示在多重比對的每一列上各種字符出現(xiàn)的概率分布Pj=(pj0,pj1,…,pj|A|)；A代表字母表，

Pjk代表字母表A中第k個(gè)字符在第j列出現(xiàn)的概率。

pj0：第0個(gè)字符是特殊的空位符號“-”。稱P為多重序列比對的特征統(tǒng)計(jì)圖譜或特征統(tǒng)計(jì)矩陣。

P反映了多重序列比對各個(gè)列的特征。表示方法2：特征統(tǒng)計(jì)圖譜Biocomputingtechnology—Multiplesequencealignment第55頁/共98頁

12345(位置)A0.80.20.20.60.0T0.00.40.60.41.0C0.20.20.20.00.0G0.00.20.00.00.0

（堿基）Biocomputingtechnology—Multiplesequencealignment

ATTAT AACTT CTTAT ACTTT AGAAT第56頁/共98頁利用保守序列或者特征統(tǒng)計(jì)圖可以判斷一個(gè)序列是否滿足一定的特征給定一個(gè)序列s=a1a2…am，定義字符a在第j位的代價(jià)為

其中，|A|代表字母表A的長度，

Ak代表A的第k個(gè)字符，

A0代表空缺字符“-”。整個(gè)序列s的代價(jià)為一條序列與特征統(tǒng)計(jì)圖相對照，如果代價(jià)值小，說明該序列具有相應(yīng)的特征，否則該序列不具備相應(yīng)的特征。Biocomputingtechnology—Multiplesequencealignment第57頁/共98頁4.2.7隱馬爾可夫模型HMMBiocomputingtechnology—MultiplesequencealignmentHMM——HiddenMarkovModel1.馬爾可夫鏈（復(fù)習(xí)）①具有馬爾可夫性的離散狀態(tài)隨機(jī)過程稱為馬爾可夫鏈。②馬爾可夫性,即無后效性,若已知現(xiàn)在的狀態(tài),將來與過去無關(guān)。③離散狀態(tài)的馬爾可夫鏈的定義:

考慮只取有限個(gè)或可數(shù)個(gè)狀態(tài)的隨機(jī)過程:{Xn,n=0,1,2,…}

假設(shè)對一切狀態(tài):i1,i2,…in-1,i,j,一切n≥0,有:P(Xn+1=j∣Xn=i,……X0=i0)=P(Xn+1=j∣Xn=i)成立,

則稱隨機(jī)過程:{Xn,n=0,1,2,…}為離散狀態(tài)的馬爾可夫鏈。第58頁/共98頁1.馬爾可夫鏈（復(fù)習(xí)）Biocomputingtechnology—Multiplesequencealignment④P(Xn+1=j∣Xn=i)稱為馬爾可夫鏈的一步轉(zhuǎn)移概率,

記為Pij(n,n+1)。⑤若轉(zhuǎn)移概率Pij(n,n+1)不依賴n,則可簡記為Pij,則稱此馬爾可夫鏈?zhǔn)菚r(shí)齊的,否則是非時(shí)齊的.⑥用P記轉(zhuǎn)移概率所對應(yīng)的矩陣,即轉(zhuǎn)移矩陣(transitionmatrix)

⑦一個(gè)馬爾可夫鏈的概率分布完全由它的轉(zhuǎn)移矩陣與初始分布決定.第59頁/共98頁2.隱馬爾可夫模型HMMBiocomputingtechnology—Multiplesequencealignment

隱馬爾可夫模型是一種概率模型,它是由馬爾可夫鏈發(fā)展擴(kuò)充而來的一種隨機(jī)模型。這種方法已經(jīng)成功地應(yīng)用于多個(gè)領(lǐng)域,如語音識別、光學(xué)字符識別等。

HMM在生物信息學(xué)領(lǐng)域中也有著重要應(yīng)用。如：序列分析，基因識別等。

HMM可以被理解為一個(gè)雙重隨機(jī)過程：

(1)系統(tǒng)狀態(tài)變化的隨機(jī)過程，

(2)由狀態(tài)決定輸出的隨機(jī)過程。

一個(gè)隱馬爾可夫鏈{Xt,Yt}包含兩部分：一個(gè)潛在的、不可觀察的有限狀態(tài)馬爾可夫鏈{Xt}

一個(gè)外顯的、可觀察的隨機(jī)過程{Yt}，

Yt的分布依賴于Xt。第60頁/共98頁2.隱馬爾可夫模型HMMBiocomputingtechnology—MultiplesequencealignmentHMM是用概率統(tǒng)計(jì)的方法來描述時(shí)變信號的過程.HMM是一個(gè)動(dòng)態(tài)的統(tǒng)計(jì)模型，可以用有限狀態(tài)機(jī)FSM來描述。

(FSM——finitestatemachine)

有限狀態(tài)機(jī)可以從一種狀態(tài)轉(zhuǎn)移到另一種狀態(tài)，每個(gè)狀態(tài)轉(zhuǎn)換有不同的概率。某狀態(tài)是否轉(zhuǎn)移到下一狀態(tài)取決于該狀態(tài)的狀態(tài)轉(zhuǎn)移概率，而在某一狀態(tài)下能看到哪一個(gè)觀測值，取決于該狀態(tài)的觀測概率。第61頁/共98頁隱馬爾可夫模型HMMBiocomputingtechnology—Multiplesequencealignment

記HMM模型為:

其中:A為狀態(tài)轉(zhuǎn)移概率矩陣，

B為觀察概率矩陣，為初始狀態(tài)分布。

HMM模型的建模工作主要為確定這三個(gè)參數(shù)。

第62頁/共98頁HMM的三個(gè)經(jīng)典問題:Biocomputingtechnology—Multiplesequencealignment問題1（評測問題，evaluation）：已知模型和輸出序列O，求由生成O的概率。問題2（譯解問題，decoding）：已知模型和輸出序列O，求最有可能生成O的狀態(tài)轉(zhuǎn)移序列。問題3（學(xué)習(xí)問題，learning）：已知模型和輸出序列O，求最有可能生成O時(shí)模型的參數(shù)。第63頁/共98頁P(yáng)rofile——概形、譜Biocomputingtechnology—Multiplesequencealignment*概形是對一組序列進(jìn)行全局多重比對時(shí)被發(fā)現(xiàn)的,是將比對中具有較高保守區(qū)域變成小的多重比對,然后得到多重比對記分矩陣.*概形由更像小的多重排列的列構(gòu)成,可以包括:匹配、失配、插入、缺失.*概形一旦生成,就可用于尋找一個(gè)可能與之匹配的目標(biāo)序列,

它利用表中記分來評價(jià)每個(gè)位置的可能性.例:25個(gè)氨基酸長的概形表格,有25列,每列將有20個(gè)記分值.

一列中每個(gè)匹配氨基酸的記分都在概形中對應(yīng)的位置上.缺點(diǎn)：偏向性第64頁/共98頁P(yáng)rofileHMM

(1)模型結(jié)構(gòu)Biocomputingtechnology—Multiplesequencealignment*對于生物序列而言,HMM的字符當(dāng)然是20個(gè)字母的氨基酸或4個(gè)字母的核苷酸。但依據(jù)不同的問題,其它的一些字符也可以使用,如64個(gè)密碼子的三聯(lián)體字母,3個(gè)字母(α,β,coil)的二級結(jié)構(gòu)等.*編碼蛋白質(zhì)的原始DNA序列,在生物的進(jìn)化過程中會(huì)受到自然環(huán)境和各種因素的影響,使翻譯出的蛋白質(zhì)序列經(jīng)歷突變、遺失或引入外源序列等變化,最后按不同的進(jìn)化路徑分化,形成多種功能相近的蛋白質(zhì)。所以,可以把這些蛋白質(zhì)看作由一個(gè)基本蛋白質(zhì)序列經(jīng)過插入、刪除或替換了某些氨基酸殘基而形成的。這個(gè)過程可以用HMM來表示。第65頁/共98頁圖4.9Biocomputingtechnology—Multiplesequencealignment圖中:方形代表匹配狀態(tài)(主狀態(tài))，即輸出的氨基酸和基本序列中對應(yīng)的氨基酸匹配;

圓形表示刪除或缺失狀態(tài)，即從原始蛋白質(zhì)序列中去掉一個(gè)特定的氨基酸。

菱形表示氨基酸的插入，即在原始蛋白質(zhì)序列插入一個(gè)氨基酸。各狀態(tài)間的箭頭表示狀態(tài)間的轉(zhuǎn)換途徑。

注意:①不同的參數(shù)會(huì)使模型以不同的概率產(chǎn)生新的氨基酸。

②一個(gè)訓(xùn)練好的模型可以代表有共同特征的蛋白質(zhì)序列。第66頁/共98頁圖4.10Biocomputingtechnology—Multiplesequencealignment第67頁/共98頁P(yáng)rofileHMM與標(biāo)準(zhǔn)的Profile的比較Biocomputingtechnology—Multiplesequencealignment①ProfileHMM有正規(guī)的概率作基礎(chǔ)，對于序列的刪除和插入狀態(tài)的記分也有較為可靠的理論依據(jù)。而標(biāo)準(zhǔn)的

Profile純粹是一種啟發(fā)式的方法。②

HMM用統(tǒng)計(jì)方法估計(jì)序列某一位點(diǎn)核苷酸或氨基酸殘基出現(xiàn)的真正概率，而標(biāo)準(zhǔn)的Profile卻是用自身的觀察頻率給核苷酸或氨基酸殘基指派分?jǐn)?shù)。③由于②,ProfileHMM方法從10至20個(gè)核苷酸序列構(gòu)成的比對中提取的信息，相當(dāng)于用標(biāo)準(zhǔn)的Profile從40至50個(gè)核苷酸序列構(gòu)成的比對中提取的信息。第68頁/共98頁(2)用HMM給序列打分Biocomputingtechnology—Multiplesequencealignment*因訓(xùn)練好的HMM代表一個(gè)蛋白質(zhì)家族，我們可以用它對序列進(jìn)行打分，根據(jù)分值來衡量這條序列是否屬于由此HMM

代表的蛋白質(zhì)家族。得分高于閾值，證明這條序列是這個(gè)家族的成員；得分低于閾值，說明它不是家族的成員。第69頁/共98頁HMM用于分析蛋白質(zhì)序列的原理:Biocomputingtechnology—Multiplesequencealignment

分析蛋白質(zhì)產(chǎn)生不同序列的概率。對于與模型相符合的序列，能以較大的概率產(chǎn)生該序列；若不與該模型符合的序例，則按此模型產(chǎn)生該序列的概率會(huì)較小。任何一條序列都可以由HMM中的一條路徑所呈現(xiàn)。對于給定的模型，計(jì)算產(chǎn)生某條序列的概率就是計(jì)算這條路徑上的所有輸出概率與轉(zhuǎn)移概率的乘積。第70頁/共98頁圖4.11Biocomputingtechnology—Multiplesequencealignment第71頁/共98頁Biocomputingtechnology—Multiplesequencealignment*如果已知確切的狀態(tài)路徑，這個(gè)計(jì)算是很簡單的。*在一個(gè)實(shí)際模型中，可以通過不同的狀態(tài)路徑產(chǎn)生同一條序列。因此，產(chǎn)生一條序列的正確的概率是所有可能的路徑產(chǎn)生該序列的概率的總和。*這種計(jì)算復(fù)雜并且難以實(shí)施，除非是一條很短的序列。*以下介紹兩種很好的替代方法：

①

Viterbi算法——可以計(jì)算出由模型產(chǎn)生某序列的最有可能的路徑。②前向算法(forwardalgorithm)——用于計(jì)算所有路徑的概率和。第72頁/共98頁①Viterbi算法Biocomputingtechnology—Multiplesequencealignment

作為一種優(yōu)化算法，動(dòng)態(tài)規(guī)劃算法在生物信息處理中得到廣泛的應(yīng)用。前期課程，我們已經(jīng)介紹了動(dòng)態(tài)規(guī)劃算法在序列比較中的應(yīng)用；這里，介紹應(yīng)用動(dòng)態(tài)規(guī)劃算法求解HMM模型的最優(yōu)路徑的方法。這個(gè)方法就是著名的Viterbi算法。第73頁/共98頁圖4.12D1D2D3I0I1I2I3M1M2M30.06Biocomputingtechnology—Multiplesequencealignment最優(yōu)路徑：開始I0M1M2M3結(jié)束在這條路徑上產(chǎn)生“ACCY”的概率為：5.9728*10-5第74頁/共98頁圖4.130.00150.00460.48500.2231Biocomputingtechnology—Multiplesequencealignment第75頁/共98頁②

前向算法forwardalgorithmBiocomputingtechnology—Multiplesequencealignment

給定一個(gè)HMM模型M和一個(gè)字符序列X(x1,x2,…xl),假定產(chǎn)生序列X的對應(yīng)路徑未知,要求計(jì)算模型M產(chǎn)生X的概率P(X/M).*Viterbi算法是在可以產(chǎn)生序列X的各種路徑中,選擇一條最優(yōu)π*,使得P(X/π*)最大.*現(xiàn)在的問題:既然有多條路徑可以產(chǎn)生序列X,那么,模型M

產(chǎn)生序列X總的可能性有多大?*解決的方法:類似于Viterbi算法,用求和替代求最大值.第76頁/共98頁圖4.141.8*10-45.52*10-43.0912*10-43.3849*10-76.8965*10-5Biocomputingtechnology—Multiplesequencealignment圖4.12的HHM產(chǎn)生“ACCY”序列的概率為：3.3849*10-7+6.8965*10-5=6.9303*10-5第77頁/共98頁③

局部和全局打分Biocomputingtechnology—Multiplesequencealignment

有時(shí)用局部打分會(huì)更合理，即序列的分值由子序列的最高分值所替代。例：圖4.11中ACCY序列，將概率轉(zhuǎn)換為對數(shù)形式，序列的全局打分就是以下四個(gè)分值之和：-13.25。

Score(Ａ)=ln（0.4）+ln（0.3）=-2.12Score(Ｃ)=ln（0.46）+ln（0.6）=-1.29Score(Ｃ)=ln（0.97）+ln（0.5）=-0.72Score(Ｙ)=ln（0.015）+ln（0.73）+ln（0.01）=-9.12(-2.12)+(-1.29)+(-0.72)+(-9.12)=-13.25分值太低了，以至于不能判定ACCY是否為被模擬的蛋白質(zhì)家族的成員。第78頁/共98頁Biocomputingtechnology—Multiplesequencealignment

假設(shè)ACCY是圖4.15中所示蛋白質(zhì)家族的成員，在這種情況下，全局分就不是一個(gè)很好的衡量標(biāo)準(zhǔn)。如果采用局部打分，最后分值就會(huì)高得多。我們發(fā)現(xiàn)最高分的子序列為CC，分值為–2.01。這個(gè)分值就足以判定ACCY為此家族的成員。這時(shí)，局部打分比全局打分更精確。第79頁/共98頁(3)模型的訓(xùn)練Biocomputingtechnology—Multiplesequencealignment①根據(jù)訓(xùn)練用蛋白質(zhì)序列的平均長度確定模型中各個(gè)狀態(tài)的數(shù)量而建立HMM。②選用訓(xùn)練算法對HMM進(jìn)行訓(xùn)練，即調(diào)整該模型的參數(shù)（轉(zhuǎn)移概率和輸出概率），使得該模型適用于訓(xùn)練所用的序列，并能夠以最大的可能性產(chǎn)生參與訓(xùn)練的觀察序列。第80頁/共98頁(3)模型的訓(xùn)練-最大似然法MLBiocomputingtechnology—Multiplesequencealignment給定在訓(xùn)練數(shù)據(jù)集中的一組序列：

s(1)，s(2)，…，s(ｎ)，一個(gè)模型究竟在多大程度上適合該數(shù)據(jù)集，可由根據(jù)該模型觀測到每一個(gè)序列的聯(lián)合概率來表征：其中，是模型產(chǎn)生第j個(gè)序列的概率。上式稱為該模型的似然，最大似然（maximumlikelihood，ML）其原理可用于測量模型的性能，即選取模型參數(shù)使上式值最大。第81頁/共98頁(3)模型的訓(xùn)練-最大后驗(yàn)概率MAPBiocomputingtechnology—Multiplesequencealignment基本思想：在給定序列下，使模型的后驗(yàn)概率最大。假設(shè)：對所有可能的模型參數(shù)存在優(yōu)先概率分布，這種概率分布體現(xiàn)模型的特點(diǎn)。為模型優(yōu)先概率分布，為歸一化常數(shù)。其中：第82頁/共98頁(3)模型的訓(xùn)練-Baum-Welch算法Biocomputingtechnology—Multiplesequencealignment方法:Baum-Welch算法目的:

給定觀測值序列,通過計(jì)算確定一個(gè)模型λ,使得

P(O/λ)最大.Baum-Welch算法的基本思路:①初始模型λ0(待訓(xùn)練模型)

②

基于λ0以及觀察值序列,訓(xùn)練新模型λ③如果logP(X/λ)-logP(X/λ0)＜δ,說明訓(xùn)練已經(jīng)達(dá)到預(yù)期效果,算法結(jié)束.

④

否則,令λ0=λ,繼續(xù)第②步工作第83頁/共98頁(4)序列加權(quán)Biocomputingtechnology—Multiplesequencealignment

如果訓(xùn)練集中有一小組序列具有高度相似性，那么訓(xùn)練出來的HMM就會(huì)過于專門化。如果HMM過于專門化,那么在對其它序列進(jìn)行打分時(shí),有可能將那些本該屬于此