6、从四重奏到系统发生树：生物进化树构建方法解析

从四重奏到系统发生树：生物进化树构建方法解析

1. 引言

在生物学领域，依据生物数据构建系统发生树（或进化树）是一个经典问题，至今仍是重大挑战。多数考虑误差的现实问题表述会引发复杂的计算难题。本文聚焦于基于四重奏（quartet）的树重建方法，先简述一般的树重建问题，探讨使用四重奏重建的动机，接着分析该重建任务的计算复杂度，最后概述从四重奏构建树的一些算法和启发式方法。

系统发生重建的目标是为一组相关生物物种（分类单元）构建一棵树，以最佳呈现其随时间的进化历程。树的叶子标注为现存分类单元，内部节点代表假设的灭绝分类单元。由于分类单元分歧事件被认为较为罕见，理想的树是二叉树，内部节点度数为 3（在数据模糊时可能需采用信息较少的多叉树）。早期主要利用形态特征研究进化，如今分子数据成为系统发生分析的主要依据，但古生物学、解剖学和形态学等其他信息源也在使用。为简化讨论，本文主要关注分子序列数据。

构建树的第一步是从更新的数据库中收集所有研究分类单元的 DNA（通常是基因）、RNA 或氨基酸序列（为简便，本文主要讨论蛋白质）。来自不同分类单元的同源序列（通过相似性或低编辑距离检测）会被分组。不同分类单元的同源序列通常具有相同功能，被认为是共同祖先序列的后代，其相似程度可指示两个分类单元分歧的时间。由于突变过程本质上是概率性的且局部发生，我们预期分歧时间越长，积累的突变越多。然而，不同蛋白质的进化速率可能差异很大，单个蛋白质单独考虑时可能对进化历史给出相互矛盾的指示。为克服这种“噪声效应”，建议使用通过连接多个同源序列获得的更长序列。

一般来说，系统发生重建方法分为基于特征和基于距离的方法。基于特征的方法直接处理代表各种生物特征的特征数据，试图生成一棵使树边总变化数最小的树；基于距离的方法首先计算分类单元对之间的“进化距离”，然后寻找一棵带权边的树，其成对树距离近似进化距离。

2. 基于特征和基于距离的重建方法

2.1 基于特征的方法

基于特征的方法考虑输入分类单元的定性特征。每个特征根据分类单元的值对输入集进行划分，每个由此定义的等价类称为特征状态。例如，DNA 序列由 4 种核苷酸字符 A、C、T、G 组成，RNA 序列由 A、C、U、G 组成，蛋白质序列由 20 种氨基酸字符 A、C、D、E、F、G、H、I、K、L、M、N、P、Q、R、S、T、V、W 和 Y 组成。来自不同物种的序列通过多序列比对过程进行比对，必要时在序列中插入空位以最大化序列之间的相似度。因此，比对后的氨基酸序列的每个位置（称为位点）是一个具有 21 个状态（20 种氨基酸和空位符号）的特征。

以下是胰岛素样生长因子 II 对应四种物种的多序列比对结果：
| 物种 | 序列 |
| — | — |
| (1) Human (PRI) | MGIPMGKSMLVLLTFLAFASCCIAAYRPSETLCGGELVDTLQFVCGDRGFYFSRPASRVS |
| (2) Mouse (MUR) | MGIPVGKSMLVLLISLAFALCCIAAYGPGETLCGGELVDTLQFVCSDRGFYFSRPSSRAN |
| (3) Sheep (RUM) | MGITAGKSMLALLAFLAFASCCYAAYRPSETLCGGELVDTLQFVCGDRGFYFSRPSSRIN |
| (4) Chicken (OUT) | MC - AARQILLLLLAFLAYALDSAAAYGTAETLCGGELVDTLQFVCGDRGFYFSRPVGRNN |
| | * . . ** . . *** * * * |
| (1) | RRS - RGIVEECCFRSCDLALLETYCATPAKSERDVSTP-----PTVLPDNFPRYPVGKFF |
| (2) | RRS - RGIVEECCFRSCDLALLETYCATPAKSERDVSTS-----QAVLPDDFPRYPVGKFF |
| (3) | RRS - RGIVEECCFRSCDLALLETYCAAPAKSERDVSAS-----TTVLPDDFTAYPVGKFF |
| (4) | RRINRGIVEECCFRSCDLALLETYCAKSVKSERDLSATSLAGLPALNKESFQKPSHAKYS |
| | ** * * * *. . .. . * . |
| (1) | QYDTWK - QSTQRLRRGLPALLRARRGHVLAKELEAFREA - KRHRPLIALPTQDPA - HGGA |
| (2) | QYDTWR - QSAGRLRRGLPALLRARRGRMLAKELKEFREA - KRHRPLIVLPPKDPA - HGGA |
| (3) | QSDTWK - QSTQRLRRGLPAFLRARRGRTLAKELEALREA - KSHRPLIALPTQDPATHGGA |
| (4) | KYNVWQKKSSQRLQREVPGILRARRYRWQAEGLQAAEEARAMHRPLISLPSQRPP - APRA |
| | . . .*. . . * . * * ** *** ** . * * |
| (1) | PPEMASNRK |
| (2) | SSEMSSNHQ |
| (3) | SSEASSD– |
| (4) | SPEATGPQE |
| | * . |

