使用示例
需要使用PHYLIP格式的数据格式

# 注：bsub是集群提交命令的方式，本地服务器使用不需要添加
# 提交作业
bsub -n 8 mpirun -n 8 ./pb -cat -mtart -d ~/phylobayes3.3b/combined_aa.phy ~/phylobayes3.3b/test/aaa1 &
bsub -n 8 mpirun -n 8 ./pb -cat -mtart -d ~/phylobayes3.3b/combined_aa.phy ~/phylobayes3.3b/test/aaa2 &
# 比较两个链
bpcomp aaa1 aaa2
# 输出结果
initialising random
seed was : 791961

aaa1.treelist : 1630 trees
aaa2.treelist : 3129 trees

maxdiff     : 0.253139
meandiff    : 0.0146816

bipartition list in : bpcomp.bplist
consensus in        : bpcomp.con.tre

Maxdiff<0.1时最好。

# 停止计算
stoppb aaa1
stoppb aaa2

# 继续计算
bsub -n 8 -m "node44" mpirun -n 8 pb ~/phylobayes3.3b/test/aaa1 &
bsub -n 8 -m "node44" mpirun -n 8 pb ~/phylobayes3.3b/test/aaa2 &

示例2，自动停止计算的方法，运行4个链

# 模型选择，混合模型，有效防止组成异质性位点造成的假象
DNA或者SR4 recode：iqtree(C60SR4);phylobayes(CAT+GTR)
protein:iqtree(LG+pmsf(C60)+F+G4);phylobayse(CAT+GTR+G4/CAT+LG+G4)

# 提交作业，10个核心运行，-x 1 100 每个树都保存，直到有100个数就停止
nohup mpirun -np 10 ~/pb_mpi -d sequence.phy -x 1 100 -gtr -cat -dgam 4 tree1 > tree1.log &
nohup mpirun -np 10 ~/pb_mpi -d sequence.phy -x 1 100 -gtr -cat -dgam 4 tree2 > tree2.log &
nohup mpirun -np 10 ~/pb_mpi -d sequence.phy -x 1 100 -gtr -cat -dgam 4 tree3 > tree3.log &
nohup mpirun -np 10 ~/pb_mpi -d sequence.phy -x 1 100 -gtr -cat -dgam 4 tree4 > tree4.log &

# 停止运行
echo 0 > chainname.run

## 测试是否收敛和有效样本大小
# bpcomp:用于评估树空间中的收敛性
# -x 1000 10：burn-in老化，去除前1000个不稳定的树；10:每隔10树取一个样本树
bpcomp -x 1000 10 <chain1> <chain2>  
# 输出结果
1. bpcomp程序将输出在所有双分割中观察到的最大(maxdiff)和平均(meandiff)差异。
2. 生成一致树文件bpcomp.con.tre
maxdiff < 0.1: good run.
maxdiff < 0.3: acceptable, 给出了一个很好的一致性树。
0.3 < maxdiff < 1:样本量不够，链不收敛，但是模型和数据集没有错，可以继续运行
maxdiff = 1：即使在10000个循环后仍是该结果，表明至少有一个运行卡在了局部最大值中。

# tracecomp:用于检验模型连续参数的收敛性
tracecomp -x 1000 <chain1> <chain2>
# 输出结果
1. 总结差异和样本的有效大小
rel diff < 0.1 and minimum effective size > 300: good run;
rel diff < 0.3 and minimum effective size > 50: acceptable run.

# 注意
1. 获得一个好的结果，通常需要10000到30000个循环，与数据集大小有关
2. 客观的评价标准： 有效样本大小(effective sample size)和结果再现性(reproducibility of the result)
3. 当使用CAT或者CAT-GTR等复杂混合模型时，应确保后验一致树在独立运行中也能恢复，同时，记录在跟踪文件中的汇总统计数据的跟踪图捕获模型的各个子组件(树长度、alpha参数、无限混合物的占用组件数、混合物的熵、可交换性熵)似乎是平稳的，并且在运行中是可重复的。

# 停止参数
-x   every until 
指定保存频率和(可选的)链应停止的点数。如果未指定此数字，则链将“永远”运行。根据定义，-x 1对应于默认的保存频率。在某些情况下，样本可能是强相关的，在这种情况下，如果磁盘空间或访问受限，那么减少保存点的频率(比如减少10倍)是有意义的:为此，可以使用-x 10选项。

官方说明文档部分

3.3 获得后验一致树和参数估计

