第四章隐马尔科夫模型和CpG岛预测

第四章 隐马尔科夫模型和CpG岛预测

隐马尔科夫模型（HMM）

u 

语音处理

l  l  l  l 

Step1:语音分帧：将信号分为多个10-20毫秒长的片段（帧） Step2:矢量量子化：将每帧分配到预先设定的类别 Step3:语音识别：类别标签序列中所蕴含的音素（单词）序列 HMM: 首先用于语音识别领域 物体（如手掌）的自动识别 人体连续动作的识别 序列分析 高通量数据建模：微阵列数据、质谱数据

u 

图像处理

l  l 

u 

生物信息学

l  l 

12-12-6

生物信息学算法

基因组中核苷酸的出现频率

u 

基因组4种核苷酸C、G、A、T的出现频率 大致相等 ~0.25 双核苷酸的出现频率？

l  l 

u 

DNA序列中双核苷酸的频率差异很大 一种特殊的双核苷酸：CpG的出现频率通常低 于1/16

12-12-6

生物信息学算法

CpG岛

u 

CpG中的C非常容易被甲 基化，然后突变成T 但是甲基化作用常常被 一个区域周围的基因所 抑制，这个区域就叫做Cp G岛，也就是CpG出现相 对频繁的区域

u 

12-12-6

生物信息学算法

CpG岛的生物学意义

12-12-6

生物信息学算法

发现CpG岛的问题

u 

CpG岛的两个问题:

l 

l 

给出一条短基因组序列, 我们将如何判断它 是否来自CpG岛? 给出一条长序列, 如果它含有CpG岛, 我们将 如何发现?

12-12-6

生物信息学算法

CpG岛的Markov链

u 

两个核苷酸的分布不独立

l 

l  l 

统计模型中一个符号的概率依赖于它前面的 一个符号 最简单的模型：Markov链 状态转移概率:

ast = P( xi = t | xi −1 = s)

12-12-6

生物信息学算法

CpG岛的Markov链

u 

对于一条观测序列 X = x1 x2 ...xL 其出现的概率

P( X ) = P( x1 , x2 ,..., xL )

= P( xL | x1, x2 ,...xL−1 ) * P( xL−1 | x1, x2 ,...xL−2 ) *...P( x2 | x1 )P( x1 )

u 

由Markov链的性质

P( xi | x1 , x2 ,... xi −1 ) = P( xi | xi −1 ) = a xi−1xi

u 

因此有 P( X ) = P( xL | xL −1 ) P( xL −1 | xL −2 )...P( x2 | x1 ) P( x1 )

L

= P( x1 )∏ a xi−1xi

i =2

12-12-6

生物信息学算法

CpG岛的Markov链

系统中增加状态起始（B ）和结束状态（E）

l 

l 

传统上, Markov 链并不 对序列结尾建模 明确地加入结束状态是 为了对序列的长度分布 建模

P( x1 = s) = aBs

P( E | xL = t ) = atE

12-12-6

生物信息学算法

CpG岛的预测

u 

解决问题的思路：

l  l  l  l 

用已知CpG岛的序列构建Markov链模型（+） 用已知非CpG岛的序列构建Markov链模型（-） (+)和（-）进行比较，判断 实验： 一组人类DNA 序列中提取总共48 个推测 的CpG 岛并建立2个Markov 链模型,其中一个 针对标记为CpG 岛的区域(“+”模型), 而另一 个则针对序列的其余部分(“-”模

型)

生物信息学算法

12-12-6

参数的估计

u 

Markov链参数的最大似然（ML）估计：

+ c + st ast = + c ∑t ' st '

12-12-6

生物信息学算法

基于模型的判别

u 

对数几率比(log-odds ratio)：

P( x | model+) S ( x) = log = P( x | model−)

∑ log a

i =1

L

+ ax i −1 xi − xi−1 xi

L

= ∑ β xi−1xi

i =1

12-12-6

生物信息学算法

CpG岛的问题（二）

u 

如何在没有注释的长序列中找到CpG 岛?

l 

策略1：Markov链模型 使用滑动窗（如长度100bp）计算分值 问题： CpG 岛有明确的边界且长度不一 策略2：为整条序列建立一个整合（+）（-） Markov链的模型

