BEAST2简介

背景

在深入进行使用BEAST2进行复杂分析之前，需要了解基本的工作流程和概念。虽然BEAST2尽可能地用户友好，但可能性的数量可能令人不知所措。

在这个简单的教程中，您将了解BEAST2的基本工作流程以及用于解释分析结果的最常用软件工具。请记住，此教程仅旨在帮助您开始使用BEAST2。本教程不会详细讨论所有选择和概念，因为它们在进一步的教程中会有所讨论。教程中穿插着一些讨论话题。这些讨论话题是可选的，但如果您逐个完成它们，您将更好地理解本教程中讨论的概念。随时跳过讨论话题，稍后再回来，可以在运行分析文件时或完成整个教程后再回来。

本练习中使用的程序

BEAST2（http://www.beast2.org）是一个免费的软件包，用于使用MCMC进行分子序列的贝叶斯进化分析，并严格面向使用根据时间测量的系统发育树进行推断。本教程是为BEAST v2.7.x编写的（[Bouckaert et al., 2014; Bouckaert et al., 2019)](#Bouckaert2014）。

BEAUti2 - 贝叶斯进化分析实用工具

BEAUti2 是用于生成 BEAST2 XML 配置文件的图形用户界面工具。

BEAST2和BEAUti2都是Java程序，这意味着完全相同的代码可以在所有平台上运行。对我们来说，这意味着界面在所有平台上都是相同的。本教程中使用的屏幕截图是在Mac OS X计算机上拍摄的；然而，这两个程序在Windows和Linux上的布局和功能都是相同的。BEAUti2作为BEAST2软件包的一部分提供，因此您无需单独安装它。

树注释器

TreeAnnotator 用于总结树的后验样本，生成最大支系可信度树。它还可以用于总结和可视化其他树参数的后验估计（例如节点高度）。

TreeAnnotator是BEAST2软件包的一部分，因此您无需单独安装它。

追踪器