MCMC采样器在平衡状态下采样的所有树的一致性(即后验一致性树的MCMC估计)通常被用作系统发生树的点估计。这种(多数规则)共识树是由bpcomp程序自动生成的(见上文)。请注意，bpcomp可以在单个链上运行(在这种情况下，它将在burn-in后简单地生成所有树的一致性树)，而不一定在多个链上运行(在这种情况下，如上所述，它将在所有链上汇集所有树的一致性树，并计算跨链的差异度量)。在多链上使用bpcomp通常会得到更稳定的后验共识树MCMC估计

readpb_mpi程序用于估计模型的一些关键参数，执行后验预测分析，计算后验平均似然和交叉验证分数。

默认情况下，readpb_mpi只计算树的总长度(整个系统发育中每个位点的替换总数)和跨位点率的离散伽玛分布的alpha参数的简单估计(平均值和95%可信区间)。readpb_mpi执行的所有其他任务都可以通过specic选项访问(请参阅下一节中readpb_mpi的详细选项)。

5. 后验概率模型检验--Posterior predictiv model checking

简述

检验模型的充分性是贝叶斯数据分析的总要步骤。一个适当的模型应该能够很好地预测在实际数据中通常观察到的模式。在目前的情况下，该模型应该正确地预测那些被怀疑与系统发育重建特别相关的模式。特别是，可接受的核苷酸或氨基酸的位点特异性限制，或物种之间的组成差异，如果没有正确建模，可能会导致树重建中的系统错误。因此，检验我们的模型是否正确地预测了跨位点和跨分类群的典型组成变化模式是很重要的。

在贝叶斯推理中，这种类型的模型检查可以使用后验预测模拟来完成。后验预测检查可以看作是参数自举的贝叶斯模拟：一旦模型以经验数据为条件，则使用由此估计的参数配置来重新模拟数据(多次)。然后，在模拟的复制上计算一些感兴趣的汇总统计值(用于捕获模型应该正确捕获和复制的关键模式)的值，从而产生模型下统计值的零分布。然后将在原始数据上计算的统计值与此零分布进行比较。通常，计算p-value，将其定义为测试统计量的值至少与观测值一样极端的模拟重复的百分比。作为后验预测p-value的一个很好的补充，z-score也可以被考虑：观察值与零分布的平均值之间的偏差，相对于零分布的标准差。在MCMC对p-value的估计等于0的情况下，z-score是有用的。除了模型检查，后验预测模拟也可以被认为是一种有原则的方法，可以根据经验数据进行校准模拟。

5.1 应用示例

后验预测模拟和检查可以基于MCMC采样器的输出来完成——对于平稳采样的一系列点中的每一个，根据相应的参数配置重新模拟新的数据。这里以ef2.ali数据集为例，在ef2.tree文件中给出的固定树拓扑下。将比较两种模型，GTR和CAT-Poisson，用于多样性统计(测量在每个位点观察到的不同氨基酸的平均数量)。

第一步是在GTR和CAT模型下运行mcmc链:

pb_mpi -d ef2.ali -T ef2.tree -dp -poisson -x 3 1100 catef2
pb_mpi -d ef2.ali -T ef2.tree -gtr -ncat 1 -x 1 1100 gtref2

然后，一旦运行完成，我们可以使用以下命令对这两条链运行后验预测测试:

readpb_mpi -x 100 1 -ppred -div catef2
readpb_mpi -x 100 1 -ppred -div gtref2

结果将写入扩展名为.div:的文件中。catef2.div和gtref2.div
在目前的情况下，根据GTR模式，这应该是:

diversity test
obs div : 3.6571
mean div: 4.46541 +/- 0.0788579
z-score : 10.2503
pp : 0

CAT模型下:

diversity test
obs div : 3.6571
mean div: 3.75094 +/- 0.0579987
z-score : 1.61798
pp : 0.053

可以解释如下:

• 原始经验序列比对平均每列3.66个氨基酸;
• 在GTR模型下模拟的比对，以及在参数值下EF2序列估计的比对，平均每列有更高的氨基酸数(4.47);与原始数据集的差异非常显著(z-score为10.25,pp≈0)。绝对差异可能看起来很小，但存在潜在的异质性：与进化较慢的位点相比，进化速度较高的地点在经验和模拟之间往往表现出更大的差异。
• 相比之下，在CAT模型下模拟的比对结果平均每个位点有3.75个氨基酸，因此更接近于经验序列比对结果。CAT诱导的多样性与观测的差异不显著(pp = 0.053)。