山大生物信息学讲义第9章 系统发育分析

1、第9章系统发育分析系统发育学研究的是进化关系，系统发育分析就是要推断或者评估这些进化关系。通过系统发育分析所推断出来的进化关系一般用分枝图表（进化树）来描述，这个进化树就描述了同一谱系的进化关系，包括了分子进化（基因树）、物种进化以及分子进化和物种进化的综合。因为”clade”这个词（拥有共同祖先的同一谱系）在希腊文中的本意是分支，所以系统发育学有时被称为遗传分类学(cladistics)。在现代系统发育学研究中，研究的重点已经不再是生物的形态学特征或者其他特性，而是生物大分子尤其是序列。尽管本章的目的是想描述一个万能的系统发育分析方法，但是这个目标过于天真，是不可能实现的(Hills et

2、al.,1993)。虽然人们已经设计了大量的用于系统发育的算法、过程和计算机程序，这些方法的可靠性和实用性还是依赖于数据的结构和大小，因此关于这些方法的优点和缺点的激烈争论成了分类学和系统发育学期刊的热门课题（Avise, 1994:Saitou, 1996; Li, 1997; Swofford et al., 1996a）。Hillis等人（1993）曾经极其简单地介绍过系统发育学。比起许多其它学科而言，用计算的方法进行系统发育分析很容易得到错误的结果，而且出错的危险几乎是不可避免的；其它学科一般都会有实验基础，而系统发育分析不太可能会拥有实验基础，至多也就是一些模拟实验或者病毒实验（Hi

3、llis et al., 1994）；实际上，系统发育的发生过程都是已经完成的历史，只能去推断或者评估，而无法再现了。尽管现在的系统发育分析的程序存在着很多缺陷，许多生物学著作中仍然大量引用系统发育分析得到的结论，这些结论的源数据只是简单地应用了某一个系统发育的分析程序。只有在很偶然的情况下，在应用某些 HYPERLINK /soft/ 软件处理数据时受到的困难会少一点，但是在概念上理解这些程序是如何处理数据仍然是个难题。本文中关于系统发育分析方法的指导有如下一些目的。首先，我们要介绍一个概念性的步骤，这个步骤包含了当前在对生物序列进行系统发育分析时最广泛使用的方法中应用到的一些最重要的原理。

4、这个步骤认为，所有的系统发育方法都刻画了一个模型去模拟真实的进化过程，然后假定程序发现的序列之间的差异是模拟的进化结果，而这个模型本身就包含了方法中固有的假定。对于一个好的使用的系统发育的分析方法，刻画一个正确的模型同描绘一个正确的进化树是同等重要的。另外，我们还将讨论一些比较通用的方法所刻画的模型，并且讨论一下这些模型是如何影响对数据组的分析的。最后，我们会针对系统发育模型和进化史的推论，提供一些具体方法应用的例子。系统发育模型的组成系统发育的建树方法都会预先假定一个进化模型（Penny et al., 1994）。比如，所有广泛使用的方法都假定进化的分歧是严格分枝的，因此我们可以用树状拓扑

5、发生图来描述已知的数据。但是在一个给定的数据组中，因为存在着物种的杂交以及物种之间遗传物质的传递，这个假定很可能会被推翻。因此，如果所观察的序列并非是严格遗传的话，大多数系统发育方法就会得到错误的结果。系统发育方法中固有的模型会导致一些额外的假定：序列有指定的来源并且正确无误（Helbig and Seibold, 1996; Hershkovitz and Lewis, 1996; Soltis et al., 1997）。 序列是同源的（也就是说，所有的序列都起源于同一祖先序列）；这些序列不是“paralog“（paralog指的是一个祖先序列通过复制等方法在 HYPERLINK /bio

6、logy/postgenomics/ t _blank 基因组中产生的歧化序列）的混合物。 序列比对中，不同序列的同一个位点都是同源的。 在接受分析的一个序列组中，序列之间的系统发育史是相同的（比如，序列组中不存在核序列与细胞器序列的混合）。 样本足以解决感兴趣的问题。 样本序列之间的差异代表了感兴趣的宽组。 样本序列之间的差异包含了足以解决感兴趣的问题的系统发育信号。 还有一些假定，在某些方法中是默认的，但是在另外一些方法中是可以修正的：样本序列是随机进化的。 序列中的所有位点的进化都是随机的。 序列中的每一个位点的进化都是独立的。 许多公开发表的系统发育分析结果的错误通常都是因为违反了上述

7、假定，特别令人恼火的是现有的方法并不能察觉对基本假定的违反。每一个序列组都必须根据上述基本假定进行评估，并且同分析过程相协调。下面我们将会描述分析过程。系统发育数据分析：比对，建立取代模型，建立进化树以及进化树评估对DNA序列进行系统发育分析的四个主要步骤是比对，建立取代模型，建立进化树以及进化树评估。一般来说，执行这四个步骤的计算机程序都是相互独立的，但是他们都是进行系统发育分析的有机部分。当前的讨论主要集中在分析DNA序列的方法上，这些方法基本上（不全是）仍然是通过碱基和密码子的替代来考察序列的差异；这个方法同样应用于对蛋白质序列的分析，但是由于氨基酸的生物化学多样性，我们必须引入更多的数

