PAML — 用最大似然法分析序列分子进化的经典工具¶

# 方法1：conda 安装（推荐）
conda install -c bioconda paml  # 安装 PAML（当前版本 4.10.x）
paml --help  # 确认安装

# 方法2：从源码编译
wget http://abacus.gene.ucl.ac.uk/software/paml4.10.8.tgz  # 下载源码
tar -xzf paml4.10.8.tgz          # 解压
cd paml4.10.8/src/               # 进入源码目录
make -f Makefile                 # 编译

# PAML 包含多个程序：
# codeml  — 蛋白质/密码子进化（最常用，计算 dN/dS）
# baseml  — DNA 碱基进化
# mcmctree — 贝叶斯分歧时间估算
# yn00    — Yang & Nielsen 2000 dN/dS 计算方法

# 1. 准备输入文件
#    aln.phy   — Phylip 格式的密码子多序列比对
#    tree.nwk  — Newick 格式的进化树（需事先建好）

# 2. 创建 codeml 控制文件 codeml.ctl
cat > codeml.ctl << 'EOF'
      seqfile = aln.phy          * 序列文件（Phylip 格式密码子比对）
     treefile = tree.nwk         * 进化树文件
      outfile = codeml_output.txt * 输出结果文件

        noisy = 9                * 屏幕输出详细程度（0-9）
      verbose = 1                * 详细输出模式
      runmode = 0                * 0=给定拓扑树，-2=成对比较

      seqtype = 1                * 1=密码子，2=氨基酸
    CodonFreq = 2                * 密码子频率模型（2=F3x4）
        model = 0                * 分支模型（0=所有分支相同 ω）
      NSsites = 0                * 位点模型（0=关闭；7 8=M7 vs M8 检测正选择）

        icode = 0                * 遗传密码（0=标准密码）
        Mgene = 0                * 多基因模型（0=单基因）
    fix_kappa = 0                * kappa 是否固定（0=估算）
        kappa = 2                * 转换/颠换率初始值
    fix_omega = 0                * ω 是否固定（0=估算）
        omega = 0.4              * dN/dS 初始值
EOF

# 3. 运行 codeml
codeml codeml.ctl  # 运行，结果写入 codeml_output.txt

# --- 运行 M7 模型（零假设，无正选择）---
cat > m7.ctl << 'EOF'
      seqfile = gene_aln.phy
     treefile = gene_tree.nwk
      outfile = m7_result.txt
      seqtype = 1
    CodonFreq = 2
        model = 0
      NSsites = 7              * M7 模型
        omega = 0.5
EOF
codeml m7.ctl  # 运行 M7

# --- 运行 M8 模型（备择假设，含正选择）---
cat > m8.ctl << 'EOF'
      seqfile = gene_aln.phy
     treefile = gene_tree.nwk
      outfile = m8_result.txt
      seqtype = 1
    CodonFreq = 2
        model = 0
      NSsites = 8              * M8 模型
        omega = 0.5
EOF
codeml m8.ctl  # 运行 M8

# --- 用 Python 提取对数似然值并做 LRT ---
python3 << 'PYEOF'
import re, scipy.stats

def extract_lnL(outfile):
    """从 codeml 输出文件中提取对数似然值"""
    with open(outfile) as f:
        content = f.read()
    # 匹配 "lnL(ntime:  ... np: ...) = -XXXXX.XXXX"
    match = re.search(r'lnL\(.*?\)\s*=\s*([-\d.]+)', content)
    return float(match.group(1)) if match else None

lnL_M7 = extract_lnL("m7_result.txt")   # 提取 M7 对数似然
lnL_M8 = extract_lnL("m8_result.txt")   # 提取 M8 对数似然

# 似然比统计量 = 2 * (lnL_M8 - lnL_M7)
LRT = 2 * (lnL_M8 - lnL_M7)
# M7 vs M8 自由度差为 2
p_value = scipy.stats.chi2.sf(LRT, df=2)  # 卡方检验 p 值

print(f"M7 lnL: {lnL_M7:.4f}")
print(f"M8 lnL: {lnL_M8:.4f}")
print(f"LRT 统计量: {LRT:.4f}")
print(f"P 值: {p_value:.6f}")
print("结论:", "存在正选择 (p<0.05)" if p_value < 0.05 else "无显著正选择")
PYEOF

# 安装
conda install -c bioconda paml

# 主程序
codeml   # 密码子/蛋白质进化分析（最常用）
baseml   # DNA 碱基替换分析
mcmctree # 贝叶斯分歧时间估算（需化石校正点）
yn00     # 快速 dN/dS 成对计算

# 常用 NSsites 组合（LRT 检验正选择）
NSsites = 1 2   # M1a vs M2a
NSsites = 7 8   # M7 vs M8（推荐）
NSsites = 8 8a  # M8 vs M8a（验证）

# 读取 BEB 正选择位点（M8 输出）
grep -A 50 "Bayes Empirical Bayes" m8_result.txt | head -60
# 后验概率 > 0.95 的位点为正选择位点

# dN/dS 解读
# ω < 1 → 纯化选择（保守进化）
# ω = 1 → 中性进化
# ω > 1 → 正选择（加速进化）