l 

12-12-6

生物信息学算法

构建新的模型

u 

新模型：

（+ ）

（-­‐）

u 

问题：

l  l 

如何从数据得到的模型参数？ 进一步用训练得到的模型检测CpG岛？

生物信息学算法

12-12-6

HMM的正式定义

u 

状态序列：

l 

本身就是一条简单的Markov链，被称为路 径（path） π akl = P(π i = l | π i −1 = k ) 实际不可见 通常一个状态可以按照某种概率分布产生所 有的可能符号 发射(emission)概率：ek (b) = P( xi = b | π i = k )

生物信息学算法

l 

u 

符号序列（观察序列）：

l 

l 

12-12-6

作弊的赌场问题

l 

在一个赌场里, 他们大部分情况下使用公平 的骰子, 但有时换成作弊的骰子

12-12-6

生物信息学算法

计算联合概率

u 

观测序列x和状态序列 π 的联合概率:

L

P( x, π ) = a0π1 ∏ eπ i ( xi )aπ iπ i+1

i =1

u 

由状态序列(C+, G-, C-, G+) 发射序列CGCG 的概率为

12-12-6

生物信息学算法

解码问题

u 

虽然通过观察相应的符号已经无法再判断系统 处于哪个状态, 但是通常我们感兴趣的是潜在 的状态序列 CGCG

(C+, G+, C+, G+) (C+, G-, C+, G-) (C-, G-, C-, G-)

12-12-6

生物信息学算法

HMM中的解码算法

u 

思路：所有可能的状态序列中，选择概 率最大的路径作为预测结果,即

π = arg max P( x, π )

π

u 

问题：可能的状态序列组合非常多，不 可能采用穷尽的方法

l  l 

在语音识别中的术语称为解码(decoding)问题 基于动态规划的Viterbi算法

生物信息学算法

12-12-6

HMM中的解码算法

u  u 

考虑曼哈顿网格模型 每一个状态序列对应于图中的一条路径

比对问题

解码问题

12-12-6

生物信息学算法

HMM中的解码算法

u 

HMM的解码问题转化为DAG中的寻找最长 路径问题

l  l 

路径长度定义为所有边的权重乘积（非求和） 图中每条路径对应于概率P(x|π)

u 

Viterbi算法的复杂度为O(n|Q|2).

12-12-6

生物信息学算法

HMM中的解码算法

u 

解码问题: 边的权重

wi,k,l

（π i = k , i）

L

（π i +1 = l , i + 1）

P( x, π

) = a0π1 ∏ eπ i ( xi )aπ iπ i+1

i =1

wi ,k ,l = eπ i +1 ( xi +1 )aπ iπ i+1 = el ( xi +1 )akl

12-12-6

生物信息学算法

HMM中的解码算法

u 

vk(i): 满足观察序列x1…xi以及πi=k条 件下最长路径对应的概率

v0 (0) = 1

vl (i + 1) = max (vk (i) wi ,k ,l )

k

（ k , i） （l , i + 1）

= max (vk (i)el ( xi +1 )akl )

k

= el ( xi +1 ) max (vk (i)akl )

k

12-12-6

生物信息学算法

HMM中的解码算法

12-12-6

生物信息学算法

HMM中的解码算法

CpG 岛模型和序列CGCG所对应v的结果

12-12-6

生物信息学算法

HMM中的解码算法

u 

作弊的赌场问题

12-12-6

生物信息学算法

HMM中的前后向算法

u 

问题：计算观察序列x的概率

l 

不同状态路径均可以生成同一条序列x, 必须 将所有可能路径的概率相加

P( x) = ∑ P( x, π )

π

l  l  l 

思路1：穷举所有路径 思路2： P( x) ∝ P( x, π * ) 思路3：基于动态规划的前向算法

12-12-6

生物信息学算法

前向算法

u 

前向算法

l 

思路与Viterbi算法相似,只需要将取最大值 的步骤替换为求和

前向因子： f k (i )

= P( x1...xi , π i = k )

l 

递归方程为：

f l (i + 1) = ∑ ( f k (i)wi ,k ,l )

k

（ k , i） （l , i + 1）

= ∑ ( f k (i)el ( xi +1 )akl )

k

