東南大學chapter5-1DNA序列分析

主講人：孫嘯

制作人：劉志華東南大學吳健雄實驗室第五章DNA序列分析第五章DNA序列分析

DNA序列分析

——基因序列

——基因表達調(diào)控信息

尋找基因牽涉到兩個方面的工作：識別與基因相關(guān)的特殊序列信號預測基因的編碼區(qū)域結(jié)合兩個方面的結(jié)果確定基因的位置和結(jié)構(gòu)

基因表達調(diào)控信息隱藏在基因的上游區(qū)域，在組成上具有一定的特征，可以通過序列分析識別這些特征。

第一節(jié)DNA序列分析步驟和分析結(jié)果評價在DNA序列中，除了基因之外，還包含許多其它信息，這些信息大部分與核酸的結(jié)構(gòu)特征相關(guān)聯(lián)，通常決定了DNA與蛋白質(zhì)或者DNA與RNA的相互作用。存放這些信息的DNA片段稱為功能位點如啟動子（Promoter）、基因終止序列（Terminatorsequence）、剪切位點（Splicesite）等。發(fā)現(xiàn)重復元素數(shù)據(jù)庫搜索分析功能位點序列組成統(tǒng)計分析綜合分析一個基本的DNA序列分析方案功能序列分析的準確性來自于對“功能序列”和“非功能序列”的辨別能力。兩個集合：訓練集（trainingset）用于建立完成識別任務的數(shù)學模型。 測試集或控制集（controlset）用于檢驗所建模型的正確性。用訓練集中實例對預測模型進行訓練，使之通過學習后具有正確處理和辨別能力。然后，用模型對測試集中的實例進行“功能”與“非功能”的判斷，根據(jù)判斷結(jié)果計算模識別的準確性。收集已知的功能序列和非功能序列實例（這些序列之間是非相關(guān)的）訓練集（trainingset）測試集或控制集（controlset）建立完成識別任務的模型檢驗所建模型的正確性對預測模型進行訓練，使之通過學習后具有正確處理和辨別能力。進行“功能”與“非功能”的判斷，根據(jù)判斷結(jié)果計算模識別的準確性。識別“功能序列”和“非功能序列”的過程

Sn

——敏感性Sp——特異性Tp是正確識別的功能序列數(shù)，Tn為正確識別的非功能序列數(shù)，F(xiàn)n是被錯誤識別為非功能序列的功能序列數(shù)，F(xiàn)p是被錯誤識別為功能序列的非功能序列數(shù)。敏感性和特異性的權(quán)衡對于一個實用程序，既要求有較高的敏感性，也要求有較高的特異性。如果敏感性很高，但特異性比較低，則在實際應用中會產(chǎn)生高比率的假陽性；相反，如果特異性很高，而敏感性比較低，則會產(chǎn)生高比率的假陰性。對于敏感性和特異性需要進行權(quán)衡，給出綜合評價指標。對于一個識別程序準確性可按下式進行綜合評價：另一個綜合評介指標為相關(guān)系數(shù)，其計算計算公式為：選擇訓練集和測試集在檢測算法的可行性時，需要從已知的數(shù)據(jù)中按照不同的方式選擇訓練集和測試集測試集的構(gòu)成非常關(guān)鍵在不同的測試集上進行測試可能會得到不同的準確性結(jié)果，甚至準確性相差很大。建立標準的功能序列測試集合。如基因轉(zhuǎn)錄剪切位點的測試集合、編碼區(qū)域的測試集合等。第二節(jié)核苷酸關(guān)聯(lián)分析對于一個給定的基因組，最簡單的計算就是統(tǒng)計DNA序列中各類核苷酸出現(xiàn)的頻率。對于隨機分布的DNA序列，每種核苷酸的出現(xiàn)是均勻分布的出現(xiàn)頻率各為0.25。而真實基因組的核苷酸分布則是非均勻的核苷酸

頻率

