分子进化分析 — PAML dN/dS 选择压力分析¶

一句话概述：PAML 的 codeml 程序通过计算非同义替换率(dN)与同义替换率(dS)的比值 ω=dN/dS，来检测蛋白编码基因是否经历了正向选择。

核心知识点速查表¶

概念	白话解释
dN（非同义替换率）	导致氨基酸改变的突变速率（改变了蛋白质）
dS（同义替换率）	不改变氨基酸的突变速率（沉默突变）
ω = dN/dS	ω<1 纯化选择；ω=1 中性进化；ω>1 正向选择
codeml	PAML 中分析密码子替代模型的核心程序
位点模型（Site Model）	假设不同密码子位点有不同的 ω
分支模型（Branch Model）	假设不同谱系（分支）有不同的 ω
分支-位点模型	同时考虑分支和位点的差异（最灵活）
LRT	似然比检验，比较两个模型谁更好
BEB	贝叶斯经验贝叶斯，识别具体哪些位点受正向选择

一、dN/dS 的原理（白话版）¶

比喻：想象 DNA 密码子是一个三个字母的密码锁。

同义突变（dS）：锁的密码改了但开的还是同一个门（氨基酸没变）→ 这种突变通常是"中性"的
非同义突变（dN）：锁的密码改了，门也换了（氨基酸变了）→ 这种突变会影响蛋白功能

ω = dN/dS 的含义： - ω < 1：非同义突变比同义突变少 → 蛋白功能重要，突变被淘汰（纯化选择） - ω = 1：两种突变速率相同 → 蛋白不受选择约束（中性进化） - ω > 1：非同义突变比同义突变多 → 蛋白功能在快速变化（正向选择）

二、分析流程¶

整体流程¶

编码序列 → 多序列比对（蛋白指导） → 系统发育树 → codeml分析 → LRT检验 → BEB识别位点

第 1 步：安装 PAML¶

# 方法1：conda安装（推荐）
conda install -c bioconda paml  # 安装PAML

# 方法2：源码编译
git clone https://github.com/abacus-gene/paml.git  # 克隆仓库
cd paml/src
make -f Makefile                # 编译
export PATH=$PATH:$(pwd)        # 添加到PATH

# 验证安装
codeml  # 会提示需要控制文件，说明安装成功

第 2 步：准备输入文件¶

# PAML 需要三个文件：
# 1. 序列比对文件（PHYLIP格式）
# 2. 树文件（Newick格式）
# 3. 控制文件（.ctl文件）

# 比对编码序列（关键：先比对蛋白质，再反推到核苷酸）
# 用 MAFFT 比对蛋白质
mafft --auto protein.fasta > protein_aligned.fasta  # 比对蛋白序列

# 用 PAL2NAL 反推到密码子比对
# PAL2NAL 确保核苷酸比对按密码子对齐（每3个碱基是一个单位）
pal2nal.pl protein_aligned.fasta \   # 蛋白比对
           nucleotide.fasta \         # 对应的核苷酸序列
           -output paml \            # 输出PAML格式
           > aligned_codon.phy       # 输出文件

# 检查PHYLIP格式
head aligned_codon.phy
# 第一行：序列数  位点数
# 后面每行：序列名  序列（密码子对齐的核苷酸）

# 准备树文件（Newick格式）
# 可以用 IQ-TREE 先建一棵树
iqtree2 -s aligned_codon.phy -m GTR+G4 -bb 1000 -nt AUTO
# 然后把 .treefile 复制为 tree.nwk

# 重要：PAML 的树文件第一行是 "物种数  树数"
echo "10  1" > tree.nwk       # 假设10个物种，1棵树
cat output.treefile >> tree.nwk  # 追加Newick格式的树

# 对于分支模型和分支-位点模型，需要在树上标记前景分支
# 用 #1 标记要检测的分支（前景分支）
# 例如：((A,B),(C #1,D));  ← C 这个分支被标记为前景

第 3 步：位点模型（检测哪些位点受正向选择）¶

# 位点模型比较：M1a(近中性) vs M2a(正选择)；M7(beta) vs M8(beta+ω>1)

