系统发育与进化分析

引言

系统发育分析（phylogenetics）研究「物种或基因之间的进化关系」，用树状结构表示共同祖先与分化历史。它是进化生物学、微生物分类、病毒溯源、物种鉴定的核心方法。2020 年新冠疫情期间的病毒溯源分析，就是系统发育方法的经典应用——通过比较病毒基因组的差异，推断传播链条与共同祖先。

与变异检测、表达分析不同，系统发育的难点在于「推断历史」而非「测量当下」。我们只有现存的序列，要反推几百万年前的进化事件。这需要统计模型（核苷酸替换模型、分子钟模型）来「倒推时间」，而模型的假设是否正确，直接决定结论的可信度。

工程上的挑战有三个：一是多序列比对（MSA）的质量，比对错了，建树必然错，且 MSA 本身是 NP-hard 问题，只能近似求解；二是计算复杂度，最大似然建树是组合优化问题，序列多、物种多时计算量爆炸；三是模型选择，核苷酸替换模型有几十种，选错了会得到错误的树。这些都是需要理解原理才能做好的。

本文按「树的基本概念 → 多序列比对 → 距离法与模型法 → 最大似然与贝叶斯 → 建树工具 → 树评估 → 分子钟 → 选择压力 → 宏基因组」的顺序展开。工具基于 IQ-TREE 2.2、MAFFT 7.5、PAML 4.10。读完你应该能独立完成一次系统发育分析，并理解「树的置信度」意味着什么。

目录

1. 系统发育树的基本概念
2. 多序列比对（MSA）
3. 距离法与核苷酸替换模型
4. 最大似然与贝叶斯推断
5. 建树工具：IQ-TREE、RAxML、MrBayes
6. 树的评估：bootstrap 与后验概率
7. 分子钟与时间树
8. 选择压力分析：dN/dS
9. 宏基因组与微生物系统发育

1. 系统发育树的基本概念

系统发育树（phylogenetic tree）由节点（node）与分支（branch）构成：

        ┌── A
     ┌──┤
     │  └── B
  ───┤
     │  ┌── C
     └──┤
        └── D

  · 叶节点（leaf）：现存物种/序列（A、B、C、D）
  · 内部节点（internal node）：推断的共同祖先
  · 分支长度（branch length）：进化距离（替换数/位点）
  · 根（root）：共同祖先的位置

几个核心概念：

• 拓扑结构（topology）：树的分支模式（谁和谁更近），这是建树的主要目标；
• 分支长度：通常表示「每位点的替换数」，反映进化量；
• 单系群（monophyletic/clade）：一个祖先及其所有后代，是分类学的目标；
• 根（root）：多数建树方法产生「无根树」（只知关系不知方向），定根需要外部信息（外群 outgroup）。

定根（rooting） 是常见难点。无根树只告诉你「A 和 B 更近」，但不告诉你「谁更古老」。定根方法：

1. 外群法：加入一个已知在目标群之外的序列，根位于外群与目标群之间；
2. 分子钟法：假设进化速率恒定，用 LSD2、TreeTime 定根；
3. 中点定根（midpoint rooting）：把根放在树的最长路径中点，简单但假设强。

理解「无根 vs 有根」很重要：很多软件（如 IQ-TREE 默认）输出无根树，如果直接拿来解读「谁先分化」就错了。

2. 多序列比对（MSA）

建树的前提是多序列比对：把同源序列对齐，使「同一列」代表「同一个进化位置」。MSA 的质量直接决定树的正确性——垃圾进，垃圾出。

MSA 是 NP-hard 问题，实用工具都用启发式。两大流派：

渐进式比对（progressive）：先比对最相似的序列，再逐步加入其他序列。快，但对「比对顺序」敏感。代表：Clustal、MAFFT（默认模式）。

迭代精炼（iterative refinement）：在初始比对基础上反复调整，优化打分。慢但更准。代表：MAFFT 的 L-INS-i 模式、PRANK。

# MAFFT 常用模式
mafft --auto input.fasta > aligned.fasta           # 自动选择策略（推荐起点）
mafft --maxiterate 1000 --localpair input.fasta > out.fasta  # L-INS-i，高精度
mafft --thread 8 --auto input.fasta > aligned.fasta

# 其他工具
muscle -in input.fasta -out aligned.fasta          # MUSCLE v5，快
clustalo -i input.fasta -o aligned.fasta --outfmt=fasta  # Clustal Omega

比对后的修剪（trimming） 是常被忽略但重要的一步：比对两端或中间常有「gap 密集区」，这些区域对齐不可靠，会引入噪声。用 trimAl 或 Gblocks 修剪：

