用贝叶斯方法重建基因进化历史

实验3 用贝叶斯方法重建基因进化历史传统的系统进化学研究一般采用的要么是表型的数据，要么是化石的证据。化石的证据依赖于考古学的发现，而表型数据往往极难量化，所以往往会得到许多极具争议的结论。如今，现代分子生物学尤其是测序技术的发展为重建进化史提供了大量的数据，如多态性数据（如SNPs或微卫星）、基因序列、蛋白序列等等。常规的做法一般都是利用某一个或者几个基因来构建物种树（species tree），但是一个基因的进化史能不能完全代表所有被研究物种的进化史呢？这是非常值得讨论的问题，但这不是我们本次实验的重点，在这里就不多赘述了。所以，我们这里所指的进化树如非特别说明，指的都是基因树（gene tree）。

经典的研究系统进化的方法主要有距离法、最大简约法（maximum parsimony，MP）、最大似然法（maximum likelihood，ML）等等。这些方法各有各的优点，也分别有其局限性，例如距离法胜在简单快速、容易理解，但是其模糊化了状态变量，将其简化为距离，也就不可避免的丧失了许多序列本身所提供的信息。而最大简约法虽然用的是原始数据，但也只是原始数据的一小部分。特别是在信息位点比较小的情况下，其计算能力还不如距离法。相对来说，最大似然法虽然考虑问题更加全面，但带来的另一个结果是其计算量大大增加，因此常常需要采用启发式（heuristic）方法推断模型参数，重建进化模型。

本实验利用的是贝叶斯方法来重建基因进化史。

1.贝叶斯方法概述

不可免俗的，我们还是要来看看贝叶斯模型，并分别对模型内部的一系列内容一一进行简单的介绍。

Bayes模型将模型参数视作随机变量（r.v.），并在不考虑序列的同时为参数假设先验分布（prior distribution）。所谓先验分布，是对参数分布的初始化估计。根据Bayes定理，可以不断对参数进行改进：

f(θ|D)=f(D|θ)f(θ)f(D)(1) 其中f(θ|D)为后验概率分布（posterior probability distribution），而f(θ)是先验概率分布（prior probability distribution），而f(D|θ)为似然值。此外

f(D)=∫f(D|θ)f(θ)

Ωdθ (2)

其中参数空间Ω=(Ψ,Φ)包含了所有可能树的集合Ψ和所有似然模型参数的集合Φ。一个树实例ψ=(τ,β)可用其拓扑结构（topology）τ和枝长参数（branch length）β表示。而似然模型则包含了其他所有参数φ。对于给定的似然模型f(D|θ)和多序列比对数据D，以及参数φ，每一组θ=(ψ,φ)代表特定的拓扑结构、枝长以及模型参数。于是，公式(1)又可写作

p(τ,β,φ|D)=L(D|τ,β,φ)p(τ,β,φ)

∑∫∫L(D|τ,β,φ)p(τ,β,φ)dφdβ

ΦBτ (3)

公式(3)中采用了∑⋅的原因是拓扑结构是离散变量，其中B为所有可能枝长的集合。如果要计算某一拓扑结构的后验概率分布，那么可以写出其边际概率分布（marginal probability distribution）

p(τ|D)=∫∫L(D|τ,β,φ)p(τ,β,φ)dφdβ

ΦB∑∫∫L(D|τ,β,φ)p(τ,β,φ)dφdβ

ΦBτ (4)

用同样的方法我们可以获得其他诸如枝长、似然模型参数的后验概率分布。贝叶斯推断依赖于后验概率，但在大部分情形下后验概率的归一化常数无法直接得到，因此需要采用数值方法如马尔可夫蒙特卡洛（MCMC）来估计参数的后验概率分布。

1.1马尔可夫链蒙特卡洛算法（Markov chain Monte Carlo, MCMC）

Metropolis-Hasting算法在参数空间Ω中进行连续的先后依赖的采样，获得一系列采样点�θ(0),θ(1),θ(2),…�使得从特定的某一采样点之后的所有采样点�θ(i+1),θ(i+2),θ(i+3),…�达到平稳分布且近似于θ的后验分布。因此根据马尔可夫链大数定律（Markov chain law of large numbers），只要保证足够的模拟代数，去掉�θ(0),θ(1),…,θ(i)�后采样拓扑结构的长期频率就近似于拓扑结构的后验概率。在参数空间Ω的马尔可夫链上，从状态θ1→θ2的概率密度函数为q(θ1,θ2)。因此Metropolis-Hasting算法的关键是构造符合上述条件的转移矩阵，使得到的稳态分布为我们想要得到的后验概率分布。

在现有状态θ下，Metropolis-Hasting算法接受（accept）新状态θ∗的概率为R=min�1,p(θ∗|D)q(θ∗,θ)

p(θ|D)q(θ,θ)� (5) 通常q都是对称的，也就是说q(θ∗,θ)q(θ,θ)=1。因此公式(5)可写作

R=min�1,L(D|θ∗)p(θ∗)

L(D|θ)p(θ)� (6) 如果R=1，接受新状态；如果R<1，则用runif(1)产生一个随机数α，如

果R>α则接受新状态，否则拒绝。

表1. MCMC算法

收敛速率可能特别慢，也就是实际上不可行。另外，MCMC还可能陷入局部最优。因此，评估收敛在MCMC算法中是一个非常重要和具有挑战性的问题。

1.2MCMC中产生新树的方法

在表1的MCMC算法2-(1)步骤中，树的更新是一个关键的问题，下面是一般的方法：

（1）随机选择一个树作为初始树（starting tree）；

（2）随机选择一个edge，采用NNI（nearest-neighbor interchange）；或者选择一个internal node，进行SPR（subtree pruning and regrafting）；再或者用TBR（tree bisection and reconnection）策略…….

（3）是否接受新的拓扑结构？

而改变枝长的方法也有很多，其中最为简单的一种策略是：

（1）随机选择一个edge，用一个在原枝长左右对称的一种分布随机产生一个数值作为候选的新枝长；

（2）是否接受新枝长？

1.3MCMC中更新模型参数的方法

对于取值范围在[0,M]的参数，如HKY模型中的转换/颠换比κ，可以随机产