8、学参数。因此，下面讨论的适用于核苷酸碱基的方法同样适用于氨基酸和密码子，而关于蛋白质的一些特殊问题和程序会在别的地方得到讨论（Felsenstein, 1996）。因为建树的标准在一定程度上依赖比对和取代模型，因此一开始就应该介绍这些方法，这些介绍非常必要。三种主要的建树方法分别是距离、最大节约（maximum parsimony, MP）和最大似然（maximum likelihood，ML）。距离树考察数据组中所有序列的两两比对结果，通过序列两两之间的差异决定进化树的拓扑结构和树枝长度。最大节约方法考察数据组中序列的多重比对结果，优化出的进化树能够利用最少的离散步骤去解释多重比对中的碱基差

9、异。最大似然方法考察数据组中序列的多重比对结果，优化出拥有一定拓扑结构和树枝长度的进化树，这个进化树能够以最大的概率导致考察的多重比对结果。这些方面将会在后面详细讨论。建立数据模型（比对）系统发育的序列数据通常会进行多重比对，也就是说，在系统发育的有关文献中通常会把那些单独的比对好的碱基所占据的位置称作位点；在进行系统发育的理论讨论时这些位点等同于特征符，而占据这些位点的真实的碱基或者空位被称为特征符状态；但是也有例外，STATALIGN程序（Thorne and Kishino, 1992）在进行系统发育分析时不会进行多重比对，它使用原始的未比对过的序列。多重比对在第七章中有详细讨论，在这里

10、我们只讨论与系统发育分析有关的多重比对方法。对准的序列位点将会被应用于系统发育分析，但是这只能描述一个先验的系统发育的结论，因为位点本身（不是真实的碱基）都是假定为同源的（Mindell, 1991; Wheeler, 1994）。因此，出于系统发育的目的，比对程序只是系统发育分析的一部分。建立一个比对模型的基本步骤包括：选择合适的比对程序；然后从比对结果中提取系统发育的数据集，至于如何提取有效数据，取决于所选择的建树程序如何处理容易引起歧义的比对区域和插入删除序列（即所谓的indel状态或者空位状态）。多重比对程序多重比对程序有如下一些属性：计算机依赖性：无依赖；部分依赖；完全依赖。系统发育

11、调用：无调用；先验调用；递归调用。比对参数评估：先验评估；动态评估；递归评估。比对特征：基本结构（比如序列）；高级结构。数学优化：统计优化；非统计优化。在系统发育分析研究中，一个典型的比对过程包括：首先应用CLUSTAL W程序，然后进行手工比对，最后提交给一个建树程序。这个过程有如下特征选项：（1）、部分依赖于计算机（也就是说，需要手工调整）；（2）、需要一个先验的系统发育标准（也就是说需要一个前导树）；（3）、使用先验评估方法和动态评估方法（推荐）对比对参数进行评估；（4）、对基本结构（序列）进行比对（对于亲水的氨基酸，推荐引入部分的二级结构特征）；（5）、应用非统计数学优化。这些特征选项

12、的取舍都依赖于下面将要描述的具体的系统发育分析方法。计算机依赖型提倡进行完全的计算机多重比对的人常常有这样的依据，他们认为手工比对通常都是含混不清的（Gatesy et al., 1993）。MALIGN（Wheeler and Gladstein, 1994）和TreeAlign（Hein, 1990, 1994）程序都试图在根据系统发育功能而优化比对的时候实现完全计算机化，尤其是试图通过初步的多重比对而优化出一个最大节约树（MP tree）。一般来说，只有当我们假定计算机程序能够进行正确的比对的时候，我们才会实行完全的计算机处理。之所以还提倡进行手工比对（比如，Thompson et al

13、., 1994），是因为现在的比对算法和程序还不能满足系统发育分析的需要。系统发育标准许多计算机多重比对程序（比如CLUSTAL, PileUp,ALIGN in ProPack）根据明确的系统发育标准（一个前导树）进行比对，这个前导树是由双重比对得到的。但是SAM（Hughey et al., 1996）和MACAW（Lawrence et al., 1993）程序在进行多重比对时并不引入明确的系统发育标准，虽然这些程序也可以模拟系统发育过程操作参数。如果在进行系统发育分析的时候，比对中引入了前导树，那么通过这个比对推导出的进化树逻辑上应该同前导树的拓扑结构相同。由CLUSTAL比对得到的前

14、导树（如图9.1）将会被转化成PHYLIP树的文件格式，然后输入到画树程序中，这些画树程序包括TreeTool(X windows), TreeDraw(Macintosh), PHYLODENDRON(Macintosh), TREEVIEW(Macintosh, Microsoft Windows) 或者PAUP(图 9.1；Macintosh, Microsoft Windows)的画树工具。按道理，我们应该回过头来为CLUSTAL比对再指定一个前导树，但是在实际操作中我们并不会这么做。有些程序（比如TreeAlign and MALIGN）为了得到优化的比对和系统发育树，程序本身就设计

15、了交叉（同步）递归优化的算法。理论上，能够解决比对系统发育难题的同步优化算法或者配套算法应该是存在的，但是递归算法必须冒一定的风险，它很可能会导致一个错误的或者不完整的结果（Thorne and Kishino, 1992）。因此，根据比对结果建立进化树之后，必须考虑另外的可能性，也就是说，如果根据其它的比对结果得到一个并不是最优化的进化树，这个次优化的进化树是不是更能够满足研究的需要。比对参数评估在比对中会出现一些序列区域，其长度是可变的，如何处理这些区域中indel状态的位点是最重要，这取决于进化模型的所有要素（比如，包括核苷酸转换/颠换速率），而且相关的参数在前导树与比对推导的进化树中应