# trimAl 自动修剪
trimal -in aligned.fasta -out trimmed.fasta -automated1

# 或按 gap 比例修剪
trimal -in aligned.fasta -out trimmed.fasta -gt 0.5   # 保留 gap < 50% 的列

一个实用原则：用多个 MSA 工具/参数分别建树，看结论是否一致。如果不同比对给出不同的树拓扑，说明数据本身的分辨率不足，结论需谨慎。

3. 距离法与核苷酸替换模型

建树方法分三大类：距离法、最大简约法、最大似然/贝叶斯。先从最简单的距离法说起。

距离法：先计算每对序列的进化距离（矩阵），再用聚类算法建树。最常用的是邻接法（Neighbor-Joining, NJ）：

# 用 dist 计算距离矩阵，NJ 建树
# （实际中多用 MEGA 或 R 的 ape 包）

距离法的关键在「如何把观测差异转换成进化距离」。观测到的差异（p-distance）会低估真实进化距离，因为同一位置可能发生多次替换（回复突变、平行突变）。核苷酸替换模型校正这个问题：

为什么 +G 很重要？ 不同位点的进化速率不同——编码区的密码子第三位变化快，第一、二位慢；非编码区更快。如果不建模速率异质性，会低估长分支的距离（「长枝吸引」问题）。+G（Gamma 分布）是几乎所有分析的标配。

模型选择用 ModelFinder（IQ-TREE 内置）自动完成：

# IQ-TREE 自动选择最优模型
iqtree2 -s aligned.fasta -m MFP -B 1000 -T 8
#  -m MFP: ModelFinder Plus，自动选模型
#  -B 1000: 1000 次 ultrafast bootstrap

ModelFinder 会评估几十个模型，用 BIC 准则选择最优。输出的 .iqtree 文件里会报告选中的模型（如 GTR+F+G4）。不要手动猜模型——自动选择又快又准。

4. 最大似然与贝叶斯推断

现代建树的主力是最大似然（Maximum Likelihood, ML）与贝叶斯推断（Bayesian Inference, BI）。

最大似然（ML）：给定序列数据和进化模型，找到「使数据出现概率最大」的树。它搜索树空间，对每棵候选树计算似然值，取最大值。

ML 的核心：
  L(tree) = P(sequence data | tree, model)
  搜索 tree 空间，最大化 L

  搜索算法：NNI（最近邻交换）、SPR（子树修剪嫁接）
  从初始树出发，反复尝试局部修改，直到无法改进

ML 的挑战是「树空间巨大」——20 个物种就有约 10^21 种无根树。搜索只能找「局部最优」，可能陷入局部极值。缓解方法是「多次随机起点 + 更激进的搜索」。

贝叶斯推断（BI）：用 MCMC 采样树的后验分布，输出「树 + 分支的后验概率」。与 ML 的区别：

选择：常规分析用 ML（快、够用），需要时间估计或复杂模型用贝叶斯。ML 的 bootstrap 支持率与贝叶斯后验概率含义不同（bootstrap 是频率学派的重采样，后验是贝叶斯概率），不能直接比较数值。

5. 建树工具：IQ-TREE、RAxML、MrBayes

三大主流工具各有侧重：

IQ-TREE 2 是当前 ML 建树的首选：速度快、内置模型选择、支持 ultrafast bootstrap：

# 标准 ML 建树
iqtree2 -s aligned.fasta -m MFP -B 1000 -alrt 1000 -T AUTO --prefix out
#  -m MFP: 自动模型选择
#  -B 1000: ultrafast bootstrap（快，推荐）
#  -alrt 1000: SH-aLRT 检验（补充置信度）
#  -T AUTO: 自动线程数

# 分区模型（多基因/多密码子位置分别建模）
iqtree2 -s concat.fasta -p partition.txt -m MFP -B 1000

RAxML 是老牌 ML 工具，适合超大规模数据：

raxml-ng --all --msa aligned.fasta --model GTR+G --bs-trees 100 --threads 8

MrBayes 是贝叶斯建树的经典工具：

# MrBayes 需要 nexus 格式输入
# 运行后输出 .con.tre（共识树）与 .p（后验概率）
mb mrbayes_block.nex

工具选择的实用建议：

• 常规分析 → IQ-TREE（快、准、易用）；
• 超大规模（上千序列） → RAxML-NG 或 FastTree（近似法）；
• 需要后验概率/时间树 → MrBayes 或 BEAST；
• 快速探索 → FastTree（极快，但精度低，适合初步看结构）。