= el ( xi +1 )∑ ( f k (i)akl )

k

12-12-6 生物信息学算法

算法的伪代码

12-12-6

生物信息学算法

后项因子的意义

u 

问题：已知观测序列，符号xi来自于状态 k的概率？

P( x, π i = k )

= P( x1...xi , π i = k ) P( xi +1...xL | x1...xi , π i = k )

= P( x1...xi , π i = k ) P( xi +1...xL | π i = k )

后向因子：

bk (i) = P( xi +1...xL | π i = k )

12-12-6

生物信息学算法

后项算法的伪代码

12-12-6

生物信息学算法

后验概率

u 

条件概率为：

P( x, π i = k ) = f k (i)bk (i)

u 

后验概率为：

f k (i)bk (i) P(π i = k | x) = P( x)

12-12-6

生物信息学算法

HMM中的前后向算法

赌场问题中处于公平骰子状态的后验概率

12-12-6

生物信息学算法

后验解码

u 

方式1：

ˆ i = arg max P(π i = k | x) π

k

u 

方式2：

G(i | x) = ∑ P(π i = k | x) g (k )

k

12-12-6

生物信息学算法

预测性能指标

12-12-6

生物信息学算法

CpG岛的预测问题（二）

目标：“A+，T+，C+，G+”：CpG岛 u  Viterbi 算法寻找贯穿序列的最可能路径

u 

l  l  l 

41个已知的CpG岛，找到了其中的39个 预测到121个新CpG岛 结果后处理：1)间隔小于500 个碱基的预测片段连接;2)删去短于500

基的预测片段

个碱

l 

预测到67个新CpG岛 41个已知的CpG岛，找到了其中的39个 同时预测出236 个假阳性岛 à 处理后为83

u 

后验解码算法寻找CpG岛

l  l 

12-12-6

生物信息学算法

调整赌场问题的参数

参数调整后骰子公

平的后验概率

12-12-6

生物信息学算法

HMM的参数估计

u 

问题：对于一个HMM，如何利用一组独立 同分布的训练数据x1,x2,…,xn，对HMM中 转移概率和发射概率等参数进行最优 估计？

l 

所有数据的联合概率为

l ( x ,..., x | θ ) = log P( x ,..., x | θ ) = ∑ log P( x j | θ )

j =1 1 n 1 n n

12-12-6

生物信息学算法

HMM参数估计的两种情况

u 

存在两种情况

l 

训练数据包括路径信息：

n  n  n 

赌场允许大家观摩作弊的过程 序列数据中CpG岛信息已经通过实验方法标记 … 赌场对作弊的情况保密 序列数据中CpG岛没有相关实验验证 …

l 

训练数据不含路径信息：

n  n  n 

12-12-6

生物信息学算法

第一种情况下的参数估计

u 

路径已知时，参数可以通过最大似然估计（ML ）进行计算： Akl Ek (b) akl = ek (b) = ∑l ' Akl ' ∑b' Ek (b' ) 伪计数：

u 

12-12-6

生物信息学算法

HMM的参数估计算法

u 

路径未知时

l 

不能直接估计参数值 思路：

1) 算法先用akl和ek(b)的当前值通过考虑训练序列 的可能路径来估计Akl和Ek(b) 2) 然后用路径已知时的公式导出a和e的新值

u 

Baum-Welch算法：一种特殊的迭代算法

l 

l  l 

可以证明迭代过程增加模型的整体对数似然 算法收敛于一个局部极大值

生物信息学算法

12-12-6

参数akl更新

expected number of transition s from state k to state l expected number of transition s from state k

akl =

expected number of transition s from state k to state l

= ∑ expected number of transition from state k on point i to state l on point i + 1

i =1 n −1

= ∑ P( x, πi = k , πi + 1 = l | θ )

i =1

n −1

12-12-6

生物信息学算法

计算联合概率

P( x, π i = k , π i +1 = l | θ ) = f k (i)akl el ( xi +1 )bl (i + 1)

12-12-6

生物信息学算法

参数akl更新

expected number of transition s from state k

= ∑ expected number of transitions from state k to state l

= ∑∑ expected number of transition from state k on point i to state l on point i + 1

i =1 n −1

l

= ∑∑ P( x, πi = k , πi + 1 = l | θ ) =

