贝叶斯建树原理及优缺点

注意：该前言部分使用 ChatGPT4.0 生成并由 Juse 纠正总结。

问题：请介绍一下系统发育分析中的贝叶斯推断法（原理、优缺点等），并指出几个代表生信软件。

在系统发育学中，贝叶斯推断是一种常用的统计方法，用于根据观测数据推断出物种间的系统发育关系。

原理：

贝叶斯推断的原理基于贝叶斯定理，该定理可以将先验概率和似然度结合起来，得到后验概率。

在系统发育分析中，贝叶斯推断法被用来估计种群或物种之间的演化关系。这种方法假设所有可能的树形图都存在某种先验概率，并使用蒙特卡洛马尔科夫链（MCMC）来在所有可能的树形图和其关联参数的空间中进行采样。每一步，系统会根据特定的概率提议一种新的树形图或参数，然后根据新的树形图或参数与当前状态的拟合度决定是否接受。

优点：

  1. 贝叶斯推断提供了一种直观的方式来理解不确定性和证据的权重。因为后验概率可以直接解释为在给定数据下模型或参数的概率。
  2. 贝叶斯推断允许使用先验信息，有助于引导分析并提高准确性。

缺点：

  1. 计算复杂性高，尤其是对于大数据集。MCMC方法需要长时间的运行才能确保采样的充分性。
  2. 对先验概率的选择可能对结果有影响。不合理的先验可能会导致偏差的结果。

代表性生物信息学软件：

  1. MrBayes：可能是最广泛使用的贝叶斯系统发育分析软件。它允许使用各种进化模型，并且可以并行处理以加快计算速度。
  2. BEAST：通常用于分子时钟和历史进化的分析。它可以结合时间信息来推断系统发育树，并估计进化参数。
  3. PhyloBayes：这是一个专门针对复杂模型的贝叶斯系统发育分析软件，常用于蛋白质序列的系统发育重建。

更新日志

2023.07.17 完善了一些细节。

2023.10.07 补充了如何设定起始树以增快收敛速度。

2023.11.08 根据留言版内容添加了 slurm 运行 Mrbayes 的方法链接。

Mrbayes 使用方法

关于 Mrbayes 的安装和运行原理介绍（最详尽版）：Mrbayes Manual

该 manual 在 Mrbayes 的根目录下 share/doc/mrbayes/Manual_MrBayes_v3.2.pdf 也能找到。

很好很易懂的其他教程：贝叶斯法构建进化树 from 陈连福

注意：请安装支持多线程的 MPI 版 Mrbayes。

本文将以尽量最简洁的方式介绍 Mrbayes 的运行及注意事项。

使用 slurm 运行 Mrbayes 的方法：https://github.com/NBISweden/MrBayes/issues/130

运行方法

所需要准备的文件及可用信息：

  • nex 格式的序列文件（不是 fasta）。
  • 分区信息及对应模型（可选）。

Mrbayes 也支持其他类型（例如形态学）数据的输入，但本文将侧重于序列数据（DNA or 蛋白）的处理。

①、使用多核运行 Mrbayes：

1
$ mpirun -np 10 mb

命令中的 10 表示使用十个核，使用的核数量越多，运行速度越快。需要注意的是，如果 Mrbayes 耗尽了内存，分析会变得更慢很多。

1
2
MrBayes > log start filename = mrbayes.log;
MrBayes > execute xxx.nex

建议先输入 log start filename = mrbayes.log; 命令将 Mrbayes 的所有屏幕输出储存到文件 mrbayes.log 中以便后续查阅（该命令会直接覆盖目标文件，因此在断点重续时注意更改文件名），确定有相关 log 记录后再通过 execute xxx.nex 读取序列文件。

此后的模型选择，如果没有相关的先验理解，那么可以先通过一些方法确定自己数据所适合的模型：

  • Mrmodeltest2
  • IQTREE 的 ModelFinder
  • PartitionFinder2（仅适用于多基因串联超矩阵序列）