Tracer (http://beast.community/tracer) 用于总结马尔可夫链采样的各种参数的后验估计。该程序可用于视觉检查和评估收敛性。它有助于快速查看参数的中位数估计值和95%的最高后验密度区间，并计算参数的有效样本量（ESS）。它还可用于研究潜在的参数相关性。我们将使用 Tracer v1.7.x。

FigTree

FigTree（http://beast.community/figtree）是一个用于查看树形图并生成出版质量图像的程序。它可以解释由TreeAnnotator在汇总树上创建的节点注释，使用户能够显示基于节点的统计信息（例如后验概率）。我们将使用FigTree v1.4.x。

DensiTree

BEAST2使用贝叶斯分析提供了树空间中不确定性的估计。这个分布由一组树表示，可能相当庞大且难以解释。DensiTree是用于对树集进行定性分析的程序。DensiTree可以快速获取树集的性质，如受支持的类群、树高的分布和拓扑不确定性的区域。

DensiTree是BEAST2软件包的一部分，因此您无需单独安装它。

实践：使用BEAST2进行简单分析

本教程将指导您分析来自十二种灵长类动物物种的序列比对。本教程的目的是共同估计基因系统发育和进化速率。更一般地，本教程旨在向新用户介绍基本工作流程，并指出在贝叶斯框架中使用BEAST2进行完整测序数据分析的步骤。

完成本教程后，您应该能够：

在 BEAUti2 中设置简单 BEAST2 分析的所有组件
使用校准节点来校准分子钟
使用 Tracer、FigTree 和 DensiTree 检查收敛性并分析结果

数据

在我们开始之前，我们需要下载教程的输入数据。在本教程中，我们使用一个名为primate-mtDNA.nex的NEXUS文件，其中包含我们将要分析的十二个灵长类动物线粒体基因组的序列比对和元数据。元数据中包含了将比对分成4个区域的信息，以及其他信息。

非编码区
第一个密码子位置
第二个密码子位置
第三个密码子位置

可以从《驯服野兽》网站https://taming-the-beast.org/tutorials/Introduction-to-BEAST2/或Github下载对齐文件。

从 taming-the-beast.org 下载

在左侧面板下方的 Data 标题下有一个指向对齐文件 primate-mtDNA.nex 的链接。右键单击该链接，然后选择 “另存为…”（Firefox 和 Chrome）或 “链接另存为…”（Safari），将文件保存到本地驱动器上的便捷位置。请注意，某些浏览器会自动将文件的扩展名从 .nex 更改为 .nex.txt。如果出现这种情况，只需再次重命名文件即可。

或者，如果您左键单击该链接，大多数浏览器将显示对齐文件。然后，您可以按 文件 > 另存为 以存储文件的本地副本。请注意，某些浏览器会向文件中注入 HTML 头部，这将使其在 BEAST2 中无法使用（这是下载数据文件的次优选项）。

同样，您也可以下载本教程中分析的示例 .xml 文件，以及 预处理 输出的 .log 和 .trees 文件。我们建议仅在检查结果或遇到严重困难时才下载这些文件。

从 Github 下载

如果您导航到 Github 存储库，您可以直接从 Github 下载原始数据文件，或者克隆/下载存储库到本地驱动器。

请注意，本教程是根据 CC BY 4.0 许可协议分发的，任何人都有权利自由使用（和修改）它，只要给予适当的信用并以相同方式许可更新的材料。

使用BEAUti创建分析文件

要使用BEAST运行分析，需要准备一个XML格式的配置文件，其中包含BEAST2运行分析所需的所有内容。BEAST2 XML文件包含：

数据（通常是序列比对）
模型规范
初始值和参数约束
MCMC算法设置
输出选项

尽管可以在文本编辑器中从头开始创建这些文件，但这可能会很复杂，而且并不是很直接。BEAUti是一个用户友好的程序，旨在帮助您生成BEAST的有效配置文件。必要时，该文件可以稍后手动编辑，但建议使用BEAUti生成文件（至少用于初始分析阶段）。

首先启动 BEAUti2。

导入对齐

为了让BEAST2访问数据，必须将序列比对添加到配置文件中。

将文件 primate-mtDNA.nex 拖放到打开的 BEAUti 窗口中（应该在 Partitions 选项卡中）。

或者，使用 文件 > 导入比对 或单击窗口左下角的 +，然后定位并单击比对文件。

一旦完成这一步，数据应该会出现在 BEAUti 窗口中，看起来应该如图 1 所示。

图 1: 数据导入到 BEAUti2。

一种常见的解决站点间速率异质性（不同站点之间的替换速率变化）的方法是使用Gamma站点模型。在这个模型中，假设速率变化遵循Gamma分布。为了使分析变得可行，Gamma分布被离散化为少量的区间（通常为4-6个）。然后，每个区间的均值作为整体替换速率的乘数。然后为每个缩放替换速率计算过渡概率。为了计算一个站点的似然性，在每个Gamma速率类别下计算P(data | tree, substitution model)，然后将结果相加以对所有可能速率进行平均。如果怀疑某些站点变异速度比其他站点快，但这些站点在比对中的确切位置是未知或随机的，这是一个方便的方法。

另一种解释站点间速率异质性的方法是将比对分割为明确的分区，并为每个分区指定独立的替代模型。当人们确切知道比对中哪些位置预计会以不同速率进化时，这一点尤为重要。在我们的示例中，我们将比对分割为编码和非编码区域，并进一步将编码区域分割为第1、第2和第3密码子位置。这些信息被编码为元数据存储在.nex文件中，BEAUti会自动处理这些信息，生成Partitions选项卡中显示的四个分区，如图1所示。

双击不同的分区（在文件列下）以查看各个对齐。

图2：灵长类动物线粒体DNA序列编码区第2密码子位置的分区。

图3：灵长类哺乳动物线粒体DNA编码区第三密码子位点的分区。

通过查看第2和第3密码子位置的配对情况（图2和图3），我们可以立即看到两个密码子位置之间的明显差异。对于第2密码子位置，许多祖先关系可以从群体之间共享的替代中清晰地看出，例如大型类人猿（人类、黑猩猩、大猩猩_和_猩猩 - 人类、黑猩猩、大猩猩和猩猩）以及旧世界猴（日本猕猴、恒河猴、猴蟒猴_和_西班牙猴 - 猕猴）。小型类人猿（由_长臂猿_ - 长臂猿代表）与大型类人猿共享大多数替代，但偶尔与猕猴共享替代。对于第3密码子位置，替代更多，群体关系不太明确。

讨论话题： 您认为在.nex文件的不同分区上使用独立的替代模型有充分的理由吗？您认为这足以考虑所有位点间的速率变异吗？

您将如何考虑每个分区内部位点之间的速率变异？

由于此数据集中的所有序列均来自线粒体基因组（据信在鸟类和哺乳动物中不会发生重组），因此它们都具有相同的祖先。默认情况下，BEAST2将为每个分区恢复一个独立的时间校准树，因此我们需要确保它使用所有数据仅恢复一个共享树。为简单起见，我们还将假设分区具有相同的进化分支速率分布，因此也共享时钟模型。

为了确保分区共享相同的演化历史，我们需要在BEAUti中链接时钟模型和树。

在Partitions面板中选择所有四个数据分区（使用shift+click），然后单击Link Trees和Link Clock Models按钮。

你会看到表中的时钟模型和树列都已更改为noncoding。现在我们将重命名这两个模型，以便以下选项和生成的日志文件更易于阅读。最终的设置应如图 4所示。

点击时钟模型列中的第一个下拉菜单，将共享的时钟模型重命名为 clock。

同样地，将共享的树重命名为 tree。

图 4: 链接的时钟和树模型。

在这个分析中，我们所有的序列都来自现存物种，因此都是在当今采样的（假定为 t=0）。因此，我们不需要设置采样日期，也跳过Tip Dates面板。接下来，我们需要在Site Model选项卡中设置替代模型。

选择 站点模型 选项卡。

此面板中可用的选项取决于配准数据是核苷酸、氨基酸、二进制数据还是一般数据。加载配准后可用的设置将包含我们通常希望修改的默认值。

左侧面板显示每个分区。请记住，在前一步中我们没有为不同分区链接替代模型，因此每个分区都在不同的替代模型下演化，即我们假设对齐中的不同位置以不同方式累积替代。我们需要为对齐的每个部分单独设置站点模型，因为这些模型是不相关的。但是，我们认为所有分区都是根据相同的一般模型演化（尽管具有不同的参数值）。

确保选择了 noncoding。

选中 Substitution Rate 的 estimate 复选框。

将 Gamma Category Count 设置为 4。

选中 Shape 参数的 estimate 复选框（应该已经选中）。

在 Subst Model 下拉菜单中选择 HKY。

确保 Kappa 的 estimate 复选框已选中（应该已经选中）。

从 Frequencies 下拉菜单中选择 Empirical。

请注意，当您为替换率选中 estimate 时，窗口底部灰色的 Fix mean substitution rate 复选框也会自动选中。

图5：站点模型设置。

现在面板应该看起来像图 5。

我们正在使用带有经验频率的HKY替代模型。这将把频率固定为在分区中观察到的比例。这种方法意味着我们可以在不明确估计这些参数的情况下很好地拟合数据。为了在每个分区内对站点间速率变异进行建模，我们使用了一个具有4个类别的离散Gamma站点模型。现在我们可以为剩余分区中的每个分区重复上述步骤，或者我们可以采取捷径。

选择剩余的三个分区（使用 shift+click）。现在窗口看起来像图 6。

选择 noncoding 并单击 OK 以克隆来自 noncoding 的其他三个分区的站点模型。

图 6: 在分区之间克隆站点模型的快捷方式。

讨论主题: 你能想出在你勾选估计替代率时为什么修正均值替代率会自动被选中的原因吗？如果你想不明白也不用担心，后续教程会详细解释。

设置时钟模型

接下来，在主窗口顶部选择时钟模型选项卡。这是我们设置分子时钟模型的地方。在本练习中，我们将保持选择默认值为严格分子时钟，因为这些数据非常类似时钟，不需要在模型中包含分支间的速率变化。

点击时钟模型选项卡，查看设置（但不要进行任何操作）。

请注意，平均时钟速率旁边的估计框未被选中。

设置先验

在先验选项卡中，可以为每个模型参数指定先验分布。在站点模型和时钟模型选项卡中做出的模型选择确定了模型中包含哪些参数。对于这些参数中的每一个，需要指定一个先验分布。还可以为每个模型参数指定超先验（以及超超先验等）。我们还需要为树指定一个先验。在这个例子中，树的先验是灵长目物种多样性随时间的空模型。

在这里，我们指定使用校准的尤勒模型作为树先验。这是一种简单的物种分化模型，通常在考虑来自不同物种的序列时是合适的。

转到 Priors 选项卡，并在 Tree.t:tree 旁边的下拉菜单中选择 Calibrated Yule Model。

在校准的尤尔模型中，birthRate参数衡量物种形成的速率。我们将为该参数设置一个不特定的Gamma先验。

对于 birthRateY.t:tree，请从下拉菜单中选择 Gamma

使用左侧的箭头按钮展开 birthRateY.t:tree 的选项。

将 Alpha（形状）参数设置为 0.001，将 Beta（比例）参数设置为 1000。

请注意，BEAUti 在右侧显示了先验分布的图表，以及其中的一些分位数。这是为了方便参考，可以帮助我们决定先验是否合适。我们正在使用的不特定的 Gamma 先验定义在 (0,∞) 区间内，并且在其2.5% 和 97.5% 分位数之间包含了非常广泛的数值范围。

如果我们想在 Gamma 先验的一个参数上添加一个超先验，我们会在参数右侧的估计框中勾选。我们还可以通过点击下拉菜单旁边的框来更改模型参数的初始值或限制。在这里不要这样做，因为我们在这个分析中没有添加任何超先验或更改限制！

我们将保留先验的其余部分为默认值！ BEAUti 面板应如图 7 所示。

请注意，通常不建议使用默认先验，因为先验应该传达您对参数的先验知识。重要的是要了解先验为MCMC分析增加了什么信息，以及这是否符合您的特定情况。在我们的情况下，默认先验适用于这种特定分析，但是对于进一步、更复杂的分析，我们将需要清楚地了解先验的含义。获得这种理解是困难的，并且需要经验。有关选择先验的主题将在后续教程中更详细地讨论。

图 7: 先前设置。

添加校准节点

由于所有样本均来自单个时间点，因此在时间单位上没有关于系统发生树实际高度的信息。这意味着系统发生树高度（tMRCA）和替代率参数将无法区分，BEAST2 只能估计它们的乘积。为了让 BEAST2 分离这两个参数，我们需要输入额外信息，以帮助校准系统发生树的时间。

在贝叶斯分析中，来自外部来源的额外信息应以先验分布的形式进行编码。因此，我们将不得不向模型添加一个新的先验。

要在模型之前添加额外的先验，点击先验列表下方的**

添加先验按钮。如果这不会自动打开分类单元集编辑器窗口，请从下拉菜单中选择MRCA 先验**。

您将看到一个对话框（分类单元集编辑器），允许您从系统树中选择一组分类单元。一旦您创建了一个分类单元集，您将能够随后为其最近共同祖先（MRCA）添加校准信息。

将分类单元集标签设置为 human-chimp。

在左侧列表中找到 Homo_sapiens，然后单击**»**按钮将其添加到 human-chimp 分类单元集。

在左侧列表中找到 Pan，然后单击**»**按钮将其添加到 human-chimp 分类单元集。

分类单元集现在应该看起来像图 8。

单击确定按钮，将新定义的分类单元集添加到先前的列表中。

图 8：校准节点分类单元。

我们添加的新节点是一个校准节点，用于人类-黑猩猩分裂，与校准尤勒先验一起使用。为了使其正常工作，我们需要强制单系性。这将限制树的拓扑结构，使得在MCMC分析过程中人类-黑猩猩分组保持单系性。

在 human-chimp.prior 旁边勾选 单系群 复选框。

我们现在需要根据化石的先验知识在校准节点上指定一个先验分布，以校准我们的树。我们将使用均值为6百万年，标准差为0.5百万年的正态分布。这将给出大约5-7百万年的中心95%范围，这大致对应于人类和黑猩猩最近共同祖先的日期的当前共识估计。

在新添加的 human-chimp.prior 右侧的下拉菜单中选择 Normal。

使用左侧的箭头按钮展开分布选项。

将分布的 Mean 设置为 6。

将分布的 Sigma 设置为 0.5。

校准节点的最终设置应如图 9所示。

图 9: 校准节点设置前。

讨论主题: 在设置校准节点之后，clockrate.c:clock 的新先验神奇地出现在先验选项卡中！如果你回到时钟模型选项卡，你会看到均值时钟速率旁边的估计框现在被选中并变灰（无法取消选中）。

刚刚发生了什么？

设置 MCMC 选项

最后，MCMC 选项卡允许我们控制 MCMC 链的长度和存储样本的频率。它还允许更改输出文件名。

转到 MCMC 选项卡。

参数链长度指定 MCMC 链在完成之前将进行的步数（即接受的提议数）。此数字取决于数据集的大小、模型的复杂性和所需答案的精度。默认值为 10'000'000 是任意的，应相应调整。对于这个小数据集，我们最初将链长度设置为 1'000'000，以便这个分析在大多数现代计算机上只需几分钟（而不是几小时）。我们将存储间隔和预燃烧字段保留在它们的默认值。

将链长度设置为 1'000'000。

在这些常规设置下，您将找到日志记录设置。单击其左侧的箭头，可以详细查看每个特定选项。您可以控制日志文件的名称以及值将存储在每个文件中的频率。

首先扩展 tracelog 选项。这是您稍后将用于分析和总结运行结果的日志文件。日志文件的 Log Every 参数应相对于链的总长度进行设置。过于频繁的采样将导致非常大的文件，而在分析的准确性方面几乎没有额外的好处。过于稀疏的采样将意味着日志文件不会记录足够的参数分布信息。我们通常希望目标是存储不超过 10,000 个样本，因此这应设置为不低于链长度/10,000。对于本次分析，我们将使 BEAST2 每 200 个样本写入日志文件。

扩展 tracelog 选项。

将 Log Every 参数设置为 200。

保持文件名不变 ($(filebase).log).

接下来，扩展 screenlog 选项。屏幕输出仅用于监控程序的进度。由于它并不是那么重要，特别是当你在远程计算机或计算机集群上运行分析时，Log Every 可以设置为任何值。然而，如果设置得太小，屏幕上显示的信息量会实际减慢程序的运行速度。对于这次分析，我们将 BEAST2 设置为每 1,000 个样本记录到屏幕上，这是默认设置。

扩展 screenlog 选项。

将 Log Every 参数保持在默认值 1’000。

最后，我们还可以通过扩展 treelog.t:tree 来更改树的记录频率。对于具有许多分类群的大树，每个单独的树已经相当大，因此如果你记录很多树，树文件很容易变得极其庞大。如果树的记录频率过高，你会惊讶于 BEAST 如何迅速填满即使是最大的驱动器！因此，通常将树的记录频率设置得低于轨迹日志是个好主意（特别是对于具有许多分类群的分析）。但是，要小心，因为某些模型的后处理步骤（例如贝叶斯天际线图）要求轨迹和树的记录频率相同！

扩展 treelog.t:tree 选项。

将 文件名 设置为 primate-mtDNA.trees。

将 每次记录 参数保持在默认值 1’000。

最终的设置应如图10所示。

图10：日志选项。

生成 XML 文件

我们现在准备创建 BEAST2 XML 文件。这是 BEAST2 可以用来执行分析的最终配置文件。

将XML文件保存为primates-mtDNA.xml，使用文件 > 保存。

运行分析

现在运行 BEAST2 并提供您新创建的 XML 文件作为输入。您还可以更改 随机数种子 以进行运行。这个数字是 BEAST2 用来生成样本的伪随机数链的起始点。由于计算机无法生成真正的随机数，我们不得不生成看似随机的确定性数字序列，但当起始种子相同时，它们将是相同的。如果您的 MCMC 运行收敛到真实后验，那么无论提供哪个随机种子，您都将能够得出相同的结论。然而，如果您想准确重现一次运行的结果，您需要使用相同的随机数种子开始它。

运行 BEAST2 程序。

选择 primates-mtDNA.xml 作为 Beast XML 文件。

将 随机数种子 设置为 777（或选择你喜欢的数字）。

勾选 如果可用则使用 BEAGLE 库 复选框。如果你之前安装过 BEAGLE，这将使分析运行得更快。

图11：BEAST2分析的设置。

BEAST2 窗口应该类似于图 11。

点击 Run 按钮运行 BEAST2.

BEAST2 将运行直到达到链中指定的步骤数。在运行期间，它会将屏幕日志值打印到控制台，并将跟踪日志和树日志值存储到与配置 XML 文件位于同一文件夹中的文件中。屏幕输出将大致如下所示图 12。

图12：BEAST2分析的屏幕输出。

当 BEAST2 完成分析时，窗口将保持打开状态。当你尝试关闭它时，你可能会看到 BEAST2 提出问题：“你想保存吗？”。请注意，无论你对这个问题选择什么答案，你的日志和树文件始终会被保存。因此，这个问题仅限于保存 BEAST2 屏幕输出（其中包含一些关于硬件配置、初始值、操作员接受率和运行时间的信息，这些信息不会存储在其他输出文件中）。

讨论主题： 在分析运行时，看看你能否识别出 BEAUti 中与数据、模型和 MCMC 算法相关的设置部分。

在你喜欢的文本编辑器中打开 XML 文件。你能识别出在 BEAUti 中设置的任何值吗？你能在 XML 文件中识别出数据、模型规范和 MCMC 设置吗？

你能在 XML 文件中找到似然、先验和超先验吗？

分析结果

一旦 BEAST2 完成运行，打开 Tracer 以获取 BEAST2 输出的概览。当主窗口打开时，选择 文件 > 导入跟踪文件... 并选择 BEAST2 创建的名为 primate-mtDNA.log 的文件，或者直接将文件从文件管理器窗口拖入 Tracer。

打开 Tracer。将 BEAST2 创建的 primate-mtDNA.log 文件拖放到打开的 Tracer 窗口中。

或者，使用 文件 > 导入跟踪文件…（或按下 + 按钮在 跟踪文件 面板下），然后找到并点击 primate-mtDNA.log。

Tracer 窗口应如图 13 所示。

图13：示踪器显示了BEAST2对灵长类动物数据的运行摘要，MCMC链长度为1'000'000。

Tracer 提供了一些关于分析结果的有用摘要统计信息。在顶部窗口的左侧，它提供了当前加载到程序中的日志文件列表。下面的窗口显示了每个文件中记录的统计信息列表。对于每个统计信息，它提供了一系列摘要值，如均值、标准误、 медиана 以及其他可以从数据中计算得出的值。摘要值显示在右上角窗口中，右下角窗口中显示了统计信息的分布直方图。

日志文件包含后验的痕迹（这是树的似然和先验密度的乘积的自然对数）、先验、似然、树的似然和其他连续参数。在左侧选择一个痕迹会在右侧显示该痕迹的摘要统计信息。首次打开时，后验痕迹被选中，并在估计标签下显示该痕迹的各种统计信息。

对于每个加载的日志文件，我们可以指定一个 Burn-In，它在 Tracer 的文件列表表（左上角）中显示。Burn-in 旨在给马尔可夫链提供时间以达到其平衡分布，特别是当它从一个不好的起点开始时。不好的起点可能导致对后验分布中实际上具有非常低概率的区域进行过度采样，然后链才会稳定到平衡分布。Burn-in 允许我们简单地丢弃链的前 N 个样本，而不使用它们来计算摘要统计量。确定要丢弃的样本数量并不是一个简单的问题，它取决于数据集的大小、模型的复杂性和链的长度。一个好的经验法则是始终丢弃整个链长度的至少 10% 作为 burn-in（然而，在某些情况下，可能需要丢弃多达 50% 的 MCMC 链）。

在左侧列表中选择 TreeHeight 统计量，以查看对齐中所有分区共同估计的树高。Tracer 绘制所选统计量的（边际后验）直方图，并提供诸如均值和中位数等摘要统计量。95% HPD 代表 最高后验密度区间，表示所选统计量中包含 95% 后验密度的最紧凑区间。可以粗略地将其视为置信区间的贝叶斯类比。TreeHeight 统计量给出了整个树根的年龄的边际后验分布（即 tMRCA）。

在 Tracer 的左下角列表中选择 TreeHeight，并查看右侧的不同摘要统计信息。

您还可以在Tracer中比较不同参数的估计值。一旦跟踪文件加载到程序中，您可以例如比较与对齐中不同分区对应的不同突变率的估计值。

通过点击第一个突变率 (mutationRate.noncoding)，然后按住 Shift 并点击最后一个突变率 (mutationRate.3rdpos) 来选择所有四个突变率。

在右侧选择 边际密度 选项卡以一起查看这四个分布。

在显示下拉菜单中选择不同的选项，以不同方式显示后验分布。

您将能够在一个图中看到所有四个分布，类似于在图14中所示的内容。

图14：示踪器显示对齐中每个分区的突变率的四个边际概率分布。顶部的图显示了边际分布，使用核密度估计（KDE）绘制，中间为小提琴图，底部为箱线图。

讨论主题： 你能从4个突变率的边际密度中推断出什么？这在生物学上有意义吗？

你认为非编码DNA的突变率为什么与第1和第2密码子位置的突变率相似？（在图中显示图例以帮助你的分析。）

分析后验估计质量

我们应该关注的两个非常重要的汇总统计量是自相关时间（Auto-Correlation Time，ACT）和有效样本量（Effective Sample Size，ESS）。ACT是MCMC链中两个样本之间需要分开的平均状态数，以使它们不相关，即使它们是来自后验的独立样本。ACT是从样本中估计的（不包括预热期）。ESS是轨迹等效的独立样本数量。这是通过链长（不包括预热期）除以ACT来计算的。

ESS通常被视为结果样本序列的质量指标。目前尚不清楚如何准确确定ESS的大小，以使分析可信。一般来说，ESS为200被认为足够高，可以使分析有用。然而，这只是一个任意的数字，您应该始终根据自己的判断来决定分析是否已收敛。如您在图13中所见，ESS值低于100的部分被标记为红色，这意味着我们不应信任统计值，而ESS值在100到200之间的部分则被标记为黄色。

如果很多统计数据的 ESS 值为红色或黄色，我们尚未充分探索后验空间。这很可能是链运行时间不够长的结果。尝试进行与之前相同的分析，但使用更长的链。

首先通过按 文件 > 加载 再次将 XML 配置文件加载到 BEAUti 中，并选择 primate-mtDNA.xml 文件。在 BEAUti 中，将 MCMC 链长度参数更改为 10’000’000，并将 tracelog 频率更改为 1’000。

更改追踪和树日志文件名，以免覆盖先前分析的结果。您可以在文件名后添加类似 _long 的内容，以获得日志文件 primate-mtDNA_long.log 和树文件 primate-mtDNA_long.trees。

使用更新的配置文件和相同的种子 777 再次运行 BEAST2。

这将需要更多时间。图 15 显示了更长运行的估计。在这种情况下，所有参数的 ESS 值都大于 200。请记住，MCMC 是一种随机算法，因此如果您设置了不同的种子，实际数字将与图中所示的数字不完全相同。

图15：示踪器显示了BEAST2运行的摘要，MCMC链长度为10'000'000。

Tracer 还允许我们在 联合边际 选项卡下查找参数之间的相关性，如图 16 所示。当两个参数高度相关时，这可能导致 MCMC 链的收敛性差（在后面的教程中会详细介绍这一点）。

图16：树高与时钟频率估计之间的相关性。

我们还可以查看多个参数之间的相关性。

再次选择所有 4 个突变率

导航到 联合边际 标签

勾选 显示点

面板应该喜欢图 17。椭圆表示参数对之间的协方差，便于识别哪些参数对是相关的或反相关的。我们的某些突变率参数之间是否存在强相关或反相关？

图17：突变率参数之间的相关性。

讨论主题： 我们还没有探索 Tracer 中的 Trace 选项卡！

Trace 选项卡主要是一个诊断工具，用于检查后验的收敛性，评估烧入期的长度以及链是否良好混合。可以说这是 Tracer 程序中 最重要 的选项卡，也是用户应该首先查看的选项卡。

查看 Trace 选项卡中各个参数的轨迹，包括短日志和长日志文件。你能找出为什么某些参数的 ESS 值高于其他参数吗？

你认为 10% 的烧入期对于这个分析来说足够吗？

分析树估计

除了生成参数估计的样本，BEAST2 还生成系统发育时间树的后验样本。在对后验估计的质量得出任何结论之前，这些也需要进行总结。

总结树的一种方法是使用程序 TreeAnnotator。该程序将获取树的集合并找到支持最好的树。然后，它将用所有节点的平均年龄和相应的 95% HPD 范围对这个代表性摘要树进行注释。它还将计算每个节点的后验分支概率。这样的树被称为 最大分支可信度 树。

打开 TreeAnnotator.

将 Burnin percentage 设置为 10%，以丢弃日志文件中前 10% 的树。

下一个选项，后验概率限制，指定一个限制，以便如果在树的样本中发现一个节点的频率低于此限制（即后验概率低于此限制），则不会进行注释。例如，将其设置为0.5意味着只有在大多数（超过50%）树中看到的节点才会被注释。默认值为0，我们将保持不变，这意味着TreeAnnotator将注释所有节点。

将 后验概率限制 保持在默认值 0。

对于 目标树类型 选项，您可以选择从文件中选择特定的树，或者要求 TreeAnnotator 在您的样本中查找一棵树。我们将保留的默认选项 最大类群可信树 找到所有节点后验概率乘积最高的树。

将 目标树类型 保持为默认值 最大分支可信树。

接下来，选择 平均高度 作为 节点高度。这将树中每个节点的高度（年龄）设置为该类群所有树的平均高度。

选择 平均高度 在 节点高度 下拉菜单中。

最后，我们必须选择输入树日志文件并设置输出文件。

点击 选择文件 在 输入树文件 旁边，并选择 primate-mtDNA.trees。

将 输出文件 设置为 Primates.MCC.tree.

设置应如图 18 所示。您现在可以运行程序。

图18：TreeAnnotator 设置

可视化树估计

最后，我们可以使用可用的软件之一，如 FigTree，来可视化树。

打开 FigTree。使用 文件 > 打开，然后找到并点击 Primates.MCC.tree。

展开树选项，勾选 排序节点，并从下拉菜单中选择递减。

展开 提示标签 选项，增加 字体大小 直到可读。

勾选 节点条 复选框，展开选项，从显示下拉菜单中选择 height_95%_HPD。

勾选 节点标签 复选框，展开选项，从显示下拉菜单中选择 posterior。

增加 字体大小 直到可读。

取消勾选 比例条 复选框。

勾选 比例轴 复选框，展开选项，勾选 反向轴 并增加 字体大小。

图19：FigTree对估计树的可视化。

您的树现在应该看起来像图19。我们首先对树节点进行了排序。因为有很多方法可以绘制相同的树，排序节点使我们更容易相互比较不同的树。我们添加的比例尺表示树中每个节点年龄的95% HPD区间，这是通过BEAST2分析估算的。我们添加的节点标签给出了后验树集中一个节点的后验概率（即，在丢弃烧入后，记录在树日志文件中的树）。我们还可以使用FigTree显示其他统计信息，例如分支长度、节点的95% HPD区间等。可用的确切统计信息将取决于使用的模型。

讨论主题： 后验概率告诉我们哪些类群得到了高度支持，比例尺告诉我们对它们的分歧时间有多自信。

所有类群都得到良好支持吗？它们的年龄如何？

查看猿类（Hylobates, Pongo, Gorilla, Pan 和 Homo sapiens）年龄的95% HPD区间。估计的年龄是否与您之前的知识一致？

旧世界猴和新世界猴之间的分歧时间呢？（Saimiri sciureus 是这个数据集中唯一的新世界猴。）

可视化树后验（可选）

MCC树是一种将树的后验分布总结为单棵树的方法，附加了一些节点的额外信息以表示树估计中的不确定性。就像将连续参数的后验分布总结为中位数和置信区间会丢失很多信息（例如分布的形状），将一组树总结为MCC树时也会丢失很多信息。然而，视觉化后验树的集合要困难得多。

一种可能性是使用程序 DensiTree。DensiTree 不需要摘要树（因此我们在使用 DensiTree 之前不需要运行 TreeAnnotator）就能够可视化估计。

打开 DensiTree。使用 文件 > 加载，然后找到并点击 primate-mtDNA.trees。

展开显示选项并勾选 共识树 复选框。

您现在应该看到许多行，对应于您的 MCMC 链所采样的所有单个树。您还可以清楚地看到所有后验树之间的模式。

为了查看拓扑的支持情况，请选择中央视图模式。

现在展开类群菜单，勾选 显示类群 复选框和支持的文本复选框。

树应该如图 20所示。

图20：树样本的DensiTree可视化。

您还可以通过选择 帮助 > 查看类群 来查看所有不同的类群及其后验概率。在这个特定的运行中，树的估计在类群分组方面几乎没有不确定性，因为几乎每个类群都有 100% 的支持。

讨论主题： Yule 物种形成模型有一个参数 (birthRateY.t:tree)，代表物种形成率。该模型假设没有灭绝，因此所有分类群都被采样。

birthRateY.t:tree 的单位是什么？根据你的分析，你能算出物种形成事件发生的平均时间吗？（查看追踪日志）。

Yule 模型在这里是否合适？

在数据集中，大猿类（4/8 种现存物种）的采样比例远大于狐猴、懒猴和新世界猴（各一种）。你认为不均等的采样比例是个问题吗？

致谢

本教程的内容基于 [与 BEAST 2.0 的分歧约会教程]（https://github.com/CompEvol/beast2/blob/master/doc/tutorials/DivergenceDating/DivergenceDatingTutorialv2.0.3.pdf?raw=true）由 Drummond、Rambaut 和 Bouckaert 编写。

背景