隐马尔可夫模型（HMM）

引言

1. 马尔可夫模型的基本概念

来对2段氨基酸序列x和y进行残基比对，认为存在3种比对关系的状态：

  1. M：残基能够比对上但不一定相等
  2. X：序列x的残基比对到1个空位，或x上发生了1次插入
  3. Y：序列y的残基比对到1个空位，或y上发生了1次插入

序列比对就是在上述3个状态中不断转换的过程：

  1. \( M(i,j) \) : \( x_i \)比对到\( y_j \)时，序列x从1到\( i \)和序列\( y \)从1到\( j \)最好的比对分数
  2. \( X(i,j) \) : \( x_i \)比对到空位时，序列x从1到\( i \)最好的比对分数
  3. \( Y(i,j) \) : \( y_j \)比对到空位时，序列y从1到\( j \)最好的比对分数

转移（从一个状态到另外一个状态）概率：

$$ a_{kl} = P (X_t=S_l|X_{t-1}=S_k) $$

$$ a_{lk} = P (X_t=S_k|X_{t-1}=S_l) $$

转移矩阵：

设定：

  1. \( \delta \) ：Gap open（d）的概率
  2. \( \epsilon \) ：Gap extension（e）的概率

马尔可夫链

先根据转移概率得到一个转移概率矩阵：

转移矩阵

假设匹配状态是XMMY：

匹配状态

计算匹配状态的概率： $$ P(XMMY)=\alpha_{XM}\alpha_{MM}\alpha_{MY} = (1-\epsilon)(1-2\delta)\delta $$

通过上面的例子，下面介绍马尔可夫模型的一些概念：

  • 马尔可夫过程：

    马尔可夫过程是一类随机过程，该过程的“将来”仅依赖“现在”，而不依赖“过去”，把一个总随机过程看作是状态的不断转移，表达式为： $$ x(t+1)=f(x(t)) $$

  • 马尔可夫链：

    时间和状态都离散的马尔可夫过程。

  • 马尔可夫状态链的状态空间： $$ S = \begin{Bmatrix} S_1,& S_2,& S_3,… \end{Bmatrix}S_i\in R $$

  • 马尔可夫链在时刻\( t \)处于状态\( S_i \)条件下，在时刻\( t+1 \)转移到状态\( S_j \)的转移概率是条件概率： $$ P_{ij}(t, t+1)=P\begin{Bmatrix} x_{t+1}=S_j|x_t=S_i \end{Bmatrix} $$

  • 由于马氏链在时刻 \( t \) 从任何一个状态\( S_i \)出发，到下一时刻 \( t+1 \)，必然转移到\( S_j \)，\( j=1,2,… \)，诸状

    态中的某一个，所以有：

    $$ \forall _{i^{\forall}}\sum_{j}P_{ij}(t,t+1)=1 $$

  • 当\( P_{ij}(t,t+1) \)与\( t \)无关时，为齐次马尔可夫链，通常说到的马尔可夫链都是指齐次的。马尔可夫链除了仅依赖于上一次观测值的一阶马尔可夫链，还可以依赖多个连续观测数据，后者称之为高阶马尔可夫链。

2. 马尔可夫模型的组成

  • 随机序列变量 $$ X = \begin{Bmatrix} x_1,& x_2,& x_3,…,x_t \end{Bmatrix} $$

  • 状态空间 $$ S = \begin{Bmatrix} S_1,& S_2,& S_3,… \end{Bmatrix}\space\space S_i\in R \space\space\space i=1,2,…,n $$

  • 转移概率矩阵 $$ P=\begin{Bmatrix} P_{ij}=P(x_{t+1}=S_j|x_t=S_i) \end{Bmatrix} $$

  • 初始状态向量 $$ \prod =\begin{Bmatrix} \pi=P(x_0=S_i) \end{Bmatrix} , \space\space\space i=1,2,…,n $$

3. 隐马尔可夫模型

隐马尔可夫模型是结构最简单的动态贝叶斯网，可以用五元组来描述分别为：

  • 状态空间\( S \)
  • 状态对应的观测空间\( X \)
  • 状态转移矩阵\( A \)
  • 每个状态下观测事件的概率矩阵\( B \)
  • 初始状态概率分布\( \pi \)

给出这五个参数就能够确定一个隐马尔可夫模型，通常用参数\( \lambda=\begin{Bmatrix}A,B,\pi\end{Bmatrix} \)来指代。

举一个新的例子：基因预测——给定序列预测编码区

