621_分子进化选择压力分析¶

# === PAML安装 ===
conda install -c bioconda paml  # 安装PAML
codeml  # 验证（会报错缺少控制文件，说明安装成功）

# === HyPhy安装 ===
conda install -c bioconda hyphy  # 安装HyPhy
hyphy --version  # 查看版本

# === ETE Toolkit（PAML的Python封装，更易用）===
pip install ete3  # 安装ETE3

# === 第1步：准备密码子比对 ===
# 需要：蛋白序列比对 + CDS核苷酸序列 → 密码子比对

# 先用MAFFT比对蛋白序列
mafft --auto proteins.fasta > proteins_aligned.fasta

# 用PAL2NAL将蛋白比对映射回密码子比对
pal2nal.pl \
    proteins_aligned.fasta \  # 蛋白比对
    cds.fasta \               # CDS核苷酸序列
    -output paml \            # 输出PAML格式
    > codon_alignment.phy     # 输出文件

# === 第2步：准备系统发育树 ===
# 用IQ-TREE建树（用密码子模型）
iqtree2 -s codon_alignment.phy -m MFP -bb 1000 -nt AUTO

# 将树转为无分支长度的NEWICK格式（PAML需要）
# 树格式示例：((human,chimp),mouse,rat);

# === 准备CODEML控制文件 ===
cat > site_model.ctl << 'EOF'
      seqfile = codon_alignment.phy   * 密码子比对文件
     treefile = tree.nwk              * 系统发育树文件
      outfile = site_model_results    * 输出文件

        noisy = 3       * 输出详细程度
      verbose = 1       * 详细输出
      runmode = 0       * 用户提供的树

      seqtype = 1       * 密码子序列
    CodonFreq = 2       * F3x4密码子频率模型
        clock = 0       * 无分子钟
       aaDist = 0       * 默认氨基酸距离

        model = 0       * 一个ω比值给所有分支
      NSsites = 0 1 2 7 8   * 测试的位点模型
                            * 0=M0(单ω), 1=M1a(近中性)
                            * 2=M2a(正选择), 7=M7(beta分布)
                            * 8=M8(beta+正选择)

        icode = 0       * 标准遗传密码
    fix_kappa = 0       * 估计kappa
        kappa = 2       * kappa初始值
    fix_omega = 0       * 估计omega
        omega = 0.5     * omega初始值

        getSE = 0       * 不计算标准误
 RateAncestor = 0       * 不重建祖先序列
   Small_Diff = .5e-6   * 收敛精度
    cleandata = 1       * 删除有gap的位置
  fix_blength = 0       * 估计分支长度
EOF

# === 运行CODEML ===
codeml site_model.ctl  # 运行（可能需要几分钟到几小时）

# 查看结果文件
grep "lnL" site_model_results  # 提取各模型的似然值

# === 似然比检验（LRT）===
# 比较M1a vs M2a（检测正选择）
# 比较M7 vs M8（检测正选择，更严格）

# === R语言做LRT ===
# M1a vs M2a
lnL_M1a <- -5123.45  # 从结果文件中获取
lnL_M2a <- -5118.90

LRT <- 2 * (lnL_M2a - lnL_M1a)  # 似然比统计量
df <- 2  # 自由度差（M2a比M1a多2个参数）
p_value <- 1 - pchisq(LRT, df)  # p值

cat("M1a vs M2a:\n")
cat("  LRT =", LRT, "\n")
cat("  p-value =", p_value, "\n")
if (p_value < 0.05) {
    cat("  结论: 存在正选择位点！\n")
} else {
    cat("  结论: 无显著正选择证据\n")
}

# 如果M2a显著优于M1a，查看BEB（Bayes Empirical Bayes）结果
# BEB后验概率>0.95的位点被认为是正选择位点

# 在树文件中标记前景分支（在分支上加 #1）
# 例如：((human #1,chimp),mouse,rat);

cat > branch_site.ctl << 'EOF'
      seqfile = codon_alignment.phy
     treefile = tree_foreground.nwk   * 标记了前景分支的树
      outfile = branch_site_results

      seqtype = 1
    CodonFreq = 2
        model = 2       * 分支-位点模型
      NSsites = 2       * 模型A
    fix_omega = 0       * 备择假设：估计omega
        omega = 1.5     * omega初始值
EOF

codeml branch_site.ctl  # 运行备择模型（Model A）

# 然后修改控制文件运行零假设
# fix_omega = 1, omega = 1（固定omega=1）
# 比较两个模型的LRT

# BUSTED: 检测比对中是否存在正选择的证据
hyphy busted \
    --alignment codon_alignment.fasta \  # 密码子比对
    --tree tree.nwk \                     # 系统发育树
    --output busted_results.json          # JSON输出

# MEME: 混合效应模型，允许ω在分支间和位点间都变化
hyphy meme \
    --alignment codon_alignment.fasta \
    --tree tree.nwk \
    --output meme_results.json

# MEME比PAML的位点模型更灵活
# 不需要预先指定前景分支

# FUBAR: 比FEL更快更敏感
hyphy fubar \
    --alignment codon_alignment.fasta \
    --tree tree.nwk \
    --output fubar_results.json

# 输出：每个位点的dN和dS后验均值
# 后验概率>0.9的位点被认为受选择

# RELAX: 检测特定分支是否经历了选择放松或增强
hyphy relax \
    --alignment codon_alignment.fasta \
    --tree tree_labeled.nwk \
    --test Foreground \      # 检测前景分支
    --output relax_results.json

# K>1: 选择增强（更强的纯化或正选择）
# K<1: 选择松弛（约束减弱）