用隐马尔可夫模型做基因预测

Posted 2020-11-29 klausage

什么是隐马尔可夫模型

隐马尔可夫模型（Hidden Markov Model，HMM） 是统计模型，它用来描述一个含有隐含未知参数的马尔可夫过程。其难点是从可观察的参数中确定该过程的隐含参数。然后利用这些参数来作进一步的分析，例如模式识别，特别是我们今天要讲的基因预测。是在被建模的系统被认为是一个马尔可夫过程【一段组装好的序列】与未观测到的（隐藏的）的状态【哪些是编码区哪些不是】的统计马尔可夫模型。

下面用一个简单的例子来阐述：
假设我手里有两个颜色不同的骰子，一个是橘色（Coding，C）的另一个是蓝色（Noncoding,N）的。但是和平常的骰子不同的是，他们稳定下来只要有四种可能，也就是上下是固定的（你可以理解为他们只会平行旋转）。这样每个骰子出现ATCG的概率都是1/4.

 

两条链，在一起

假设我们开始投骰子，我们先从两种颜色选一个，挑到每种骰子的概率都是1/2。然后我们掷骰子，我们得到ATCG中的一个。不停地重复以上过程，我们将会得到一串序列，每个字符都是ATCG中的一个。例如CGAAAAAATCG

这串序列就叫做可见状态链。但是在隐马尔可夫模型中，我们不仅仅有这么一串可见状态链，还有一串隐含状态链。在这个例子里，这串隐含状态链就是你用的骰子的序列。比如，隐含状态链有可能是：C C N N N N N N N C C C。

一般来说，HMM中说到的马尔可夫链其实是指隐含状态链，因为隐含状态（骰子）之间存在转换概率（transition probability）。在我们这个例子里，C的下一个状态是N。C，N的概率都是1/2。这样设定是为了最开始容易说清楚，但是我们其实是可以随意设定转换概率的。比如，我们可以这样定义，N后面不能接两个C，或C 的概率是0.1。这样就是一个新的HMM。

同样的，尽管可见状态之间没有转换概率，但是隐含状态和可见状态之间有一个概率叫做输出概率（emission probability）。就我们的例子来说，Coding（C）产生A的概率是1/4，Noncoding（N）产生A的概率是1/4，当然这些概率是我人为定义的，你可以定义为其它的值。

隐含状态转换关系示意图

其实对于HMM来说，如果提前知道所有隐含状态之间的转换概率和所有隐含状态到所有可见状态之间的生成概率，做模拟是相当容易的。对，我们就捡容易的来。

用隐马尔可夫模型做基因预测

接下来，我们来做一个最简单的基因预测。给定一段基因组DNA序列，我们来预测其中的编码区。按照之前说道的隐马尔可夫模型，我们先要区分不能直接观测到的隐状态和可以直接观测的显符号。

两条链在一起

在这个例子里，我们可以很容易看出，给定的基因组DNA序列是可以观测到的符号串。而编码/非编码是不能直接观测的隐状态。因此，我们可以画出状态转换图。首先我们有编码和非编码两个状态。因为基因组会同时包好编码和非编码区域，因此这两个状态之间可以转换。当然，每个状态也可以转换到自己，表示连续的编码或非编码区域。这样，我们就有了一个2*2的转移矩阵。

转移概率

接下来，我们需要写生成概率。这个也很直观，无论是编码还是非编码状态，都有可能产生A,C,G,T四个碱基，因此我们由可以分别有两个矩阵。

生成概率

现在，我们需要一个训练集（Training set）,来把这三个矩阵的格子中填上具体的数值。具体来说，我们需要一个已经事先注释过得——也就是说正确标记了编码，非编码区域的DNA序列，通常这个序列要比较长，以便有充足的数据统计。

假定我们经过训练集的分析，分别填好了转移矩阵概率和生产概率矩阵。我们需要根据这些数据，来对一个未知的给定基因组序列反推出最可能的状态路径，也就是概率最大的那个状态路径。因此，我们还是和之前一样利用动态规划的算法，写出迭代公式，以及最后的终止点公式（Termination equation）

训练的结果

从公式里面，我们看到，我们需要做大量测乘法。这个不仅比较慢，而且利用计算机操作时，随着连乘次数的增加，很容易数值过小而出现下溢（underflow）的问题。因此，我们通常会引入对数计算，从而将乘法转换成加法。具体来说，就是对转移和生成概率都预先取log10。

取log

好，我们正式开始，假设我组装了一段序列（咦，怎么这么短？为了简单^_）：

CGAAAAAATCG

首先，让我们和之前一样，画出动态规划的迭代矩阵，其中包含两个状态，非编码状态N与编码状态C。接下来，我们需要设定边界条件（boundary condition）,也就是这两个状态默认的分布比例。为了计算方便，我们分别设为0.8和0.2，经过log10转换后，分别为-0.097和-0.699.接下来我们逐步填格子

之后，我们碰到的第一个碱基C,由生成概率可知C在非编码状态下的log10生成的概率是-0.523，将之与-0.097相加，就可以得到-0.62。类似的，在编码状态下，这个数是-0.699+（-0.699）=-1.4。

接下来，我们要前进一个碱基，就需要进行状态转移，我们先来看第一种情形，也就是从非编码状态到非编码状态的转换。从转移举证可以看到，这里的转移概率是-0.097，再加上非编码状态下下一个碱基G的生产概率是-0.523.我们能就可以得到（-0.699）+（-0.097 ）+（-0.523）=-1.24。类似的我们来计算，在这个位点从编码状态到非编码状态的转换，也就是-1.40+（-0.398）+（-0.523）=-2.32.这个值比从非编码状态转移得到的-1.24小，因此不会被保留。（舍去概率小的可能路径）

类似的，我们可以继续完成后续的迭代，把后面所有的格子都一个个填满。如下：

移步换景

接下来，我们来做回溯。首先，选出最终概率值最大的那个值。以它为起点，依次来回溯，那么我们得到的回溯路径就可以得到最终的结果。在回溯路径中，如果下一步有两种可能，就走向概率大的那一家。

回溯路径

把这一路上走过的NC标记下来，就可以得到最后的结果：

NNCCCCCCNNN

也就是说，我们把输入的序列CGAAAAAATCG分为了非编码区N和编码区C.

结果

由于时间有限，我们的MSGP（The Most Simple Gene Predictor）非常简单。但它很容易被扩展，只需要你引入更多的状态，唯一的限制是，不同的状态对应的生成概率--在这里也就是碱基的组分--必须存在显性的差异。这样，我们才可能由你的观测序列反推出状态来。

比如说，Chris Burge 1996年提出基因预测算法GenScan针对外显子，内含子以及UTR等设定了独立的状态，从而大大提高了预测的准确度，是最成功的基因预测工具之一。但在基本原理上，它与我们刚刚讲的，最简单的MSFP并没有区别。类似的，我们还以可以用类似的方法去做5’剪切位点的预测等等。

事实上，通过将状态和可观测的符号分离开，隐马尔可夫模型为生物信息学的数据分析提供了一个有效的概率框架，是当代生物信息学最常用的算法模型之一。

本文基本上是对下面两篇博文的复述，对作者表示敬意和谢意。
另外，也参考了吴军老师的《数学之美》一书。