隐马尔科夫链模型

1、第三节 后验解码问题 前面我们介绍了在给定一个符号序列时，在不同的可能状态序列中找出概率（可能性）最大的状态序列（或路径）。在实际情况中可能会碰到一个问题：不同的路径其概率相差不大，这为我们选择哪一条路径增加了一定的困难。但我们知道一条状态序列，它是由各个位置的状态组成的，比如图6-11中，从第1个位置到第45个位置均由状态F组成，紧接后面的21个是L状态，依此类推。因此克服上述困难我们可以每个位置各个不同状态的概率，这里面我们必须分清一点：某一条序列对应于某个状态序列的概率含义与状态序列中某个位置上具体某个状态的概率是不同的，我们将它们的区分总结于如图6-14：图6-14 某个DNA序列对应

2、的不同可能的CpG岛分布（状态序列）的概率与各个位置上某个状态的概率之间的区别（“+”代表CpG岛区域；“-”代表非CpG岛区域）。我们在上节的“Viterbi”算法中重点介绍的是如何计算图6-14中右边部分，并将其中最大的概率对应的路径找出来。而我们这一节要介绍的是如何计算图6-14中最下面的概率，这就是通常我们所说的“后验解码问题”：计算概率为多少？为了计算这个概率，首先得介绍所给定序列对应的概率，由于不同的状态序列可以产生同一条序列，因此，我们有图6-15的描述：图6-15 隐马尔克夫链中符号序列的概率与路径的关系（这里假设符号序列集只有两条序列与，对应的状态序列（路径）中也只有两个路径

3、与。 图6-15中的与就是要求的某个符号序列的概率，由于此时路径与已知，因此它与符号序列的联合概率表达式可写为：（6-11）公式（6-11）罗列了一大堆符号，这对学生物学的人来说要比较快地接受它们可能会要多花点时间，为使我们有一个清晰的直观印象，我们这里以两条DNA序列及相应的CpG岛分布图来说明。图6-16 DNA序列中CpG岛分布的隐马尔科夫链模型中在公式（6-6）中的示意图如果我们以代表中的某个序列，相应的状态序列有，则有： (6-12)在实际情况中，可能的路径（状态序列）很多，而且随着符号序列的长度的增加，姿态态序列的个数N呈指数函数的方式往上增长，以CpG岛为例：当序列长度为3时，其

4、可能的状态序列个数为，当序列长度为4时，则有。因此如何计算（6-11）则需要一定的技巧。借鉴Viterbi算法，我们也不妨来个沿着符号序列逐步计算，这种逐步的方式有两种：一种是从符号序列的左端向右端（DNA序列的5端到3端，蛋白质序列的N端到C端），相应的算法叫前向算法（forward algorithm）；另一种是从右端向左端，相应的算法为后向算法(backward algorithm)。我们首先介绍前向算法。一、前向算法(Forward Algorithm) 在介绍前向算法之前，我们先看一下式（6-7）的第一个公式即： （6-13）从中我们可发现它取的是每个位置最大值概率，因为对符号序列，

5、其每个位置的概率等其各个状态概率之和，因此，整个符号序列的概率是每个位置概率之积，沿着符号序列的方向，我们有： (6-14)整个算法可写成如下：初始化： (6-15a)迭代过程： (6-15b)终止： (6-15c)我们可将前向算法以图示方式表示：图6-17 前向算法的基本框架图图6-18 前向算法的上体计算图。这里我们仅选状态1的计算过程，目的是使图形显得简洁明了，因为其它状的算法与此是相同的。计算例子：符号序列“3151”。我们仍以Viterbi算法中的例子即符号序列“3151”来说明：计算第一个点：初始条件如Viterbi算法所示，同时设定。首先，我们计算第一个位置即，我们有：当时，我们

