用于beast并行化运行的beagle的安装

此处的beagle和用于填充vcf文件缺失基因型的beagle不是同一个软件。

# 下载源码
git clone https://github.com/beagle-dev/beagle-lib.git
cd beagle-lib
mkdir build && cd build

# 配置编译（安装到路径/share/home/XXX/soft/beast2/beagle）
# 清理之前的配置
cd ~/soft/beast2/beagle-lib/build
rm -rf *
##编译使用cpu版本的beagle
cmake .. \
  -DCMAKE_INSTALL_PREFIX=/share/home/XXX/soft/beast2/beagle \
  -DBUILD_SHARED_LIBS=ON \
  -DBUILD_CUDA=OFF \
  -DBUILD_JNI=ON \
  -DBUILD_OPENCL=OFF \
  -DCMAKE_CXX_COMPILER=/share/home/XXX/soft/gcc/8.3.0/bin/g++ \
  -DCMAKE_C_COMPILER=/share/home/XXX/soft/gcc/8.3.0/bin/gcc 
# 编译和安装
make -j8
make install
# 把下面的一行内容添加到~/.bashrc
export LD_LIBRARY_PATH=/share/home/XXX/soft/beast2/beagle/lib:$LD_LIBRARY_PATH
source ~/.bashrc

上面安装的是cpu版本的beagle,我在此处指定了编译的gcc和g++的绝对路径。
如果你的服务器有GPU的话，可以使用下面的命令编译支持显卡的beagle

#编译使用GPU版本的beagle 
cmake .. \
  -DCMAKE_INSTALL_PREFIX=/share/home/XXX/soft/beast2/beagle \
  -DBUILD_SHARED_LIBS=ON \
  -DBUILD_CUDA=ON \
  -DBUILD_JNI=ON \
  -DCUDA_TOOLKIT_ROOT_DIR=/share/home/XXX/CUDA11.8 \
  -DBUILD_OPENCL=OFF \
  -DCMAKE_CXX_COMPILER=/share/home/XXX/soft/gcc/8.3.0/bin/g++ \
  -DCMAKE_C_COMPILER=/share/home/XXX/soft/gcc/8.3.0/bin/gcc

编译cuda版本的beagle报错fatal error: libhmsbeagle/GPU/kernels/BeagleCUDA_kernels.h: No such file or directory，官方github说是gcc的版本不对，我此处修改为8.3.0版本可以成功编译。
注意cmake指定参数后的所有路径都要使用绝对路径，不要使用相对路径，例如~,这个会报错。

准备xml文件(在windows上使用BEAUti)

现在假设你已经拥有比对好的多个物种的单拷贝序列的fa文件。
如果你的序列很长，例如超过10K，而且物种很多，则使用BEAUti会很慢。
所有我使用的魔法工具来完成。

截取前1K的序列

#提取序列的前1-1000个碱基用于测试数据
seqkit subseq -r 1:1000 All.min4.fasta|seqkit seq -w 0  >test.fa

test.fa导入BEAUti，进行设置，获得xml文件。
打开 BEAUti2（通常位于 BEAST2 安装目录）。
导入序列文件：
File → Import Alignment → 选择 test.fasta。
设置分区（Partitions）：
如果数据有多个基因/位点，可定义分区模型，此处选择 nucleotide（DNA序列）或aminoacid （氨基酸序列）

选择进化模型(Model site)：
Site Model → 选择 Substitution Model（如 HKY、GTR）。
设置 Gamma 或 Invariant Sites（如 +G 或 +I）。
此处可以使用iqtree估计的最佳模型也可以使用beast的拓展包bModelTest做核酸模型的自动选择。

设置分子钟（Clock Model）：
Strict Clock（严格分子钟）或 Relaxed Clock（松弛分子钟，如 Log Normal）。实际使用Log Normal对数松弛分子钟

设置树先验（Tree Prior）：
Coalescent（种群历史）或 Yule（物种分化）。

设置 MCMC 参数：
一个粗略的计算方法是总代数=3000*(序列条数的平方)。
这里有100条序列，就是30,000,000代
Chain Length（链长，通常 10,000,000 次迭代）。

样本数量=Chain Length /采样频率，要求样本数量至少有10,000个。
所以此处采样频率最大为30,000,000/10,000=3000，此处实际设置为1000
Log Parameters（采样频率，如每 1000 代记录一次）。
保存 .xml 文件：
File → Save → 生成 test1K.xml。

使用我自定义的脚本fa_to_xml.py
对于大的fa文件，可以使用我自己编写的脚本fa_to_xml.py,把fa序列转为需要的xml格式，使用这个脚本后可以节省几天的时间。

python3 fa_to_xml.py All.min4.fasta test1K.xml test.xml

fa_to_xml.py脚本需要3个参数，

