HMM 不負責講座

未經原作者同意，請勿任意轉載，謝謝。
--

一、前言

不才斗膽，準備以這篇文章替大家介紹 Hidden Markov Model (HMM)。我會根據下列這本書的內容來改寫：

Biological sequence analysis: Probabilistic models of proteins and nucleic acids
by R. Durbin, S. Eddy, A. Krogh, and G. Mitchison
Cambridge University press

當然我會儘量寫清楚一點。我會假設各位已經具備機率的基礎知識。

在此先釐清一個觀念：當我們說，可以用 HMM 來描述或解決某問題時，不代表這個問題就一定要用 HMM 來解；
也不表示 HMM 就是描述這個問題最好最適當的模型。

此外，各位可能會發現，有些生物問題用 HMM 來解可能效果不錯，但是從生物學的觀點來看，
此 HMM 的行為卻不見得直觀或合理。
原因可能是這個 HMM 已經藉由統計學習的方式將我們未知的各種參數變化一網打盡了。

二、CpG islands

先為大家介紹 CpG islands。

這個東西在學習 HMM 的過程中很重要嗎？其實不是。只不過我們可以用 HMM 來解這個問題，
而此問題本身的定義很簡單，因此我決定在之後的文章中儘量以 CpG islands 為主角來舉例說明。
當然有些情況下我也會舉些與 CpG islands 無關且更簡單的例子來說明 HMM。

CpG islands
在 DNA 的單股上 (DNA 序列為雙股螺旋結構)，相鄰的核苷酸 C 與 G 就稱為一組 CpG。
我們不直接寫 CG 的原因是：通常在雙股的 DNA 序列中，
一股上的 C 會與另一股上的 G 透過氫鍵連結而配對在一起，我們通常將這類配對稱為 CG base pair。
為了避免搞混，所以我們將單股上相鄰的 C 與 G 稱為 CpG。

就人類的基因體來說，因為化學特性的關係，CpG 中的 C 會被甲基化 (methylation)，
而此甲基化過後的 C 常常會異變成 T，因此，CpG 其實數量很少。
但是在很多基因的啟動子 (promoter) 附近，CpG 的數量會明顯地增多 (應該是甲基化反應被抑制)，
並且 C 與 G 的總體數量也會上升。我們就將這種 CpG 分佈頻繁的區段稱為 CpG islands；
通常 CpG islands 的長度約數百至數千個鹼基對 (base pairs)。很顯然地，
如果我們有一套很好的方法可以找出 CpG islands 的話，那麼我們就可能找到啟動子，進而找出基因所在。

關於 CpG islands，通常我們好奇的問題是：第一、給定一小段 DNA 序列，它是不是 CpG islands？
第二、給定一大段 DNA 序列，其內部有沒有 CpG islands？如果有的話，在哪個區段呢？
接下來會陸續為各位提出關於這兩個問題的解決方案
(請注意：所提出的解決方案不見得是唯一且最好的喔)。

三、馬可夫鏈 (Markov Chains)

先來看看第一個問題：給定一小段 DNA 序列 X，它是不是 CpG islands？

如果我們有一個可以產生 CpG islands 序列的數學模型的話，
那麼我們就可以看看這個模型生出 X 的機率有多高。如果這個模型有很高的機率可以產生 X 的話，
我們就可以預測 X 的確是 CpG islands；反之 X 就不是 CpG island。因此，
我們的目標就是造出這樣的數學模型。

喔，那這樣的模型應該有什麼特性呢？我們其實可以換個角度來看，我們關心的就是在序列中，
G 出現在 C 後面的機率高不高？換言之，這個模型所產生出來的序列，
其中 G 接在 C 後面的機率應該是挺高的。
因此，為了簡化起見，我們的模型只要具備以下特性應該就可以達成這個目標了：
某個字元出現的機率唯一取決於它的前一個字元。
這樣的模型可以很簡單地訂出 G 出現在 C 後面的機率，而不會受到其他參數的影響。
有了這樣的模型之後，我們只要計算出 training set 內 G 出現在 C 後面的機率，
然後把這個值設定在模型內就好了。

有這樣的數學模型嗎？有的，就是馬可夫鏈 (Markov Chain)。