多序列比对的输出是一个 (|S|\times|C|) 矩阵，其中 (S) 是分类单元集，(C) 是特征集（在上述示例中，(|S| = 4)，(|C| = 189)）。每个条目表示特定分类单元对于给定特征所呈现的状态，该矩阵作为基于特征方法的输入。给定每个分类单元的比对序列（插入空位后），一种自然的方法是构建一棵树，其内部节点和外部节点都用序列标注。叶子的标签作为输入，树的拓扑结构和内部标签由算法计算得出。目标是最小化相邻树边为解释叶子数据所需的“突变”数量，这导致一个最小化问题，目标函数是树边相邻序列之间汉明距离的总和。更精细的方法会为不同变化赋予不同的“价格”，该“价格”衡量局部变化的对数似然性，因为并非所有点突变的可能性都相同。对齐位置的成本总和决定了任意一对序列的成本函数，最大简约性的全局优化标准是最小化相邻序列之间的成本总和。一般的最大简约性问题是 NP 难的，但对于给定拓扑结构找到内部标签可以高效完成，这促使算法扫描许多不同的拓扑结构并输出最简约的树。

另一种方法试图构建一棵树，使得对于每个特征，对应任何特征状态的节点集形成一个连通子图。最大化满足此条件的特征数量的问题称为最大兼容性问题，该问题也是 NP 难的。同时满足最大简约性和最大兼容性的系统发生树称为完美系统发生树，确定数据是否支持完美系统发生树已被证明是 NP 难的（通过与三角化彩色图问题的等价性以及从介于性问题的归约）。尽管在一般情况下很难，但在某些输入参数（状态数量或特征数量）固定的情况下，存在多项式时间的完美系统发生树算法。

2.2 基于距离的方法

基于距离的系统发生方法通常首先对输入序列进行所有成对比对，计算每对序列之间的编辑距离，得到一个对称的零对角 (|S|\times|S|) 距离矩阵。目标是生成一棵其诱导度量能最好表示输入数据的系统发生树。如果输入距离矩阵 (M) 可以由一棵树及其诱导路径长度实现，则称 (M) 是可加的。所有叶子到根的距离相同的特殊情况称为超度量，超度量树对应于不同物种的替换事件以相同速率发生的生物学理论（即“分子钟”假设，如今该假设普遍受到质疑）。

给定一个可加度量，构建树相对容易。但现实输入存在误差，部分误差源于假设的进化模型本身。因此，我们寻求其诱导度量在某种标准下近似“最佳可能”树度量的系统发生树。同样，该领域的大多数问题都是 NP 难的。

一种流行的启发式方法是凝聚聚类，它迭代工作，每次迭代将两个节点合并为一个父节点，例如流行的邻接法（Neighbor Joining method）、Fitch - Margoliash 方法和 BIONJ 方法。

事实证明，在 (L_1) 准则下近似超度量可以在多项式时间内完成，因此存在使用近似超度量来生成可加度量近似的算法。

3. 基于四重奏重建的动机

以哺乳动物为例，人类是研究最广泛的物种，其次是大鼠、小鼠等模式物种以及牛、羊等具有经济/农业重要性的物种。而犰狳等物种在可预见的未来不太可能获得政府或行业资助进行基因组项目，因此不同分类单元之间可用序列数据存在显著差异的情况将持续多年。