②、确定最佳模型后，输入相关参数，以 GTR substitution model with gamma-distributed rate variation across sites and a proportion of invariable sites 为例（多基因串联情况见后文）：

1
MrBayes > lset nst=6 rates=invgamma

接着设置运行参数，可以以下列参数作为参考：

1
MrBayes > mcmcp append=no ngen=100000 printfreq=1000 samplefreq=1000 nchains=10 nruns=2 savebrlens=yes checkpoint=yes checkfreq=5000;

各个参数的含义及设置参考：

  • append=no，表示新开始一个分析，如果之前已经运行过但因为某种原因中断，可以将其调整为 append=yes 从上一个检查点继续运行，需要搭配 checkpoint=yes。
  • ngen=100000，表示分析总共迭代十万代，对于小数据集一般只需几万代就可以收敛，而大数据集可能则要至少数百万代才能收敛。可以根据自己的数据集调整，如果迭代完依旧未收敛可以再添加该值继续运行（Mrbayes 会询问）。
  • printfreq=1000，表示每一千代打印一次链状况到屏幕输出上。
  • samplefreq=1000，表示每一千代采一次样，对于需要长时间才能收敛的大数据集而言可以将该值调高以降低采样频率，避免最终文件包含的树和参数信息过多。采样量为 ngen/samplefreq。
  • nchains=10，表示每个运行使用十条链，总链数为 nchains*nruns。
  • nruns=2，表示进行两个独立的分析，用于判断收敛情况。有时，提升该值（例如到 3 或 4）会对有效样本大小（ESS）增加有奇效。

链数跟 nchains 和 nruns 有关，所有的链会被平均分到各个核中，当每个分析中链数量大于 1 时（nchains > 1）会设置一条冷链（其他均为热链），热链的存在在某些情况下对于收敛来说是必要的。在一些大数据集中热链越多收敛越快。因此设置合理的链数和内核数对于分析来说也很重要。

在 Mrbayes 的屏幕输出中，以圆括号框起来的即为热链，方括号框起来的为冷链，星号分离不同的分析。

1
0 -- [-8849.255] (-8762.724) (-8987.103) (-8764.951) * [-8802.317] ...

  • savebrlens=yes，表示记录分支长度。如果只关注树的拓扑结构且不需分支长度信息可以设置为 no 节省空间。

  • checkfreq=5000，表示每五千代进行一次收敛诊断，输出 ASDSF，例：

    1
    Average standard deviation of split frequencies: 0.178972

设定好所有参数以后，输入 mcmc; 开始分析。

在运行中及结束以后，有几个判断分析收敛程度的指标：

  • 平均标准分歧（ASDSF），低于 0.01 时表明收敛良好。
  • 潜在似然（PSRF），在 sump 输出中出现，理想情况应该接近 1，若不是则说明可能仍未收敛。
  • 有效样本大小（ESS），在 sump 输出中出现，通常情况下以大于 200 作为样本量足够的标准（PhyloBayes 则为 300）。

一般来说，三个指标都达标时我们才有充分的信心相信贝叶斯推断的结果足够可靠，以下是得到后两种信息的方法：

1
MrBayes > sump relburnin=yes burninfrac=0.25;

其中 relburnin=yes 表示使用比例，此处意为丢弃冷链中前 25% 的样本。

输出的结果中，min ESS 和 PSRF 两列即为上述要关注的两个指标。

③、确保分析收敛后，可使用以下命令得到贝叶斯推断树：

1
MrBayes > sumt conformat=Simple contype=Halfcompat relburnin=yes burninfrac=0.25;

其中 conformat=Simple 时输出的树较为整洁，能被大多数的软件识别接受，另一种选项为 Figtree，会输出一个适用于 FigTree 程序的更丰富的树格式。contype=Halfcompat 表示输出 majority rule consensus tree，该参数指定为 Allcompat 时输出 strict consensus tree。关于这两者的描述可见：Consensus tree program

输出的树文件以 .con.tre 结尾，建议使用 Figtree 打开并转换为其他所需格式。