基本上，馬可夫鏈是由兩種元素所組成：一是可能的狀態，一是不同狀態間的轉換機率。
我們現在來看看如何用馬可夫鏈來描述「CpG islands 產生器」。
首先，這個馬可夫鏈必須有四個狀態：A、G、C 與 T，分別對應到 DNA 序列上可能出現的四種核苷酸。
其次，它必須有十六個狀態轉換機率，分別對應到 A to A、A to G、A to C、A to T、G to A、
G to G、G to C、G to T、C to A、C to G、C to C、C to T、T to A、T to G、T to C 以及
T to T 的機率。至於這些機率怎麼來的，我們可以針對 training set 中的序列分別來計算，
當然這些序列必須是生物學家所提供的已驗證的 CpG islands 序列 (見圖一與表一)。


圖一 「CpG islands 產生器」的馬可夫鏈 (以狀態轉換圖來表示)


表一 CpG islands 內 (+) 與 CpG islands 外 (-) 的狀態轉換機率


從上表中可以看出，CpG 在 CpG islands 內出現的機率比在 CpG islands 外高很多 (0.274 vs. 0.078)。
有了這個模型之後，我們就可以算 X 出現的機率啦。假設 X 是由 L 個元素所組成的：

X = X(1)X(2)...X(L)

則 X 出現的機率為：

P(X) = P(X(L), X(L-1), ..., X(1)) = P(X(L) | X(L-1), ..., X(1))P(X(L-1) | X(L-2), ..., X(1))...P(X(1))

因為馬可夫鏈中的元素只受前一個元素的影響，所以

P(X(i) | X(i-1), ..., X(1)) = P(X(i) | X(i-1)) = a_X(i-1)X(i)

其中 a_ij
是從狀態 i 轉到狀態 j 的機率。 因此原式變成

P(X) = P(X(1))a_X(1)X(2)a_X(2)X(3)...a_X(L-1)X(L)

換句話說，就是把序列上所有的字元轉換的機率都乘起來就對啦。
第一個字元出現的機率我們可以很簡單地假設成四分之一。例如，CGTAG 出現的機率就是
0.25 x (C to G 的機率) x (G to T 的機率) x (T to A 的機率) x (A to G 的機率)。

當然啦，其機率值會隨著序列長短不同而有變化。因此，為了方便比較，
我們通常會以下式來給定一個相對分數，藉此作為判斷的標準：

S(X) = log(P(X | +) / P(X | -))

也就是 X 出現在 CpG islands 內的機率除以 X 出現在 CpG islands 之外的機率再取對數。
取對數可以將乘法降為加法，在計算上簡單且快速許多。

這樣一來，我們就可以根據 S(X) 來判斷 X 是否為 CpG islands 了。

四、HMM 的動機 (以 CpG islands 問題而言)

讓我們複習一下第二個問題：給定一大段 DNA 序列 X，裡面有哪些區段是 CpG islands？
我們之前看過的馬可夫鏈可不可以解這個問題呢？

可以說是「可以」；也可以說是「不可以」。

我們可以每次看 X 的一小段區段，利用馬可夫鏈來判定這個區段是不是 CpG island，
藉此把所有的 CpG islands 找出來。例如，先看 X(1) 到 X(k)，看看這一段是不是 CpG island、
再看 X(2) 到 X(k+1)，看看這一段是不是 CpG island、...、最後看 X(n-k+1) 到 X(n)，
看看這一段是不是 CpG island。在這裡我們假設 X 的長度是 n，並且 X(i) 指的是 X 的第 i 個元素。
相信大家很快地就可以發現這個方法的限制，就是我們如何決定 k？
其實 CpG islands 的長度不見得固定, 這樣做感覺上效果有限。並且 k 取得太大或太小的話，
都會影響 CpG islands 的判斷準確度。

因此，我們需要另一個數學模型來幫助我們判斷兩件事：

1. CpG islands 可能的所在；
2. CpG islands 的邊界。
   

我們可以利用 Hidden Markov Model (HMM) 來替我們解決這兩件事。
接下來我們會替大家介紹 HMM 的定義，並且舉一個例子 (喔，謝天謝地，這個例子不是 CpG islands 了)。

再次強調，HMM 不見得是這個問題的最佳解或是唯一解。它只是解決 CpG islands 的一個可行方案罷了。

五、HMM 的定義與實例

讓我們先來看個小問題。