16、该保持一致。比对参数应该随着进化的分叉动态变化（Thompson et al., 1994），只有这样才能保证碱基错配的几率能够满足序列趋异的需要；比对参数应该随时调整（Thompson et al., 1994, Hughey et al., 1996），以防止引入过多的近似序列而导致比对序列的信息量不足，可以通过降低近似序列的比对分值权重来防止这种情况。CULSTAL程序兼顾了这两种情况（参数动态变化），而SAM程序引入了序列权重。利用基本结构或者高级结构进行比对根据二级或者三级序列结构进行比对，比起直接利用一级序列进行比对的可信度要好，因为在同源性评估中，人们一直认为复杂结构的保守性高于

17、简单特征（核苷酸，氨基酸）的同源保守性，而且，立足于复杂结构的比对程序还可以 HYPERLINK / 搜索到一些特殊的关联位点，这些位点是进化的功能区域。实际上，基于系统发育的结构多重比对并没有将问题简化，也就是说，序列比对必须服从结构进化，而结构进化则同系统发育保持一致。有一个探索式的手工程序（如图9.2），是用来对核糖体DNA进行结构比对的（Gutell et al., 1994），这个程序要考察相关取代的样式，但是相关性必须通过系统发育树中的多个独立的补偿性突变推导得到（cf. Harvey and Pagel, 1991）。数学优化有些比对程序（比如，MACAW, SAM）根据一个统计

18、模型进行优化，但是这些统计同系统发育模型的关系并不清楚。仅仅根据一个系统发育模型是没有办法比较多重比对方法的优劣的。总结：对于系统发育分析，最好的比对程序是什么呢？对于一个给定的系统发育问题，如果预先不知道其中的系统发育关系，就很难知道最适用的比对方法（cf. Morrison and Ellis, 1997）。一般来说，我们并不赞成直接把计算机比对结果提交给建树程序，因为建树程序不能发现比对的错误；尤其是那些包含在比对程序包中（比如，CLUSTAL and TREE in ProPack）的建树程序，特别要注意这一点（Feng and Doolittle, 1996），因为在这些程序包中的建

19、树程序更加不严格。我们必须通过分子结构功能和碱基取代过程作出一些假定，并且结合另外一些独立的系统发育证据，对整个比对进行考察。从比对中提取系统发育数据集在某些比对中，比对长度是可变的，这时，系统发育数据集同比对就不会完全吻合；即使 在一些长度不变的比对中，数据集也可能同比对结果不一致-举一个很简单的例子，有时候我们只需要处理第一个和第二个密码子位点，就不需要全部的比对结果，这个话题我们在后面讨论取代模型的时候还会涉及到。如果比对中出现可变长度，我们通常会根据比对的不确定性程度和处理indel状态的原则这两个标准对比对结果进行取舍，从中选择所需的系统发育数据集；其中针对indel状态的处理方法取

20、决于建树方法以及从比对结果中发掘出的系统发育信息，最极端的方法是把包括空位在内的所有indel位点从比对中清除出去，在分析时不加考虑（cf. Swofford et al., 1996a），这个方法的好处是可以把序列的变化包容在取代模型中，而不需要特别的模型来处理indel状态，但是它的缺点也很明显：indel区域的系统发育信息完全被忽略了。在提取数据集时保留indel区域但是忽略所有的空位分值，将会保留包括空位在内的位点碱基变化信息。某些长度可变区域在部分序列或者全部序列中很难对准，在这种情况下，这些难以对准的碱基的分值应该清零；这个方法存在很大的缺陷，MP和ML建树方法会不加考虑地把这些清

21、零的或者被忽略的分值理解为零分歧，但是实际上隐藏在这些分值下面的实际的数据（不管是空位还是难以对准的碱基），一般来说，反映出的分歧度都很大。PAUP 4.0中的距离建树方法（将在下文描述）允许通过非空位区域外推得到空位区域的距离。最大节约（MP；见下）是允许把可比对的空位合并，并将其视为特征符的唯一方法；可以通过两种方式达到目的：作为一个附加的特征符状态（第五种核酸碱基或者第二十一种氨基酸），或者作为一套独立于碱基取代的特征符集。当空位占据了不止一个位点的时候，前一个方法行不通，因为每一个空位位点都会被统计为一次独立的特征符状态变化。当比对的序列的局部出现很好的可比对的空位的时候，后一种方法非

22、常有用。我们可以把一套空位特征符附加到比对序列数据集中，也可以用额外的碱基程序在适当的位置对空位计分，但是在计分的时候，空位位点中只有一个作为空位计分，其余的将会被忽略。PAUP将会执行这个方法。对于某些比对而言，比对程序会忽略所有的空位分值或者忽略所有低于预设值的空位分值；但是，还没有任何一种程序会忽略单个序列的单个位点。如果比对在序列组内部相当明确，但是处身其中时却不太清楚，此时必须对比对做“手术”，确保同序列组相关的明确的信息被保留，而除去模糊的信息。图9.3给出了一个比对“手术”的例子。在空位区域，我们必须作出决定：在可供选择的比对中，哪一个更加合理，尤其重要的是，哪一个更加适合于建立

23、进化树分析。如果手工解决比对的不确定性，就必须考虑系统发育关系、取代过程（比如，转换和颠换）和碱基组成；在这个阶段，用系统发育证据解决不确定性非常合理。在倾向于变长的序列区域，关系非常疏远的序列和序列组的比对就可以侧向展开（就是说，引入人工空位，并且忽略分值），最终的结果使得只有关系很近的序列区域对准在一起。某些序列中的某些位点虽然对准了，但是并不确定，他们的分值可以在计分时忽略；这个方法的优点是可以保留同这些序列相关的明确信息，缺点是最大节约和最大似然的建树方法会把这些“缺失”的分值看做是零分歧。由MALIGN（Wheeler and Gladstein, 1994）和TreeAlign得到