A0.3248693727808C0.1751306272192G0.1751306272192T0.3248693727808酵母基因組核苷酸出現(xiàn)頻率在統(tǒng)計過程中，如果同時計算DNA的正反兩條鏈，則根據(jù)堿基配對原則，A和T、C和G的出現(xiàn)頻率相同。如果僅統(tǒng)計一條鏈，則雖然A和T、C和G的出現(xiàn)頻率不同，但是非常接近。核苷酸

頻率

A0.344C0.155G0.157T0.343單鏈核苷酸出現(xiàn)頻率

基因和其它功能區(qū)域在正反兩條鏈上出現(xiàn)的可能性通常一樣核苷酸出現(xiàn)頻率也不應該有偏差正反兩條鏈在信息的組織結(jié)構(gòu)方面不應該有差別單鏈上A和T、C和G的出現(xiàn)頻率相近。正反兩條鏈堿基互補的原則

單鏈上A和T、C和G的出現(xiàn)頻率相近的解釋兩聯(lián)核苷酸頻率不同基因組中兩個連續(xù)核苷酸出現(xiàn)的頻率也是不相同的4種核苷酸可以組合成16種兩聯(lián)核苷酸酵母基因組兩聯(lián)核苷酸頻率表對酵母基因組兩聯(lián)核苷酸的統(tǒng)計結(jié)果其中核苷酸對出現(xiàn)頻率最高的達到0.119而出現(xiàn)頻率最低的只有0.028令:

Pij

——代表兩聯(lián)核苷酸（i，j）的出現(xiàn)頻率

Pi——代表核苷酸i的出現(xiàn)頻率則：

Pij’=Pij/(PiPj)

的值反應核苷酸i和j的關(guān)聯(lián)關(guān)系如果Pij’=1，則在兩個連續(xù)的位置上，核苷酸i和j的出現(xiàn)是相對獨立的。關(guān)聯(lián)性分析

對于酵母基因組

PA=0.3248PAA=0.1193PAA’=0.1193/（0.3248*0.3248） =1.131>1

表明在兩個連續(xù)位置上“A”的出現(xiàn)不是獨立的，而是相關(guān)的。關(guān)聯(lián)性分析

同樣，對于相隔一定距離k（k代表核苷酸個數(shù)）的兩個核苷酸，也可能具有一定的相關(guān)性。假設(shè)Pij(k)代表核苷酸j出現(xiàn)在核苷酸i之后第k個位置的頻率，則可定義一個反應統(tǒng)計相關(guān)性的互信息I(k)I(k)值得大小實際上反應了距離為k的兩個核苷酸之間的相關(guān)性的程度三聯(lián)核苷酸——基因密碼子在進行編碼區(qū)域識別時，常常需要對三聯(lián)核苷酸進行統(tǒng)計分析，這實際上是分析密碼子的使用偏性。由于密碼子的簡并性（degeneracy），每個氨基酸至少對應1種密碼子，最多有6種對應的密碼子。在基因中，同義密碼子的使用并不是完全一致的。不同物種、不同生物體的基因密碼子使用存在著很大的差異基因密碼子的使用與基因編碼的蛋白的結(jié)構(gòu)和功能有關(guān)，與基因表達的生理功能有著密切的聯(lián)系蛋白的三級結(jié)構(gòu)與密碼子使用概率有密切的關(guān)系通過對密碼子的聚類分析，可以很清晰地將具有不同三級結(jié)構(gòu)蛋白質(zhì)的編碼基因分成不同的類，而具有相似三級結(jié)構(gòu)蛋白的編碼基因則大致聚在同一類中，從而證明基因密碼子的使用偏性與蛋白質(zhì)三級結(jié)構(gòu)具有密切的相關(guān)性。在不同物種中，類型相同的基因具有相近的同義密碼子使用偏性對于同一類型的基因由物種引起的同義密碼子使用偏性的差異較小針對酵母第一染色體的分析結(jié)果第三節(jié)功能位點分析功能位點（functionalsite）與特定功能相關(guān)的位點，是生物分子序列上的一個功能單元，或者是生物分子序列上一個較短的片段。功能位點又稱為功能序列（functionalsequence）、序列模式（motif）、信號（signal）等。核酸序列中的功能位點包括轉(zhuǎn)錄因子結(jié)合位點、轉(zhuǎn)錄剪切位點、翻譯起始位點等。在蛋白質(zhì)序列分析中，常使用序列模式這個名詞，蛋白質(zhì)的序列模式往往與蛋白質(zhì)結(jié)構(gòu)域或者作用部位有關(guān)。功能位點示意基因組序列中若干個相鄰的功能位點組合形成功能區(qū)域（functionalregion）。功能位點分析的任務發(fā)現(xiàn)功能位點特征識別功能位點1、利用共有序列搜索功能位點共有序列（consensus）又稱一致性片段共有序列是關(guān)于功能位點特征的描述，它描述了功能位點每個位置上核苷酸進化的保守性