分享多基因串联处理方法前的一些建议：

  • 官方 manual 中有更多细节，包括在收敛困难时可以用于改善收敛情况的手段以及让程序的运行效率变得更高的方法等，遇到这些问题的时候可以在里面寻找到相应的解决方案。

    • Mrbayes 可以从一个已有的树开始进行分析，因此可以先使用最大似然法等其他方法完成建树后，再提交给 Mrbayes 进行分析（有助于收敛）。详情可见 Manual 的第 94 页，TREE block 格式见 77 页。

      假设此时你已经有一个使用其他方法得到的系统发育树，那么你可以将其添加到 execute xxx.nex 中 xxx.nex 的结尾：

      1
      2
      3
      begin trees;
      tree usertree = 这里输入你的树
      end;

      需要注意的几点：

      ①、输入的树若为有根树，则需要在其前面加上 [&R] 进行标记。

      ②、输入的树可以有枝长，但一定不能有 bootstrap 值，可以通过正则表达式将其替换掉。

      之后在执行时需输入以下命令：

      1
      2
      3
      > mcmcp nperts=3 append=no ngen=100000 printfreq=1000 samplefreq=1000 nchains=xx nruns=x savebrlens=yes checkpoint=yes checkfreq=5000;
      > startvals tau=usertree;
      > mcmc;

      因为在设定起始树后再重新设置 nruns 和 nchains 会使其失效，因此在 mcmcp 设置参数后再设定起始树，mcmcp 中新增的 nperts 会对起始树进行扰动从而让不同的链起始树具有略微差异，避免不同运行难以检测收敛问题。

  • Mrbayes 有运行示例，可以跟着做一遍。

  • 关于最终收敛的判断，还可以使用 Tracer R 等软件辅助进行。

多基因串联时的分区方法

需准备的软件：ModelFinder（IQTREE）。

首先通过下列命令寻找最佳的分区方案和对应进化模型：

1
2
3
4
$ iqtree2 -s concatenation_ortho.fasta \
-spp IQ_partition.txt \
-m TESTMERGEONLY -mset mrbayes \
-nt AUTO -pre mybayes_aa

此处，-s 指定串联得到的序列文件，-spp 指定对应分区文件，-m 和 -mset 将模型寻找限制在 Mrbayes 所具有的模型中并仅寻找最佳模型不进行建树。如果已有多个基因的比对文件但不知道怎么串联可见之前的文章：

运行结束后，最佳的分区方案和进化模型会出现在以 .best_scheme.nex 结尾的文件中，以上述代码示例为例即为 mybayes_aa.best_scheme.nex。

需要注意的是，该结果并没法直接套用于 Mrbayes 中，因此还需要进一步的转换，对此，我基于 Phylosuite 的源码（详情可见脚本中所给链接）编制了一个专门用于转换的脚本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
# mf2mrbayes.py
# Author: Juse
# Transform result from modelfinder to format used in Mrbayes.
# Usage: python mf2mrbayes.py xxx.best_scheme.nex

import sys
import os

mf_result = open(sys.argv[1],"r")
abs_path = os.path.abspath(sys.argv[1])
op_path = os.path.dirname(abs_path)

numPar = 0
numMol = 0
conPar = ''
conPar_set = []
conMol = ''

# Reference: https://github.com/dongzhang0725/PhyloSuite/blob/master/PhyloSuite/src/Lg_Mrbayes.py
# Line73 - Line124
dict_models = {"JC": "1", "F81": "1", "K80": "2", "HKY": "2", "TrNef": "6", "TrN": "6", "K81": "6",
               "K81uf": "6", "K2P": "2", "JC69": "1", "HKY85": "2", "K3P": "6",
               "TIMef": "6", "TIM": "6", "TVMef": "6", "TVM": "6", "SYM": "6", "GTR": "6", "TPM2": "6",
               "TPM2uf": "6", "TPM3": "6", "TPM3uf": "6", "TIM2ef": "6", "TIM2": "6", "TIM3ef": "6",
               "TIM3": "6"}