i =1 l

n −1

l

∑f

i =1

n −1

k

(i )bk (i )

n −1 i =1

akl =

∑∑ P(π = k , π

i i =1 l

i =1 n −1

∑ P(π = k , π

i

n −1

i +1

= l | x, θ )

=

∑f

k

(i )akl el ( xi +1 )bl (i + 1)

i +1

= l | x, θ )

∑f

i =1

n −1

k

(i )bk (i )

12-12-6

生物信息学算法

参数ek(b)更新

expected number of times in state k and observing symbol b ek (b) = expected number of times in state k

ek (b) =

{i| xi =b} n −1

∑f

k

(i )bk (i )

∑f

i =1

k

(i )bk (i )

12-12-6

生物信息学算法

算法伪代码

Initialization: Pick the best-guess for model parameters (or arbitrary) Iteration:

1.  2.  3.  4.  5. 

Forward for each x Backward for each x Calculate expected numb

ers Calculate new akl, ek(b) Calculate new log-likelihood

Until log-likelihood does not change much

12-12-6 生物信息学算法

Viterbi训练算法

u 

思路：

1）使用Viterbi算法得到训练序列的最可能路径 2）估计参数并再次迭代

u 

特点：

1）算法精确收敛 2）通常性能不如Baum-Welch算法，最大化

P( x1 ,..., x n , π * ( x1 ),..., π * ( x n ) | θ )

12-12-6

生物信息学算法

HMM的参数估计算法

u 

作弊的赌场问题

实际模型

300个数据训练得到的模型

30000个数据训练得到的模型

12-12-6 生物信息学算法

HMM的扩展问题

模型拓扑的选择 u  时长建模 u  HMM算法的数值稳定性

u 

l  l 

对数变换 概率缩放

12-12-6

生物信息学算法

HMM的扩展问题(一)

u 

模型拓扑的选择：允许任何状态间转移？ 拥有充足的训练数据？

模型的限制越少, 局部极大的问题就越严重！

u 

l 

u 

l 

构建HMM需要基于被研究问题的知识

CpG岛建模：CG双核苷酸的概率不同

12-12-6

生物信息学算法

HMM的扩展问题(二)

u 

时长建模

l 

序列中连续出现l个状态的概率按照指数衰减

P(l residues) = (1- p) p l −1

l 

有时长度分布重要且显著不同于指数分布 需要设计特殊的HMM拓扑结构

l 

12-12-6

生物信息学算法

HMM的扩展问题(二)

u 

一些特殊的HMM拓扑结构

⎛ l − 1 ⎞ l − n n P(l ) = ⎜ ⎜ n − 1⎟ ⎟ p (1 − p) ⎝ ⎠

u 

时长建模:

l  l 

直接对长度分布建模 算法速度慢

生物信息学算法

12-12-6

HMM的扩展问题(二)

u 

一些特殊的HMM拓扑结构

12-12-6

生物信息学算法

HMM的扩展问题(三)

u 

下溢错误

l  l 

~ ( x ) + max(V (i) + a ~ ) 策略1：对数变换 Vl (i + 1) = e l i +1 k kl k 策略2：概率的缩放

n 

对每个i 定义一个缩放变量si

~ f l (i) = f l (i )

∏

i j =1

sj

n 

令 ∑ f (i) = 1

k k

~

~ ~ 1 f l (i + 1) = el ( xi +1 )∑ f k (i )akl si +1 k

~ si +1 = ∑ el ( xi +1 )∑ f k (i)akl

l k

n 

对于变量b

~ 1 bk (i ) = si ~ a b ∑ kl l (i + 1)el ( xi+1 )

l

12-12-6

生物信息学算法

总结

u 

优点：

l  l  l  l 

坚实的统计学基础和有效的学习算法 可以直接从原始数据中学习 兼容不同长度的输人序列 结构灵活、适用面广 经常含有大量无结构的参数 受其1阶马尔可夫性质的限制

u 

缺点：

l  l 

12-12-6

生物信息学算法

HMM的其他应用

GPHMM: an integrated hidden Markov model for identification of copy number alteration and loss of heterozygosity in complex tumor samples using whole genome SNP arrays

12-12-6

生物信息学算法