24、的比对不需要在比对后用这些方法中的建树方法进行数据修饰，即使这些比对中仍然有一些同样类型的不确定性，这些不确定性在另外一个程序进行分析时需要修正。如前所述，这些程序会根据由比对得到的最好的MP系统发育进化树，对比对参数进行递归优化。MALIGN还会利用一套空位为代价，对以连接的可供选择的比对为基础的建树方法进行优化；在这个方法中，在最有可能的几种比对中出现的比对特征将会被加权。这就提供了一种方法，可以捕获序列分歧的数量（在取消不确定的比对区域的分值的时候这些分歧是被忽略掉的），因为这些区域的所有可能的比对方式都将显示这些区域的最大的序列分歧。处于不确定的比对区域中的位点很可能不是同源的，因此在

25、进化树中需要加入一些噪声干扰或者偏向。决定取代模型对于取代模型，应该给予同比对和建树同样的重视。就像前面暗示的那样，取代模型既影响比对，也影响建树；因此需要采用递归方法。现在，对于核酸数据而言，可以通过取代模型中的两个要素进行计算机评估（Swofford, 1997），但是对于氨基酸和密码子数据而言，没有什么评估方案（Felsenstein, 1996）。其中一个要素是碱基之间相互取代的模型；另外一个要素是序列中不同位点的所有取代的相对速率。还没有一种简单的计算机程序可以对较复杂的变量（比如，位点特异性或者系统特异性取代模型）进行评估，同样，现有的建树 HYPERLINK /soft/ 软件也

26、不可能理解这些复杂变量。碱基取代速率模型一般而言，生物化学性质相近的碱基之间的取代频率较高；在DNA中，四种转换（A G, G A, C T, T C）的频率比八种颠换（A C, A T, C G, G T, 以及前四种的反向取代）的频率要高；这些偏向会影响两个序列之间的预计的分歧。各个残基之间的相对取代速率一般由方阵形式列出；对于碱基而言，行数和列数都是4，对于氨基酸而言，行数和列数都是20（比如PAM方阵），对于密码子而言，行数和列数都是61（除去了中止密码子）。非对角线元素对应于一个碱基变为另一个碱基的相对代价，而对角线元素则代表不同序列拥有同一个碱基的代价。这些代价值可以固定为先验的代

27、价表，以确保建树方法在计分时对每一种取代都使用确定的代价值。固定的代价方阵是典型的静态权重方阵，MP建树方法（如图9.4）使用的就是这种方阵。如果使用这种权重，那么这个方法就会被称为“加权节约”。又如，ML建树方法，代价值是由即时的速率方阵得到的，这个方阵（如图9.5）代表了各种取代可能会发生的概率的ML估计值。MP权重方阵只涉及简单的算术，而应用距离和ML速率方阵则可以引入复杂的代数。为了避免盲目使用不适当的方法，建议大家熟悉其内部的基本原理（见Li, 1997, and / or Swofford et al., 1996a）。实际上，“前进”和“反向”取代速率被认为是相同的；这个取代模型

28、被称为是“时间可逆”；这个模型拥有“静态”的性质，因为在所有的碱基频率中没有预知的变化。在系统发育的特殊历史中，不同序列中的碱基频率不同表明，前进和反向速率实际上可能会不同；而传统的取代权重或者速率方阵不能包容这个“非静态”环境；本节的结尾将讨论一个基于非静态取代模型（“logdet”）的建树方法，这个方法将会提供一个可供选择的计算方法。通常，特征符状态的权重方阵都会或多或少地通过观察进行过估值，当然也可以从速率矩阵衍生得到。比如，如果假定两个转化的其中一个，发生的频率是每个颠换的两倍，那么据此就可以确定一个权重方阵，比如，AG的转换代价为1，而AT的颠换代价为2（图9.4）。（节约方法规定对

29、角线元素值，或者说是不同序列中拥有相同碱基的代价值为零。这是节约方法的一个缺点详见下文）在随后的建树步骤中，这套假定会把颠换的总数降至最低值，而力求把那些主要差异是转换的序列集中在一起。任何一种“时间可逆”的核苷酸取代模型都可以用图9.5所示的方阵刻画，只是其中一个速率和其它速率的差异；在任意组合中，最多可以达到只有六个参数，其中每一个速率参数都是独立的（Swofford et al., 1996a；Li, 1997）。如果平衡的碱基频率不相等，则需要额外的参数；如果平衡的碱基频率不相等，但是却假定这些频率相等，那么系统发育进化树的最终结果将会出错（Li, 1997）。侧线（paralinea

30、r）（Lake, 1994）和“logdet”（Lockhart et al., 1994）做了一些修正（见Swofford et al., 1996a）来满足非静态环境的需要；这个方法只适用于距离进化树的建立；在这个方法中，对于每一个序列匹配，各种类型和变化方向的原始取代的数目都会计算在一个44的方阵中（如图9.6）。每个方阵都会有一个代数行列式，这个行列式的log值是评估序列差异性的一个要素，因此被称为“logdet”。对那些拥有各种各样的碱基频率的序列进行双重比较，就会得到各种各样的方阵，也就会得到各种各样的行列式值；因此，在评估序列两两之间的距离的时候，就要受到序列两两之间的行列式值的