隐马尔可夫模型在状态的基础上，增加符号的概念，每个状态可以以不同的概率产生可观测到的符号。在例子中，给定的基因组序列为观测到的符号串，编码和非编码为2种隐状态，即编码与否是未知的，需要通过已知的符号来推测。

示例：马尔可夫链

转移矩阵：

基因组会同时包含编码和非编码区域，自我转移的箭头表示状态的连续。

转移矩阵

转移概率矩阵：

转移概率矩阵

$$ a_{kl} = P (X_t=S_l|X_{t-1}=S_k) $$

生成概率：

生成概率

这里的状态路径无法进行观测，需要根据符号路径来推测状态，引入生成概率，状态\( S_k \)时产生符号\( b \)的概率： $$ e_{k}(b) = P (y_i=b|X_{i}=S_k) $$ 生成概率矩阵（在这里有两个）：

2个生成概率矩阵

训练集：正确标记好了编码和非编码区域的DNA序列。

根据训练集填好上述转移概率矩阵和生成概率矩阵，来对未知的给定基因组序列反推最可能的基因组状态路径。

Logarithmic transformation：引入对数计算将乘法变成加法，对转移/生成矩阵概率取\( log_{10} \)（在计算机运算中，很容易因为连乘的次数的增加，很容音因为数值过小，出现下溢的问题）

测试序列：CGAAAAAATCG

根据训练集得到转移概率矩阵和生成概率矩阵：

取log10的转移概率矩阵

取log10的生成概率矩阵

动态规划进行序列比对：

序列比对

假设n状态和c状态默认的分布比例分别是\( log_{10}(0.8) \)和\( log_{10}(0.2) \)，使用序列比对的方法，找到概率最大的值为起点进行回溯。最终找到的状态链为：NNCCCCCCNNN。

根据示例我们可以得到一个大概的步骤来对隐马尔可夫模型进行应用：

  • 首先要确定具体问题中什么是状态，什么是符号
  • 其次根据状态列出转移概率矩阵，根据状态和符号列出生成概率矩阵
  • 根据训练集填上矩阵的具体数值
  • 根据具体问题使用相应的解决方法(本例子中根据数据来对未知的给定基因组序列反推出最可能的基因组状态路径，迭代使用动态规划进行序列比对)

实际应用中关注隐马尔可夫模型的3个基本问题以及对应的解决方法：

  • 估算问题：

    给定模型\( \lambda=\begin{Bmatrix}A,B,\pi\end{Bmatrix} \)，如何有效计算观测序列\( x=\begin{Bmatrix} x_1,x_2,…,x_n\end{Bmatrix} \)出现的概率？也就是模型和观测序列之间的匹配程度。

    解决：foreward和backward算法

  • 解码问题：

    给定模型\( \lambda=\begin{Bmatrix}A,B,\pi\end{Bmatrix} \)和观测序列\( x=\begin{Bmatrix} x_1,x_2,…,x_n\end{Bmatrix} \)，如何找到和观测序列最匹配的状态序列？也就是找到隐藏的模型状态。

    解决：Viterbi算法

  • 学习问题或训练问题：

    给定观测序列\( x=\begin{Bmatrix} x_1,x_2,…,x_n\end{Bmatrix} \)，如何调整模型参数\( \lambda=\begin{Bmatrix}A,B,\pi\end{Bmatrix} \)，使得该观测序列发生的概率最大？也就是根据训练样本学得最优的参数模型。

    解决：Baum-Welch算法

4. 隐马尔可夫模型判断CNV

以PennCNV为例介绍HMM在实际判断CNV中的应用。PennCNV是一个免费的检测SNP分型阵列芯片的拷贝数变异的工具，目前可以处理Illumina和Affymetrix array的数据来得到信号强度，若使用其他类型的SNP芯片数据和寡核苷酸芯片需要预先处理文件的格式。PennCNV使用隐马尔可夫模型整合多个来源的信息来对单个样本推断CNV。它与基于分割的算法不同，除了单独考虑信号强度外，还考虑了SNP等位基因比例分布等因素。

PennCNV通过基因型来检测拷贝数变异的算法流程（Wang K et al. 2007）

在PennCNV中使用的是一阶HMM，隐状态是人为设定的离散值1，2，3，4，5，6，各自对应的总拷贝数数值和CNV基因型如上表所示。注意这里最大的拷贝数状态设定为拷贝数为4，因为4个和4个以上的拷贝无法进行区分。符号是信号强度值的两种形式：BAF和LRR。结合BAF和LRR可以判断不同的拷贝数以及区分出拷贝中性的LOH（文章中没有具体说怎么结合BAF和LRR数据的，但是通过公示我的理解是BAF和LRR同时发生的联合概率作为HMM的生成概率）。

