什么是隱馬爾可夫模型
隱馬爾可夫模型（Hidden Markov Model，HMM） 是統(tǒng)計模型，它用來描述一個含有隱含未知參數(shù)的馬爾可夫過程。其難點是從可觀察的參數(shù)中確定該過程的隱含參數(shù)。然后利用這些參數(shù)來作進一步的分析，例如模式識別，特別是我們今天要講的基因預(yù)測。是在被建模的系統(tǒng)被認為是一個馬爾可夫過程【一段組裝好的序列】與未觀測到的（隱藏的）的狀態(tài)【哪些是編碼區(qū)哪些不是】的統(tǒng)計馬爾可夫模型。
下面用一個簡單的例子來闡述：
假設(shè)我手里有兩個顏色不同的骰子，一個是橘色（Coding，C）的另一個是藍色（Noncoding,N）的。但是和平常的骰子不同的是，他們穩(wěn)定下來只要有四種可能，也就是上下是固定的（你可以理解為他們只會平行旋轉(zhuǎn)）。這樣每個骰子出現(xiàn)ATCG的概率都是1/4.
假設(shè)我們開始投骰子，我們先從兩種顏色選一個，挑到每種骰子的概率都是1/2。然后我們擲骰子，我們得到ATCG中的一個。不停地重復(fù)以上過程，我們將會得到一串序列，每個字符都是ATCG中的一個。例如CGAAAAAATCG
這串序列就叫做可見狀態(tài)鏈。但是在隱馬爾可夫模型中，我們不僅僅有這么一串可見狀態(tài)鏈，還有一串隱含狀態(tài)鏈。在這個例子里，這串隱含狀態(tài)鏈就是你用的骰子的序列。比如，隱含狀態(tài)鏈有可能是：C C N N N N N N N C C C。
一般來說，HMM中說到的馬爾可夫鏈其實是指隱含狀態(tài)鏈，因為隱含狀態(tài)（骰子）之間存在轉(zhuǎn)換概率（transition probability）。在我們這個例子里，C的下一個狀態(tài)是N。C，N的概率都是1/2。這樣設(shè)定是為了最開始容易說清楚，但是我們其實是可以隨意設(shè)定轉(zhuǎn)換概率的。比如，我們可以這樣定義，N后面不能接兩個C，或C 的概率是0.1。這樣就是一個新的HMM。
同樣的，盡管可見狀態(tài)之間沒有轉(zhuǎn)換概率，但是隱含狀態(tài)和可見狀態(tài)之間有一個概率叫做輸出概率（emission probability）。就我們的例子來說，Coding（C）產(chǎn)生A的概率是1/4，Noncoding（N）產(chǎn)生A的概率是1/4，當(dāng)然這些概率是我人為定義的，你可以定義為其它的值。
其實對于HMM來說，如果提前知道所有隱含狀態(tài)之間的轉(zhuǎn)換概率和所有隱含狀態(tài)到所有可見狀態(tài)之間的生成概率，做模擬是相當(dāng)容易的。對，我們就撿容易的來。
用隱馬爾可夫模型做基因預(yù)測
接下來，我們來做一個最簡單的基因預(yù)測。給定一段基因組DNA序列，我們來預(yù)測其中的編碼區(qū)。按照之前說道的隱馬爾可夫模型，我們先要區(qū)分不能直接觀測到的隱狀態(tài)和可以直接觀測的顯符號。
在這個例子里，我們可以很容易看出，給定的基因組DNA序列是可以觀測到的符號串。而編碼/非編碼是不能直接觀測的隱狀態(tài)。因此，我們可以畫出狀態(tài)轉(zhuǎn)換圖。首先我們有編碼和非編碼兩個狀態(tài)。因為基因組會同時包好編碼和非編碼區(qū)域，因此這兩個狀態(tài)之間可以轉(zhuǎn)換。當(dāng)然，每個狀態(tài)也可以轉(zhuǎn)換到自己，表示連續(xù)的編碼或非編碼區(qū)域。這樣，我們就有了一個2*2的轉(zhuǎn)移矩陣。
接下來，我們需要寫生成概率。這個也很直觀，無論是編碼還是非編碼狀態(tài)，都有可能產(chǎn)生A,C,G,T四個堿基，因此我們由可以分別有兩個矩陣。
現(xiàn)在，我們需要一個訓(xùn)練集（Training set）,來把這三個矩陣的格子中填上具體的數(shù)值。具體來說，我們需要一個已經(jīng)事先注釋過得——也就是說正確標記了編碼，非編碼區(qū)域的DNA序列，通常這個序列要比較長，以便有充足的數(shù)據(jù)統(tǒng)計。
假定我們經(jīng)過訓(xùn)練集的分析，分別填好了轉(zhuǎn)移矩陣概率和生產(chǎn)概率矩陣。我們需要根據(jù)這些數(shù)據(jù)，來對一個未知的給定基因組序列反推出最可能的狀態(tài)路徑，也就是概率最大的那個狀態(tài)路徑。因此，我們還是和之前一樣利用動態(tài)規(guī)劃的算法，寫出迭代公式，以及最后的終止點公式（Termination equation）
從公式里面，我們看到，我們需要做大量測乘法。這個不僅比較慢，而且利用計算機操作時，隨著連乘次數(shù)的增加，很容易數(shù)值過小而出現(xiàn)下溢（underflow）的問題。因此，我們通常會引入對數(shù)計算，從而將乘法轉(zhuǎn)換成加法。具體來說，就是對轉(zhuǎn)移和生成概率都預(yù)先取log10。
好，我們正式開始，假設(shè)我組裝了一段序列（咦，怎么這么短？為了簡單_）：
CGAAAAAATCG
首先，讓我們和之前一樣，畫出動態(tài)規(guī)劃的迭代矩陣，其中包含兩個狀態(tài)，非編碼狀態(tài)N與編碼狀態(tài)C。接下來，我們需要設(shè)定邊界條件（boundary condition）,也就是這兩個狀態(tài)默認的分布比例。為了計算方便，我們分別設(shè)為0.8和0.2，經(jīng)過log10轉(zhuǎn)換后，分別為-0.097和-0.699.接下來我們逐步填格子
之后，我們碰到的第一個堿基C,由生成概率可知C在非編碼狀態(tài)下的log10生成的概率是-0.523，將之與-0.097相加，就可以得到-0.62。類似的，在編碼狀態(tài)下，這個數(shù)是-0.699+（-0.699）=-1.4。
接下來，我們要前進一個堿基，就需要進行狀態(tài)轉(zhuǎn)移，我們先來看第一種情形，也就是從非編碼狀態(tài)到非編碼狀態(tài)的轉(zhuǎn)換。從轉(zhuǎn)移舉證可以看到，這里的轉(zhuǎn)移概率是-0.097，再加上非編碼狀態(tài)下下一個堿基G的生產(chǎn)概率是-0.523.我們能就可以得到（-0.699）+（-0.097 ）+（-0.523）=-1.24。類似的我們來計算，在這個位點從編碼狀態(tài)到非編碼狀態(tài)的轉(zhuǎn)換，也就是-1.40+（-0.398）+（-0.523）=-2.32.這個值比從非編碼狀態(tài)轉(zhuǎn)移得到的-1.24小，因此不會被保留。（舍去概率小的可能路徑）
類似的，我們可以繼續(xù)完成后續(xù)的迭代，把后面所有的格子都一個個填滿。如下：
接下來，我們來做回溯。首先，選出最終概率值最大的那個值。以它為起點，依次來回溯，那么我們得到的回溯路徑就可以得到最終的結(jié)果。在回溯路徑中，如果下一步有兩種可能，就走向概率大的那一家。
把這一路上走過的NC標記下來，就可以得到最后的結(jié)果：
NNCCCCCCNNN
也就是說，我們把輸入的序列CGAAAAAATCG分為了非編碼區(qū)N和編碼區(qū)C.
由于時間有限，我們的MSGP（The Most Simple Gene Predictor）非常簡單。但它很容易被擴展，只需要你引入更多的狀態(tài)，唯一的限制是，不同的狀態(tài)對應(yīng)的生成概率--在這里也就是堿基的組分--必須存在顯性的差異。這樣，我們才可能由你的觀測序列反推出狀態(tài)來。
比如說，Chris Burge 1996年提出基因預(yù)測算法GenScan針對外顯子，內(nèi)含子以及UTR等設(shè)定了獨立的狀態(tài)，從而大大提高了預(yù)測的準確度，是最成功的基因預(yù)測工具之一。但在基本原理上，它與我們剛剛講的，最簡單的MSFP并沒有區(qū)別。類似的，我們還以可以用類似的方法去做5’剪切位點的預(yù)測等等。
事實上，通過將狀態(tài)和可觀測的符號分離開，隱馬爾可夫模型為生物信息學(xué)的數(shù)據(jù)分析提供了一個有效的概率框架，是當(dāng)代生物信息學(xué)最常用的算法模型之一。
本文基本上是對下面兩篇博文的復(fù)述，對作者表示敬意和謝意。
另外，也參考了吳軍老師的《數(shù)學(xué)之美》一書。
用隱馬爾可夫模型建立預(yù)測模型
一文搞懂HMM（隱馬爾可夫模型）
