系统发育分析中的模型比较 MrBayes

系统发育分析时，有太多的参数和模型需要选择，实际分析中总是被一系列复杂的参数和模型搞得头晕眼花。到底该怎样判断使用哪些模型、修改哪些参数，貌似是系统发育分析中最棘手的问题了。 MrBayes 里提供了两种可以尝试的方法。这里，以构建系统树时应该使用哪一种树分枝先验模型为例，介绍这些模型比较方法。

构建系统树时可选的分枝模型大致分三种：non-clock, Strict-clock, 以及 Relaxed-clock。对于一组目标序列来说，哪一种分枝模型最适于构建系统树？大致有两种方法来判断。

其一，是根据各自分枝模型的谐波均值（harmonic mean）大小，粗略代表各模型的边际似然性（marginal likelihoods），进而使用贝叶斯因子（Bayes factor）找出最佳模型。其二，是根据更准确的踏脚石取样（stepping-stone sampling method）。

本例使用的是田螺科部分物种的 16S、28S、COI、H3 序列文件，序列总长度为 2067 ，共计 56 个个体。其中第 34 号和第 35 号是相对古老的分布于欧洲的 Viviparus 属，拟作为系统树的外群。物种图片参考 e.g. ^[Sengupta, M. E., T. K. Kristensen, H. Madsen, and A. Jørgensen. 2009. Molecular phylogenetic investigations of the Viviparidae (Gastropoda: Caenogastropoda) in the lakes of the Rift Valley area of Africa. Mol Phylogenet Evol 52: 797-805.] ^[Van Bocxlaer, B., and E. E. Strong. 2016. Anatomy, functional morphology, evolutionary ecology and systematics of the invasive gastropod Cipangopaludina japonica (Viviparidae: Bellamyinae). Contributions to Zoology 85: 235-263.] ^[Cianfanelli, S. et al. 2017. First European record of Sinotaia quadrata (Benson, 1842), an alien invasive freshwater species: accidental or voluntary introduction? (Caenogastropoda: Viviparidae). Bollettino malacologico 53: 150-160.] ^[Stelbrink, B., von Rintelen, T., Albrecht, C. et al. Forgotten for decades: Lake Lanao and the genetic assessment of its mollusc diversity. Hydrobiologia 843, 31–49 (2019) doi:10.1007/s10750-018-3666-0]。

谐波均值法（Harmonic Mean）

首先运行 non-clock 模型，对该模型来说，定义外群是重要的一个步骤。在 MrBayes 中仅能定义一个 taxa 为外群。

MrBayes > execute viviparidae_65ind.nex
MrBayes > outgroup 34
MrBayes > lset nst=6 rates=invgamma
          # set autoclose=yes nowarnings=yes
MrBayes > mcmc ngen=100000
MrBayes > sumt
MrBayes > sump

接着运行 strict-clock 模型。对于 strict 和 relaxed 模型来说，定义外群显得不那么重要，但是可以在树结构中强制定义树的根（MrBayes Manual v3.2 p56 ）。这里强制定义树的根为 34 号和 35 号，此两种模型下提供的根位置的信息量很弱，并不像 non-clock 中定义的外群那样有强烈的作用。这种根，有助于模型正确地推断树基部的变异率，或者避免人为定义树根造成的误差，还能加速运行中的收敛过程（MrBayes Manual v3.2 p59 ）。

MrBayes > lset nst=6 rates=invgamma
MrBayes > constraint ingroup partial = 1-33 36-. : 34 35
MrBayes > prset topologypr=constraints(ingroup)
MrBayes > prset brlenspr = clock:uniform
MrBayes > mcmc ngen=100000
MrBayes > sumt
MrBayes > sump

再运行 relaxed-clock 模型。同样强制定义树的根为 34 号和 35 号，并使用 birthdeath 模型作为 relaxed 模型。

MrBayes > lset nst=6 rates=invgamma
MrBayes > constraint ingroup partial = 1-33 36-. : 34 35
MrBayes > prset topologypr=constraints(ingroup)
MrBayes > prset brlenspr = clock:birthdeath
MrBayes > mcmc ngen=100000
MrBayes > sumt
MrBayes > sump

