619_系统发育分析BEAST2¶

# === 下载安装BEAST2 ===
# 从 https://www.beast2.org/ 下载对应系统的安装包
# Linux:
wget https://github.com/CompEvol/beast2/releases/download/v2.7.8/BEAST.v2.7.8.Linux.x86.tgz
tar -xzf BEAST.v2.7.8.Linux.x86.tgz  # 解压
export PATH=$PATH:$(pwd)/beast/bin  # 添加到PATH

# 验证
beast -version  # 查看BEAST版本
beauti           # 打开BEAUti图形界面

# 安装辅助工具
# Tracer: https://github.com/beast-dev/tracer/releases
# FigTree: https://github.com/rambaut/figtree/releases

FASTA比对 → BEAUti配置XML → BEAST运行MCMC → Tracer检查收敛 → TreeAnnotator总结树 → FigTree可视化

1. 打开BEAUti
2. File → Import Alignment → 选择比对好的FASTA/NEXUS文件
3. 在Partitions标签页确认数据正确

Clock Model标签页:
- Strict Clock（严格钟）: 所有分支速率相同，简单但不现实
- Relaxed Clock Log Normal（松弛钟-对数正态）: 推荐！每条分支速率从对数正态分布中独立抽取
- Relaxed Clock Exponential: 速率从指数分布抽取

Priors标签页:
1. 点击 "+ Add Prior" 添加校准先验
2. 选择树上的节点（MRCA）
3. 设置时间先验分布（如正态分布：均值=10Ma，标准差=2Ma）
4. 至少需要一个校准点来锚定时间尺度

MCMC标签页:
- Chain Length: 链长度（通常1000万-1亿）
- Store Every: 每N步记录一次（通常1000-5000）
- Log Every: 每N步记录日志
- 预估：简单数据10M步，复杂数据100M步

# === 运行BEAST ===
beast \
    -threads 8 \               # 使用8个线程
    -beagle_SSE \              # 启用SSE加速
    -overwrite \               # 覆盖已有结果
    analysis.xml               # BEAUti生成的XML配置文件

# === 后台运行（推荐，因为很耗时）===
nohup beast -threads 16 analysis.xml > beast.log 2>&1 &

# === 恢复中断的运行 ===
beast -resume analysis.xml     # 从上次中断处继续

# 打开Tracer查看.log文件
tracer analysis.log

# 关键检查项：
# 1. ESS (有效样本大小) — 所有参数的ESS必须>200
# 2. Trace plot — 链应该像"毛毛虫"一样上下波动，不应有趋势
# 3. Posterior — 后验概率应平稳
# 4. TreeLikelihood — 似然值应平稳

# === 用Python快速检查ESS ===
import re  # 正则表达式

# 读取BEAST日志的最后几行
with open("analysis.log", 'r') as f:
    lines = f.readlines()

# 跳过注释行
header = None
for i, line in enumerate(lines):
    if not line.startswith('#'):
        header = line.strip().split('\t')
        data_start = i + 1
        break

print("参数列表:")
for j, param in enumerate(header):
    print(f"  {j}: {param}")

# 注意：完整ESS计算需要用Tracer
# 这里只是快速检查链是否在运行
print(f"\n总步数: {lines[-1].split()[0]}")

# === 用TreeAnnotator生成最大可信度树（MCC tree）===
treeannotator \
    -burnin 10 \               # burn-in比例（去掉前10%的样本）
    -heights median \          # 节点高度用中位数
    analysis.trees \           # BEAST输出的树文件
    mcc_tree.nex               # 输出MCC树

# heights选项：
# median — 节点高度取中位数（推荐）
# mean — 取均值
# ca — 共同祖先高度（推荐用于分化时间估计）

library(ggtree)  # 加载ggtree
library(treeio)  # 读取BEAST树

# 读取BEAST的MCC树
tree <- read.beast("mcc_tree.nex")  # 读取

# 画时间标尺树
ggtree(tree, mrsd = "2024-01-01") +  # mrsd: 最近样本日期
    geom_tiplab(size = 3) +  # 叶节点标签
    geom_range("height_0.95_HPD",  # 95% HPD区间
               color = "steelblue", 
               alpha = 0.3,
               size = 2) +
    theme_tree2() +  # 时间轴主题
    scale_x_continuous(
        name = "Time (years ago)"  # x轴标签
    )

ggsave("beast_time_tree.pdf", width = 12, height = 8)

<!-- 松弛分子钟设置示例 -->
<branchRateModel id="RelaxedClock" 
    spec="beast.evolution.branchratemodel.UCRelaxedClockModel"
    rateCategories="@rateCategories" 
    tree="@Tree">
    <LogNormal id="LogNormalDistributionModel" 
        name="distr" 
        meanInRealSpace="true">
        <parameter id="ucldMean" name="M" value="1.0"/>
        <parameter id="ucldStdev" name="S" value="0.3"/>
    </LogNormal>
</branchRateModel>