# === M1a 模型（零假设：没有正向选择） ===
cat << 'EOF' > codeml_M1a.ctl
      seqfile = aligned_codon.phy   * 序列文件
     treefile = tree.nwk            * 树文件
      outfile = M1a_results.txt     * 输出文件

        noisy = 9                   * 输出详细程度（0-9）
      verbose = 1                   * 输出详细结果
      runmode = 0                   * 0=使用用户提供的树

      seqtype = 1                   * 1=密码子序列
    CodonFreq = 2                   * 2=F3x4模型（从数据估计）
        clock = 0                   * 0=无分子钟（分支速率可不同）
        model = 0                   * 0=所有分支相同的ω
      NSsites = 1                   * 1=M1a模型（近中性）

        icode = 0                   * 0=标准遗传密码
    fix_kappa = 0                   * 0=估计κ（转换/颠换比）
        kappa = 2                   * κ的初始值
    fix_omega = 0                   * 0=估计ω
        omega = 1                   * ω的初始值

    fix_alpha = 1                   * 1=固定alpha
        alpha = 0                   * 0=不使用gamma分布
       Malpha = 0
        ncatG = 4                   * gamma分类数

        getSE = 0                   * 0=不计算标准误
 RateAncestor = 0                   * 0=不做祖先重建
   Small_Diff = .5e-6
    cleandata = 1                   * 1=删除含gap的位点
EOF

codeml codeml_M1a.ctl  # 运行M1a模型

# === M2a 模型（备择假设：存在正向选择） ===
# 只需要修改 NSsites 和 outfile
sed 's/NSsites = 1/NSsites = 2/' codeml_M1a.ctl > codeml_M2a.ctl  # M2a模型
sed -i 's/M1a_results/M2a_results/' codeml_M2a.ctl  # 修改输出文件名
codeml codeml_M2a.ctl  # 运行M2a模型

# === M7 vs M8 比较（更敏感的检测） ===
sed 's/NSsites = 1/NSsites = 7/' codeml_M1a.ctl > codeml_M7.ctl
sed -i 's/M1a_results/M7_results/' codeml_M7.ctl
codeml codeml_M7.ctl

sed 's/NSsites = 1/NSsites = 8/' codeml_M1a.ctl > codeml_M8.ctl
sed -i 's/M1a_results/M8_results/' codeml_M8.ctl
codeml codeml_M8.ctl

第 4 步：似然比检验（LRT）¶

# 提取两个模型的对数似然值 (lnL)
grep "lnL" M1a_results.txt  # 例如：lnL = -3456.789
grep "lnL" M2a_results.txt  # 例如：lnL = -3450.123

# 计算 LRT 统计量
# LRT = 2 × (lnL_M2a - lnL_M1a)
# 自由度 df = M2a参数数 - M1a参数数 = 2（M2a比M1a多2个参数：p2和ω2）

# 用R计算p值
cat << 'EOF' > lrt_test.R
lnL_M1a <- -3456.789          # M1a的对数似然值（替换为实际值）
lnL_M2a <- -3450.123          # M2a的对数似然值（替换为实际值）

LRT <- 2 * (lnL_M2a - lnL_M1a)  # 计算LRT统计量
df <- 2                          # 自由度（M1a vs M2a 为2）
p_value <- pchisq(LRT, df, lower.tail=FALSE)  # 卡方检验p值

cat("LRT statistic:", LRT, "\n")  # 打印LRT值
cat("df:", df, "\n")               # 打印自由度
cat("p-value:", p_value, "\n")     # 打印p值

# p < 0.05 → M2a 显著优于 M1a → 存在正向选择
# p >= 0.05 → 不能拒绝M1a → 没有检测到正向选择
EOF
Rscript lrt_test.R

第 5 步：BEB 识别正选择位点¶

# 如果 LRT 显著（p < 0.05），查看 M2a 或 M8 结果中的 BEB 分析
grep -A 50 "Bayes Empirical Bayes" M2a_results.txt

# 输出格式：
# 位点编号  氨基酸  后验概率  ω值±SE
#   45   T   0.987*  3.456 +- 1.234
#  123   K   0.956*  2.789 +- 0.987
# * 表示后验概率 > 0.95
# ** 表示后验概率 > 0.99

# 后验概率 > 0.95 的位点被认为是正选择位点
# 这些位点的蛋白位置可以映射到3D结构上进行功能分析

第 6 步：分支-位点模型（检测特定谱系的正向选择）¶

# 分支-位点模型（Branch-Site Model）
# 检测某个特定分支上是否有位点经历了正向选择