31、影响，而且序列两两之间的比较允许适用不同的取代模型，因此沿着系统发育进化树的不同树枝，将会产生多元化。Logdet尤其对位点之间的速率差异（见下文）敏感，因此，碱基频率的偏向可能只存在于那些承受变化的位点。位点内速率差异模型除了取代模型的多元化以外，一个序列中各个不同位点之间取代速率的差异也会对建立进化树的结果产生很深远的影响（Swofford et al., 1996a）；关于位点之间的速率差异（或者叫做位点异质性），有一个最明显的例子，就是在一个编码序列中，三联体编码的位点差异：在三联体编码中，第三个编码位点比另外两个位点更加容易发生变化；正是出于这个原因，许多系统发育分析方法在分析编码序

32、列时，都会把第三个编码位点排除在外；但是在某些情况下，速率差异模型会更加敏锐（比如，对应于蛋白质或者rRNA的保守序列）。对位点差异的取代速率进行估值的方法有非参数化模型（W.M. Yang et al., 1996），不变式模型和gamma分布模型（Swofford et al., 1996a）。非参数化方法源于特异位点的相对速率的范畴；这个方法可以在MP建树方法中使用，只要根据相对的变异频率对特异位点进行简单加权就可以了，当然进行加权时需要有关于真实进化树的预备知识；这个方法同样也可以适用于ML建树方法，但是在计算上被认为是不切实际的（W.M. Yang et al., 1996）。不变式

33、模型对一定比例的位点进行估值，这些位点并不能自由变化；剩余的位点假定为等概率变化。至于gamma模型方法，它假定一个给定的序列变化的概率是遵守gamma分布规律的，据此指定位点的取代概率；gamma分布的形状（有形状参数描述）描述了一个序列中各个位点的取代频率的分布（Swofford et al., 1996a, p. 444, Figure 13; cf. Li, 1997, p. 76, Figure 3.10; 注意尺度差异）。在一个混合方法中，可以假定一部分位点是不变的，而剩余的位点则是按照gamma分布变化的。实际上，gamma修正可以是连续的，离散的或者自离散的（W.M. Yang

34、 et al., 1996）。连续gamma的意思是各个位点沿一条连续的概率曲线变化；目前，这个方法在绝大多数情况下无法计算。离散gamma逼近方法指定各个位点的概率，使得这些（大量的）概率值逼近gamma曲线。自离散模型假定相邻的位点的变化速率是相关联的；许多组位点被分为许多类，其中每一类中的位点的变化速率可能被假定为常量或者异类值。进化树建立程序使用各种各样的位点速率差异修正方法。对于核酸数据，PAUP 4.0在单独或者混合使用时间可逆的距离建树方法和最大似然建树方法时，既使用不变式gamma模型，也使用离散的gamma模型；在使用logdet距离建树方法时，使用不变式模型（见下）。对于核

35、酸，氨基酸和编码子数据，PAML使用连续的，离散的和自离散的gamma模型。对于核酸和氨基酸数据，PHYLIP使用一种离散的gamma模型。究竟使用哪一种取代模型呢？在一个给定的序列集中，即使一个取代模型中的任何参数都可以证明是很有判断力的，但是最好的模型并不一定总是拥有最多的蚕室；相反，参数越少，模型越好，这是因为对每一个参数进行估值都会引入一个相关的变量，而每引入一个附加的参数维度，都会使得整体的变数增加，有时候甚至会对模型起抑制作用（见Li, 1997: p. 84, 表4.1）。对于一个给定的序列比较，如果模型只有两个参数，那么碱基差异的总和要被划分为两个类别；如果模型有六个参数，那么

36、碱基差异的总和要被划分为六个类别；很明显，如果分为六个类别，那么每一个类别中的位点样本数目将会很小，很可能小得无法进行合理地估值。PAUP的“描述树”特征是对DNA序列的取代模型进行规范的一个较好的策略，它使用似然方法同时评估六个可逆的取代速率、gamma分布的形状参数和不变的位点的比例（图9.7）。这些参数可以通过相等的或者指定的碱基频率进行估值。通常，任何一个合理的系统发育进化树（比如，很容易就可以得到一个相邻连接的进化树）都适用于这个程序，因为很明显，对参数的估值在很大程度上受到特征符模式的影响，而不是进化树的拓扑结构（Swofford et al., 1996b）。这个估值程序对于50

37、个序列而言，并不会耗费太多的时间。如果序列较多，或者时间较紧，可以对试验的进化树进行精简，在保留全部的系统发育范围和结构的同时，减少分类数目。通过这些估算的取代参数，我们可以通过比较由较多参数和较少参数分别评估得到的似然分值，决定一个简化的模型是否合理（比如，六个取代类别是否可以减少到两个）。有时候，参数和不变位点的比例可以相互替换，所以我们应该比较每一个单独使用时得到的似然分值和两个同时使用时得到的似然分值。注意，和MP以及ME不同，用不同的参数值得到ML分值可以直接比较（Swofford et al., 1996b）。对于编码蛋白质的DNA序列，根据样本的分歧程度，有时候很明显地，有用的变

