隱馬爾科夫鏈模型

1、第三節(jié) 后驗解碼問題 前面我們介紹了在給定一個符號序列時，在不同的可能狀態(tài)序列中找出概率（可能性）最大的狀態(tài)序列（或路徑）。在實際情況中可能會碰到一個問題：不同的路徑其概率相差不大，這為我們選擇哪一條路徑增加了一定的困難。但我們知道一條狀態(tài)序列，它是由各個位置的狀態(tài)組成的，比如圖6-11中，從第1個位置到第45個位置均由狀態(tài)F組成，緊接后面的21個是L狀態(tài)，依此類推。因此克服上述困難我們可以每個位置各個不同狀態(tài)的概率，這里面我們必須分清一點：某一條序列對應(yīng)于某個狀態(tài)序列的概率含義與狀態(tài)序列中某個位置上具體某個狀態(tài)的概率是不同的，我們將它們的區(qū)分總結(jié)于如圖6-14：圖6-14 某個DNA序列對應(yīng)

2、的不同可能的CpG島分布（狀態(tài)序列）的概率與各個位置上某個狀態(tài)的概率之間的區(qū)別（“+”代表CpG島區(qū)域；“-”代表非CpG島區(qū)域）。我們在上節(jié)的“Viterbi”算法中重點介紹的是如何計算圖6-14中右邊部分，并將其中最大的概率對應(yīng)的路徑找出來。而我們這一節(jié)要介紹的是如何計算圖6-14中最下面的概率，這就是通常我們所說的“后驗解碼問題”：計算概率為多少？為了計算這個概率，首先得介紹所給定序列對應(yīng)的概率，由于不同的狀態(tài)序列可以產(chǎn)生同一條序列，因此，我們有圖6-15的描述：圖6-15 隱馬爾克夫鏈中符號序列的概率與路徑的關(guān)系（這里假設(shè)符號序列集只有兩條序列與，對應(yīng)的狀態(tài)序列（路徑）中也只有兩個路徑

3、與。 圖6-15中的與就是要求的某個符號序列的概率，由于此時路徑與已知，因此它與符號序列的聯(lián)合概率表達式可寫為：（6-11）公式（6-11）羅列了一大堆符號，這對學(xué)生物學(xué)的人來說要比較快地接受它們可能會要多花點時間，為使我們有一個清晰的直觀印象，我們這里以兩條DNA序列及相應(yīng)的CpG島分布圖來說明。圖6-16 DNA序列中CpG島分布的隱馬爾科夫鏈模型中在公式（6-6）中的示意圖如果我們以代表中的某個序列，相應(yīng)的狀態(tài)序列有，則有： (6-12)在實際情況中，可能的路徑（狀態(tài)序列）很多，而且隨著符號序列的長度的增加，姿態(tài)態(tài)序列的個數(shù)N呈指數(shù)函數(shù)的方式往上增長，以CpG島為例：當(dāng)序列長度為3時，其

4、可能的狀態(tài)序列個數(shù)為，當(dāng)序列長度為4時，則有。因此如何計算（6-11）則需要一定的技巧。借鑒Viterbi算法，我們也不妨來個沿著符號序列逐步計算，這種逐步的方式有兩種：一種是從符號序列的左端向右端（DNA序列的5端到3端，蛋白質(zhì)序列的N端到C端），相應(yīng)的算法叫前向算法（forward algorithm）；另一種是從右端向左端，相應(yīng)的算法為后向算法(backward algorithm)。我們首先介紹前向算法。一、前向算法(Forward Algorithm) 在介紹前向算法之前，我們先看一下式（6-7）的第一個公式即： （6-13）從中我們可發(fā)現(xiàn)它取的是每個位置最大值概率，因為對符號序列，