6、有：当，我们有：当时，我们有：在终止阶段，我们有：与Viterbi算法相似，前向算法也存着“因计算机浮点数的限制导致后面的结果失真”，因此同样的需要用对数运算来代替相应的乘法运算，有文献曾试图通过如下的转换：克服这一点。但我们认为这样做其实没有实质性的意义，因为中间一个“exp”的运算，如果log(p)小到一定程度，exp(log(p)仍旧受到浮点数的影响。因此，我们并不认为这是一个好的处理方式，甚至认为这种做法除了增加运算量外，没有其它的实际意义。有人曾引进一个标度因子，我们认为这种方法相比较理想，而且，我们也自编了相应的程序，发现它的确可以克服计算机浮点数有限的限制。为此，我们这里详细介绍

7、这个方法。该方法的基本原理是将每一步（或每一时刻）进行归一化计算，具体为：每一步的前向概率进行如下转换： (6-16a)即有 (6-16b)通过这样的转换，我们可以得到式（6-15b）的迭代公式为： (6-17)通过这样的转换，要求转换后的符合如下条件： (6-18)根据公式（6-18），我们可以得到式（6-17）中的归一化因子的表达式为： (6-19)应用这种方法计算前向概率与后向概率显然可以在基本上克服计算机浮点数带来的不足，因为计算得到的其值肯定在状态数倒数值附近浮动，与序列长度相比，值显然要小得多，如在蛋白质二级结构中有，在CpG岛中有。应用这种校正方法有一个很显著的优点就是计算序列的

8、概率值即，它与归一化因子的关系为： (6-20a)因此很自然地就可以将它转化为相应的对数即： (6-20b)这里提前说明一下：有关文献在介绍这种归一化因子法是，将我们刚才前面介绍的前向算法的归一化因子也用在我们接下来要介绍的后向算法中。我们通过深入的研究发现：在后向算法中直接应用前向算法中的归一化因子不是很妥当，在计算后验解码问题中“计算机浮点数有限”的问题仍未得到有效的解决。为此，我们在后向算法中重新计算归一化因子，具体地将在后向算法中介绍。虽然这种改变是微小的，但它在后验解码问题及后面我们要介绍的HMM模型的概率参数即的求解算法之一的Bauch-Welch算法中会发现其具有非常独特的优越性

9、。二、后向算法（Backward Algorithm）后向算法基本上与前向算法基本类似，只不过计算的方向不同：从符号序列的右边向左边运算（对DNA序列是从3端到5端；对蛋白质序列则是从C端到N端）。下面我们直接将其算法叙述如下：初始化： 对所有状态，有 (6-20a)迭代过程： (6-20b)终止： (6-20c)计算例子：符号序列“3151”我们仍以符号序列“3151”来说明后向算法的计算过程：计算第一个点：设定，根据初始条件的设置我们有：首先，我们计算倒数第二个位置即，此时，我们有： 其次，计算倒数第三个位置即时，我们有： 最后，我们计算倒数第四个位置也就是第一个位置即时，我们有：这样我们

10、计算得到的有： 我们将后向算法的结果与前向算法的结果相比，它们基本一致，其差异来自于计算过程中的舍入误差。 我们在前一小节中提到“如果将前向算法的归一化因子用于后向算法，则在后验解码问题及概率参数的求解仍未摆脱计算机浮点数有限的束缚”，并说明可以重新对后向算法各迭代步骤中的概率进行归一化，具体的方法很简单，其归一化因子为： (6-21)这里的代表后向算法中的归一化因子，转化后的后向概率即与其实际的概率即的关系为： (6-22)这样我们就有： (6-23)三 后验解码问题的算法本节一开头我们就提到过：计算后验解码问题需要前向算法与后向算法。现在我们已经介绍完它们，自然地，接下来就要计算后验解码问

11、题即计算： 根据条件概率公式： (6-24)我们有： (6-25)而上式右边的分子则为：于是式（6-20）可转化为： (6-26)当然可以用前向算法或后向算法得到。从式（6-25）我们可以看出的具体意义就是对某个符号序列，当第个位置的符号为状态时的概率。在CpG岛问题中，如果为“+”，则说明这个位置的碱基在CpG岛区域内；如果为“-”，则说明这个位置的碱基不在CpG岛区域内。现在我们根据式（6-26）来说明针对后向算法我们自己提出的归一化因子在后验解码中的意义。将实际的前向概率与后向概率即与分别应用归一化后的概率即与来表示，即将式(6-16b)，（6-20a）及（6-23）代入式（6-26），

12、我们有： （6-27）由于一般会在“1”附近，因此上式中它们的乘积部分显然不会太小或也不会太大，即可以将它们的乘积控制在计算机浮点数允许的范围内。自然就会摆脱了计算机浮点数有限的束缚。按老办法，我们现在仍以前面的例子即符号序列“来说明“后验解码”的计算方式。后验解码的计算例子显然我们有，根据我们前面计算得到的一系列与，我们很容易就可以计算出各个位置中各个状态的概率即。当时，我们有：当时，我们有：当时，我们有：当时，我们有：以上我们仅选了前4个点的计算，为了使读者对全部数据序列有一个比较直观的认识，我们应用C语言编制了一个计算机程序，将应用前向算法和后向算法计算得到的各位置上的概率即列于图6-1

13、5，同时将也将各个位置上对应的各状态概率即总结在图6-16中。同时，应用前向算法与后向算法计算得到该序列的概率均为：。 读者如果对隐马尔科夫链方法感兴趣，可以根据我们提供的例子手算几个位置，这样可以加深对隐马尔科法的印象。而且，正如我们前面所说的，一旦我们熟悉了具体的算法，则可以很“自然”将它们应用到我们熟悉的其的领域中去，显然对扩大隐马尔科夫链方法的应用以及生物信息学新方法的建立是很有帮助的。此外，根据图6-20，我们显然也可以据此判断每个位置是“F态”还是“L态”。具体的方法是如果位置的某个状态概率大于另一个状态概率，则可以确定张三用的是这个状态所代表的骰子，具体为： 如果，则该位置的状态

14、为F，反之则为L我们将上面的结果图示如下： -3.00 -5.32 -7.58 -9.76 -11.86 -12.27 -14.57 -16.81 -17.35 -19.68 -21.96 -22.55 -23.29 -25.67 -28.03 -30.35 -32.60 -34.79 -36.88 -38.90 -40.85 -42.76 -44.64 -46.50 -46.74 -49.00 -51.20 -53.30 -53.72 -56.02 -58.26 -60.41 -62.49 -64.48 -66.42 -68.32 -68.58 -70.85 -73.05 -75.16 -7

15、7.19 -79.16 -81.07 -82.95 -84.82 -85.06 -87.32 -89.52 -90.02 -92.33 -2.48 -4.27 -6.07 -7.89 -9.72 -11.55 -13.34 -15.16 -16.98 -18.75 -20.56 -22.37 -24.13 -25.76 -27.49 -29.27 -31.08 -32.90 -34.73 -36.56 -38.39 -40.23 -42.06 -43.90 -45.73 -47.54 -49.36 -51.18 -53.01 -54.81 -56.62 -58.44 -60.27 -62.10

16、 -63.93 -65.77 -67.60 -69.41 -71.23 -73.05 -74.88 -76.72 -78.55 -80.38 -82.22 -84.05 -85.86 -87.68 -89.50 -91.29-516.29 -514.19 -512.00 -509.75 -507.43 -506.68 -504.42 -502.10 -501.35 -499.08 -496.76 -496.01 -495.36 -493.52 -491.68 -489.83 -487.97 -486.10 -484.20 -482.27 -480.27 -478.20 -476.04 -473

17、.80 -473.11 -471.08 -468.97 -466.77 -466.11 -464.25 -462.36 -460.44 -458.47 -456.43 -454.31 -452.10 -451.43 -449.52 -447.55 -445.52 -443.41 -441.21 -438.94 -436.61 -434.25 -433.47 -431.11 -428.73 -427.94 -425.55-515.35 -513.52 -511.70 -509.89 -508.11 -506.38 -504.57 -502.79 -501.06 -499.25 -497.47 -

18、495.74 -493.93 -492.10 -490.26 -488.43 -486.59 -484.76 -482.92 -481.09 -479.26 -477.43 -475.61 -473.80 -472.00 -470.17 -468.35 -466.53 -464.72 -462.89 -461.05 -459.22 -457.39 -455.56 -453.73 -451.92 -450.11 -448.28 -446.45 -444.62 -442.79 -440.97 -439.17 -437.39 -435.67 -434.08 -432.34 -430.71 -429.

19、37 -427.78图6-19 前50个位置应用前向算法和后向算法的概率：显然，如果将换成，换成则“掷骰子问题”即刻变成蛋白质二级结构预测问题了。结合前面最大路径的求解，我们马上可以得到蛋白质二级结构预测的两个方法。足见将数学方法彻底搞懂的好处。图6-20 300个位置上各种状态（F与L）的概率四、后验解码问题的应用后验解码问题的应用之一就是我们前面刚提到的选取某个位置概率最大的状态，这些状态同样也可以组成一条路径即： （6-28）我们将式（6-18）与最大几率状态序列（路径）相比： （6-28）可以发现：后验解码问题将各个位置的状态概率分开来考虑，其路径是由每个位置概率最大的状态组成，但由此路径与符