例如:NTATN利用共有序列進行功能位點分析牽涉到兩個方面的問題，如何構(gòu)造共有序列如何利用共有序列在給定的核酸序列上搜索尋找功能位點，并計算所找到的功能位點的可靠性共有序列具有以下幾個方面的特征：（1）共有序列中既有保守的位置，也有可變的位置；（2）任何位置上的核苷酸可以用15種類型之一來表示：核苷酸表示符號符

號含

義說

明GG腺嘌呤AA鳥嘌呤TT胸腺嘧啶CC胞嘧啶RGorA嘌呤YTorC嘧啶MAorC氨基KGorT羧基SGorC強氫鍵(3個氫鍵)WAorT弱氫鍵(2個氫鍵)HAorCorT非GBGorTorC非AVGorCorA非T(非U)DGorAorT非CNGorAorTorC任意堿基共有序列構(gòu)造過程:(1)初始化共有序列為一系列可變位置，以“N”代表；(2)在可變位置尋找出現(xiàn)次數(shù)最多的核苷酸，并將該位置轉(zhuǎn)化為保守位置；(3)對當前所得到的共有序列進行特異性檢查，若通過檢查，轉(zhuǎn)（5），否則轉(zhuǎn)（4）；(4)形成與當前共有序列一致的位點子集，轉(zhuǎn)（2）；(5)從原位點集合中刪除與當前共有序列一致的位點，若還有剩余位點，則轉(zhuǎn)（1），構(gòu)造另外的共有序列。TTATGATATATACGCTTGTC

TCCAC

TTATGATATATACGCTTGTC

TCCAC

TNNNNtTATG

tACGC

tTGTC

tCCAC

tTATG

tACGC

tTGTC

tCCAC

TNNNC

[1]

[2]

[3]

[4]

[2]

[3]

NNNNNTNNNN非特異

TNNNC非特異

tACGc

tTGTc

tCCAc

[4]

[2]

tACGc

tTGTc

tCCAc

[3]

TNSNC特異

[5]

Consensus1：

TNSNC剩余位點：

TTATG

ATATA

[5]

Consensus2：

NTATNTNNSC在給定的序列中搜索與共有序列一致的序列片段數(shù)據(jù)庫搜索共有序列表示方法的缺點：是關(guān)于序列特征的一種定性描述，對于DNA序列，它能夠說明序列每個位置可能出現(xiàn)的堿基類型，但是不能準確地說明各位置上不同類型堿基出現(xiàn)的可能性大小。2、用感知矩陣分析功能位點用權(quán)系數(shù)描述功能位點各位置上每種核苷酸的相對重要性感知矩陣（或加權(quán)矩陣）根據(jù)一系列功能位點的多重對比排列結(jié)果而建立的其大小為4n4代表堿基的種類數(shù)目，n代表功能位點的長度矩陣的每一個元素M(a,j)的值代表第a種核苷酸在功能位點第j個位置上出現(xiàn)的得分,a{A,T,G,C}。123456A18227-319T26142-10G3110-50-19C5-916880感知矩陣示例對于一個序列s=a1a2…an，根據(jù)對應位置上核苷酸的類型，取感知矩陣中對應的權(quán)值，加和以后得到該序列的得分設(shè)S=ATTGCA，則