隐状态，拷贝数数值及其对应的描述（Wang K et al. 2007）

下面介绍一些定义：

  • 原始的信号强度数据需要经过“5步标准化”（http://icom.illumina.com/iom/software.ilmn），每个SNP的X值和Y值分别代表实验得到经过标准化的等位基因A和B的信号强度值，\( R \)为总信号强度： $$ R = X + Y $$

  • \( \theta \)为相对的等位信号强度比率： $$ \theta = \frac{arctan(Y/X)}{\pi/2} $$ 该计算分母为\( \pi/2 \)可以将数值\( \theta \)压缩在-1～1之间。

  • B Allele Frequency （BAF）通常指的是标准化的等位基因B和A的相对信号强度：

    BAF公式（Wang K et al. 2007）

    0代表只检测到了A这个allele对应的荧光信号，分型结果为AA；1代表只检测到了B这个allele对应的荧光信号，分型结果为BB；0.5代表A和B这两个allele的荧光信号强度相等，分型结果为AB。\( \theta_{AA} \),\( \theta_{AB} \),\( \theta_{BB} \)都是根据正常样本得到的值，荧光信号存在一定程度的扰动，因此BAF的取值是在0-1范围波动的值。

  • \( R_{observed} \)是观测值，\( R_{expected} \)是通过算法拟合得到的数值代表了正常样本中的信号强度，每个SNP的\( log \) R 比率（LRR）：

$$ LRR={log_{2}}(R_{observed}/R_{expected}) $$

​ LRR = 0 拷贝数为2；LRR > 0 拷贝数增加；LRR < 0 拷贝数减少

介绍完定义之后，来确定转移概率和生成概率以及对应矩阵。

  • LRR的生成概率：

    是混合的均匀和正态分布模型， $$ P(r|z)=\pi_r+(1-\pi_r)\phi(r;\mu_{r,z},S_{r,z}) $$ \( (\phi·；·) \)是均值为\( \mu_{r,z} \)，标准差为\( S_{r,z} \)的正态分布，均匀分布用于对于芯片中随机的信号波动和可能的基因组错误注释和错误assmbly的建模。

  • BAF的生成概率：

    对于每一个隐状态（除了状态1），有不同可能的基因型也就有不同模式的BAF。

    BAF生成概率（Wang K et al. 2007）

  • 对于染色体X进行特殊处理：

    将染色体X中所有的SNP的LRR值减去1个常数，使得针对女性，平均的LRR值不是0，针对男性，LRR不是单拷贝删除的LRR的参考值，因为正常的男性染色体X的状态就是state2不应该是state3。

  • 隐状态的转移概率：

    即两个相邻SNP发生拷贝数状态改变的概率，一般来讲，临近的SNPs拷贝数状态不太可能发生变化，但是较远的SNPs的拷贝数状态更容易发生改变。因此需要找到转移概率和距离的关系，那么转移概率为：

    $$ p(z_i=l|z_{i-1}=j)=\left\{\begin{matrix} 1-\sum_{k=2}^{b}p_{j,k-1}(1-e^{d_i/D}),if \space l=j \\ p_{j,l-1}(1-e^{d_i/D}),if \space l\neq j \end{matrix}\right. $$
    \( D \)是一个常数，当状态为4时设为100Mb，其他状态为100kb，p是未知参数，用Baum-Welch算法进行估计，得到最优的参数后，再使用Viterbi算法进行隐状态的判断。在PennCNV的分析中，排除了所有包含小于等于2个SNP的CNV，因为这些CNV中可能存在较高的假阳性比例。

介绍算法之前先定义几个概念：

  • \( O = O_1,O_2,…,O_t \) 输出的观察序列符号
  • \( P(O|\lambda) \) 给定模型参数时，输出符号序列\( O \)的概率
  • \( a_{ij} \) 从状态\( S_i \)到状态\( S_j \)的转移概率
  • \( b_j(O_t) \) 在状态\( S_j \)时，输出\( O_t \)的概率
  • \( a_t(j) \) 输出部分符号序号\( O = O_1,O_2,…,O_t \)达到状态\( S_j \)时的前向概率

前向算法：

前向算法

前向算法步骤

Vertibi算法

介绍完定义，明确完我们的五元组，我们再明确一下我们的问题是解码问题，解码问题对应的解法是Viterbi算法，介绍一下该算法具体过程：

Vertibi算法步骤

参考

目 录
comments powered by Disqus