欢迎回到北京大学生物信息学:导论与方法网上课程。
我是北京大学生物信息中心高歌,
我们来继续课程,Markov Model
在上两个单元里,我们引入了Markov Chain,并介绍了隐马尔科夫模型在序列比对中的应用
在这个单元里,我们会进一步介绍隐马尔科夫模型在预测以及特征识别中的应用。
我们之前介绍了隐马尔科夫模型中状态与观测相分离的特点,这使得同一个观察到的符号序列,可能来自多个不同的状态路径。
每个状态都有自己的生成概率分布,按照不同的概率产生一组可以被观测到的符号,从而使从观测到的符号序列去反推其对应的状态路径成为可能。
具体来说,我们可以对每个可能的状态路径,计算其产生观测符号串的可能性,而其中概率最大的那个状态路径,也就是最可能产生这个字串的路径
接下来,我们用一个例子来具体说明一下
我们来做一个最简单的基因预测:
给定一段基因组DNA序列,我们来预测其中的编码区。
按照上节讲过的隐马尔科夫模型,我们先要区分不能直接观测的隐状态和可以直接观测的显符号。
在这个例子里,可以很容易的看出,给定的基因组DNA序列是可以观测到的符号串。
而编码/非编码则是不能直接观测的隐状态。
因此,我们可以画出状态转换图。
首先,我们有编码和非编码两个状态。
因为基因组会同时包含编码和非编码区域,因此这两个状态之间可以转换。
当然,每个状态也可以转换到自己,表示连续的编码或非编码区域
这样,我们就有了一个2x2的转移矩阵,也就是alpha k,l
接下来,我们需要写生成概率,也就是e,k,b。
这个也很直观,无论是编码或非编码状态,都有可能产生A, C, G, T四个碱基,因此我们又可以分别有两个矩阵。
接下来,我们需要一个训练集(Training set),来把这三个矩阵的格子中填上具体的数值。
具体来说,我们需要一个已经事先注释过的——也就是说正确标记了编码、非编码区域的DNA序列,
通常这个序列要比较长,以便有充足的数据供统计
好,假定我们经过对训练集的分析,分别填好了转移概率矩阵和生成概率矩阵。
我们需要根据这些数据,来对一个未知的给定基因组序列反推出最可能的状态路径,也就是概率最大的那个状态路径。
因此,我们还是和之前一样利用动态规划算法,写出迭代公式
以及最后的终止点公式(Termination equation)
好,那么现在数据和公式都已经写好了,让我们开始算吧!
别着急,从公式中可以看到,我们需要做大量的乘法。
这不但比较慢,而且利用计算机操作时,随着连乘次数的增加,很容易数值过小而出现下溢(underflow)的问题。
因此,我们通常会引入对数计算,从而将乘法转换成加法。
具体来说,也就是对转移和生成概率都预先取log10。
好,接下来让我们用一个测序的序列正式开始!
首先,让我们和之前一样,画出动态规划的迭代矩阵,其中包含两个状态:
非编码状态n与编码状态c。
接下来,我们需要设定边界条件(boundary condition),也就是这两个状态默认的分布比例,
为了计算方便,我们分别设成0.8和0.2 经过log10转换后,分别为-0.097和-0.699
接下来,我们开始逐步填格子、迭代。
首先,我们碰到了第一个符号:C,
那么,生成概率矩阵可以知道,它在非编码状态下的log10生成概率是-0.523
将之与-0.097相加,就可以得到-0.62
类似的,对于编码状态下,这个数是-0.699 + (-0.699) = -1.40
接下来,要前进一个碱基,就需要进行状态的转移,我们先来看第一种情形,也就是从非编码状态到非编码状态的转换,
从转移矩阵可以看到,这里的转移概率是-0.097,再加上非编码状态下下一个碱基G的生成概率-0.523,
我们就可以得到-1.24
类似的,我们来计算在这个位点从编码状态到非编码状态的转换,也就是-1.40 + (-0.398) + (-0.523) = -2.32
值得注意的是,这个值比从非编码状态转移来得到的-1.24小,因此不会被保留
类似的,我们可以继续完成后续的迭代,把所有的格子都一个个填满
好
接下来,我们来做回溯。首先,选出最终概率值最大的那个值。
以它为起点,依次来回溯……
好,那么我们拿到回溯路径,就可以得到最终的结果。
也就是说,我们把输入序列CGAAAAAATCG分成了非编码区N和编码区C
由于时间有限,我们的MSGP非常简单,
但它可以被很容易的被扩展,只需要你引入更多的状态。
唯一的限制是,不同的状态对应的生成概率——在这里也就是碱基的组分——必须存在显著的差异,
这样我们才可能由你的观测序列反推出状态来。
例如,Chris Burge 1996年提出的基因预测算法GenScan针对外显子、内含子以及UTR等都设定了独立的状态,
从而大大提高了预测的准确度,是最成功的基因预测工具之一。
但在基本原理上,它与我们刚刚讲过、最简单的MSGP并没有区别。
类似的,你还可以用类似的方法去做5’剪切位点的预测等等……
事实上,通过将状态与可观测的符号分离开,
隐马尔科夫模型为生物信息学的数据分析提供了一个有效的概率框架,是当代生物信息学研究最常用的算法模型之一。
有着非常广泛的应用。
我们也会通过补充材料的形式,为感兴趣的同学提供进一步的资料。