# 备择假设（Model A）：前景分支有正选择位点
cat << 'EOF' > codeml_bsA.ctl
      seqfile = aligned_codon.phy
     treefile = tree_marked.nwk    * 树文件中前景分支用#1标记
      outfile = bsA_results.txt

        model = 2                   * 2=多分支模型
      NSsites = 2                   * 2=分支-位点模型
    fix_omega = 0                   * 0=估计ω
        omega = 1.5                 * ω的初始值

        noisy = 9
      verbose = 1
      runmode = 0
      seqtype = 1
    CodonFreq = 2
        clock = 0
        icode = 0
    fix_kappa = 0
        kappa = 2
    fix_alpha = 1
        alpha = 0
       Malpha = 0
        ncatG = 4
        getSE = 0
 RateAncestor = 0
   Small_Diff = .5e-6
    cleandata = 1
EOF
codeml codeml_bsA.ctl  # 运行备择假设

# 零假设（Model A null）：前景分支的ω固定为1
cat << 'EOF' > codeml_bsA_null.ctl
      seqfile = aligned_codon.phy
     treefile = tree_marked.nwk
      outfile = bsA_null_results.txt
        model = 2
      NSsites = 2
    fix_omega = 1                   * 1=固定ω
        omega = 1                   * ω固定为1（零假设）
        noisy = 9
      verbose = 1
      runmode = 0
      seqtype = 1
    CodonFreq = 2
        clock = 0
        icode = 0
    fix_kappa = 0
        kappa = 2
    fix_alpha = 1
        alpha = 0
       Malpha = 0
        ncatG = 4
        getSE = 0
 RateAncestor = 0
   Small_Diff = .5e-6
    cleandata = 1
EOF
codeml codeml_bsA_null.ctl  # 运行零假设

# LRT 比较（自由度 = 1）
grep "lnL" bsA_results.txt
grep "lnL" bsA_null_results.txt

三、常见报错与解决¶

报错信息	原因	解决方法
`Error: seq length not multiple of 3`	序列长度不是3的倍数	检查密码子比对，确保无移码
`stop codon found`	序列中有终止密码子	删除终止密码子或检查比对
`too many parameters`	物种数太少但模型太复杂	增加物种或用更简单的模型
收敛到局部最优	ω 初始值不好	多次运行，用不同初始值(ω=0.5,1,2)
`lnL = nan`	数据质量差或模型不适合	检查比对质量，尝试 `cleandata = 1`
BEB 没有显著位点	LRT 显著但无单个位点显著	正常，可能是多个位点共同贡献
`tree topology different`	树文件格式不正确	检查物种数和 Newick 格式

四、面试高频问题¶

Q1：什么情况下 ω 会大于 1？

免疫基因（如 MHC）、病毒受体基因、生殖隔离基因等经常检测到 ω>1。这些基因面临"军备竞赛"式的选择压力（如宿主与病原体共进化），需要不断改变蛋白功能。

Q2：位点模型和分支-位点模型的区别？

位点模型假设正向选择影响所有分支上的某些位点，适合检测长期的适应性进化。分支-位点模型只检测特定分支（前景分支）上的选择，适合检测特定谱系的适应事件（如人类祖先特异的基因加速进化）。

Q3：为什么不能直接用整个基因的 ω 来判断正向选择？

因为大多数蛋白质中只有少数关键位点受正向选择（ω>1），而大多数位点受纯化选择（ω<<1）。平均下来 ω 几乎总是 <1，即使存在正向选择也检测不到。所以需要位点模型来区分不同位点的 ω。

五、速查表¶

# === PAML codeml 速查 ===

# 模型设置关键参数：
# model = 0, NSsites = 0  → M0（一个ω）
# model = 0, NSsites = 1  → M1a（近中性：ω0<1, ω1=1）
# model = 0, NSsites = 2  → M2a（正选择：加ω2>1）
# model = 0, NSsites = 7  → M7（beta分布）
# model = 0, NSsites = 8  → M8（beta + ω>1）
# model = 2, NSsites = 2  → 分支-位点模型

# 标准比较对：
# M1a vs M2a (df=2)  — 检测位点水平正向选择
# M7 vs M8 (df=2)    — 更灵敏的位点检测
# Model A vs null (df=1) — 分支特异的正向选择

# LRT 计算：
# LRT = 2 × (lnL_alternative - lnL_null)
# 查卡方分布表得p值

# BEB 结果解读：
# 后验概率 > 0.95 (*) → 正选择位点
# 后验概率 > 0.99 (**) → 强正选择位点

参考资料：Beginner's Guide to PAML (2023, MBE)、PAML Tutorial on GitHub (abacus-gene/paml-tutorial)、PAML Demo 2024