HMMs and genome annotation¶
约 5934 个字 预计阅读时间 20 分钟
在上一节中,我们介绍了获取DNA信息的一系列方法,例如通过表观遗传学捕获基因活动信息,来识别该段基因可能代表的生物学功能和意义。
为了分析未注释的DNA序列,即尚未被标注为特定基因或功能区域的序列,在生物信息学中,我们通过统计和建立模型的方式,尝试用计算机方法对这些序列进行预测,步骤如下:
查找匹配序列(比对):
- 使用 比对工具(alignment) 将该序列与已知基因组或数据库中的序列进行比较,以发现可能的相似序列。这样可以帮助推测序列的功能或起源。
观察核苷酸或k-mer频率:
- k-mer 是指长度为k的连续核苷酸片段。通过观察某些片段的频率分布,可以发现序列中是否存在某些功能片段(如基因、调控元件等)。
建模推测序列的潜在状态:
- 通过计算或生物信息学方法,推测序列的功能状态,例如它可能是外显子、内含子、启动子,或是包含 CpG岛(常见于启动子附近,未甲基化的区域)。
Markov chains¶
我们首先介绍马尔可夫模型,并进一步推导出隐马尔可夫模型在序列预测方面的应用。
马尔可夫模型(Markov Model)是一种随机过程模型,用于描述系统从一个状态到另一个状态的转移。它基于马尔可夫性质,意味着当前状态的未来行为只依赖于当前状态,与过去状态无关。
马尔可夫模型通过状态(States)和转移概率(Transition Probabilities)来描述系统。假设一个系统在任意时间点只能处于有限个状态之一,然后根据某种转移规则(通常是概率性规则)在这些状态之间移动。
马尔可夫性质的基本定义是无记忆性,即:
这意味着系统未来的状态只依赖于当前的状态,而与过去的任何状态无关。
状态空间(State Space) 是马尔可夫系统中可能存在的所有状态的集合。状态空间可以是有限的,也可以是无限的,但通常我们讨论有限状态空间的马尔可夫模型。
在生物序列分析中,状态可能代表核苷酸(A、T、C、G),或外显子和内含子的标识。
转移概率矩阵(Transition Probability Matrix) 是描述从一个状态转移到另一个状态的概率。这些概率通常表示为矩阵形式,每个元素表示从状态 \(i\) 转移到状态 \(j\) 的概率:
- 矩阵的每一行的和为1,因为对于一个给定的状态,必须转移到另一个状态。
初始状态分布(Initial State Distribution) 定义了系统最开始处于某一状态的概率。通常表示为向量,每个元素表示系统一开始处于状态 \(i\) 的概率:
因而,这种最简单的马尔可夫模型描述的就是通过 当前状态 + 状态转移矩阵 = 预测下一个状态。例如,假设当前的核苷酸为 A,我们通过对该类核苷酸序列进行总体统计,得到了状态转移矩阵,则可以预测出下一个为 A、T、G、C 的概率分别为 0.1 0.2 0.2 0.5。而对于一整个核苷酸序列的出现概率的预测则为状态转移概率的累乘,即
其中 \(a_{x_{i-1}x_i}\) 表示了从 \(x_{i-1}\) 状态到 \(x_i\) 状态的概率。
一个具体的应用是比较一个序列更可能来自哪一个模型,对应于基因中不同的生物条件下产生的序列或不同的基因区域。
我们使用 对数优势比(log odds ratio),来比较该序列由模型 1 或模型 2 生成的概率。
其中:\(P(x|model +)\) 和 \(P(x|model -)\) 分别表示序列在模型 1 和模型 2 下的生成概率。\(a_{x_{i-1}x_i}^+\) 和 \(a_{x_{i-1}x_i}^-\) 是两个模型中状态转移矩阵中的转移概率,表示从前一个状态转移到当前状态的概率。
这个对数比例的和可以帮助我们判断这个序列是更有可能来自模型 1 还是模型 2。
在基因序列中,当我们处理长序列时,可能并不是整个序列都来自同一个模型。因此,我们需要找出序列中的哪些部分更有可能来自某个模型。解决这个问题,我们可以采用隐马尔可夫模型。
Hidden Markov model¶
相比于马尔可夫模型,隐马尔可夫模型在它的基础上,将状态改编为了隐状态 \(\pi_i\),用于表示当先观测到的序列 \(x_i\) 对应的最可能的隐状态。
具体而言,HMM 由两个主要部分组成:状态序列和观察序列。区别于传统马尔可夫模型的显式状态,HMM 的状态是隐藏的,即不可直接观察到的。我们只能通过观察序列来推断隐藏状态。
HMM 的各个组成部件如下:
状态 (States):
-
HMM 假设系统在任意时刻处于某个状态。系统可以在不同状态之间跳转,并且状态转移遵循马尔可夫性质,即系统的下一个状态只取决于当前状态,与过去的状态无关。
-
状态是隐藏的,即我们无法直接看到这些状态。
观测 (Observations):
-
每个隐藏状态会发出一个可观察的信号,这个信号就是观测序列中的元素。
-
观测值可以是离散的(如字母、单词)或连续的(如频率、测量值等)。
状态转移概率 (State Transition Probabilities):
- 定义从一个隐藏状态转移到另一个隐藏状态的概率。常用状态转移矩阵 \(A\) 来表示,元素 \(a_{ij}\) 表示从状态 \(i\) 转移到状态 \(j\) 的概率。
发射概率 (Emission Probabilities):
- 每个隐藏状态发射一个观察值的概率。发射概率描述了在某个隐藏状态下,观测到某个特定符号或值的可能性。发射概率可以用矩阵 \(B\) 表示,元素 \(b_i(o_t)\) 表示在状态 \(i\) 下观测到 \(o_t\) 的概率。
初始状态分布 (Initial State Distribution):
- 初始状态分布定义了系统在时间 \(t=1\) 时处于每个可能状态的概率。用向量 \(\pi\) 表示,\(\pi_i\) 表示系统开始时处于状态 \(i\) 的概率。
因此,一个序列出现的概率公式可以表示为
其中 \(a\) 表示状态间转移概率,\(e\) 表示该状态下发生 \(x\) 的概率。
这种概率的计算需要我们知道隐藏状态序列,然后去推断这种隐藏状态下序列出现的概率,这在实际应用中往往是不现实的。通常情况下,我们是通过观测序列来推断隐藏状态,也就是解码问题。
问题可以具体定义为:给定 HMM 和观察序列,找出最可能的隐藏状态序列。即已知观测序列,推测隐藏状态序列。
Viterbi algorithm¶
维特比算法的目的是找到最可能的整个隐藏状态序列。即给定观测序列 \(x\) 和模型参数,我们希望找到一条完整的隐藏状态路径 \(\pi\),使得它的概率最大化。
公式表示为:
这意味着我们在所有可能的隐藏状态序列中,找到使得联合概率 \(P(x,π)\) 最大的路径。
显然,一种最直接暴力的办法就是枚举出所有隐状态序列,然后判断哪个的联合概率最大,通过之前的序列概率出现公式得到。但更显然的是,我们无法做到阶乘的时间复杂度的枚举。
维特比算法通过动态规划来高效地计算这个最优路径。算法的核心是递归地计算每一步中最可能的路径,避免暴力枚举所有可能的状态序列。算法的具体步骤如下:
定义维特比表:
- 定义一个维特比表 \(V[t][j]\),表示在时间 \(t\) 时刻,以隐藏状态 \(j\) 结束的最可能的路径的最大概率。
- \(V[t][j] = \max_{\pi_{1:t-1}} P(\pi_1, \pi_2, \dots, \pi_t = j, x_1, x_2, \dots, x_t)\)
初始化(时间 \(t = 1\) ): - 对每个可能的隐藏状态 \(j\):
也就是说,初始时刻以状态 \(j\) 开始的路径概率是初始状态概率 \(\pi_0(j)\) 乘以在该状态下生成观测值 \(x_1\) 的发射概率 \(b_j(x_1)\)。
递归(时间 \(t = 2, 3, \dots, T\) ):
- 对于每一个时刻 \(t\) 和每一个状态 \(j\):
其中:
-
\(V[t-1][i]\) 是前一时刻以状态 \(i\) 结束的最优路径的概率。
-
\(a_{ij}\) 是从状态 \(i\) 转移到状态 \(j\) 的概率。
-
\(b_j(x_t)\) 是在状态 \(j\) 下生成观测值 \(x_t\) 的发射概率。
-
同时记录每个时刻 \(t\) 的最佳路径来自于哪个状态 \(i\):
终止:
- 在时间 \(T\) 时刻,找到所有可能的状态中最大概率的状态:
- 这个状态 \(\pi_T^*\) 是隐藏状态序列的最后一个状态,它的路径概率是最大的。
回溯:
- 利用
backpointer
表,从 \(\pi_T^*\) 开始,回溯到 \(\pi_1^*\),重构最可能的隐藏状态序列:
维特比算法的时间复杂度为 \(O(T \times N^2)\),其中 \(T\) 是观测序列的长度,\(N\) 是隐藏状态的个数。
Forward-backward algorithm (posterior decoding)¶
维特比算法帮助我们得到了观测序列对应的完整的最可能的隐藏序列的形式。然而,有时候我们只想得到某个时刻观测状态下对应的最可能隐藏状态,而不是整个最可能的隐藏序列。这时候我们使用后验解码(Posterior Decoding)来获得单独位置上的最优解。
后验解码想要找到最可能的隐状态:
根据链式法则,我们可以展开得到:
其中 \(P(x)=\sum_\pi P(x,\pi)\),代表了在 HMM 中观测到特定观测序列的总概率。
为了计算得到分子分母对应的概率从而得到 \(P(\pi_i=k|x)\),我们需要进一步展开得到递推形式:
我们设定 \(f_k(i)=P(x_1,\dots,x_i,\pi_i=k)\) 为前向变量,\(b_k(i)=P(x_{i+1}\dots x_L|\pi_i=k)\) 为后向变量,通过递推的方式得到对应的值。
前向算法 (Forward Algorithm)
-
目的是计算给定观测序列下,以每个隐藏状态 \(k\) 结束的部分路径的概率,即前向变量 \(f_k(i)\)。
-
初始化: 当 \(i = 0\) 时,设定:
- 递推: 对于每个位置 \(i = 1, 2, \dots, L\),计算:
其中,\(e_k(x_i)\) 是第 \(k\) 个隐藏状态下发射第 \(i\) 个观测值的概率,\(a_{lk}\) 是状态从 \(l\) 转移到 \(k\) 的概率。
- 终止: 在序列的最后位置 \(i = L\),计算:
这给出了完整的观测序列 \(x\) 的总概率 \(P(x)\)。
后向算法 (Backward Algorithm)
-
后向算法的目的是计算给定每个位置 \(i\) 后面的观测值的条件下,隐藏状态 \(k\) 的后向变量 \(b_k(i)\)。
-
初始化: 在序列的最后一个位置 \(i = L\) 时,设定:
- 递推: 从 \(i = L - 1, L - 2, \dots, 1\) 递归计算:
- 终止: 计算:
结合前向和后向算法的结果,后验解码用于计算在给定整个观测序列 \(x\) 下,每个位置 \(i\) 的隐藏状态 \(k\) 的后验概率 \(P(\pi_i = k \mid x)\):
- \(f_k(i)\) 表示观测序列前 \(i\) 个观测值的前向概率,\(b_k(i)\) 表示从位置 \(i\) 开始的观测值的后向概率。
通过前向和后向算法分别计算出每个位置 \(i\) 的前向概率 \(f_k(i)\) 和后向概率 \(b_k(i)\),并结合公式 \(P(\pi_i = k \mid x) = \frac{f_k(i) b_k(i)}{P(x)}\),可以得到给定观测序列 \(x\) 下的最优隐藏状态 \(\pi_i\),整体时间复杂度为 \(O(LK^2)\)。
find parameters for an HMM¶
前面的方法假设我们知道所有的模型参数,也就是状态转移概率 \(a_{kl}=P(\pi_i=l|\pi_{i-1}=k)\) 以及状态表现概率 \(e_k(b)=P(x_i=b|\pi_i=k)\)。现在我们将介绍如何获取这些模型参数。
我们可以将学习方法分为两种:有监督学习和无监督学习。其中有监督学习提供标记了隐状态的核苷酸序列。
目标(Goals for Parameter Estimation):
我们希望找到能最大化某个参数集 \(\theta\) 下观测序列集合的似然参数。给定我们的似然函数如下:
给定 \(n\) 个独立的观测序列 \(x^{(1)}, x^{(2)}, \dots, x^{(n)}\),它们的联合概率(即总的似然函数)可以表示为这些序列各自概率的乘积:
这里,\(P(x^{(j)} | \theta)\) 表示在给定参数集 \(\theta\) 下生成第 \(j\) 个序列 \(x^{(j)}\) 的概率。
为了找到最优参数,我们需要最大化这个似然函数。我们通过最大化似然估计来寻找最佳参数集:
在有监督学习的条件下,问题变得十分简单,序列经过标注了隐状态之后,我们可以统计训练序列中每个隐藏状态之间的转移次数,以及从某个隐藏状态到观测状态的发射次数。
转移概率矩阵的估计:
其中,\(A_{kl}\)表示在标注的训练序列中,从状态 \(k\) 转移到状态 \(l\) 的次数。
发射概率矩阵的估计:
其中,\(E_k(b)\):表示在标注的训练序列中,从状态 \(k\) 发射到观测值 \(b\) 的次数。
例如:
有时候,为了避免统计概率为 \(0\) 的情况出现,我们会添加我们定义的先验信息,也就是 \(\hat a_{kl}=\frac{A_{kl}+r_{kl}}{\sum_{l'} A_{kl'}+r_{kl}}\)。
下面介绍无监督学习下的情况。
无监督学习可以被概括为上面这张图,我们通过设置初始概率(发射概率,转移概率),从而找到目前初始概率下的隐藏状态序列,再通过这个隐藏状态序列作为 label,去更新新的初始概率。
Viterbi training¶
第一种无监督学习的方法是使用维特比训练进行 HMM 参数估计,通过找到最可能的隐藏状态序列来进行参数更新。
在初始化参数后,使用维特比算法来解码,找到给定观测序列时最有可能的隐藏状态序列。之后,根据该序列重新估计模型的参数:
-
通过计算最可能的隐藏状态序列中每个状态之间的转移次数来更新转移概率矩阵。
-
通过计算最可能的隐藏状态序列中每个状态生成的观测符号频率来更新发射概率矩阵。
以上两个步骤(估计隐藏状态和更新参数)会反复循环迭代,直到参数不再发生明显变化,即达到收敛。
维特比训练是贪婪的,每次都寻找最可能的隐藏状态序列。这通常会导致模型快速收敛,但可能卡在局部最优解中。
Baum-Welch algorithm¶
Baum-Welch算法 是一种基于EM算法(期望最大化算法)的隐马尔科夫模型(HMM)训练方法,用来估计HMM模型的参数。
Baum-Welch 算法的核心思想
不同于维特比训练只考虑最可能的隐藏状态序列,Baum-Welch算法会计算所有可能的隐藏状态序列的期望值。也就是说,它不是通过硬解码来计算从一个状态到另一个状态的转移次数,而是通过对隐藏状态概率进行加权平均来计算期望的转移和发射次数。
计算期望转移次数 \(A_{kl}\)
- \(A_{kl}\) 表示:从状态 \(k\) 转移到状态 \(l\) 的期望次数,这是在给定当前参数的情况下计算得到的期望值。
Baum-Welch通过前向-后向算法,计算每个时刻可能处于某个隐藏状态的概率,然后利用这些概率计算出从某个状态转移到另一个状态的期望转移次数。
这里,我们对每一个时间点 \(i\) 进行求和,计算在时刻 \(i\) 时隐藏状态是 \(k\),在时刻 \(i+1\) 时隐藏状态是 \(l\) 的联合概率。
根据贝叶斯公式,联合概率可以写成:
其中:
-
\(P(x, \pi_i = k, \pi_{i+1} = l)\) 是观测序列 \(x\) 和隐藏状态 \(k\) 转移到 \(l\) 的联合概率。
-
\(P(x)\) 是观测序列的总概率。
我们可以将联合概率 \(P(x, \pi_i = k, \pi_{i+1} = l)\) 分解为两部分:
其中:
-
\(P(x_1, \dots, x_i, \pi_i = k)\) 是在时刻 \(i\) 前,观测到序列 \(x_1, \dots, x_i\) 并且隐藏状态为 \(k\) 的概率(通过前向算法计算,即 \(f_k(i)\))。
-
\(P(x_{i+1}, \dots, x_L, \pi_{i+1} = l | x_1, \dots, x_i, \pi_i = k)\) 是在时刻 \(i+1\) 到 \(L\) 时,观测到剩余序列并且隐藏状态从 \(k\) 转移到 \(l\) 的概率。
应用条件独立性,这一项可以进一步分解为:
-
\(P(\pi_{i+1} = l | \pi_i = k)\) 就是从状态 \(k\) 转移到状态 \(l\) 的转移概率 \(a_{kl}\)。
-
\(P(x_{i+1} | \pi_{i+1} = l)\) 是在状态 \(l\) 时观测到 \(x_{i+1}\) 的概率,即发射概率 \(e_l(x_{i+1})\)。
-
\(P(x_{i+2}, \dots, x_L | \pi_{i+1} = l)\) 是后续观测的概率,可以通过后向算法计算,即 \(b_l(i+1)\)。
因此,联合概率可以表示为:
代入到 \(A_{kl}\) 的公式中:
这就表示在所有时间点上,从状态 \(k\) 转移到状态 \(l\) 的期望次数。
多个训练序列的推广
对于多个训练序列 \(j\),我们可以将公式推广为:
计算期望发射次数 \(E_k(b)\)
- \(E_k(b)\) 表示:在状态 \(k\) 下生成观测符号 \(b\) 的期望次数。这也是在给定当前模型参数的情况下,通过所有可能的隐藏状态概率进行加权计算。
类似地,Baum-Welch通过对所有可能的隐藏状态进行加权来计算生成每个观测符号的期望发射次数。
即,对于所有观测值为 \(b\) 的时间点 \(i\),我们计算在该时间点隐藏状态为 \(k\) 的概率 \(P(\pi_i = k | x)\),并求和。
根据贝叶斯公式,状态 \(k\) 在时刻 \(i\) 的条件概率可以表示为:
其中:
-
\(P(x, \pi_i = k)\) 是观测序列 \(x\) 和在时刻 \(i\) 隐藏状态为 \(k\) 的联合概率。
-
\(P(x)\) 是观测序列 \(x\) 的总概率。
为了进一步计算 \(P(x, \pi_i = k)\),我们可以将其分解为两部分:
这两部分分别通过前向算法和后向算法来计算:
-
前向概率 \(f_k(i)\):表示在时刻 \(i\) 前,观测到部分序列 \(x_1, \dots, x_i\) 并且隐藏状态为 \(k\) 的概率。
-
后向概率 \(b_k(i)\):表示在时刻 \(i\) 后,观测到剩余序列 \(x_{i+1}, \dots, x_L\) 并且当前隐藏状态为 \(k\) 的概率。
因此:
而总的观测序列的概率 \(P(x)\) 可以通过前向算法最后一个时间步的结果求和得到:
于是,我们可以将条件概率表示为:
代入到期望发射次数 \(E_k(b)\) 的公式中:
这就是在状态 \(k\) 下,观测到符号 \(b\) 的期望次数。
对于多个训练序列 \(j\),我们可以将其推广为:
即对所有训练序列进行求和,其中每个序列都会有相应的前向概率 \(f_k^j(i)\) 和后向概率 \(b_k^j(i)\)。
更新转移概率 \(a_{kl}\)
- 转移概率 \(a_{kl}\) 的更新公式为:
即,从状态 \(k\) 转移到状态 \(l\) 的转移概率等于从状态 \(k\) 转移到状态 \(l\) 的期望次数 \(A_{kl}\) 除以从状态 \(k\) 转移到任何状态的总期望次数。
更新发射概率 \(e_k(b)\)
发射概率 \(e_k(b)\) 的更新公式为:
即,在状态 \(k\) 下生成观测符号 \(b\) 的发射概率等于在状态 \(k\) 下生成符号 \(b\) 的期望次数 \(E_k(b)\) 除以在状态 \(k\) 下生成任意符号的总期望次数。