一个常见误区是「用 FastTree 做最终结果」。FastTree 是近似算法，速度快但精度有限，适合探索性分析；正式发表的树应该用 ML（IQ-TREE）或贝叶斯。

6. 树的评估：bootstrap 与后验概率

树是「推断」的结果，必须给出置信度。评估方法：

Bootstrap（自举）：从比对中「有放回地重采样」列，重建树，重复数百次，统计每个分支在多少比例的树中出现。这个比例就是 bootstrap 支持率：

分支支持率 95% → 该分支在 95% 的重采样树中出现
  ≥ 95%: 强支持
  70-95%: 中等支持
  < 70%: 弱支持（不可信）

Ultrafast Bootstrap（UFBoot）：IQ-TREE 的创新，用近似方法把 bootstrap 加速 100 倍以上。1000 次 UFBoot 只需几分钟（传统 bootstrap 要几小时）。当前标准做法。

SH-aLRT：另一种分支检验，常与 UFBoot 联合报告。一个分支若同时满足 UFBoot ≥ 95% 且 SH-aLRT ≥ 80%，可信度高。

后验概率（贝叶斯）：MCMC 采样中该分支出现的比例。≥ 0.95 通常认为可信。注意后验概率通常比 bootstrap 数值高（贝叶斯方法倾向于给出更乐观的置信度），不要直接对比。

# IQ-TREE 输出中的支持率
# 查看 .treefile，分支上的数字就是 UFBoot/SH-aLRT 支持率
# 用 FigTree 或 iTOL 可视化

一个重要的心态：低支持率不代表「树错了」，而代表「数据不足以分辨」。如果关键分支支持率低，可能需要更多数据（更多基因、更长序列）或更好的模型，而非换工具。

7. 分子钟与时间树

系统发育树只给拓扑，不给时间。要估计「分化发生在多少百万年前」，需要**分子钟（molecular clock）**假设：序列以大致恒定的速率进化，所以分支长度可以换算成时间。

严格分子钟：假设所有分支速率相同。简单，但常被违反（不同物种/基因的进化速率差异大）。

松弛分子钟（relaxed clock）：允许速率在不同分支上变化。更符合实际，是现代主流。

校准（calibration）：分子钟需要「时间锚点」，即已知的分化时间（来自化石记录、地理事件、已知的采样时间）。校准方法：

工具：BEAST 2 是时间树的标准工具（贝叶斯 + 松弛分子钟），TreeTime 更快但更简单。

# BEAST 需要 XML 配置（通常用 BEAUti 图形界面生成）
beast -threads 8 analysis.xml

# TreeTime：快速时间树（适合病毒等有采样日期的场景）
treetime --aln aligned.fasta --tree ml_tree.nwk --dates dates.csv

病毒溯源是时间树的典型应用：病毒的采样日期已知，序列差异随时间累积，用「采样时间校准」的分子钟可以估计「最近共同祖先（TMRCA）」的时间——这就是「新冠病毒某谱系起源于某年某月」这类结论的来源。

一个警示：分子钟的结论依赖校准的准确性。化石校准的年龄本身有误差（可能 ±几千万年），传导到时间树估计上就是巨大的不确定性。报告中必须给出置信区间，而非单点估计。

8. 选择压力分析：dN/dS

系统发育不只是「看关系」，还能分析「哪些基因受选择」。核心指标是 dN/dS（ω）：

dN = 非同义替换率（改变氨基酸，通常有害）
dS = 同义替换率（不改变氨基酸，近似中性）

ω = dN / dS
  ω < 1: 纯化选择（purifying，有害突变被淘汰）—— 多数基因
  ω = 1: 中性进化（neutral）
  ω > 1: 正选择（positive，有利突变被保留）—— 少数，重要

为什么 dS 可以当「中性参照」？ 因为同义突变不改变蛋白，大多不受选择，其替换率反映「背景突变率」。用 dN 与之比较，就能剥离背景突变率，看出选择压力。

工具是 PAML 的 codeml（需要系统树 + 密码子比对）：

# codeml 需要配置文件（.ctl），指定树、序列、模型
codeml codeml.ctl

# 常用模型：
#   M0: 单一 ω
#   M1a vs M2a: 检验正选择
#   M7 vs M8: 检验正选择（更严格）
#   branch-site: 检验特定分支的正选择

位点模型（site model） 分析「哪些位点受选择」，分支模型（branch model） 分析「哪些分支受选择」，分支-位点模型结合两者。一个经典应用是「检测免疫基因的正选择位点」——被病原体驱动的基因往往在特定位点有 ω > 1 的信号。

