系统演化树重建的贝叶斯方法：基本原理

演化
目录

在分子系统演化学（molecular phylogenetics）领域，贝叶斯系统演化学（Bayesian phylogenetics）是在20世纪90年代末，借助于计算机技术和算法的进步而发展起来，并且得到广泛应用的一种方法论1。在此，我对该方法的基本原理做一个简明扼要的总结。作为初学者，本文难免存在谬误，欢迎读者指正。对系统生成树重建感兴趣的读者，不妨阅读文末的参考文献。

1. 基本属性

1.1. 贝叶斯概率

一言以蔽之，贝叶斯概率论研究的是如何放马后炮，为事后诸葛亮服务。简单地说，某事件的贝叶斯概率就是：你原本以为这件事发生的可能性（先验概率 prior probability）× 如果它确实发生了，产生已知结果的可能性 ÷ 观察到已知结果的可能性 = 已知结果的情况下，该事件真实发生过的可能性（后验概率 posterior probability）。

在演化树重建分析中，虽然贝叶斯方法或者说贝叶斯推断（Bayesian inference）与最大似然法（maximum likelihood）均基于个体特征（character-based），但二者采用截然相反的视角1, 5。根据贝叶斯概率理论，第i棵演化树（包括其拓扑结构τi和参数集θi）和数据D（例如一组序列排列）存在如下关系1： P(τi,θi|D)=P(τi,θ)P(D|τi,θi)P(D)

其中，

  • P(τi,θ|D)：在获知结果（数据）D的情况下，真实存在该演化树的概率（后验概率）。
  • P(τi,θ)：本演化树（拓扑与参数集）为真实系统演化树的先验概率；
  • P(D|τi,θ)：在本演化树为真实的系统演化树时，观察到特定数据D的（条件）概率——似然2；
  • P(D)：不考虑其它任何因素，单纯观察到特定数据D的概率。

假设对于任意连续参数τ∈Tc和θ∈Θc，演化树均包含某特定分支（clade，或者split，对节点的一种特定划分），则该分支C的后验概率为包含该分支的所有演化树的后验概率积分： P(C|D)=∫Tc∫ΘcP(τ,θ|D)dθdτ

假设在演化树上，共有n个末端节点（tips，包括个体、基因组或物种等对象），那么令B(n)为树空间的大小（也就是所有可能的演化树的数量），则对于有根树（rooted tree）2， B(n)=(2n−3)!2n−2(n−2)!

对于无根树（unrooted tree）4， B(n)=(2n−5)!2n−3(n−3)!
然后，利用全概率公式（formula of total probability），我们可以把观测到数据的概率P(D)展开为2： P(D)=B(s)∑i=1P(D|τi,θi)P(τi,θi)
演化树的参数（parameters）包括各个分支长度（branch length）v与突变参数（parameters of character change）s1。假设这些参数均为连续实数，则对于特定的演化树T(τi,θi)，条件概率P(D|τi,θi)本身可以用全概率公式展开为一个似然函数（likelihood function）f(D|τi,v,s)在参数空间V与S上的二元积分2：
P(D|τi,θi)=∫V∫Sf(D|τi,v,s)f(v,s)dsdv

1.2. 二叉树的节点和分支数

假设二叉树（binary tree）具有n个末端节点，则利用二叉树内部节点逐级递增的特点，我们可以计算出如下拓扑特征。

有根树

  • n−1个内部节点（internal nodes或ancestral nodes）——这很好理解，因为从根节点出发，每添加一个内部节点，末端节点也会相应地增加1个，而最简单的有根树只包括1个根节点和2个末端节点。显然，每个有根树具有n−1+n=2n个节点。
  • 2n−2条分支（因为从最简单的、只有0个内部分支和2个末端分支的有根树出发，每增加一个内部分支，也就是增加一个内部节点，末端分支的数量也会增加一个，所以内部分支与末端分支的数量始终相差2），包括n−2条内部分支和n条末端分支4。

无根树

  • n−2+n=2n−2个节点（因为最简单的无根树只有2个末端节点，而每添加一个内部节点，末端节点同样会增加1个），包括n−2个内部节点。
  • 2n−3个分支，包括n−3个内部分支（因为从最简单的、只具有1个末端分支的无根树出发，需要增加至少2个末端分支才能开始添加内部节点和分支，所以内部分支与末端分支永远保持3的差值）与n个末端分支。

2. 后验概率

2.1. 计算

