隐马尔可夫模型在序列比对和基因预测中的应用

马尔可夫模型由马尔可夫所建立。此模型中，某个事件的发生与其前n个事件的发生有关(比如根据前几天的天气来预测后面的天气)，并且在状态转移间是存在概率的，这种概率可以根据事件发生情况被推算出来。用公式可以表示为P(x_t | x₁…x_t-1) = P(x_i | x_t-n…x_t-1)，其中n为阶数，一阶则代表每个事件都只与上一个事件的发生有关。

隐马尔可夫模型之所以称为“隐”，是因为在我们可以观测到的符号（在序列比对中即为各种碱基或者氨基酸等）背后还有一种不能被观察到的状态（例如序列比对中的匹配状态和对上空位的状态）。每个状态都有自己的生成概率及状态间的转换概率（匹配与空位间的转换），这些概率可以从观察到的符号序列推算出来并作为其它用途（如基因预测）。

---

1.在序列比对中的应用

由于在马尔可夫模型中能够根据某事件之前的事件预测当前发生的各种事件的概率，将此一阶模型引入序列比对时，引入事件发生概率如图：
其中z表示从匹配状态到引入空位的概率，w表示某个链进行空位延伸的概率，根据这两个概率即可计算剩下各事件的概率及序列比对的概率。如P(XMYM) = P(X→M)×P(M→Y)×P(Y→M) = (1-w)×z×(1-w)，但是此式子只考虑了各种状态间的转换概率，并未考虑各状态出现的概率（即各种状态在基因链中出现的概率），因此若加上生成概率的话，就应该变为P(XMYM) = P_(X)×P_(X→M)×P_(M)×P_(M→Y)×P_(Y)×P_(Y→M)×P_(M)，其中P_(X)、P_(M)、P_(Y)即为各种状态出现的生成概率。因此根据上图还可以将序列比对结果的概率公式列为： $$P_M(i,j)=P(x,y)\times max\{\begin{array}{c}(1-2z)\times P_M(i-1,j-1)\\ (1-w)\times P_X(i-1,j-1)\\ (1-w)\times P_Y(i-1,j-1)\end{array}$$$$P_X(i,j)=Qx\times max\{\begin{array}{c}z\times P_M(i-1,j)\\ w\times P_Y(i-1,j)\end{array}$$$$P_Y(i,j)=Qy\times max\{\begin{array}{c}z\times P_M(i,j-1)\\ w\times P_Y(i,j-1)\end{array}$$
其中P_M(i,j)、P_X(i,j)、P_Y(i,j)分别表示为(x_i,y_j)位置上为匹配，x_i对到空位或者y_j对到空位的整条链的形成概率。Qx、Qy、P(x,y)则分别表示出现这三种状态即x匹配空位、y匹配空位或者匹配的概率，即为这三种状态单独的生成概率。

将隐马尔可夫模型引入序列比对，能够给予序列比对以概率学解释，并且若将所有可能的比对结果的概率相加即将上面的概率公式中的求最大值变为求和，那么还能够得到两序列间的最大相似性。

2.在基因预测中的作用

隐马尔可夫模型在生物信息学中更重要的应用还是在基因预测方面，基因预测即根据生物序列的普遍特点从而对基因的编码区及非编码区进行判断，较常用的基因预测工具有Glimmer、Genescan、Genemark等。

要使用隐马尔可夫模型方法就需要先得到在编码区和非编码区中的各种碱基的生成概率和是否编码的状态间的转换概率。这就需要预先根据已知链中的编码区及非编码区来产生训练集，以足够大的训练集来得到相对准确的各种生成概率和转换概率即可用于基因预测。Genescan的原理就与此模型类似，但基因预测所使用的模型一般不是一阶能处理的，有的基因预测甚至已经达到了35阶。可以想象阶数越高，预测的结果会越准确。因为其更注重整条完整链的协同性和整体性。不管是通过多少阶的模型进行预测，都需要通过计算整条链的状态形成概率从而寻找到概率最大的状态组合，从而得出较为可靠的编码区和非编码区预测。

一阶隐马尔可夫模型预测基因的概率计算公式：

其中P表示链各种状态的形成概率，e_coding及e_noncoding分别表示编码区和非编码区的生成概率，a_k→coding或a_{k→noncoding}则分别表示状态间的转换概率，最终从算到的生成概率中取最大概率作为结果即可得到预测结果[^1]。

题外话：其实自己发的内容别人很多都已经发过相关的了，但是自己做的记录总是会印象更深刻一些，当然也有自己的一点理解在里面。另外下次准备分享一下生信最基本的本地blast的操作。冲冲冲！！！

ps：若有侵权情况或者有错误，请通知我，我立刻修改！

[^1]本文内容主要来自高歌老师的 北京大学生物信息学:学习方法(完整带课件）系列视频。