Ws=1+6+14-5+8+19=43T——功能位點閾值T‘——非功能位點閾值如果WsT，則S是功能位點；如果WsT'，則S是非功能位點。感知矩陣M的構(gòu)造算法令A+代表功能位點集合

A-代表非功能位點集合過程如下：（1）初始化M為零矩陣；（2）執(zhí)行過程（3）-（6）的循環(huán)；（3）逐步取訓練集合中的每個實例Si，如果SiA+，轉(zhuǎn) 過程（4）；如果SiA-，轉(zhuǎn)過程（5）；（4）如果W（Si）T，M不變，否則根據(jù)Si的核苷酸分布 將M中所有對應元素的值加1；轉(zhuǎn)（6）；（5）如果W（Si）T‘，M不變，否則根據(jù)Si的核苷酸分 布將M中所有對應元素的值減1；轉(zhuǎn)（6）；（6）若訓練集合中的所有實例都處理過，則循環(huán)結(jié)束， 轉(zhuǎn)（7），否則繼續(xù)執(zhí)行循環(huán)體，直到處理完所有實 例；（7）如果M穩(wěn)定，則結(jié)束；否則轉(zhuǎn)（2）。上述算法反復調(diào)整感知矩陣M的元素值，直到M矩陣能夠正確識別訓練集中的所有功能位點和非功能位點。對于最終得到的感知矩陣，要求其具有敏感性和特異性，每一列上的元素值應該盡可能地有明顯的差別，以便反應功能位點各個位置上的特點。與感知矩陣類似，如果令矩陣每一個元素M(a,j)的值代表第a種核苷酸在功能位點第j個位置上出現(xiàn)的概率，則M是一個概率矩陣。假設(shè)各個位置上出現(xiàn)的堿基是相互獨立的，即任何兩個位置上的堿基是不相關(guān)的，那么對于給定一個序列s=a1a2…an，可以計算出功能位點序列為s的概率：如果分別統(tǒng)計功能位點和非功能位點，通過計算可以形成兩個矩陣M和M’，進一步計算可以判斷一個給定的序列究竟屬于功能位點，還是屬于非功能位點。給定一個序列s=a1a2…an，定義似然比LR(M,M’,s)：在進行功能位點檢測時，計算LR(M,M’,s)，并與給定的閾值L比較，如果LR(M,M’,s)>L，則序列s可能是一個功能位點。概率矩陣M和M’的每個元素是一個0和1之間的正數(shù)。如果令一個4n新矩陣U的元素(a,j)的值為

log2(M(a,j)/M’(a,j))則矩陣U的每個元素值可能是正值，也可能是負值。實際上，矩陣U就是感知矩陣。第四節(jié)隱馬爾柯夫模型1、馬爾柯夫鏈（Markovchain）考慮一個具有多個狀態(tài)的系統(tǒng)S，S={s1,s2,…,s|s|},令S0、S1、…、St為一系列在各個時刻系統(tǒng)狀態(tài)的變量，即狀態(tài)鏈。對于每個1到|S|的整數(shù)，它們分別與狀態(tài)鏈中的一個狀態(tài)相聯(lián)系，并且在任何時刻，這條鏈都處于一個特殊的狀態(tài)。當且僅當對于任何t有

則St形成一條馬爾柯夫鏈。簡單地說，就是系統(tǒng)未來的狀態(tài)僅依賴于當前狀態(tài)。St稱為在時刻t系統(tǒng)鏈的狀態(tài)。一條馬爾柯夫鏈完全決定于初始分布P(S0)和轉(zhuǎn)換概率Pt=P(St+1|St)。令狀態(tài)轉(zhuǎn)換矩陣為

F=(fij)

