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

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

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

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

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

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

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

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

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

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

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

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

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

14、為F，反之則為L我們將上面的結果圖示如下： -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個位置應用前向算法和后向算法的概率：顯然，如果將換成，換成則“擲骰子問題”即刻變成蛋白質二級結構預測問題了。結合前面最大路徑的求解，我們馬上可以得到蛋白質二級結構預測的兩個方法。足見將數學方法徹底搞懂的好處。圖6-20 300個位置上各種狀態（F與L）的概率四、后驗解碼問題的應用后驗解碼問題的應用之一就是我們前面剛提到的選取某個位置概率最大的狀態，這些狀態同樣也可以組成一條路徑即： （6-28）我們將式（6-18）與最大幾率狀態序列（路徑）相比： （6-28）可以發現：后驗解碼問題將各個位置的狀態概率分開來考慮，其路徑是由每個位置概率最大的狀態組成，但由此路徑與符