38、化基本上都是第一位和第二位的编码位点，而在整个数据集中，第三位点通常都是随机的，或者第三位点变化而第一位点和第二位点不变。尽管除去“无用”位点可以提高剩余位点的不同速率的估值精确度，上面所述的程序还是要对这个速率差异进行修正。对于节约进化树的建立，我们有时候会把随机的第三位点从分析中除去，因为这些位点只会引入噪声，如果碱基频率不相等，这些位点还会引入错误。怎么样才能确定数据集中的非静态因素是否会成为一个问题呢？最简单的方法可能就是去比较PAUP中两种通过不同方法得到的建树结果和进化树评估结果，一种方法是使用时间可逆的方法，另一种方法是logdet距离建树方法。下面的章节中将会涉及到这些程序。还

39、没有什么好的计算方法，能够直接从序列数据中评估非静态因素的影响。PAUP中有一个命令，会列出所有序列的碱基频率。这个程序应该使用排除不变位点（Exclude Constant Sites）的选项。序列中的碱基频率可以很直观地比较出来。数据文件应该指定 gapmode=missing，或者PAUP要把一个空位特征符计算为一个碱基变化。碱基频率的命令还要对数据执行一个“chi平方（chisquare）”测试，但是这个测试并不切实际，因为它假定数据是从一个随机样本中提取的，而结果就会假定所观察到的不相等是相互独立的而不是系统发育结构的结果。一个毫无意义的chisquare分值并不能除去非静态因素，而

40、一个有意义的分值很可能会进一步证实非静态因素。PAUP中的碱基组成命令已经被用来证明：在被子植物和绿 HYPERLINK /biology/Special/zaolei/Index.shtml t _blank 藻的5.8S rDNA序列中，两者之间变化比较大的位点集中存在中碱基偏好，但是并不一定在某个序列上产生变化（Hershkovitz and Lewis, 1996）。建树方法现有的软件中所采用的建树方法已经在某些著作（Saitou, 1996; Swofford et al., 1996a; Li, 1997）中有很详尽的讨论；这部分只是简单地描述一些最常用的方法。建树方法可以分为两类

41、，每类有两种不同的方法：基于算法的和基于标准的。基于算法的建树方法根据一系列的步骤得到一个进化树；而基于标准的建树方法则是根据一些优化的功能对可选的进化树进行评估。相邻连接方法（NJ）是一个纯粹的基于算法的建树方法，这个方法只得到一个进化树，这个进化树拥有令人满意的性质，它的距离附加值接近或者是非常优化的（见下）；而一个基于标准的距离建树方法将会根据进化树的附加值优化的标准，对所有可能的进化树（不管这些进化树是否产生）进行评估。 基于距离的和基于特征符的。历史上的和现有的许多关于系统发育的讨论描述了各种各样的基于距离的和基于特征符的建树方法的效用（比如，Saitou, 1996: Li, 19

42、97）。距离建树方法根据一些尺度计算出双重序列的距离，然后抛开真实数据，只是根据固定的距离建立进化树；而基于特征符的建树方法在建立进化树时，优化了每一个特征符的真实数据模式的分布，于是双重序列的距离不再固定，而是取决于进化树的拓扑结构。最常用的基于特征符的建树方法包括MP和ML。 距离建树方法距离建树方法根据双重序列比对的差异程度（距离）建立进化树。如果所有的起源分歧事件都很精确地记录在序列中（Swofford et al., 1996a），那么距离建树方法将会重构真实的进化树。然而，当序列突变达到饱和时，分歧程度就会达到上限。一对分歧序列的其中之一在某个特定位点发生突变后，后续的突变不管发生

43、在哪一条序列，都不会再造成更多的差异。实际上，后续的变异很可能会使得前一个位点恢复成相等的状态，从而屏蔽掉前一个突变，因此，绝大多数基于距离的建树方法都会针对这样“不可见”的替换进行修正。在实际操作中，使用的速率方阵都会很有效地假定：在所观察的相等的碱基对中，实际上存在一定比例的位点，这些位点经受了多次突变，而且，随着整个序列差异程度的增加，这些位点的比例也在增加。有些程序（至少是随意地）会计算出不正确的距离进行，比方说，MEGA程序（Kumar et al., 1994）只针对密码子和氨基酸数据会计算出错误的距离；除非序列的整体差异程度很小，这个程序实质上将会保证给出错误的结果。双重序列差异

44、使用最大似然方法的取代速率计算得到。最常用的距离建树程序使用的模型局限于时间可逆模型，只包含很有限的几个取代模型；但是PAUP 4.0 基本上会测试时间可逆模型的所有变化，包括用最大似然方法从数据中评估得到的实际模型，以及针对非静态数据的logdet距离建树方法所得到的模型。与最大似然方法相比，距离建树方法的计算强度很小，但是可以使用序列进化的相同模型，这是它们最主要的优点。这个方法的缺点是屏蔽了真实的特征符数据。最常用的距离建树方法是“不加权配对组算术方法”（UPGMA，unweighted pair group method with arithmetic mean）、相邻连接方法（NJ，

45、neighbor joining）和对距离进化树的附加值进行优化的方法，包括最小进化方法（ME，minimum evolution）。有些方法被不止一个系统发育软件包引用，但是在使用时不一定会使用相同的指定参数和（或）进化树优化特征（比如，树枝交换见下）。不加权配对组算术方法（UPGMA，unweighted pair group method with arithmetic mean）UPGMA是一种聚类或者说是分类方法它按照配对序列的最大相似性和连接配对的平均值的标准将进化树的树枝连接起来。它还不是一种严格的进化距离建树方法（Li, 1997）。只有当序列分歧是基于一个分子钟或者近似等于原

