生物信息中的Markov链

Markov链

Markov链是状态的集合，也可以理解为一条有时序性的状态链，它的一个基本假设是下一状态只与当前状态有关，与过去的状态无关，用数学语言描述便是：
$$P(x_{i}|x_{i-1},x_{i-2},......,x_1)=P(x_i|x_{i-1})~~~~~~（1）$$
$x_{i-1}$到$x_{i}$状态之间转转移的概率定义为：
$$a_{x_{i-1}x_i}$$
对于一条序列$x=(x_i,x_{i-1},......,x_1)$,其发生的概率为：
$$P(x)=P(x_1)\prod_{i=2}^{L}a_{x_{i-1}x_i}~~~~~~(2)$$

利用Markov链对CpG岛识别问题建模

Markov链在生物信息中的一个重要应用就是CpG岛序列的识别问题。人类基因组中只要出现双核苷酸序列CG，核苷酸C一般就会被甲基化修饰，使得这种C转化为T的几率高，通常情况下基因组中CpG双核苷酸序列出现的概率通常小于C和G独立概率出发所得到的期望。而又因为某些重要的生物学原因，在诸如基因启动子或“起始”区附近，这种甲基化过程被抑制，导致这些区域一般会看到比其他地方更多的CpG双核苷，这些区域在生物学上有重要意义，在学术上称作CpG岛，它们的长度通常有几百到几千个碱基。那么问题来了：
给定一段短序列，如何判断它是否属于CpG岛上的序列呢？
首先考虑第一个问题，CpG岛序列的判别问题。我们把短序列$x=(x_i,x_{i-1},......,x_1)$作为一条有时序性的观测序列，然后假设CpG岛序列是由一个（+）模型生成的，非CpG岛序列是由一个（-）模型生成的，用Markov链建立这两个建模。序列的每一时刻的观测都有A，C，G，T四种可能，每个状态都能相互转移，当然也可以原地转移（例如产生AAAAAAA），下图能够清晰地描述这个转移关系：

CpG序列可以由上面的模型生成，非CpG岛序列也能够由上面的模型生成，但是生成CpG岛序列的Markov链（（+）模型）与生成非CpG岛的Markov链（（-）模型）的转移概率$a_{x_{i-1}x_i}$是不一样的。直观地理解，CpG岛区域CG双核苷酸出现的概率会比较高，所以由C与G之间的转移概率会比基因组中其他区域高。实际上也确实如此，下面是从人类基因组序列中提取48个CpG岛并用最大似然估计后的转移概率，其中的最大似然估计方法：
$$a_{st}^{+}=\frac{c_{st}^{+}}{\sum _{l}c_{st}^{+}}$$
估计后（+）模型与（-）模型的转移概率为：

Markov模型的参数确定好之后，就可以用贝叶斯模型去做判别了。我们这里用对数概率比去做贝叶斯判别：
$$S(x)=log\frac{P(x|model+)}{P(x|model-)}~~~~~~（3）$$
它直观理解就是如果$\frac{P(x|model+)}{P(x|model-)}$很大，那可以判定x来自（+）模型，如果$\frac{P(x|model+)}{P(x|model-)}$很小，则判定来自（-）模型。
我们把（2）代入（3）式有：
$$S(x)=log\frac{P(x|model+)}{P(x|model-)}=\sum _{i=1}^{L}log\frac{a_{x_{i-1}x_i}^{+}}{a_{x_{i-1}x_i}^{-}}$$
最后代入模型的转移概率，就能够做判别了。

Reference：Durbin R, Eddy S R, Krogh A, et al. Biological sequence analysis: probabilistic models of proteins and nucleic acids[M]. Cambridge university press, 1998.