HMM简介

HMM是一种概率图模型

image.png

即有马尔可夫链这个随机过程,在马尔科夫链对应状态点上,又有相应的观测点,状态点之间转移满足马尔可夫链的转移矩阵,状态点与观测点之间的传递通过发射矩阵

例子

假设说我有多条序列:

image.png

对于每一个postion来说,可能出现的的情况如下:

image.png

特别的,第4个位置:

第4个位置

image.png

引入HMM可以看到:

image.png

解释下这幅图:
箭头表示转移方向,上面的数字表示走这条路径的概率,框里的柱形图表示先验的碱基出现的个数

这个随机过程可以理解为:第一个位置往第二个位置转移(计算碱基出现的频率)

以第三个位置为例:
A : 0.8, C : 0.2
那么第四个位置: 从第三个位置向上转移的概率为0.6;向右侧转移的概率为0.4

假设我要计算ACACATC这条序列出现的概率:

image.png

其中(0.8;0.8;0.8;0.4;1;0.8;0.8为碱基出现概率)
其中(1;1;0.6;0.6;1;1为转移概率)

当然也有另一个算法:

image.png

P(S)是碱基发生概率,L是sequence长度

image.png

image.png

例子

当然呢,我也有看到大佬这么写的,利用HMM预测非编码序列

image.png

两个颜色不同的骰子，一个是橘色（Coding，C）的另一个是蓝色（Noncoding,N）的
黑色部分即是马尔可夫链隐含部分---隐含层,蓝色部分是可见部分(即你看到的部分)---状态层

事实上对于编码序列和非编码序列,我们选择概率为1/2,然后ATCG四个序列每个碱基选择概率为1/4

那么我们分开两层:隐含层用于表示是编码还是非编码序列,C C N N N N N N N C C C(C为coding ,N 为nocoding)
状态层用于表示是何种类型碱基,CGAAAAAATCG

利用马尔可夫模型预测基因

我们来做一个最简单的基因预测。给定一段基因组DNA序列，我们来预测其中的编码区。
按照我们之前说的,隐含层用于判断是否为编码/非编码;状态层用于表示观察到底是什么序列

image.png

根据之前分析的,我们在隐含层的转移用转移矩阵(图1),从隐含层到对应状态层的转移用发射矩阵(图2)

图1

图2

根据一些先验信息和统计数据来制定相应的概率值

image.png

image.png

假设我有一段基因组序列: CGAAAAAATCG
非编码状态N的分布设为0.8,编码状态C的分布设为0.2，经过log10转换后，分别为-0.097和-0.699

image.png

我们先假设一种情形:
我们说过HMM模型最关键的两个矩阵，一个叫状态转移矩阵，另一个叫发射矩阵
那么我们将编码区和非编码区看成两个状态，状态的转换通过状态转移矩阵实现；每个状态到观察到的序列通过发射矩阵连接

我们先看状态转移矩阵：

这个矩阵表示对于非编码区，当前非编码状态向下一个非编码状态转移的概率为-0.097，当前非编码状态向下一个编码状态转移的概率为-0.699，当前编码状态向下一个非编码状态转移的概率为-0.398，当前编码状态向下一个编码状态转移的概率为-0.222

对于发射矩阵：

发射矩阵分两种情况：1.若当前状态为非编码状态时，观察到A的概率为-0.699；观察到C的概率为-0.523；观察到G的概率为-0.523；观察到T的概率为-0.699
2.若当前状态为编码状态时，观察到A的概率为-0.398；观察到C的概率为-0.699；观察到G的概率为-0.699；观察到T的概率为-0.699。

假如只考虑状态转移矩阵：
假设说第一个位置C,来自的是非编码序列,那么该事件发生概率是-0.523,那么它向下一个状态:

1. 非编码区域 : -0.523 + (-0.097) = -0.62

2.编码区: -0.523+（-0.699）= -1.222
显然我们选择最大值,-0.62,即下一个任然是非编码区,且出现G的概率为-0.523

如果既考虑状态转移矩阵，又考虑发射矩阵：
1.第一个序列为非编码取的C碱基,第二个位置为非编码的G碱基的总概率和为: -0.523 + (-0.097) + (-0.523) = -1.143

2.第一个序列为非编码取的C碱基,第二个位置为编码的G碱基的总概率和为: -0.523 + (-0.699) + (-0.699) = -1.921
显然选择第一种可能（值最大），接下来的序列就根据这个原理往下推

那么假设第二种情形也是一样的:
假设说第一个位置C,来自的是编码序列,那么该事件发生概率是 -0.699,仿照第一种情形往下推,选取概率最大值
原理:

原理

然后回溯回去,找每个状态的概率最大值回溯

image.png

然后就可以寻找到最终的结果了:

image.png

参考:
<Hidden Markov Models for Biological Sequences>
https://www.jianshu.com/p/3b367bc14147