fij代表從狀態(tài)si移動到狀態(tài)sj的概率。生物序列可以被描述為一個隨機過程的輸出，其中對于一個給定的核酸在位置p出現(xiàn)的概率依賴于已占據(jù)前面k個位置的核酸，這樣一種表示稱為k階馬爾柯夫模型。ATCGTAGCAT…….一個序列具有不同的統(tǒng)計性質(zhì) （如二目頻率或三目周期性）不同的功能區(qū)域（如編碼區(qū)域、非編碼區(qū)域）對應于不同的馬爾柯夫模型。馬爾柯夫鏈在識別CpG島中的應用CpG島是一類長度在幾百bp的特殊DNA序列，其中CG核苷酸對出現(xiàn)的頻率非常高。

ACGCGCGTACGCGAATCpG島在基因組中有重要的生物學意義，而識別CpG島有助于在基因組序列中確定我們感興趣的區(qū)域。CpG島的識別問題表述為：給定一段DNA序列

X=(x1，x2，…,xL)，確定X是否是一個CpG島。設(shè)字母表A={a,t,c,g}，對于字母表中的任何兩個字符s、t，定義轉(zhuǎn)換概率為fst=p(xi=t|xi-1=s)，即字符s后面出現(xiàn)字符t的概率。假設(shè){xi}是一個隨機過程，隨機變量xi的取值僅依賴于xi-1，即對于所有x1,x2,…,xiA，整個序列X的發(fā)生概率為為了處理方便，添加兩個特殊的字符‘B’（begin）和‘E’（end），使得x0=‘B’,xL+1=‘E’，則上述公式簡化為：令fst+為CpG

島內(nèi)的字符轉(zhuǎn)換概率

fst-為CpG

島外的字符轉(zhuǎn)換概率 則X的對數(shù)似然得分為上述計算值越大，則X越可能是CpG島。

CpG島內(nèi)部和外部的轉(zhuǎn)換概率

另外一個待解決的問題是：給定DNA序列，確定CpG島的位置。 直接的方法：對窗口內(nèi)的子序列計算得分Score(Xk)，具有正值的Xk

就是可能的CpG島子序列起始位置為k+1，長度為l問題： 事先不知道CpG島的長度 但是假設(shè)CpG島的長度為l如果l比較大，而真實的CpG島又比較小，則上述概率計算值不足以證實CpG島；如果l取值比較小，則難以找出整個CpG島。這是該算法的最大不足之處，需要考慮其他的算法。HMM2、隱馬爾柯夫模型（HMM）功能位點的正則表達式來表示相當于一致性序列這里的正則表達式描述了一個功能位點的構(gòu)成規(guī)律，或者說描述了功能位點各個位置上核苷酸的組成。TGCC—AGG ???ACAC—ATC問題： 對于每個位置，僅僅說明可能的取值，而沒有說明各種取值出現(xiàn)的可能性大小例如，用這樣的方法無法區(qū)分下面兩條序列究竟哪一個更可能屬于功能位點：

TGCC--AGG ACAC—ATC第一個序列中，假設(shè)所有位置上都是取已知出現(xiàn)次數(shù)最少的字符而對于第二個序列，所有位置上都是取已知出現(xiàn)次數(shù)最多的字符。顯然，第一個序列幾乎肯定不是功能位點，而第二個序列幾乎可以肯定是功能位點，但是用正則表達式表達卻無法將兩種極端的情況分開。隱馬爾柯夫模型可以用于生物序列分析，該模型在生物信息分析方面有重要的應用。一階隱馬爾柯夫模型包括有限數(shù)目的系統(tǒng)狀態(tài)、離散的字母表、狀態(tài)轉(zhuǎn)換矩陣和字符釋放概率。一個HMM模型是一個三元組M=（A，S，Θ）A是字母表S是有限狀態(tài)集合，每個狀態(tài)可以釋放字母表中的字符。Θ為概率集合，包括兩個部分：狀態(tài)轉(zhuǎn)換概率fkl

k,lS，表示從狀態(tài)k轉(zhuǎn)換到狀態(tài)l的概率；字符釋放概率，記為ek(b)

kS,bA

表示在狀態(tài)k下釋放出字符b的概率。令路徑

=（1,2,…,L）是一個相繼狀態(tài)序列

X=(x1,x2,…,xL)是一個字符序列按下述方式定義狀態(tài)轉(zhuǎn)換概率和字符釋放概率：對于給定的路徑，可以按下面的公式計算出產(chǎn)生序列X的概率：