参数1是完整的输入的fa文件
参数2是输入上面使用BEAUti生成的xml的模板文件
参数3是输出xml的文件名
经过此处的修改后，输出的test.xml文件即是下一步真正分析的xml文件。

使用beast2进行分析

如果你只安装了cpu版本的beagle

beast -beagle -threads 24 test.xml

强制使用cpu的beagle

beast -beagle_CPU -threads 24 test.xml

强制使用GPU的beagle，放心增大使用的核心数，GPU的核心数比cpu多的多

beast -beagle_GPU -beagle_order 1,2 -threads 400 -overwrite  All.xml

beast常用的参数

-overwrite参数用于覆盖之前断掉的输出文件
-threads 40 指定使用的cpu数量
-beagle/-beagle_CPU/-beagle_GPU 指定使用beagle进行加速，如果不指定，则只能使用1个cpu，会非常慢
-D chainLength=1000000 指定新的chainlength,替换xml中的值
-prefix 字符串指定输出文件的前缀
-seed 9527 指定随机数，不同批次的分析指定不同的随机数种子和输出前缀，这样后面可以把多次运行的log和tree合并到一起。
注意：我这里使用了参数-beagle_order 1,2来指定使用编号为1和2的GPU显卡。
beast -beagle_info命令会显示，你当前环境中可用的beagle的资源，我的是显示如下：说明有2个GPU,所以可以指定使用order 1,2.

0 : CPU (x86_64)
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_SSE VECTOR_NONE THREADING_CPP THREADING_NONE PROCESSOR_CPU FRAMEWORK_CPU

1 : NVIDIA GeForce RTX 4090
    Global memory (MB): 24210
    Clock speed (Ghz): 2.52
    Number of cores: 16384
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_CPP THREADING_NONE PROCESSOR_GPU FRAMEWORK_CUDA PARALLELOPS_STREAMS PARALLELOPS_GRID

2 : NVIDIA GeForce RTX 4090
    Global memory (MB): 24210
    Clock speed (Ghz): 2.52
    Number of cores: 16384
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_CPP THREADING_NONE PROCESSOR_GPU FRAMEWORK_CUDA PARALLELOPS_STREAMS PARALLELOPS_GRID

如果你也遇到报错CUDA error: "Out of memory" (2) from file </beagle-lib/libhmsbeagle/GPU/GPUInterfaceCUDA.cpp>, line 596. 说明你的xml的数据太大了，显卡资源不足。

会输出.log文件，

使用Tracer评估上面的MCMC的结果

使用Tracer github评估结果是否收敛。
这个是在windows端的软件。

image.png

File > import 导入上面输出的log文件即可
我此处的示例因为跑的mcmc次数太低，所以ESS的值很小，并不是可以使用的结果。ESS >200, 才表示采样充分。需要增大chain length

合并后验树（TreeAnnotator）

打开 TreeAnnotator设置：
Burnin（丢弃前 10% 的树）。
Posterior Probability Limit（如 0.5，仅保留高支持率节点）。
选择 Maximum Clade Credibility Tree（最大分支可信树）。
输出 .tree 文件（可用于可视化）。
或者是linux命令行操作
treeannotator -burnin 30 -limit 0.5 test-test.trees test.trees

绘制得到的树文件可视化，可以使用Figtree或beast2自带的DensiTree或Y叔的ggtree包。

调参优化ESS值

LogCombiner 是一个用于合并多个独立运行的MCMC链（如重复运行或不同种子下的结果）或调整日志文件采样频率的工具，可有效提高ESS值（Effective Sample Size）。
如果你的ESS值始终不理想，可以使用Beast2自带的一个程序LogCombiner来优化log文件和tree文件。
注意：前提条件是你生成所有的log的xml参数除了随机数不同之外，其他的都必须一样。否则无法合并。

合并多个log文件

logcombiner -log run1.log -log run2.log -log run3.log -o combined.log -b 10 -resample 5000

参数说明：

-b 丢弃的burn-in比例默认是10，即10%。
-log 指定输入log文件
-resample 重新采样频率（步数）合并后如果文件太大了，可以设置稍大的采样数，比如从1000调整到5000
-o 输出log文件名

合并tree文件

logcombiner -log run1.trees -log run2.trees -log run3.trees  -o combined.trees -b 10 -resample 5000

合并后再使用Tracer查看log文件，看ESS值是否大于200。

beast的插件安装

beast支持win,linux系统。
BEAUti是其中一个程序，如果linux没有X11，没有的话没法用窗口模式操作。这一步可以在windows上进行操作。

beast的windows和linux的插件是相同的。
例如：安装使用用于SNP数据分析SNAPP插件
windows端点击左侧的File > manage package 就可以看到下面的界面，点击下面的?就可以出现你的默认的包的安装位置了。

image.png

上面显示我的包的位置是C:\Users\XXX\BEAST\2.7