for line in mf_result:
    # 处理分区集信息
    if line.strip().startswith('charset'):
        numPar += 1
        ParPos = line.split('=')[1].strip()
        conPar += f'charset subset{numPar} = {ParPos}\n'
        conPar_set.append(f'subset{numPar}')
    # 处理分区模型信息
    elif ':' in line:
        model_name = line.split(':')[0].strip()
        numMol += 1
        if '+G' in model_name and '+I' in model_name:
            rates = " rates=invgamma"
        elif '+G' in model_name:
            rates = " rates=gamma"
        elif '+I' in model_name:
            rates = " rates=propinv"
        else:
            rates = ''
        model_used = model_name.split('+')[0]
        if model_used in dict_models:
            # nt序列
            conMol += f'lset applyto=({numMol}) nst={dict_models[model_used]}{rates};\n'
        else:
            # aa序列
            conMol += f"lset applyto=({numMol}){rates};\n"
            conMol += f"prset applyto=({numMol}) aamodelpr=fixed({model_used.lower()});\n"
            if '+F' in model_name:
                conMol += f"prset applyto=({numMol}) statefreqpr=fixed(empirical);\n"

# 生成分区名
parnam = f'partition Names = {len(conPar_set)}:{", ".join(conPar_set)};\n' \
         f'set partition=Names;\n'

# 汇总输出
total_output = conPar + parnam + conMol
total_output = total_output.replace("amodelpr=fixed(jtt)", "amodelpr=fixed(jones)")
total_output = total_output.replace("  ", " ")

with open(f'{op_path}/Mrbayes_par.txt', 'w') as t:
    t.write(total_output)

print('Finished!')

用法：

1
$ python mf2mrbayes.py mybayes_aa.best_scheme.nex

请将脚本和最佳分区文件替换成自己的路径，脚本在 GitHub 可下载。

Juse 注：

注意！上述脚本中关于模型的选择不完全正确，例如 SYM 模型虽然与 GTR 模型基本相同，但它假设所有核苷酸的频率相等，而 GTR 模型则允许这些频率不等，因此事实上在选用纯粹的 SYM 模型时理应多加一项 statefreqpr=fixed(equal) 进行指定。

虽然说 GTR 相比之下更灵活更复杂，往往能更好的拟合数据，但是在一些情况下这也有可能导致过拟合问题，如果推测出来的最佳模型是 SYM 而不是 GTR，那么我们就有理由去怀疑以上述标准做的正确性。

幸运的是，这个问题仅在核苷酸数据集中出现，如果你使用的是蛋白序列数据，那么就不用担心模型选择的问题，因为这些已经在相对应的参数中直接进行了选择。

运行后，在最佳分区文件的同一目录下会出现名为 Mrbayes_par.txt 的文件，里面记录了 Mrbayes 指定分区和进化模型所需键入的命令，将其替换 运行方法 ② 部分 中第一部分的代码即可。

关于 Mrbayes 的使用体验

Mrbayes 是一个操作较为复杂的系统发育推断软件，此外贝叶斯法本身也需要消耗大量计算资源，所以资源有限的情况下可以有限考虑更高效同时功能强大的最大似然法或其他方法。

此外，对于不同的数据集 Mrbayes 也会存在不同的表现，以个人经验来讲：

  • 对于小数据集，Mrbayes 会以极快的速度收敛，但是相比最大似然法依旧需要消耗更多时间。
  • 就算收敛速度很快（ASDSF 低于 0.01 且 PSRF 接近 1），有效样本大小可能依旧不达标（可能这时所得到的树已经基本正确）。如果追求所有指标的达成可能会不得不添加一定的迭代次数。
  • 对于大数据集，Mrbayes 可能会收敛的非常慢，这种现象随着数据缺失（gap 比例）增多而更明显。

所以后面可能还会出一篇关于 PhyloBayes 的文章，它相比之下就操作更简单许多并且也有一些独特的优势。