這里，令0為起始狀態(tài)，L+1為終止狀態(tài)。例如，對于前面給出的兩個序列ACACATC和TGCTAGG，它們的得分分別為：P（ACACATC）=0.81.00.81.00.80.60.40.61.01.00.81.00.8=4.710-2P(TGCTAGG)=0.21.00.21.00.20.60.20.61.01.00.21.00.2=0.002310-2從上述計算結(jié)果可以看出，兩個序列差別非常大

一個功能位點的HMM模型是通過對一系列的功能位點實例進行機器學習而形成的用這樣的模型可以定量的計算一個序列片段是功能位點的可能性計算方法是從模型的第一個狀態(tài)出發(fā)，根據(jù)序列的核苷酸組成，將相應的狀態(tài)值與狀態(tài)轉(zhuǎn)換值連乘，結(jié)束于最后一個狀態(tài)一個檢測CpG島的HMM模型有8個狀態(tài)，狀態(tài)名稱和釋放的字符為：狀態(tài)：A+C+G+T+A-C-G-T-

釋放字符：ACGTACGT

其中，帶有“+”號的狀態(tài)表示在CpG島內(nèi)部，用“-”號標記的狀態(tài)代表CpG島外部。假設(shè)字符處于CpG島內(nèi)的概率是p

處于CpG島外的概率是q

可以得到狀態(tài)轉(zhuǎn)換概率CpG島HMM模型中的狀態(tài)轉(zhuǎn)換概率解碼問題： 給定一個隱馬爾柯夫模型M=（A，S，Θ）和一個字符序列X，在M中為X尋找一條最優(yōu)路徑*，在路徑中的每一個狀態(tài)都選擇釋放一個字符，要求使得P(X|*)最大，記為：

在處理CpG島問題中，最優(yōu)路徑可以幫助我們尋找CpG島所在的位置。如果找到最優(yōu)路徑*，則這條路徑穿過的“+”狀態(tài)將對應于CpG島。3、Viterbi算法求解HMM模型的最優(yōu)路徑基本思想： 動態(tài)規(guī)劃算法給定一個字符序列X=(x1,x2,…,xL)

，以vk(i)代表序列前綴（x1,x2,…,xi）終止于狀態(tài)k（kS,1≤i≤L）的最可能路徑的概率。求解過程如下：（1）初始化

（2）對于每個i=0,…,L-1及每個lS，按下式進行遞歸計算：

（3）最后，計算序列X終止于狀態(tài)“end”最可能的路徑概率，即P(X|*)的值在正向的遞歸計算過程中，保持向前推進的反向指針，這樣，在正向計算完成后，根據(jù)反向指針重構(gòu)最優(yōu)路徑*。算法的時間復雜度為O(L|S|2),空間復雜度為O(L|S|)。在概率的計算過程中，需要使用大量的乘法運算，在有限計算精度的情況下，會產(chǎn)生誤差。如果使用對數(shù)值，可以解決這個問題。因此，以vk(i)代表序列前綴（x1,x2,…,xi）終止于狀態(tài)k（kS,1≤i≤L）的最可能路徑的對數(shù)得分值，則初值按如下方式設(shè)置遞歸計算及最終得分計算改為（5-26）4、前向概率和反向概率給定一個隱馬爾可夫模型M=（A，S，Θ） 一個字符序列X=(x1,x2,…,xL)

要求計算模型M產(chǎn)生X的概率P(X|M)與最優(yōu)路徑問題不一樣 前面的問題是在可以產(chǎn)生序列X的各種路徑中，選擇一條最優(yōu)路徑*，使得P(X|*)最大。 而現(xiàn)在的問題是： 既然有多條路徑可以產(chǎn)生序列X，那么模型M產(chǎn)生序列X總的可能性有多大？如果有一條從狀態(tài)“begin”出發(fā)，終止于狀態(tài)“end”的路徑=（0,1,2,…,L,L+1），其中0=“begin”，L=1=“end”,該路徑中各狀態(tài)所釋放的字符組成的序列與X相同，則模型M產(chǎn)生X的概率為這里代表所有那些從狀態(tài)“begin”出發(fā)、終止于狀態(tài)“end”的路徑。（5-27）