46、始的序列差异性的时候，我们才会期望UPGMA会产生一个拥有真实的树枝长度的准确的拓扑结构（ultrametric; Swofford et al., 1996a）。正如前面讨论的那样，在实际问题中很少会遇到这种情况。相邻连接方法（NJ，neighbor joining）相邻连接算法在距离建树中经常会用到，而不会理会使用什么样的优化标准。完全解析出的进化树是通过对完全没有解析出的“星型”进化树进行“分解”得到的，分解的步骤是连续不断地在最接近（实际上，是最孤立的）的序列对中插入树枝，而保留进化树的终端（图9.8）。于是，最接近的序列对被巩固了，而“星型”进化树被改善了，这个过程将不断重复。这个方

47、法相对而言很快，也就是说，对于一个50个序列的进化树，只需要若干秒甚至更少。FitchMargoliash（FM）FitchMargoliash（FM）方法设法通过把所有可能观察到的距离相对于进化树中所有可能的路径长度的偏差的平方极小化，将进化树中观察到的双重距离的合适度极大化（Felsenstein, 1997）。有一些变量，在怎样对错误进行加权的处理上不太一样。对变化的估值并非是完全独立的，因为所有的进化树内部的树枝中存在的错误都至少被计算了两次（Rzhetsky and Nei, 1992）。最小进化方法（ME,Minimum Evolution）最小进化方法先使用同FM相同的方式计算出

48、路径长度，然后根据路径长度优化出最短的进化树；也就是说，它要求将观察到的距离相对于基于进化树的距离的偏差的平方最小化（Rzhetsky and Nei, 1992; Swofford et al., 1996a; Felsenstein, 1997）。同FM方法不同，ME方法并不使用所有可能的双重序列距离和所有可能的相关的进化树路径长度，而是先根据到外层节点的距离固定进化树内部节点的位置，然后根据这些观察点之间的最小计算误差，对内部的树枝长度进行优化。于是这个方法声称它消除了FM计算的依赖性。究竟哪一个基于距离的建树程序是最好的ME和FM似乎是最好的程序程序，它们在模拟研究中所取得的成绩几乎相

49、同（Huelsenbeck, 1995）。ME在计算机程序中的应用越来越广泛了，包括METREE（Rzhetsky and Nei, 1994）和PAUP都在使用ME。对于蛋白质数据，PHYLIP中的FM程序提供了最多的时间可逆取代模型，但是没有对位点内部的取代速率差异进行修正。MEGA（Kumar et al., 1994）和METREE软件包包括一个针对蛋白质的gamma修正，但是只有同一个原始（“p距离”）的分歧模型（没有距离修正或者偏好修正）结合时才进行修正，这种方法只有当分歧很小时才合理（Rzhetsky and Nei, 1994）。MEGA也计算同义的和异义的位点的分离距离，但是

50、只有当没有取代或者碱基频率偏好的时候，以及没有对位点内部的速率差异进行修正的时候，这个方法才合理。因此，对于绝大多数数据集而言，对核苷酸数据应用一个更加理想的模型可能会比MEGA方法要好。模拟研究指出，对于一个大范围的进化树形状空间，UPGMA的可操作性很差（Huelsenbeck, 1995）。我们并不赞成使用这种方法，而之所以在这里会提及这个方法，是因为在现在的出版物上经常会出现UPGMA的应用，这一点可以由当前的出版物中出现的UPGMA“基因进化树”证明（Huelsenbeck, 1995）。很显然，NJ是最快的程序，并且所产生的进化树同ME进化树相比，虽不能说一样，但也已经非常相近了（

51、Rzhetsky and Nei, 1992; Li, 1997）。但是，NJ只产生一个进化树。根据数据结构，有大量不同的进化树可能和NJ进化树一样好，甚至比NJ进化树好得多（Swofford et al.1996a）。我们中的一个（MAH）和D.Hillis（个人通讯）已经各自独立地得到了一些（虽然很少）数据集，对于这些数据集，NJ给出的进化树很糟糕；但是使用PAUP方法中的“最接近（closest）”选项后，进行逐步加总，得到了一个较好的距离进化树（Swofford， 1997）。特征符建树方法除了在所有的分析步骤中都会用到特征符数据以外，基于特征符的建树方法毫无共同之处。这就使得我们能够

52、评估一个比对中每一个碱基位点对其它所有的碱基位点的依赖性。最常用的基于特征符的建树方法是最大节约方法和最大似然方法。最大节约方法（MP，Maximum Parsimony）最大节约方法是一种优化标准，这个标准遵循“奥卡姆剃刀原则（Occams razor）”：对数据最好的解释也是最简单的，而最简单的所需要的特别假定也最少。在实际应用中，MP进化树是最短的也是变化最少的进化树，根据定义，这个进化树的平行变化最少，或者说是同形性最低。MP中有一些变量与特征符状态改变的可行方向不尽相符（Swofford et al., 1996a）。为了能够包容取代偏好，MP必须加权；比如，颠换的变换相对于转换（见

53、上文）被加权。进行加权的最简单的方法就是建立一个加权步骤方阵，在这个方阵里，权重用速率的倒数，这个速率是由上面描述的ML方法评估得到的。步骤方阵加权可能会极大地减慢MP的计算速度。如果比对内部的位点确实存在这速率差异，那么MP方法就会执行得很不好（Huelsenbeck, 1995）。对于这个问题，还没有什么很好的解决方法。有一个方法，就是对数据集进行修改，使得最终分析的数据集中只包括那些内部差异性很小的位点，这些位点是由前面提到的似然分析方法判断得到的。更常见的情况是，MP分析只简单地屏蔽掉那些高度同形的可疑位点（比如说，某些序列比对中的第三个编码位点）。另外一个方法是根据引导树中所能观察到