以下是 6 个分类单元（5 个哺乳动物分类单元和 1 个非哺乳动物“外类群”）与 15 种蛋白质序列的部分可用数据矩阵示例：
| 蛋白质 | PRI | MUR | RUM | LAG | INS | OUT |
| — | — | — | — | — | — | — |
| insulin - like growth factor II | + | + | + | - | - | + |
| insulin | + | + | + | + | - | + |
| angiotensin - converting enzyme | + | + | - | + | - | + |
| HCG beta | + | + | + | + | - | + |
| acetylcholine receptor | + | + | + | - | + | + |
| GABA - A receptor alpha - 1 subunit | + | + | + | - | - | + |
| guanine nucleotide - binding protein | + | + | + | - | - | + |
| growth hormone | + | + | + | + | - | + |
| prolactin | + | - | + | + | - | + |
| endothelin 2 | + | + | + | + | - | - |
| TGF - beta2 | + | + | - | - | - | + |
| histamine H2 receptor | + | + | - | - | - | - |
| serotonin receptor | + | + | - | + | - | - |
| PAI - 1 | + | + | + | - | - | - |
| alpha globin | + | + | + | + | + | + |

在这个小表格中，所有 15 种蛋白质序列都已知于灵长目（PRI），10 种已知于兔形目（LAG），而食虫目（INS）仅已知 2 种。实际情况中，“ - ”条目比表格所示更多。从 HOVERGEN 数据库提取的关于 23 个哺乳动物分类单元和 1 个外类群分类单元的最新数据显示，测序蛋白质的数量从灵长目的 621 种到金毛鼹科的 3 种不等，每个分类单元已知序列的蛋白质平均数量为 135 种，中位数为 19 种。特定蛋白质测序的分类单元数量也有类似变化，大多数蛋白质（共 363 种）仅在 4 个分类单元中测序，只有 3 种在多达 20 个分类单元（共 24 个分类单元）中测序。

由于数据差异问题，分子系统发生学家经常需要在分类单元数量和研究中使用的分子数据量之间进行权衡。如果只考虑所有研究分类单元都有的共同序列，可能会忽略绝大多数数据，这样构建的树会强烈偏向少数过度代表的蛋白质的进化。另一方面，如果坚持选择有大量蛋白质序列已知的分类单元，最终得到的分类单元数量会很少（这称为“分类单元采样”）。基于特征和基于距离的方法都以所有研究分类单元都已知的序列列表为起点，因此这两种方法都会受到数据差异问题的影响。这意味着为了利用所有可用数据，需要一种不同的方法。

一种尝试绕过数据差异问题的方法是为每个分类单元创建一个由单个序列连接而成的长序列，在序列信息未知的地方连接适当数量的“缺失数据”符号，这些长序列可以作为基于特征或基于距离方法的基础。但问题是，这些长连接序列大多主要由“缺失数据”符号组成，导致生成的树质量往往较差。

一种试图利用所有可用数据同时避免分类单元采样问题的方法是文献[24]提出并在[26, 26]中使用的四分类单元方法。关键思想是一次考虑一小部分分类单元（例如大小为 (l)），对于每个这样的子集，选取所有该子集中所有分类单元都已知的蛋白质。使用这些共同序列的集合，应用最大简约性或基于距离的方法来推断每个子集的系统发生树。由于小的系统发生树比大的更容易推断，这一步在计算上是可行的。这种方法的优点是每个至少在 (l) 个分类单元中共同的蛋白质都会影响树的构建，不必在所有 (n) 个分类单元中都共同，这意味着更多的序列将被利用，更少的序列会“被浪费”。

我们能期望得到什么样的信息呢？不同子集通常会共享不同的蛋白质，不同蛋白质的进化速率可能差异很大，因此我们不能以统一的方式跨不同子集使用任何度量（或距离）信息，只能得到拓扑信息——每个 (l) 个叶子的无根树。显然，较小的 (l) 能更好地利用序列信息，但 (l) 能小到什么程度呢？有 (l = 2) 个叶子的树只是连接两个叶子的一条边，对于 (l = 3)，只有一种无根树（星形树），因此 (l = 2) 和 (l = 3) 都不能提供有信息的数据。能产生依赖于输入的有信息的树拓扑的 (l) 的最小值是 (l = 4)。对于每个 4 个分类单元的子集，有三种可能的无根二叉系统发生树，这种具有相关二叉拓扑的四个分类单元的集合称为四重奏。