（5-28）

由于一個HMM模型中可能的路徑非常多，窮舉每條路徑顯然是不合適的。下面介紹解決該問題的前向算法（forwardalgorithm）與反向算法(backwardalgorithm)。算法的根本任務是對于每個1≤i≤L及kS,計算概率P(i=k|X，M)。定一個序列X=(x1,x2,…,xL)，令k(i)為釋放前綴（x1,x2,…,xi）后到達狀態(tài)i=k的概率。前向算法初始值的設(shè)置與Viterbi算法一樣：遞歸計算過程和最終計算如下：（5-29）

（5-30）

（5-31）

與前向算法相對應，給定一個序列X=(x1,x2,…,xL)，令k(i)為在給定狀態(tài)i=k下后綴（xi+1,xi+2,…,xL）的概率。反向算法初始化如下：遞歸計算和終止計算如下：（5-32）

（5-33）

（5-35）

（5-34）

利用正向和反向概率，可以計算出P(i=k|X)。由于HMM的階數(shù)為1，當前的狀態(tài)僅依賴于前一個狀態(tài)，則根據(jù)條件概率的定義，我們得到解（5-36）（5-37）5、HMM模型的參數(shù)估計應用中假設(shè)有一個HMM模型，其中的狀態(tài)轉(zhuǎn)換概率和字符釋放概率都是已知的。然而在實際中，情況并非如此。 我們所知道的僅僅是一些實例問題是要根據(jù)給定的n個字符串重構(gòu)M，使得M產(chǎn)生這n個字符串具有最大的概率。由于各個字符串是獨立產(chǎn)生的，則若使用對數(shù)表示，則目標就是尋找一個*，使得其中這里的n個字符串X(1),X(2),…,

X(n)通常被稱為“訓練序列”。（5-38）（5-39）（5-40）特殊情況：假設(shè)已知與字符串序列X(1),X(2),…,X(n)

相對應的狀態(tài)序列(1),(2),…,(n)，可以計算從狀態(tài)k到狀態(tài)l的轉(zhuǎn)換數(shù)Fkl和在狀態(tài)k下釋放字符b的次數(shù)Ek(b)。則關(guān)于最大似然估計值為：為了避免零概率，當處理數(shù)量較少的樣本時，需要對Fkl和Ek(b)進行修正：rkl、rk(b)為拉普拉斯修正項，通常情況下為1，可以解釋為預先假設(shè)的均勻分布。但是在某些情況下，這些修正項可能取其他的值，例如已知狀態(tài)轉(zhuǎn)換或字符釋放的信息，或已有的先驗知識。在一般情況，不知道狀態(tài)序列(1),(2),…,(n)，這時，尋找最優(yōu)參數(shù)集在數(shù)學上是一個NP-完全問題，可以用Baum-Welch算法或期望最大（EM）算法解決這個問題。具體的求解算法如下：（1）初始化，給中的參數(shù)賦予初值；（2）計算從狀態(tài)k到狀態(tài)l轉(zhuǎn)換的期望次數(shù)，使用與計算P(X,i=k)時相同的參數(shù)（見公式5-36），則（5-45）這樣，對所有訓練序列X(j)（j=1,…,n）的所有位置i（i=1,…,L(j)，L(j)為序列X(j)的長度）進行求和運算，按下式計算期望值Fkl：其中k(j)(i)是針對序列X(j)的正向計算結(jié)果，k(j)(i+1)是反向計算結(jié)果。接下來計算在狀態(tài)k釋放字符b的期望次數(shù)：（5-47）

（5-46）