5、其每個位置的概率等其各個狀態(tài)概率之和，因此，整個符號序列的概率是每個位置概率之積，沿著符號序列的方向，我們有： (6-14)整個算法可寫成如下：初始化： (6-15a)迭代過程： (6-15b)終止： (6-15c)我們可將前向算法以圖示方式表示：圖6-17 前向算法的基本框架圖圖6-18 前向算法的上體計算圖。這里我們僅選狀態(tài)1的計算過程，目的是使圖形顯得簡潔明了，因為其它狀的算法與此是相同的。計算例子：符號序列“3151”。我們?nèi)砸訴iterbi算法中的例子即符號序列“3151”來說明：計算第一個點：初始條件如Viterbi算法所示，同時設(shè)定。首先，我們計算第一個位置即，我們有：當(dāng)時，我們

6、有：當(dāng)，我們有：當(dāng)時，我們有：在終止階段，我們有：與Viterbi算法相似，前向算法也存著“因計算機浮點數(shù)的限制導(dǎo)致后面的結(jié)果失真”，因此同樣的需要用對數(shù)運算來代替相應(yīng)的乘法運算，有文獻曾試圖通過如下的轉(zhuǎn)換：克服這一點。但我們認(rèn)為這樣做其實沒有實質(zhì)性的意義，因為中間一個“exp”的運算，如果log(p)小到一定程度，exp(log(p)仍舊受到浮點數(shù)的影響。因此，我們并不認(rèn)為這是一個好的處理方式，甚至認(rèn)為這種做法除了增加運算量外，沒有其它的實際意義。有人曾引進一個標(biāo)度因子，我們認(rèn)為這種方法相比較理想，而且，我們也自編了相應(yīng)的程序，發(fā)現(xiàn)它的確可以克服計算機浮點數(shù)有限的限制。為此，我們這里詳細介紹

7、這個方法。該方法的基本原理是將每一步（或每一時刻）進行歸一化計算，具體為：每一步的前向概率進行如下轉(zhuǎn)換： (6-16a)即有 (6-16b)通過這樣的轉(zhuǎn)換，我們可以得到式（6-15b）的迭代公式為： (6-17)通過這樣的轉(zhuǎn)換，要求轉(zhuǎn)換后的符合如下條件： (6-18)根據(jù)公式（6-18），我們可以得到式（6-17）中的歸一化因子的表達式為： (6-19)應(yīng)用這種方法計算前向概率與后向概率顯然可以在基本上克服計算機浮點數(shù)帶來的不足，因為計算得到的其值肯定在狀態(tài)數(shù)倒數(shù)值附近浮動，與序列長度相比，值顯然要小得多，如在蛋白質(zhì)二級結(jié)構(gòu)中有，在CpG島中有。應(yīng)用這種校正方法有一個很顯著的優(yōu)點就是計算序列的

8、概率值即，它與歸一化因子的關(guān)系為： (6-20a)因此很自然地就可以將它轉(zhuǎn)化為相應(yīng)的對數(shù)即： (6-20b)這里提前說明一下：有關(guān)文獻在介紹這種歸一化因子法是，將我們剛才前面介紹的前向算法的歸一化因子也用在我們接下來要介紹的后向算法中。我們通過深入的研究發(fā)現(xiàn)：在后向算法中直接應(yīng)用前向算法中的歸一化因子不是很妥當(dāng)，在計算后驗解碼問題中“計算機浮點數(shù)有限”的問題仍未得到有效的解決。為此，我們在后向算法中重新計算歸一化因子，具體地將在后向算法中介紹。雖然這種改變是微小的，但它在后驗解碼問題及后面我們要介紹的HMM模型的概率參數(shù)即的求解算法之一的Bauch-Welch算法中會發(fā)現(xiàn)其具有非常獨特的優(yōu)越性

9、。二、后向算法（Backward Algorithm）后向算法基本上與前向算法基本類似，只不過計算的方向不同：從符號序列的右邊向左邊運算（對DNA序列是從3端到5端；對蛋白質(zhì)序列則是從C端到N端）。下面我們直接將其算法敘述如下：初始化： 對所有狀態(tài)，有 (6-20a)迭代過程： (6-20b)終止： (6-20c)計算例子：符號序列“3151”我們?nèi)砸苑栃蛄小?151”來說明后向算法的計算過程：計算第一個點：設(shè)定，根據(jù)初始條件的設(shè)置我們有：首先，我們計算倒數(shù)第二個位置即，此時，我們有： 其次，計算倒數(shù)第三個位置即時，我們有： 最后，我們計算倒數(shù)第四個位置也就是第一個位置即時，我們有：這樣我們

