Mrbayes构建贝叶斯树复习

最近又重新构建多片段的贝叶斯树，对这些参数又有了新的理解，在此记录，以备查询。

多线程安装

以前倒是写过多线程的安装，不过都比较繁琐，参考贝叶斯法构建进化树：MrBayes | 陈连福的生信博客 (chenlianfu.com)和MrBayes安装linux系统下2021-07-29 - 简书 (jianshu.com)，并进行调整，更新了一下多线程的安装方法：

git clone --depth=1 https://github.com/NBISweden/MrBayes.git
cd MrBayes
./configure  --with-beagle=no --with-mpi=yes
make -j 8

这个安装就轻松多了。注意3.2.7a版本是--with-mpi=yes不是--enable-mpi了

参数调试

这里以崔师姐的矩阵为例，一般来说，从fasta转换到nexus我都习惯使用网站ALTER (ALignment Transformation EnviRonment) (sing-group.org)，这个转格式非常稳定，一直都可以用。然后需要在矩阵的后方填入自己需要的指令，以如下所示代码为例：

begin mrbayes;
log start file=4.txt;
Outgroup 150;
Outgroup 149;
Outgroup 148;
charset ITS=1-774;
charset LSU=775-1689;
partition favored = 2: ITS, LSU;
set partition = favored;
Lset applyto=(1,2) nst=6 rates=invgamma;
unlink statefreq=(all) revmat=(all) shape=(all) pinvar=(all);
Prset applyto=(all) ratepr=variable;
mcmc ngen=50000000 samplefreq=10000 printfreq=10000 nchains=4 stoprule=no;
log stop;
end;

开始运行和log文件命名

begin mrbayes;
log start file=cyy.txt

外群设置

Outgroup 150;
Outgroup 149;
Outgroup 148;

每行只能设置一个外群，如需多个外群需要多行设置

定义基因分区

charset ITS=1-774;
charset LSU=775-1689;
partition favored = 2: ITS, LSU;
set partition = favored;

模型参数选择

此部分摘选自：好好先生 —— MrBayes 操作说明 | Bin YE (bin-ye.com)
模型一般在MrModeltest等软件中进行计算筛选，这里就不展开讲了。
[站外图片上传中...(image-c82c19-1696821241551)]
跑完模型，把这里的两行粘贴到分区的下边就行，

Lset  nst=6  rates=invgamma;
Prset statefreqpr=fixed(equal);

需要注意的是，在多个基因片段中需要在Lset和Prset后填写该模型所属的基因顺序，例如ITS的模型就 Lset applyto=(1)...，若二者都适合这个模型，就填Lset applyto=(1,2)，就比如我两个模型可以填成：

lset applyto=(1,2) nst=6 rates=invgamma;
Prset applyto=(1) statefreqpr=fixed(equal);
Prset applyto=(2) statefreqpr=dirichlet(1,1,1,1);

• lset 定义似然模型参数，要计算似然性，就必须假设碱基替换模型，在 MrBayes 中输入命令行 help lset 即可查看相应帮助文件。
• applyto 定义分区的范围，即将模型应用到对应的分区，在括号 () 中列出应用相同模型的分区。
• nst 定义位点替换的模型（实际可以理解为 AC、AT、AG、GC、GT、TC 每对碱基的替换），用不同数字代替不同模型：1 代表所有替换率相同（如模型 JC69 或 F81）； 2 代表转换和颠换可以有不同的替换率（如模型 K80 或 HKY85）； 6 代表所有替换率都不同（如模型 GTR）。
• rates 定义位点之间的替换率，有以下几种选择： equal（位点的替换率无差异）、 gamma（位点的替换率呈 gamma 分布）、 adgamma（位点的替换率自相关，边缘位点替换率呈 gamma 分布，相邻位点有相关的替换率）、 propinv（一定比例位点的替换率是恒定的）、 invgamma（一定比例位点的替换率是恒定的，剩下位点的替换率呈 gamma 分布）。
• prset 命令可定义一系列系统发育模型的先验，设置先验参数分布是贝叶斯分析的必不可少的一步，先验数据代表观测数据之前的预先判断。详细帮助文件可在 MrBayes 中输入命令行 help prset 查看。本文也具体对此作了中文翻译说明，详见 附章 2。常用的先验设置主要有 revmatpr, shapepr, pinvarpr, topologypr, brlenspr, popvarpr, clockratepr, clockvarpr, ratepr。

