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

实验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模型中的转换/颠换比κ，可以随机产生一个符合U~Uniform(−δκ,δκ)，从而产生新值κ∗=κ+U，接受与否完全看是否在取值范围内。

而对于一组和必须满足条件的数值参数来说，如四种核苷酸的分布，则可在原值的基础上构造Dirichlet分布，构造新值。

例如当前值如果是z= (z1,z2,…,z k)，其中∑z i=c i。

则我们可以构造新值z∗=cY，其中Y是从参数为(αz1,αz2,…,αz k)的Dirichlet分布中随机抽样得到的数值，这里的α为调整参数（tuning parameter），α越大，得到的新值的值就越是接近于现值。

新值接受与否的Hasting比率值是两个Dirichlet分布密度的比值。

2.设定先验概率分布对于Bayes方法来说，选择一个好的先验概率分布是非常关键的。

如果原始数据包含了大量的参数的信息，那么先验概率对后验概率分布的影响并不大，这时候也就接近于最大似然法，我们甚至可以采用毫无信息的先验概率分布（例如，Uniform分布）。

MrBayes为用户提供了大量的先验概率以供选择，用户可以采用不同的先验，并比较其结果，看看结果是否对先验概率敏感。

2.1碱基稳态频率所谓碱基稳态频率，指的是当进化达到稳态时，某一个位置上四种碱基的出现频率{p A,p C,p T,p G}，这种频率不再随着碱基取代的发生而改变。

由于p A+p C+ p T+ p G=1。

所以我们在设定参数时，常用Dirichlet分布来设定其先验概率。

2.2碱基取代速率矩阵对于DNA序列来说，在不考虑gap的情形下，其进化的历史实际上就是碱基替换的演变历史。

很自然的，我们可以用一个4×4的矩阵来表示碱基取代模型，这样就有了16个参数，但是其中有4个参数是可以通过其他参数来确定的，因此就有了16-4=12个自由参数。

这个速率矩阵具有以下性质：（1）非对角元素非负；（2）对角元素是本行其他非对角元素和的相反数；（3）因此，每行的和为0这是一个简单的例子：一般都写作Q=�a b c d e f g h i j k lm n o p �如今常用的碱基取代模型主要有：（1）JC69（Jukes-Cantor，1969）假设b=c=d=e=g=⋯=n=o，则a=f=k=p=−3b，且满足p(A)=p(C)=p(G)=p(T)=0.25；所以参数个数为1，对应MrBayes的设置lset nst=1；（2）K2P（Kimura，1980）：假定四种碱基的稳态比例均为0.25；转换（transition，嘌呤<->嘌呤或者嘧啶<->嘧啶）和颠换（transversion，嘌呤<->嘧啶）具有不同的速率。

所以，参数个数为2，对应MrBayes的设置为lset nst=2；（3）HKY模型（Hasegawa-Kishino-Yano 1985）：其他设置与K2P相同，唯有碱基的稳态比例可以自由设定；（4）GTR模型：四种碱基的组成是自由的，取代矩阵是对称的，所以取代矩阵中有6个参数，对应于lset nst=6。

2.3碱基取代速率差异模型大部分方法认为整条序列的碱基取代速率是均匀的，但这个假设并不符合我们的常识。

因此，MrBayes还引入了碱基取代速率差异模型，用Gamma分布来作为这个模型参数的先验概率。

之所以采用Gamma分布，乃是因为Gamma分布是一种非常灵活多变且强大的概率分布模型，只要设定Gamma分布的shape参数就可以定义一个Gamma分布。

2.4进化树模型进化树模型有2类参数：一是拓扑结构（topology）参数，另一个是枝长（branch length）参数。

附录一、MrBayes教程注：< >内为需要输入的内容，但不包括括号。

所有命令都需要在MrBayes >的提示下才能输入。

1.输入文件格式：Nexus文件输入格式为Nexus文件（ASCII，一种格式化的文本文件，如图）：这个文件中（每一个完整的语句后面都是一个分号）：dimensions——表示输入多序列比对的大小，包括序列个数（ntax），多序列比对的序列长度（nchar）；format——告诉程序输入的文本的类型，这里datatype=dna宣示为DNA序列，序列中gap=-，缺失核苷酸用？表示（missing=？）；matrix——宣示接下来是序列的主体；或者还有其他信息：interleave=yes代表数据矩阵为交叉序列interleaved sequencesnexus文件可由MacClade或者Mesquite生成。