（3）重新計算的參數(shù)值Fkl和Ek(b)，正如在第一種情況所做的一樣（參見公式（5-41）和公式（5-42））；（4）反復執(zhí)行步驟（2）、（3），直到Score(X(1),X(2),…,

X(n)|)的增量小于給定的一個值很小的參數(shù)為止。EM算法保證目標函數(shù)Score(X(1),X(2),…,

X(n)|)單調(diào)增加，并且概率的對數(shù)值接近于0，保證算法收斂。需要注意，收斂的是目標函數(shù)，并非是的參數(shù)。當目標函數(shù)變化趨緩時，的參數(shù)值可能波動較大，這意味著算法所得到的結(jié)果不穩(wěn)定。Baum算法的主要問題是目標函數(shù)存在若干局部極大，算法不能保證找到全局最大點，算法收斂的點可能是局部極大點?？朔植繕O大缺陷的一種方法是執(zhí)行算法若干遍，每次給取不同的初始值。如果算法多次計算結(jié)果到達同一個極大點，則可以認為該點是全局最大點。6、基于HMM模型的序列比對可以利用HMM將一個序列與一個序列統(tǒng)計特征（profile）進行比對，從而解決多重比對問題。定義一個長度為L的序列統(tǒng)計特征P是一系列的概率集合ei(b)，ei(b)表示在第i（1≤i≤L）個位置上出現(xiàn)字母表中字符b的概率。這樣，在給定條件P下序列X=(x1,x2,…,xL)發(fā)生的概率為：如果不考慮“空位”，則X與P的比對得分為：這里，p(b)是字符b的背景出現(xiàn)頻率。（5-49）定義一個基本HMM模型，有L個“匹配”狀態(tài)M1,M2,…,ML，它們對應與統(tǒng)計特征的匹配。所有這些狀態(tài)順序連接起來，即狀態(tài)Mj連接到后繼Mj+1，如圖5.5所示。其中從狀態(tài)Mj釋放字符b的概率為ej(b)。為了在比對中允許插入“空位”的操作，在上述基本模型中加入“插入”狀態(tài)I1,I2,…,IL，并假設(shè)每個插入狀態(tài)Ij,有一個來自相應匹配狀態(tài)Mj的連接，有一個到匹配狀態(tài)Mj+1的連接，還有一個自循環(huán)連接。根據(jù)“空位”的懲罰原則，給這些狀態(tài)轉(zhuǎn)換賦予適當?shù)母怕省＃?-50）

同樣，為了允許“刪除”操作，可以進一步假如“刪除”狀態(tài)D1,D2,…,DL，這些狀態(tài)不能釋放任何字符。刪除狀態(tài)依然順序連接，同時增加從Dj到Ij的連接及從Ij到Dj+1的連接。完整的HMM模型如圖5.5所示：D1D2D3I2I3I4BeginEndM1M2M3I1

圖5.5用于序列多重比對的HMM模型

下面介紹一種Viterbi類似算法，將X=(x1,x2,…,xm)與長度等于L的統(tǒng)計特征P進行比對。對于每一個1≤j≤L和1≤i≤m，定義：（1）vjM(i)代表子序列(x1,x2,…,xi)與HMM模型P的匹配對數(shù)得分值，該匹配以狀態(tài)Mj釋放字符xi作為最后操作;

（2）vjI(i)代表子序列(x1,x2,…,xi)與HMM模型P的匹配對數(shù)得分值，該匹配以狀態(tài)Ij釋放字符xi作為最后操作;

（3）vjD(i)代表子序列(x1,x2,…,xi)與HMM模型P的匹配對數(shù)得分值，該匹配以狀態(tài)Dj結(jié)束（不釋放任何字符）。模型P中特殊狀態(tài)“begin”的初始值為：為了計算vjM(i)、vjI(i)和vjD(i)的值，使用Viterbi算法中的相同技術(shù)，但現(xiàn)在的模型有兩個特點：（1）模型中的每一個狀態(tài)最多只有3個引入連接，如上圖所示；（2）“刪除”狀態(tài)不釋放任何字符。（5-51）“匹配”狀態(tài)Mj的三個前驅(qū)同屬于上一層,即j-1層，有“插入”狀態(tài)Ij的三個前驅(qū)屬于同一層,