10、計算得到的有： 我們將后向算法的結(jié)果與前向算法的結(jié)果相比，它們基本一致，其差異來自于計算過程中的舍入誤差。 我們在前一小節(jié)中提到“如果將前向算法的歸一化因子用于后向算法，則在后驗解碼問題及概率參數(shù)的求解仍未擺脫計算機浮點數(shù)有限的束縛”，并說明可以重新對后向算法各迭代步驟中的概率進行歸一化，具體的方法很簡單，其歸一化因子為： (6-21)這里的代表后向算法中的歸一化因子，轉(zhuǎn)化后的后向概率即與其實際的概率即的關(guān)系為： (6-22)這樣我們就有： (6-23)三 后驗解碼問題的算法本節(jié)一開頭我們就提到過：計算后驗解碼問題需要前向算法與后向算法?，F(xiàn)在我們已經(jīng)介紹完它們，自然地，接下來就要計算后驗解碼問

11、題即計算： 根據(jù)條件概率公式： (6-24)我們有： (6-25)而上式右邊的分子則為：于是式（6-20）可轉(zhuǎn)化為： (6-26)當(dāng)然可以用前向算法或后向算法得到。從式（6-25）我們可以看出的具體意義就是對某個符號序列，當(dāng)?shù)趥€位置的符號為狀態(tài)時的概率。在CpG島問題中，如果為“+”，則說明這個位置的堿基在CpG島區(qū)域內(nèi)；如果為“-”，則說明這個位置的堿基不在CpG島區(qū)域內(nèi)?，F(xiàn)在我們根據(jù)式（6-26）來說明針對后向算法我們自己提出的歸一化因子在后驗解碼中的意義。將實際的前向概率與后向概率即與分別應(yīng)用歸一化后的概率即與來表示，即將式(6-16b)，（6-20a）及（6-23）代入式（6-26），

12、我們有： （6-27）由于一般會在“1”附近，因此上式中它們的乘積部分顯然不會太小或也不會太大，即可以將它們的乘積控制在計算機浮點數(shù)允許的范圍內(nèi)。自然就會擺脫了計算機浮點數(shù)有限的束縛。按老辦法，我們現(xiàn)在仍以前面的例子即符號序列“來說明“后驗解碼”的計算方式。后驗解碼的計算例子顯然我們有，根據(jù)我們前面計算得到的一系列與，我們很容易就可以計算出各個位置中各個狀態(tài)的概率即。當(dāng)時，我們有：當(dāng)時，我們有：當(dāng)時，我們有：當(dāng)時，我們有：以上我們僅選了前4個點的計算，為了使讀者對全部數(shù)據(jù)序列有一個比較直觀的認(rèn)識，我們應(yīng)用C語言編制了一個計算機程序，將應(yīng)用前向算法和后向算法計算得到的各位置上的概率即列于圖6-1

13、5，同時將也將各個位置上對應(yīng)的各狀態(tài)概率即總結(jié)在圖6-16中。同時，應(yīng)用前向算法與后向算法計算得到該序列的概率均為：。 讀者如果對隱馬爾科夫鏈方法感興趣，可以根據(jù)我們提供的例子手算幾個位置，這樣可以加深對隱馬爾科法的印象。而且，正如我們前面所說的，一旦我們熟悉了具體的算法，則可以很“自然”將它們應(yīng)用到我們熟悉的其的領(lǐng)域中去，顯然對擴大隱馬爾科夫鏈方法的應(yīng)用以及生物信息學(xué)新方法的建立是很有幫助的。此外，根據(jù)圖6-20，我們顯然也可以據(jù)此判斷每個位置是“F態(tài)”還是“L態(tài)”。具體的方法是如果位置的某個狀態(tài)概率大于另一個狀態(tài)概率，則可以確定張三用的是這個狀態(tài)所代表的骰子，具體為： 如果，則該位置的狀態(tài)