三种模型运行结果均可用 sump 命令，来查看结果参数中总体谐波均值（Harmonic mean）的大小，总体谐波均值（Harmonic mean）越大则对应模型的边际似然值（marginal likelihood）越大，即代表模型更优。

结果中，首先要看的是，每个模型下的若干 run （这里各自仅有 2 个 run）之间的谐波均值差别如何，如果差别过大，则说明模型不稳定。Non-clock 模型的 2 个 run 较稳定，的总体谐波均值为 -10158.04； strict-clock 模型的 2 个 run 不太稳定，总体谐波均值为 -10304.86；relaxed-clock 模型的 2 个 run 较稳定，总体谐波均值为 -10280.29。那么，其中最优的是 non-clock 模型，其次为 relaxed-clock 模型，最不优的是 strict-clock 模型。但是模型之间差别的显著程度该如何衡量？也就是说，两个模型的总体谐波均值相差多大，才能有效地比较出模型的优劣呢？

根据 Kass & Raftery, 1995 （如下表所示），模型之间相差 5 个 log 似然值单位时，即可视为强力证据判断模型优劣。本例中，non-clock 比 relaxed-clock 的谐波均值大了 122 个 log 似然值单位， relaxed-clock 又比 strict-clock 大了 24 个 log 似然值单位。由此可见，这三种模型的差异是极其明显的（MrBayes Manual v3.2 p57）。

10 | > 150 | Very strong

模型的 Bayes factor 之间可方便地直接比较，不需要像层级似然比检验（hierarchical likelihood ratio test）那样繁琐。两个模型之间的 Bayes factor (B12) 的对数值，就等于即为两模型的边际似然性（P(D|M1), P(D|M2)）的对数值之差（Kass & Raftery, 1995）。

边际似然值（marginal likelihood）的确切计算非常繁琐，但是在这里，当 MCMC 取样结果中出现谐波均值（Harmonic mean）时，就可以把 Harmonic mean 看做是 marginal likelihood 的对数，即 loge(P(D|Mx))。这样一来，两个 Harmonic mean 的差就可以看做是 loge(B12) 了。这个差值再乘以系数 2 ，即可得到上表中的 2loge(B12) 进行判定了（MrBayes Manual v3.2 p94）。即

就本例的结果来说，non-clock 比 relaxed-clock 的谐波均值大了 122 个 log 似然值单位，那么这两个模型之间的 2loge(B12) 即是 244，差异相当明显。应当使用 non-clock 模型构建系统树。

踏脚石取样法（Stepping-stone Sampling）

该方法可以直接估算出模型的边际似然值（marginal likelihood）的对数，即 loge(P(D|Mx))， MrBayes 中使用的命令是 ss ，而不再是常规的 mcmc 命令了。

Non-clock 模型：

MrBayes > execute viviparidae_65ind.nex
MrBayes > outgroup 34
MrBayes > lset nst=6 rates=invgamma
MrBayes > ss ngen=255000 diagnfreq=5000

Strict-clock 模型：

MrBayes > execute viviparidae_65ind.nex
          # 重新 execute ，以恢复到所有默认参数
MrBayes > lset nst=6 rates=invgamma
MrBayes > constraint ingroup partial = 1-33 36-. : 34 35
MrBayes > prset topologypr=constraints(ingroup)
MrBayes > prset brlenspr = clock:uniform
MrBayes > ss ngen=255000 diagnfreq=5000

Relaxed-clock 模型：

MrBayes > execute viviparidae_65ind.nex
MrBayes > lset nst=6 rates=invgamma
MrBayes > constraint ingroup partial = 1-33 36-. : 34 35
MrBayes > prset topologypr=constraints(ingroup)
MrBayes > prset brlenspr = clock:birthdeath
MrBayes > ss ngen=255000 diagnfreq=5000

各模型的边际似然值如下图所示，注意值的后面有 (ln)：

这时，直接将彼此的结果值相减，即可得出模型之间的 Bayes factor 对数值，即 loge(B12) 。

本例中，non-clock 比 strict-clock 大了 83 个对数边际似然值单位，strict-clock 比 relaxed-clock 大了 56 个对数边际似然值单位。乘以系数 2 后对比 Table 1 ，同样可得： non-clock 模型最适于本序列文件构建系统树。

更多请参考：MrBayes Manual v3.2