unlink revmat=(all) pinvar=(all) shape=(all) statefreq=(all);

另一个参数： unlink 命令是让每个分区序列所使用的模型都不连锁，每个分区单独应用各自的模型或参数。一般情况下，分区是根据形态或核酸等数据类型进行的，也可根据自定义的分区（如前述 set partion=<name/number>）。软件运行时，如果应用于不同分区的参数有相同的先验，那么软件将会将这些参数连锁起来，也就是说，软件将会使用同一个参数，软件默认运行就是这样“吝啬”。但是，如果你想对各个分区使用不同的参数，比如，设置不同分区有不同的转换或颠换比率，那么就需要使用这一 unlink 命令。可在 MrBayes 中输入命令行 help unlink 查看相应帮助文件及一系列参数。
更多与模型相关的信息可以去源博中查看，他写了挺多的，但是我这里都没用，也许读者可以用上。

运行参数

此部分摘选自：贝叶斯建树之 Mrbayes 篇 | Juse's Blog (biojuse.com)

mcmc ngen=2000000 samplefreq=100 printfreq=100 nchains=4 stoprule=yes stopvar=0.01;

• append=no，表示新开始一个分析（no也可以不写），如果之前已经运行过但因为某种原因中断，可以将其调整为 append=yes 从上一个检查点继续运行，需要搭配 checkpoint=yes。
• ngen=2000000: 分析总迭代次数200万次，一般需要数百万代运算甚至更多
• printfreq=1000，表示每一千代打印一次链状况到屏幕输出上。
• samplefreq=1000，表示每一千代采一次样，对于需要长时间才能收敛的大数据集而言可以将该值调高以降低采样频率，避免最终文件包含的树和参数信息过多。采样量为 ngen/samplefreq。
• nchains=10，表示每个运行使用十条链，总链数为 nchains*nruns。
• nruns=2，表示进行两个独立的分析，用于判断收敛情况。
• savebrlens=yes，表示记录分支长度（yes也可不写）。如果只关注树的拓扑结构且不需分支长度信息可以设置为 no 节省空间。
• checkfreq=5000，表示每五千代进行一次收敛诊断（默认就是5000，可以不写，或者自己调一下），输出Average standard deviation of split frequencies (ASDSF)
  **[站外图片上传中...(image-80291f-1696821241551)]

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

结束指标

此部分摘选自：贝叶斯建树之 Mrbayes 篇 | Juse's Blog (biojuse.com)
在运行中及结束以后，有几个判断分析收敛程度的指标：

• 平均标准分歧（ASDSF），低于 0.01 时表明收敛良好。
• 潜在似然（PSRF），在 sump 输出中出现，理想情况应该接近 1，若不是则说明可能仍未收敛。
• 有效样本大小（ESS），在 sump 输出中出现，通常情况下以大于 200 作为样本量足够的标准（PhyloBayes 则为 300）。这个可以将输出的***.run1.p文件拖进tracer软件中查看，ESS的大小。
  一般来说，三个指标都达标时我们才有充分的信心相信贝叶斯推断的结果足够可靠，以下是得到后两种信息的方法：

sump burninfrac=0.25
sumt burninfrac=0.25

这样树就基本上计算完了，最终树*.con.tre可以通过Figtree等软件打开。

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

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

不收敛问题

此部分摘选自：贝叶斯建树之 Mrbayes 篇 | Juse's Blog (biojuse.com)

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

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

需要注意的几点：

• 输入的树若为有根树，则需要在其前面加上 [&R] 进行标记。
• 输入的树可以有枝长，但一定不能有 bootstrap 值，可以通过正则表达式将其替换掉。

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

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 多核心并行运算效率 (douban.com)
[图片上传失败...(image-90a43b-1696821241551)]
在我看来多线程还是会和单线程的结果有一些不同，但是如果树结构非常稳定，应该是问题不大的。