还有一个快速替代工具 HyPhy（尤其 MEME、FEL 方法），比 PAML 更快且能检测「瞬时正选择」。选择压力分析的常见陷阱是「序列太相似导致 dS 无法估计」（dS ≈ 0 时 ω 趋于无穷），此时需要足够的序列分歧。

9. 宏基因组与微生物系统发育

系统发育在微生物研究中有特殊地位，因为微生物形态简单、难以用传统分类，只能靠序列（尤其是 16S rRNA 基因）划分。

16S rRNA 扩增子分析：16S rRNA 基因有「保守区」（所有细菌共有，用于设计引物）和「可变区」（V1-V9，物种特异，用于分类）。分析流程：

# DADA2：扩增子序列变异（ASV）分析
# （R 包，典型流程）
# 1. 去引物、质控 → 2. 去噪（学习错误率）→ 3. 合并双端
# 4. 去嵌合体 → 5. 生成 ASV 表 → 6. 物种注释（SILVA 数据库）

ASV vs OTU：传统方法把相似度 > 97% 的序列聚成 OTU（操作分类单元），但 97% 是人为阈值。现代方法用 ASV（精确序列变异），保留单碱基分辨率，更可复现。

系统发育在微生物组的应用：

1. 物种注释：把 ASV 比对到参考数据库（SILVA、Greengenes）的进化位置；
2. 系统发育多样性：用树计算群落多样性（Faith’s PD），比单纯的物种计数更能反映进化信息；
3. 系统发育比较：看不同群落的进化组成差异（UniFrac 距离）。

# 构建 ASV 的系统发育树（用于 UniFrac）
mafft --auto asv_seqs.fasta > asv_aligned.fasta
iqtree2 -s asv_aligned.fasta -m MFP -B 1000 -T 8
# 用得到的树做 UniFrac 分析（QIIME2 或 phyloseq）

宏基因组的系统发育比 16S 更复杂：它测的是全部基因，需要「分箱（binning）」把序列归到物种，再建树。工具如 GTDB-Tk 用 120 个单拷贝标记基因建树，是微生物基因组分类的标准。

一个趋势是「系统发育感知的分析」：不再把物种当独立单位，而是用进化树的结构来做统计（系统发育回归、系统发育信号检测）。这是因为近缘物种的性状不独立（共享祖先），忽略系统发育会导致假阳性。

权衡取舍

常见坑清单

1. MSA 质量差就建树：比对错了树必错；用 trimAl 修剪并多工具交叉验证。
2. 不做模型选择：默认模型不当导致错误的树；用 ModelFinder 自动选。
3. 忽略速率异质性：不建模 +G 导致长枝吸引；几乎所有分析都加 +G。
4. 无根树当有根解读：直接说「谁先分化」是错的；需外群或分子钟定根。
5. bootstrap 太低仍下结论：支持率 < 70% 的分支不可信；数据不足需补充。
6. 混淆 bootstrap 与后验概率：两者数值不可直接比较；分开解读。
7. 分子钟无校准：时间树没有锚点，结论无意义；必须提供校准点。
8. dS ≈ 0 时算 ω：序列太相似导致 ω 无穷大；需要足够的序列分歧。
9. 用 FastTree 出正式结果：近似算法精度有限；正式分析用 ML/贝叶斯。
10. 忽略比对修剪：gap 密集区引入噪声；用 trimAl/Gblocks 修剪不可靠列。

小结

系统发育分析是「从现存数据反推历史」的推断艺术，它的每一个结论都带着模型假设的影子。MSA 决定输入质量，替换模型决定距离校正，建树方法决定搜索策略，bootstrap/后验概率决定置信度，分子钟校准决定时间估计的可靠性。理解这条链条上的每一环，才能判断一棵树是否可信。

工程上，系统发育最该建立的认知是「不确定性是常态」。树不是「唯一正确答案」，而是「在当前数据和模型下最可能的推断」，且这个推断带着置信区间。报告结果时，分支支持率、模型选择、校准来源都应透明呈现。另一个要点是「多方法交叉验证」：用不同比对、不同模型、不同软件各建一次树，看核心拓扑是否稳定——稳定的结论才值得相信。

下一步可以看 蛋白质结构预测与 AlphaFold 了解从序列到结构的推断，或 序列比对 复习比对的基础。如果你建的树「拓扑奇怪」或「支持率普遍很低」，先检查 MSA 质量与模型选择，再考虑数据是否足够。