后验概率是贝叶斯推断的目标。具有最大后验概率的演化树可以被当作对真实演化树的最优估计8。然而在现实中，后验概率往往因为计算量太大而难以直接获得。于是，人们采用马尔科夫链蒙特卡洛方法（Markov chain Monte Carlo3，MCMC，即使用蒙特卡洛仿真产生马尔科夫链的一系列算法）对后验概率总体采样并据此估计模型参数。一般来说，这些MCMC算法包括两大步骤（初始化、迭代式链增长与采样），详细的过程可以阅读参考文献1中的BOX 2。重要的是，这些算法借助后验概率比例避免了计算P(D)，同时，在迭代足够多次以后，特定的树(τi,θi)在马尔科夫链中出现的频率将近似P(τi,θi|D)，从而实现对该树后验概率的估计1, 6。

2.2. 意义

一般地，一棵演化树的后验概率有两种解释：

  • 在给定数据和突变模型9的前提下，该树为真的概率1, 6；
  • 特定分支的后验概率P(C|D)近似该分支的参数自举法7 （parametric bootstrap）估计，但与非参数自举（nonparametric bootstrap）估计不相关6。这种相关与不相关也是后验概率的性质之一。

2.3. 性质

后验概率具有如下特点：

  • 对突变模型的不当使用（model violations）敏感1, 6。特别地，当模型参数不足（过于简化）时，后验概率会明显上浮，而在模型具有过多参数（过于复杂）时，后验概率受到的影响较小（参数估计的方差会增大）6。
  • 由真实数据得到的演化树和分支经常表现出过高的后验概率（可能接近100%）1。

3. 先验概率

3.1. 计算

假设演化树空间由任意连续参数向量τ∈T和θ∈Θ确定11，则贝叶斯演化树的先验概率分布的概率密度函数可以表示为f(D|τ,θ)10。在贝叶斯分析中，我们仅使用包含所有参数的联合概率密度函数（joint probability density）6。显然，当有参数相互独立时，联合先验概率密度函数为每个参数的先验密度函数的乘积。

3.2. 选择

在主观贝叶斯分析（subjective Bayesian analysis）中，先验概率（分布）反映人们对假设或参数的预判6。然而，这样的选择会损害研究结论的客观性。尽管先验概率应该涵盖突变模型的所有方面，但为了减少先验概率对参数估计的影响，人们通常选择模糊（通用）的先验概率，例如最常用的、与拓扑结构无关的均匀分布（uniform prior）6, 12。在实践中，研究者应尝试不同的先验概率组合以证明这些选择不会对结论产生重大影响6。

参考文献

  1. Yang, Z., & Rannala, B. (2012). Molecular phylogenetics: principles and practice. Nature Reviews Genetics, 13(5), 303–314. https://doi.org/10.1038/nrg3186. （比较不同生成树重建方法的经典文章）
  2. Huelsenbeck, J. P., & Ronquist, F. (2001). MRBAYES: Bayesian inference of phylogenetic trees . Bioinformatics, 17(8), 754–755. https://doi.org/10.1093/bioinformatics/17.8.754. （介绍MrBayes的首发版本）
  3. A Zero-Math Introduction to Markov Chain Monte Carlo Methods. （深入浅出地介绍了MCMC的原理）
  4. 在参考文献1的BOX 1中，误作“s−1 branch-length parameters”，应改为s−2。
  5. 最大似然法建立在频率统计（frequentist statistics）之上。与贝叶斯概率论相对，频率统计是经典概率论的基础。
  6. Alfaro, M. E., & Holder, M. T. (2006). The Posterior and the Prior in Bayesian Phylogenetics. Annual Review of Ecology, Evolution, and Systematics, 37(1), 19–42. https://doi.org/10.1146/annurev.ecolsys.37.091305.110021.（一篇详尽的关于贝叶斯系统演化、先验和后验概率的文献综述）
  7. 关于自举法名称的由来和基本原理，可阅读《关于Bootstrap的通俗讲解》一文。
  8. Huelsenbeck, J. P., Ronquist, F., Nielsen, R., & Bollback, J. P. (2001). Bayesian Inference of Phylogeny and Its Impact on Evolutionary Biology. Science, 294(5550), 2310 LP – 2314. http://science.sciencemag.org/content/294/5550/2310.abstract.
  9. 可以使用与最大似然演化树法相同的突变模型1。
  10. 显然，对于离散参数，演化树的先验概率为P(D|τ,θ)。在这里，我仅讨论更一般、也更符合实际的情况——先验概率密度函数。
  11. 每棵演化树的拓扑和分支长度都各由一组参数确定。比如，如我在1.2.1中所述，具有n个末端节点的有根树需要2n−2个分支长度参数。为了数学公式的简洁性，在本文中我将这两组参数用向量τ和θ表示。每个向量由若干参数元素组成。
  12. 文献2指出，一个常用的均匀先验概率密度函数为1/B(n)。