假設我們去一個賭場賭骰子。您知道的嘛，賭場哪有不作弊的？根據某消息靈通人士所言，
這個賭場使用兩個骰子，一個是正常的骰子，也就是 1 到 6 出現的機率一樣，都是六分之一；
另一個骰子是動過手腳的，6 出現的機率特別高，佔了二分之一，其餘 1 到 5 出現的機率一樣，
都是十分之一。此外，賭場會偷換這兩個骰子，當他們使用正常骰子時，
下一次偷換成動過手腳的骰子的機率是 0.05；當他們使用有問題的骰子時，
下一次偷換成正常骰子的機率是 0.1。當然囉，我們只能看到每次骰子丟出來的點數是多少，
我們並不知道賭場用的是哪一個骰子。

根據以上的描述，我們可以畫出對應的狀態轉換圖 (圖二)，這個狀態轉換圖與馬可夫鏈轉換圖略有不同：
此處的狀態圖中，一個狀態有好幾種可能的對應符號。而在馬可夫鏈中，一個狀態就只能對應一種符號，
例如在我們先前舉的 CpG islands 例子中，狀態 A 所對應的符號就只能是核苷酸 A。


圖二 「賭場用骰子」的狀態轉換圖


在這個狀態轉換圖中，我們有兩個狀態，一個代表正常的骰子 (F)，另一個代表不正常的骰子 (U)。
我們有四個狀態轉換機率：從正常骰子換到正常骰子的機率是 0.95；
從正常骰子換到不正常骰子的機率是 0.05；從不正常骰子換到不正常骰子的機率是 0.9；
而從不正常骰子換到正常骰子的機率是 0.1。

此外，每個狀態會有所謂的「符號產生機率」。就正常的骰子這個狀態而言，會有六個符號產生機率，
分別代表點數 1 到點數 6 出現的機率 (在此都是六分之一)；就不正常的骰子這個狀態而言，
也會有六個符號產生機率，也分別代表了點數 1 到點數 6 出現的機率 (6 出現的機率是二分之一，
而 1 到 5 出現的機率都是十分之一) 。

像這樣的一個狀態轉換圖，就是一個 Hidden Markov Model (HMM)。

一個 HMM 包含四項元素：可能的符號、可能的狀態、狀態間的轉換機率、以及狀態內的符號產生機率。
以上述的例子來說，可能的符號就是 1 到 6 。可能的狀態就是「正常骰子」與「不正常骰子」。
狀態間的轉換機率是「正常骰子」->「正常骰子」為 0.95、「正常骰子」–>「不正常骰子」為 0.05、
「不正常骰子」->「不正常骰子」為 0.9、「不正常骰子」->「正常骰子」為 0.1。
而狀態內的符號產生機率有兩組，分別是：就「正常骰子」來說，
出現 1、2、3、4、5、6 的機率都是六分之一，而就「不正常骰子」來說，
出現 1、2、3、4、5 的機率都是十分之一，出現 6 的機率是二分之一。

這樣的模型可以產生兩種序列：狀態序列與符號序列。就以上的例子而言，狀態序列就是記錄所使用的骰子。
例如「正常」「正常」「正常」「不正常」「正常」。這個序列表示第一次用的是正常的骰子、
第二次用的也是正常的骰子、第三次用的還是正常的骰子、第四次用的是不正常的骰子，
而最後一次 (也就是第五次) 用的是正常的骰子。而符號序列則是記錄每次骰子擲出的點數。例如 65326，
分別表示從第一次到第五次所擲出的點數。相信各位可以發現，我們可以觀察到的是符號序列，
而我們無法得知的是狀態序列，這也就是為什麼這個模型叫做 Hidden Markov Model，
因為我們的狀態序列是「Hidden (隱藏)」的。

現在我們可以想一下，該如何造一個 HMM 來描述 CpG islands 的第二個問題：
給定一長段 DNA 序列 X，X 中的哪些區段可能是 CpG islands？

我們可以造出如下的 HMM：可能的符號就是四種核苷酸 A、C、G、T。而可能的狀態有 A+、
C+、G+、T+、A-、C-、G-、與 T- 這八個。其中 A+ 表示位於 CpG islands 內產生核苷酸 A 的狀態，
而 A- 表示位於 CpG islands 外產生核苷酸 A 的狀態，以此類推。而狀態轉換機率有 A+ to A+、
A+ to C+、…、A+ to G-、… 共 64 種轉換機率。而狀態內的符號產生機率有八組，其中以狀態 A+ 來說，
產生 A 的機率是 1，而產生 C、G、與 T 的機率都是零。以狀態 G- 來說，產生 G 的機率是 1，
而產生 A、C、與 T 的機率都是零。其餘六組以此類推。在這個模型之下，
狀態序列 C+G-C-G+ 所對應的符號序列就是 CGCG；
而我們知道第一與第四個核苷酸是位於 CpG islands 之內，而第二與第三個核苷酸則位於 CpG islands 之外。
請注意，這只是其中一種可以用來描述這個問題的 HMM，您應該可以輕易造出另一個可以描述此問題的 HMM。

