MrBayes操作指南.pdf

MrBayes 教程
传统的系统进化学研究一般采用的要么是表型的数据,要么是化石的证据。
化石的证据依赖于考古学的发现,而表型数据往往极难量化,所以往往会得到许
多极具争议的结论。如今,现代分子生物学尤其是测序技术的发展为重建进化史
提供了大量的数据,如多态性数据(如 SNPs 或微卫星)、基因序列、蛋白序列等
等。常规的做法一般都是利用某一个或者几个基因来构建物种树(species tree),
但是一个基因的进化史能不能完全代表所有被研究物种的进化史呢?这是非常
值得讨论的问题,但这不是我们本次实验的重点,在这里就不多赘述了。所以,
我们这里所指的进化树如非特别说明,指的都是基因树(gene tree)。
经典的研究系统进化的方法主要有距离法、最大简约法(maximum
parsimony,MP)、最大似然法(maximum likelihood,ML)等等。这些方法各有
各的优点,也分别有其局限性,例如距离法胜在简单快速、容易理解,但是其模
糊化了状态变量,将其简化为距离,也就不可避免的丧失了许多序列本身所提供
的信息。而最大简约法虽然用的是原始数据,但也只是原始数据的一小部分。特
别是在信息位点比较小的情况下,其计算能力还不如距离法。相对来说,最大似
然法虽然考虑问题更加全面,但带来的另一个结果是其计算量大大增加,因此常
常需要采用启发式(heuristic)方法推断模型参数,重建进化模型。
本实验利用的是贝叶斯方法来重建基因进化史。
1. 贝叶斯方法概述
不可免俗的,我们还是要来看看贝叶斯模型,并分别对模型内部的一系列
内容一一进行简单的介绍。
Bayes 模型将模型参数视作随机变量(.),并在不考虑序列的同时为参数
假设先验分布(prior distribution)。所谓先验分布,是对参数分布的初始化估计。
根据 Bayes 定理,可以不断对参数进行改进:
f( |D) = f(D| )f( ) (1)
f(D)
θθ
其中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)又可写作φθψφ
L(D| , , )p( , , )
p( , , |D) = (3)
L(D| , , )p( , , )d d
B τβφτβφ
τβφ∑τ∫∫Φτβφτβφφβ
公式(3)中采用了的原因是拓扑结构是离散变量,其中B为所有可能枝长的
∑
集合。如果要计算某一拓扑结构的后验概率分布,那么可以写出其边际概率分布⋅
(marginal probability distribution)
L(D| , , )p( , , )d d
p( |D) = B (4)
L(D| , , )p( , , )d d
∫ B∫Φτβφτβφφβ
τ∑τ∫∫Φτβφτβφφβ
用同样的方法我们可以获得其他诸如枝长、似然模型参数的后验概率分布。
贝叶斯推断依赖于后验概率,但在大部分情形下后验概率的归一化常数无法直接
得到,因此需要采用数值方法如马尔可夫蒙特卡洛(MCMC)来估计参数的后验
概率分布。
马尔可夫链蒙特卡洛算法(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的概率