生物信息学导论-北大-马尔可夫模型

ref: https://www.coursera.org/learn/sheng-wu-xin-xi-xue/home

本文主要来自本课的讲义。

---

Markov Model

markov链描述了一个在连续时间段的离散随机过程，其中从一个状态向其他状态（包括自身）的转换，是由一个概率分布（transition probability）控制的。而且下一个状态完全由当前状态决定。

Hidden Markov Model

token即observable symbols，y(t)，由对应状态x(t)产生，也是由一个概率分布控制的（emission probability）。注意，每个状态都可以产生所有种类的token，但是产生每种token的概率是不一样的，因此，一条token序列可以通过不同的状态转移路径生成，每一种路径的概率应该也是不一样的。

那么给定一个HMM模型，一条token是如何产生的？

• 达到某个状态时，emission概率分布决定了这个token产生的概率
• 要进入下一个状态时，transition概率分布决定了去往哪个状态

先计算某一种状态转移路径产生的概率。假设

Transition Probability是：（$a_{kl}$$x_i$表示一个状态，$S_k{S}_l$都是状态值）

$$a_{kl}=P(x_t=S_l|x_{t?1}=S_k)$$

Emission probability是：（$e_k(b)$表示状态$S_k$下产生b的概率）

$$e_k(b)=P(y_i=b|x_i=S_k)$$

那么：（其实就是将每个位置的ae乘积再作连乘）

$$P(X,Y)=\prod _{i=1}^L(e_{x_i}(y_i)?a_{x_i{x}_{i+1}})$$

应用到序列比对问题上，两条序列的每一种比对其实都可以看作一个状态转换路径（状态M表示匹配，X或Y表示匹配了gap），那么每种比对都有概率，只要找到概率最大的就可以了，也就是求解：（ali表示一种比对，或者一个状态转移路径，S1和S2就是进行比对的序列）

$$argma{x}_{ali}(P(S1,S2,ali))$$
代入概率，如下图所示：

假设$P(X,Y,ali)$表示经过当前的状态转换路径（ali）的概率（token是核苷酸或者残基，state就是M,X,Y，表示匹配或者gap），因为对于给定的一个token序列，可能有多条路径产出，所以这条token出现的概率就是把各种状态路径的概率相加，如下面的公式：

$$P(X,Y)=\sum _{ali}P(X,Y,ali)$$

因为要求和，所以

$$P(X,Y)=P_M(n,m)+P_X(n,m)+P_Y(n,m)$$
下图是推导过程。其实就是把求最大值贬成求和。

The Most Simple Gene Predictor(MSGP)

问题：给定一条序列，如何预测编码区和非编码区？

放到HMM中，核苷酸序列就是token序列，每个核苷酸可能是n状态（非编码区），也可能是c状态（编码区），那么需要获取transition概率和emission概率，再通过上面的计算获取状态序列。

虽然我们不知道真实的transition概率和emission概率，但可以使用已标记好编码区和非编码区的序列进行估算。

计算的过程类似上面的序列比对，因为此时只有两种状态，所以更简单一点。

为了计算时更简单，还可以将概率x转化为$lo{g}_{10}(x)$，因为计算时有很多连乘，而

$$log(a?b)=log(a)+log(b)$$

例子

给定序列：CGAAAAAATCG

transition概率：（经过$lo{g}_{10}$处理了）

	n	c
n	-0.097	-0.699
c	-0.398	-0.222

emission概率：（经过$lo{g}_{10}$处理了）

	A	C	G	T
n	-0.699	-0.523	-0.523	-0.699
c	-0.398	-0.699	-0.699	-0.699

从左往右填表，到结尾后，发现得分最高的是-7.774，所以从这里往左回溯（绿色箭头），可知状态转移路径是 NNCCCCCCNNN 。这样就完成了最简单的预测（MSGP）。

当然，如果对状态做细分，比如把编码区再拆成UTR、intron、exon等，就可以更好地预测序列了。（举例：GenScan）