東南大學(xué)chapter5-1DNA序列分析

主講人：孫嘯

制作人：劉志華東南大學(xué)吳健雄實(shí)驗(yàn)室第五章DNA序列分析第五章DNA序列分析

DNA序列分析

——基因序列

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

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

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

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

Sn

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

頻率

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

頻率

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

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

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

Pij

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

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

Pij’=Pij/(PiPj)

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

對(duì)于酵母基因組

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

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

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

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

號(hào)含

義說(shuō)

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

TTATG

ATATA

[5]

Consensus2：

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

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

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

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

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

F=(fij)

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

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

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

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

fst-為CpG

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

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

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

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

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

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

這里，令0為起始狀態(tài)，L+1為終止?fàn)顟B(tài)。例如，對(duì)于前面給出的兩個(gè)序列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從上述計(jì)算結(jié)果可以看出，兩個(gè)序列差別非常大

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

釋放字符：ACGTACGT

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

處于CpG島外的概率是q

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

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

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

（2）對(duì)于每個(gè)i=0,…,L-1及每個(gè)lS，按下式進(jìn)行遞歸計(jì)算：

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

要求計(jì)算模型M產(chǎn)生X的概率P(X|M)與最優(yōu)路徑問(wèn)題不一樣 前面的問(wèn)題是在可以產(chǎn)生序列X的各種路徑中，選擇一條最優(yōu)路徑*，使得P(X|*)最大。 而現(xiàn)在的問(wè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）

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

（5-30）

（5-31）

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

（5-33）

（5-35）

（5-34）

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

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

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

（5-46）

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

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

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

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

圖5.5用于序列多重比對(duì)的HMM模型

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

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

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