graph LR
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px;
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
    A([开始]):::startend --> B(收集分类单元序列):::process
    B --> C(多序列比对):::process
    C --> D{选择重建方法}:::process
    D -->|基于特征| E(构建特征矩阵):::process
    D -->|基于距离| F(计算距离矩阵):::process
    E --> G(最大简约性或兼容性计算):::process
    F --> H(凝聚聚类等方法):::process
    G --> I(输出树):::process
    H --> I
    I --> J([结束]):::startend

以上流程图展示了从收集分类单元序列到最终输出系统发生树的一般过程，包括多序列比对、选择重建方法（基于特征或基于距离）以及相应的计算步骤。

4. 四重奏重建问题及计算复杂度

四重奏重建问题可以表述为：给定一组四重奏（即四个分类单元的子集及其对应的二叉拓扑结构），目标是构建一棵包含所有分类单元的系统发生树，使得该树与给定的四重奏集合在某种程度上最匹配。

已知的计算复杂度结果表明，四重奏重建问题是 NP 难的。这意味着在一般情况下，找到最优解是非常困难的，可能需要指数级的时间和资源。不过，这并不意味着我们无法解决这个问题，只是需要采用一些启发式算法或近似算法来在合理的时间内得到一个较好的解。

5. 算法和启发式方法

为了解决四重奏重建问题，研究人员提出了多种算法和启发式方法。以下是一些常见的方法：

5.1 自底向上的方法

• 步骤 ：
  1. 从单个分类单元开始，逐步合并四重奏。
  2. 在每一步中，选择合适的四重奏进行合并，使得合并后的树与已知四重奏的一致性尽可能高。
  3. 重复步骤 2，直到所有分类单元都被合并到一棵树上。

这种方法的优点是计算相对简单，能够在较短的时间内得到一个可行的解。但缺点是可能会陷入局部最优，无法保证得到全局最优解。

5.2 基于优化的方法

• 步骤 ：
  1. 定义一个目标函数，用于衡量树与四重奏集合的一致性。
  2. 使用优化算法（如模拟退火、遗传算法等）来搜索使目标函数值最大的树。
  3. 不断迭代优化过程，直到满足停止条件（如达到最大迭代次数或目标函数值不再改善）。

基于优化的方法能够更全面地搜索解空间，有更大的机会找到全局最优解，但计算复杂度较高，需要较长的计算时间。

5.3 启发式选择方法

• 步骤 ：
  1. 根据四重奏的一些特征（如支持度、多样性等）对四重奏进行排序。
  2. 优先选择排序靠前的四重奏进行合并，逐步构建树。
  3. 在合并过程中，根据已构建的树和剩余的四重奏动态调整选择策略。

这种方法结合了自底向上方法的计算效率和基于优化方法的搜索能力，能够在一定程度上平衡计算复杂度和结果质量。

6. 计算结果和开放问题

在实际应用中，上述算法和启发式方法都取得了一定的成果。通过对模拟数据和真实生物数据的实验，我们可以评估不同方法的性能，比较它们在准确性、计算时间和内存使用等方面的差异。

然而，四重奏重建问题仍然存在一些开放问题，例如：
- 最优解的搜索 ：如何设计更高效的算法来找到四重奏重建问题的全局最优解，尤其是在大规模数据集上。
- 数据噪声的处理 ：实际生物数据中往往存在噪声和误差，如何在重建过程中有效地处理这些噪声，提高重建结果的准确性。
- 方法的可扩展性 ：随着分类单元数量和四重奏数量的增加，现有的算法和启发式方法可能会面临计算性能和内存使用的挑战，如何设计具有更好可扩展性的方法。

graph LR
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px;
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
    A([开始四重奏重建]):::startend --> B{选择方法}:::process
    B -->|自底向上| C(逐步合并四重奏):::process
    B -->|基于优化| D(定义目标函数):::process
    B -->|启发式选择| E(四重奏排序):::process
    C --> F(构建树):::process
    D --> G(优化搜索):::process
    E --> H(动态选择合并):::process
    F --> I(评估结果):::process
    G --> I
    H --> I
    I --> J{结果满意?}:::process
    J -->|是| K([结束]):::startend
    J -->|否| B

以上流程图展示了四重奏重建的一般过程，包括选择不同的方法、构建树、评估结果以及根据结果决定是否继续迭代。

总之，基于四重奏的系统发生树重建是一种有潜力的方法，能够在一定程度上解决数据差异问题，更好地利用生物序列信息。但要实现准确、高效的重建，还需要进一步深入研究和改进算法。通过不断探索和创新，我们有望在生物进化研究中取得更深入的理解和更准确的结果。