14、為F，反之則為L我們將上面的結(jié)果圖示如下： -3.00 -5.32 -7.58 -9.76 -11.86 -12.27 -14.57 -16.81 -17.35 -19.68 -21.96 -22.55 -23.29 -25.67 -28.03 -30.35 -32.60 -34.79 -36.88 -38.90 -40.85 -42.76 -44.64 -46.50 -46.74 -49.00 -51.20 -53.30 -53.72 -56.02 -58.26 -60.41 -62.49 -64.48 -66.42 -68.32 -68.58 -70.85 -73.05 -75.16 -7

15、7.19 -79.16 -81.07 -82.95 -84.82 -85.06 -87.32 -89.52 -90.02 -92.33 -2.48 -4.27 -6.07 -7.89 -9.72 -11.55 -13.34 -15.16 -16.98 -18.75 -20.56 -22.37 -24.13 -25.76 -27.49 -29.27 -31.08 -32.90 -34.73 -36.56 -38.39 -40.23 -42.06 -43.90 -45.73 -47.54 -49.36 -51.18 -53.01 -54.81 -56.62 -58.44 -60.27 -62.10

16、 -63.93 -65.77 -67.60 -69.41 -71.23 -73.05 -74.88 -76.72 -78.55 -80.38 -82.22 -84.05 -85.86 -87.68 -89.50 -91.29-516.29 -514.19 -512.00 -509.75 -507.43 -506.68 -504.42 -502.10 -501.35 -499.08 -496.76 -496.01 -495.36 -493.52 -491.68 -489.83 -487.97 -486.10 -484.20 -482.27 -480.27 -478.20 -476.04 -473

17、.80 -473.11 -471.08 -468.97 -466.77 -466.11 -464.25 -462.36 -460.44 -458.47 -456.43 -454.31 -452.10 -451.43 -449.52 -447.55 -445.52 -443.41 -441.21 -438.94 -436.61 -434.25 -433.47 -431.11 -428.73 -427.94 -425.55-515.35 -513.52 -511.70 -509.89 -508.11 -506.38 -504.57 -502.79 -501.06 -499.25 -497.47 -

18、495.74 -493.93 -492.10 -490.26 -488.43 -486.59 -484.76 -482.92 -481.09 -479.26 -477.43 -475.61 -473.80 -472.00 -470.17 -468.35 -466.53 -464.72 -462.89 -461.05 -459.22 -457.39 -455.56 -453.73 -451.92 -450.11 -448.28 -446.45 -444.62 -442.79 -440.97 -439.17 -437.39 -435.67 -434.08 -432.34 -430.71 -429.

19、37 -427.78圖6-19 前50個位置應(yīng)用前向算法和后向算法的概率：顯然，如果將換成，換成則“擲骰子問題”即刻變成蛋白質(zhì)二級結(jié)構(gòu)預(yù)測問題了。結(jié)合前面最大路徑的求解，我們馬上可以得到蛋白質(zhì)二級結(jié)構(gòu)預(yù)測的兩個方法。足見將數(shù)學(xué)方法徹底搞懂的好處。圖6-20 300個位置上各種狀態(tài)（F與L）的概率四、后驗解碼問題的應(yīng)用后驗解碼問題的應(yīng)用之一就是我們前面剛提到的選取某個位置概率最大的狀態(tài)，這些狀態(tài)同樣也可以組成一條路徑即： （6-28）我們將式（6-18）與最大幾率狀態(tài)序列（路徑）相比： （6-28）可以發(fā)現(xiàn)：后驗解碼問題將各個位置的狀態(tài)概率分開來考慮，其路徑是由每個位置概率最大的狀態(tài)組成，但由此路徑與符