好啦，HMM 的定義講完啦。那麼，如果我們有一個 HMM，還有一堆 training 用的序列組合
(狀態序列與其對應的符號序列)，我們希望完成哪些事呢？可能有以下幾項：

1. 算出此 HMM 的狀態轉換機率與狀態內的符號產生機率 (Baum-Welch 與 Viterbi training)。
2. 給定符號序列，求出相對應的狀態序列 (Viterbi)。就骰子的例子來說，給定觀察到的點數序列，
   求出分別是用哪個骰子所擲出的。
3. 此 HMM 產生某段符號序列 X 的機率是多少 (forward)？
4. 給定某段符號序列 X，X 的第 i 個元素是由狀態 S 所產生的機率是多少 (forward 與 backward)？
   

各位應該發現了，在每個問題之後都有個括弧。括弧內就是解決這些問題的既有方法。
這些既有的方法都有相當好的效果。就以骰子的例子與第二個問題來說，我們可以看看圖三：

Rolls   315116246446644245311321631164152133625144543631656626566666 
Die     FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFUUUUUUUUUUUUUUU
Viterbi FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFUUUUUUUUUUUU

Rolls   651166453132651245636664631636663162326455236266666625151631 
Die     UUUUUUFFFFFFFFFFFFUUUUUUUUUUUUUUUUFFFUUUUUUUUUUUUUUFFFFFFFFF

Viterbi UUUUUUFFFFFFFFFFFFUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUFFFFFFFF

Rolls   222555441666566563564324364131513465146353411126414626253356
Die     FFFFFFFFUUUUUUUUUUUUUFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFUU
Viterbi FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFU


Rolls   366163666466232534413661661163252562462255265252266435353336
Die     UUUUUUUUFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
Viterbi UUUUUUUUUUUUFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

Rolls   233121625364414432335163243633665562466662632666612355245242
Die     FFFFFFFFFFFFFFFFFFFFFFFFFFFUUUUUUUUUUUUUUUUUUUUUUFFFFFFFFFFF
Viterbi FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFUUUUUUUUUUUUUUUUUUUFFFFFFFFFFF

圖三 以 Viterbi 法來猜測使用的骰子

圖三記錄了三百次的骰子結果。Rolls 是擲出的點數；Die 是實際使用的骰子，F 是指正常的骰子，
U 是指不正常的骰子；Viterbi 是我們用 Viterbi 法所猜出的骰子 (在知道 Rolls 的情況下)。
紅色就是我們感興趣的部分：賭場用不正常骰子作弊的部分。我們可以看到，
Viterbi 法其實已經猜得相當準了。HMM 其他問題的演算法也可以得到類似這樣相當不錯的效果。

因此，我們擔心的並不是「操作 HMM 的工具」，而是 HMM 本身。換言之，
當我們想要用 HMM 來解決問題時，狀態轉換圖才是我們的藝術所在，
也是模擬及預測的效果好壞最重要的因素。只要 HMM 狀態轉換圖設計得好，模擬及預測的效果就會很好。
也就是說，之前提到可以模擬 CpG islands 的 HMM 不只一種，
但是本文所介紹的 CpG islands HMM 應該會比其他的 HMM 效果好些。
利用 HMM 來做 Gene finding 的論文也很多，但是 GENESCAN 所造出的 HMM 效果就是比別人好。

這篇介紹 HMM 的文章就要在此告一段落了。我不打算進一步說明 HMM 會用到的各種方法，
因為其中牽涉了一些統計與演算法的理論，可能大家看起來會比較吃力。
並且這些東西其實都已經有免費的程式套件可以使用了，各位也實在不必多花心思在這上面。
我想傳達的理念是，HMM 的效果好壞取決於它的狀態轉換圖，而我們對問題的本質瞭解得越深入，
畫出來的狀態轉換圖就會越完備，效果也會越好。