HMM簡(jiǎn)介
HMM是一種概率圖模型
即有馬爾可夫鏈這個(gè)隨機(jī)過(guò)程,在馬爾科夫鏈對(duì)應(yīng)狀態(tài)點(diǎn)上,又有相應(yīng)的觀測(cè)點(diǎn),狀態(tài)點(diǎn)之間轉(zhuǎn)移滿足馬爾可夫鏈的轉(zhuǎn)移矩陣,狀態(tài)點(diǎn)與觀測(cè)點(diǎn)之間的傳遞通過(guò)發(fā)射矩陣
例子
假設(shè)說(shuō)我有多條序列:
對(duì)于每一個(gè)postion來(lái)說(shuō),可能出現(xiàn)的的情況如下:
這個(gè)是什么意思呢?
例如[AT],表示第一個(gè)位置可能為A堿基也可能為T堿基
特別的,第4個(gè)位置:
熟悉動(dòng)態(tài)規(guī)劃NW算法的同學(xué)都知道,這個(gè)是gap penalty導(dǎo)致的加"-"
所有第4個(gè)位置有4個(gè)堿基選擇:
分別是A,T,C,G
引入HMM可以看到:
解釋下這幅圖:
箭頭表示轉(zhuǎn)移方向,上面的數(shù)字表示走這條路徑的概率,框里的柱形圖表示先驗(yàn)的堿基出現(xiàn)的個(gè)數(shù)
這個(gè)隨機(jī)過(guò)程可以理解為:第一個(gè)位置往第二個(gè)位置轉(zhuǎn)移(計(jì)算堿基出現(xiàn)的頻率)
以第三個(gè)位置為例:
A : 0.8, C : 0.2
那么第四個(gè)位置: 從第三個(gè)位置向上轉(zhuǎn)移的概率為0.6;向右側(cè)轉(zhuǎn)移的概率為0.4
假設(shè)我要計(jì)算ACACATC這條序列出現(xiàn)的概率:
其中(0.8;0.8;0.8;0.4;1;0.8;0.8為堿基出現(xiàn)概率)
其中(1;1;0.6;0.6;1;1為轉(zhuǎn)移概率)
當(dāng)然也有另一個(gè)算法:
P(S)是堿基發(fā)生概率,L是sequence長(zhǎng)度
例子
當(dāng)然呢,我也有看到大佬這么寫的,利用HMM預(yù)測(cè)非編碼序列
兩個(gè)顏色不同的骰子，一個(gè)是橘色（Coding，C）的另一個(gè)是藍(lán)色（Noncoding,N）的
黑色部分即是馬爾可夫鏈隱含部分---隱含層,藍(lán)色部分是可見(jiàn)部分(即你看到的部分)---狀態(tài)層
事實(shí)上對(duì)于編碼序列和非編碼序列,我們選擇概率為1/2,然后ATCG四個(gè)序列每個(gè)堿基選擇概率為1/4
那么我們分開(kāi)兩層:隱含層用于表示是編碼還是非編碼序列,C C N N N N N N N C C C(C為coding ,N 為nocoding)
狀態(tài)層用于表示是何種類型堿基,CGAAAAAATCG
利用馬爾可夫模型預(yù)測(cè)基因
我們來(lái)做一個(gè)最簡(jiǎn)單的基因預(yù)測(cè)。給定一段基因組DNA序列，我們來(lái)預(yù)測(cè)其中的編碼區(qū)。
按照我們之前說(shuō)的,隱含層用于判斷是否為編碼/非編碼;狀態(tài)層用于表示觀察到底是什么序列
根據(jù)之前分析的,我們?cè)陔[含層的轉(zhuǎn)移用轉(zhuǎn)移矩陣(圖1),從隱含層到對(duì)應(yīng)狀態(tài)層的轉(zhuǎn)移用發(fā)射矩陣(圖2)
根據(jù)一些先驗(yàn)信息和統(tǒng)計(jì)數(shù)據(jù)來(lái)制定相應(yīng)的概率值
為簡(jiǎn)化計(jì)算機(jī)計(jì)算,取log,轉(zhuǎn)化為加法運(yùn)算
假設(shè)我有一段基因組序列: CGAAAAAATCG
非編碼狀態(tài)N的分布設(shè)為0.8,編碼狀態(tài)C的分布設(shè)為0.2，經(jīng)過(guò)log10轉(zhuǎn)換后，分別為-0.097和-0.699
我們先假設(shè)一種情形:
我們說(shuō)過(guò)HMM模型最關(guān)鍵的兩個(gè)矩陣，一個(gè)叫狀態(tài)轉(zhuǎn)移矩陣，另一個(gè)叫發(fā)射矩陣
那么我們將編碼區(qū)和非編碼區(qū)看成兩個(gè)狀態(tài)，狀態(tài)的轉(zhuǎn)換通過(guò)狀態(tài)轉(zhuǎn)移矩陣實(shí)現(xiàn)；每個(gè)狀態(tài)到觀察到的序列通過(guò)發(fā)射矩陣連接
我們先看狀態(tài)轉(zhuǎn)移矩陣：
這個(gè)矩陣表示對(duì)于非編碼區(qū)，當(dāng)前非編碼狀態(tài)向下一個(gè)非編碼狀態(tài)轉(zhuǎn)移的概率為-0.097，當(dāng)前非編碼狀態(tài)向下一個(gè)編碼狀態(tài)轉(zhuǎn)移的概率為-0.699，當(dāng)前編碼狀態(tài)向下一個(gè)非編碼狀態(tài)轉(zhuǎn)移的概率為-0.398，當(dāng)前編碼狀態(tài)向下一個(gè)編碼狀態(tài)轉(zhuǎn)移的概率為-0.222
對(duì)于發(fā)射矩陣：
發(fā)射矩陣分兩種情況：1.若當(dāng)前狀態(tài)為非編碼狀態(tài)時(shí)，觀察到A的概率為-0.699；觀察到C的概率為-0.523；觀察到G的概率為-0.523；觀察到T的概率為-0.699
2.若當(dāng)前狀態(tài)為編碼狀態(tài)時(shí)，觀察到A的概率為-0.398；觀察到C的概率為-0.699；觀察到G的概率為-0.699；觀察到T的概率為-0.699。
假如只考慮狀態(tài)轉(zhuǎn)移矩陣：
假設(shè)說(shuō)第一個(gè)位置C,來(lái)自的是非編碼序列,那么該事件發(fā)生概率是-0.523,那么它向下一個(gè)狀態(tài):
- 非編碼區(qū)域 : -0.523 + (-0.097) = -0.62
2.編碼區(qū): -0.523+（-0.699）= -1.222
顯然我們選擇最大值,-0.62,即下一個(gè)任然是非編碼區(qū),且出現(xiàn)G的概率為-0.523
如果既考慮狀態(tài)轉(zhuǎn)移矩陣，又考慮發(fā)射矩陣：
1.第一個(gè)序列為非編碼取的C堿基,第二個(gè)位置為非編碼的G堿基的總概率和為: -0.523 + (-0.097) + (-0.523) = -1.143
2.第一個(gè)序列為非編碼取的C堿基,第二個(gè)位置為編碼的G堿基的總概率和為: -0.523 + (-0.699) + (-0.699) = -1.921
顯然選擇第一種可能（值最大），接下來(lái)的序列就根據(jù)這個(gè)原理往下推
那么假設(shè)第二種情形也是一樣的:
假設(shè)說(shuō)第一個(gè)位置C,來(lái)自的是編碼序列,那么該事件發(fā)生概率是 -0.699,仿照第一種情形往下推,選取概率最大值
原理:
然后回溯回去,找每個(gè)狀態(tài)的概率最大值回溯
然后就可以尋找到最終的結(jié)果了:
參考:
<Hidden Markov Models for Biological Sequences>
http://www.itdecent.cn/p/3b367bc14147