windows安装SNAPP包

windows你可以直接选择对应的包，然后点击下面的Install/Update按钮即可安装。如果安装失败，可能由于网络问题，则可以直接去下载对应的zip压缩包，解压缩后放入上面的包的路径中。
linux的位置是~/.beast/2.7/
如果你有多个版本的beast,比如2.7.6和2.7.7，实际使用时会优先使用最高版本的beast.

linux安装SNAPP包

packagemanager -add SNAPP

这里的packagemanager是beast2自带的脚本，如果上面还是报错，会提示无法访问网络信息。
进入下面的路径https://raw.githubusercontent.com/CompEvol/CBAN/master/packages2.7.xml ，去里面找到对应的beast2版本的SNAPP的下载路径。

mkdir ~/.beast/2.7/SNAPP
cd ~/.beast/2.7/SNAPP
wget -c https://github.com/BEAST2-Dev/SNAPP/releases/download/v1.6.0/SNAPP.addon.v1.6.0.zip
unzip SNAPP.addon.v1.6.0.zip
rm -rf SNAPP.addon.v1.6.0.zip

windows同样把上面的zip解压到文件夹SNAPP,然后把这个文件夹放到C:\Users\XXX\BEAST\2.7路径下。

下面以使用命令行的SNAPP分析SNP数据为例

##删除原有的beauti.properties
rm -rf ~/.beast/2.7/beauti.properties
##再次运行packagemanager查看安装情况
packagemanager -list

#从github下载ruby的脚本
wget https://raw.githubusercontent.com/mmatschiner/snapp_prep/master/snapp_prep.rb
#生成xml文件
ruby snapp_prep.rb -a SNAPPER -v NC_031969.f5.sub5.vcf -t individuals.txt -c constraints.txt -m 1000 -l 100000

snapp_prep.rb的命令参数

-a 指定分析时使用的是SNAPP或SNAPPER,默认是SNAPP
-v vcf文件
-t txt文本文件，第1列是物种信息，第2列是样本信息，tab分割符号
-c 化石证据年龄信息文件
-l MCMC的代数，默认是500000
-m 最多使用多少SNP.无默认值
-x 输出的xml文件的前缀
-o 输出的tree和log文件的前缀

individuals.txt的文件格式如下（所有的样本名都需要在这个文件里，tab分割）

species individual
物种名称1   样本名1
物种名称2   样本名3
物种名称1   样本名2

constraints.txt （指定物种的分化时间，中间的19数字代表19MYA）

  lognormal(0,19,0.16)  crown   物种名称1,物种名称2,物种3

Tracer的界面详细解析

image.png

使用Tracer导入log文件之后的界面如右侧。

点击屏幕右侧上面的Trace
可以看到入上面的分布图，横坐标是训练的chain length,可以看出随着训练的次数的增加，左侧对应参数的变化。左侧有对应参数的ESS值，ESS小于100是红色，100<ESS<200 是黄色，ESS >200是正常的黑色。黑色的ESS值才表示模型可信。如果ESS较小，算出的95% HPD的区间会非常大。
2 .点击屏幕右侧上面的Estimates

image.png

此时能看到分布的柱状图和上面的表格文件
能够看到95% HPD interval和ESS,均值和中间值等。
STDEV标准差，
number of samples 样本数，这个需要达到至少10000才可信。
如果你有多次运行的log文件，也可以点击左侧的+导入进来

image.png

现在我们就可以看到combine后的分布，可以看到我的这个是咧，ESS值才51，右侧的Trace分布图，后面还出现了断层，说明还需要加大训练次数，即增加chain length。
一般要求的ESS>200,才被接受，如果某些杂志要求严格，可能要ESS>400.

增大ESS的方法

增加采样频率()
增加chain的长度（chain length）
多个log文件合并分析
如果你反复运行同一个模型，比如超过10亿次，还不收敛（ESS<200）,此时建议你换一个模型试试。不同的模型的log文件无法合并。

物种分化时间估计工具beast2

物种分化时间估计工具beast2

用于beast并行化运行的beagle的安装

准备xml文件(在windows上使用BEAUti)

fa_to_xml.py脚本需要3个参数，

使用beast2进行分析

使用Tracer评估上面的MCMC的结果

合并后验树（TreeAnnotator）

绘制得到的树文件可视化，可以使用Figtree或beast2自带的DensiTree或Y叔的ggtree包。

调参优化ESS值

合并多个log文件

合并tree文件

beast的插件安装

windows安装SNAPP包

linux安装SNAPP包

下面以使用命令行的SNAPP分析SNP数据为例

individuals.txt的文件格式如下（所有的样本名都需要在这个文件里，tab分割）

constraints.txt （指定物种的分化时间，中间的19数字代表19MYA）

Tracer的界面详细解析

增大ESS的方法

推荐阅读更多精彩内容