54、的位点变化倾向，对位点进行循环加权。这个“连续逼近”方法在PAUP中会自动简化，但是这个方法的出错的倾向程度取决于引导树的出错程度。MP分析期望能够产生大量的（有时候是成千上万个）具有相同分值的进化树。因为这些进化树之间的优化程度相同，所以只有所有进化树都严格一致的分组才算是得到数据的支持。距离建树方法和ML建树方法期望得到一个最简单的最好的进化树，是因为这些方法的算法中包括了除法和小数，而MP建树方法只是计算离散的步长。对于一个给定的数据集，存在着比最佳的进化树差不了多少的ME进化树或者ML进化树，如果这些进化树严格一致的话，很可能会产生同MP一致进化树多少有些可比性的分解。不幸的是，MP的

55、使用者传统上都会给出严格一致的进化树（有时候是进化树的一致性，其中有一个或者两个步骤不太好），而ME和ML用户通常就不会这样做。现在已经开发了大量的进化树统计量值来描述一个MP进化树。除了进化树的长度外，这些统计量包括各种各样的比例，用以描述同形的数量，比如稠度指数，或者一个进化树中一个特征符改变状态的平均次数（Swofford, 1990）。尽管在MP分析中仍然经常报告这些指数，但是已经证明这些指数受到许多人为因素的干扰，比如被分析的分类群的数量。模拟研究（Huelsenbeck, 1995）已经表明，如果从世系分歧开始的序列进化的数量远大于世系分离的分歧数量（比如，一个进化树中，终端的树枝

56、很长，而内部节间的树枝很短），MP的效果同ME一样差，比ML还要更差。这种情况会造成“长树枝效应”长树枝可能会被人工连接，因为累积的非同源的相似性的数量超过了因为真实的相关性而保留的同源相似性的数量（Swofford et al., 1996a）。在这种情况下，特征符加权会改善MP的执行效果（Huelsenbeck, 1995）。最大似然方法（ML, Maximum Likelihood）ML对系统发育问题进行了彻底搜查。ML期望能够搜寻出一种进化模型（包括对进化树本身进行 HYPERLINK / 搜索），使得这个模型所能产生的数据与观察到的数据最相似。实际上，ML是从比对的每一个碱基位点衍生

57、出来的。给定一个特定的进化树和观察到的全部的碱基频率，我们可以计算出似然值，具体方法是要计算一个位点遵循一个特定取代过程时所得到的变化模式的概率；似然值就是把在这个特定的取代过程中每一个可能的取代的再现的概率进行加和。所有位点的似然值相乘就得到了整个进化树的似然值（也就是说，数据集的概率给出了进化树和进化过程）。大家可以想象一下，对于一个特定的进化树，数据集的似然值在某些位点偏低，而另外一些位点偏高。如果进化树比较好，那么大多数位点的似然值都会较高，因此整个似然值较高；如果进化树不太好，似然值就会比较低。如果数据集中没有系统发育的信号，所有随机的进化树的似然值上都会相差无几。取代模型应该得到优

58、化，以适应观察到的数据的需要。比方说，如果存在着转化的偏好（其明显表现为有大量的位点只包含嘌呤或者只包含嘧啶），那么，如果计算数据的似然值时所采用的模型没有考虑偏好的话，其效果显然不如采用考虑了偏好的模型。同样地，如果有一部分位点确实只包含一种碱基，而另外一部分位点以相同的概率包含各种碱基，那么，如果计算数据的似然值时所采用的模型假定所有位点的进化都平等的话，其效果显然不如采用考虑了位点内部的速率差异的模型。对于一个特定的进化树，改变取代参数就意味着将改变与之相关联的数据集的似然值；因此，在某一个取代模型下，进化树可以取得很高的似然值，但是，在另一个取代模型下，进化树所取得的似然值就可能会很低

59、。因为ML要耗费大量的计算机机时，所以对于一个给定的数据集，想要在优化取代模型和进化树的同时进行完全的 HYPERLINK / 搜索几乎是不可实现的。值得推荐的是一个很经济的探索式的程序（Adachi and Hasegawa, 1996; Swofford et al., 1996a）。在这一点上，可能最好的节约时间的方法是前面所说的对取代模型的ML评估（图9.7）。这个程序可以反复进行， HYPERLINK / 搜索到较好的ML进化树，然后重新对参数进行评估，然后 HYPERLINK / 搜索更好的进化树。随着算法、计算机和对系统发育的认识的进步，ML标准在分子系统发育分析中也变得越来越流

60、行。在模拟研究中，如果对由相同的模型产生的数据进行数据分析的话，ML做得总是比ME和MP要好（Huelsenbeck, 1995）。在所有的方法中，ML的计算强度最大，所以在某些情况下，它总是无法实现；而且，同样的模拟研究表明在许多情况下，ME和MP方法同ML方法的执行效果一样好（或者一样差）。距离方法、节约方法和最大似然方法的差异距离方阵方法简单的计算两个序列的差异数量。这个数量被看作进化距离，而其准确大小依赖于进化模型的选择。然后运行一个聚类算法，从最相似（也就是说，两者之间的距离最短）的序列开始，通过距离值方阵计算出实际的进化树，或者通过将总的树枝长度最小化而优化出进化树。用最大节约方法