许多计算机多重比对程序(比如CLUSTAL, PileUp,ALIGN in ProPack)根据明确的系统发育标准(一个前导树)进行比对,这个前导树是由双重比对得到的。但是SAM(Hughey et al., 1996)和MACAW(Lawrence et al., 1993)程序在进行多重比对时并不引入明确的系统发育标准,虽然这些程序也可以模拟系统发育过程操作参数。 如果在进行系统发育分析的时候,比对中引入了前导树,那么通过这个比对推导出的进化树逻辑上应该同前导树的拓扑结构相同。由CLUSTAL比对得到的前导树(如图9.1)将会被转化成PHYLIP树的文件格式,然后输入到画树程序中,这些画树程序包括TreeTool(X windows), TreeDraw(Macintosh), PHYLODENDRON(Macintosh), TREEVIEW(Macintosh, Microsoft Windows) 或者PAUP(Macintosh, Microsoft Windows)的画树工具。按道理,我们应该回过头来为CLUSTAL比对再指定一个前导树,但是在实际操作中我们并不会这么做。有些程序(比如TreeAlign and MALIGN)为了得到优化的比对和系统发育树,程序本身就设计了交叉(同步)递归优化的算法。理论上,能够解决比对----系统发育难题的同步优化算法或者配套算法应该是存在的,但是递归算法必须冒一定的风险,它很可能会导致一个错误的或者不完整的结果(Thorne and Kishino, 1992)。因此,根据比对结果建立进化树之后,必须考虑另外的可能性,也就是说,如果根据其它的比对结果得到一个并不是最优化的进化树,这个次优化的进化树是不是更能够满足研究的需要。 比对参数评估 在比对中会出现一些序列区域,其长度是可变的,如何处理这些区域中indel状态的位点是最重要,这取决于进化模型的所有要素(比如,包括核苷酸转换/颠换速率),而且相关的参数在前导树与比对推导的进化树中应该保持一致。比对参数应该随着进化的分叉动态变化(Thompson et al., 1994),只有这样才能保证碱基错配的几率能够满足序列趋异的需要;比对参数应该随时调整(Thompson et al., 1994, Hughey et al., 1996),以防止引入过多的近似序列而导致比对序列的信息量不足,可以通过降低近似序列的比对分值权重来防止这种情况。CULSTAL程序兼顾了这两种情况(参数动态变化),而SAM程序引入了序列权重。 利用基本结构或者高级结构进行比对 根据二级或者三级序列结构进行比对,比起直接利用一级序列进行比对的可信度要好,因为在同源性评估中,人们一直认为复杂结构的保守性高于简单特征(核苷酸,氨基酸)的同源保守性,而且,立足于复杂结构的比对程序还可以搜索到一些特殊的关联位点,这些位点是进化的功能区域。实际上,基于系统发育的结构多重比对并没有将问题简化,也就是说,序列比对必须服从结构进化,而结构进化则同系统发育保持一致。有一个探索式的手工程序(如图9.2),是用来对核糖体DNA进行结构比对的(Gutell et al., 1994),这个程序要考察相关取代的样式,但是相关性必须通过系统发育树中的多个独立的补偿性突变推导得到(cf. Harvey and Pagel, 1991)。 数学优化 有些比对程序(比如,MACAW, SAM)根据一个统计模型进行优化,但是这些统计同系统发育模型的关系并不清楚。仅仅根据一个系统发育模型是没有办法比较多重比对方法的优劣的。 从比对中提取系统发育数据集 在某些比对中,比对长度是可变的,这时,系统发育数据集同比对就不会完全吻合;即使 在一些长度不变的比对中,数据集也可能同比对结果不一致--举一个很简单的例子,有时候我们只需要处理第一个和第二个密码子位点,就不需要全部的比对结果,这个话题我们在后面讨论取代模型的时候还会涉及到。 如果比对中出现可变长度,我们通常会根据比对的不确定性程度和处理indel状态的原则这两个标准对比对结果进行取舍,从中选择所需的系统发育数据集;其中针对indel状态的处理方法取决于建树方法以及从比对结果中发掘出的系统发育信息,最极端的方法是把包括空位在内的所有indel位点从比对中清除出去,在分析时不加考虑(cf. Swofford et al., 1996a),这个方法的好处是可以把序列的变化包容在取代模型中,而不需要特别的模型来处理indel状态,但是它的缺点也很明显:indel区域的系统发育信息完全被忽略了。 在提取数据集时保留indel区域但是忽略所有的空位分值,将会保留包括空位在内的位点碱基变化信息。某些长度可变区域在部分序列或者全部序列中很难对准,在这种情况下,这些难以对准的碱基的分值应该清零;这个方法存在很大的缺陷,MP和ML建树方法会不加考虑地把这些清零的或者被忽略的分值理解为零分歧,但是实际上隐藏在这些分值下面的实际的数据(不管是空位还是难以对准的碱基),一般来说,反映出的分歧度都很大。PAUP 4.0中的距离建树方法(将在下文描述)允许通过非空位区域外推得到空位区域的距离。 最大节约是允许把可比对的空位合并,并将其视为特征符的唯一方法;可以通过两种方式达到目的:作为一个附加的特征符状态(第五种核酸碱基或者第二十一种氨基酸),或者作为一套独立于碱基取代的特征符集。当空位占据了不止一个位点的时候,前一个方法行不通,因为每一个空位位点都会被统计为一次独立的特征符状态变化。当比对的序列的局部出现很好的可比对的空位的时候,后一种方法非常有用。我们可以把一套空位特征符附加到比对序列数据集中,也可以用额外的碱基程序在适当的位置对空位计分,但是在计分的时候,空位位点中只有一个作为空位计分,其余的将会被忽略。PAUP将会执行这个方法。 对于某些比对而言,比对程序会忽略所有的空位分值或者忽略所有低于预设值的空位分值;但是,还没有任何一种程序会忽略单个序列的单个位点。如果比对在序列组内部相当明确,但是处身其中时却不太清楚,此时必须对比对做“手术”,确保同序列组相关的明确的信息被保留,而除去模糊的信息。 在空位区域,我们必须作出决定:在可供选择的比对中,哪一个更加合理,尤其重要的是,哪一个更加适合于建立进化树分析。如果手工解决比对的不确定性,就必须考虑系统发育关系、取代过程(比如,转换和颠换)和碱基组成;在这个阶段,用系统发育证据解决不确定性非常合理。在倾向于变长的序列区域,关系非常疏远的序列和序列组的比对就可以侧向展开(就是说,引入人工空位,并且忽略分值),最终的结果使得只有关系很近的序列区域对准在一起。某些序列中的某些位点虽然对准了,但是并不确定,他们的分值可以在计分时忽略;这个方法的优点是可以保留同这些序列相关的明确信息,缺点是最大节约和最大似然的建树方法会把这些“缺失”的分值看做是零分歧。 由MALIGN(Wheeler and Gladstein, 1994)和TreeAlign得到的比对不需要在比对后用这些方法中的建树方法进行数据修饰,即使这些比对中仍然有一些同样类型的不确定性,这些不确定性在另外一个程序进行分析时需要修正。如前所述,这些程序会根据由比对得到的最好的MP系统发育进化树,对比对参数进行递归优化。MALIGN还会利用一套空位为代价,对以连接的可供选择的比对为基础的建树方法进行优化;在这个方法中,在最有可能的几种比对中出现的比对特征将会被加权。这就提供了一种方法,可以捕获序列分歧的数量(在取消不确定的比对区域的分值的时候这些分歧是被忽略掉的),因为这些区域的所有可能的比对方式都将显示这些区域的最大的序列分歧。处于不确定的比对区域中的位点很可能不是同源的,因此在进化树中需要加入一些噪声干扰或者偏向。
对于取代模型,应该给予同比对和建树同样的重视。就像前面暗示的那样,取代模型既影响比对,也影响建树;因此需要采用递归方法。现在,对于核酸数据而言,可以通过取代模型中的两个要素进行计算机评估(Swofford, 1997),但是对于氨基酸和密码子数据而言,没有什么评估方案(Felsenstein, 1996)。其中一个要素是碱基之间相互取代的模型;另外一个要素是序列中不同位点的所有取代的相对速率。还没有一种简单的计算机程序可以对较复杂的变量(比如,位点特异性或者系统特异性取代模型)进行评估,同样,现有的建树软件也不可能理解这些复杂变量。 碱基取代速率模型 一般而言,生物化学性质相近的碱基之间的取代频率较高;在DNA中,四种转换(A G, G A, C T, T C)的频率比八种颠换(A C, A T, C G, G T, 以及前四种的反向取代)的频率要高;这些偏向会影响两个序列之间的预计的分歧。 各个残基之间的相对取代速率一般由方阵形式列出;对于碱基而言,行数和列数都是4,对于氨基酸而言,行数和列数都是20(比如PAM方阵),对于密码子而言,行数和列数都是61(除去了中止密码子)。非对角线元素对应于一个碱基变为另一个碱基的相对代价,而对角线元素则代表不同序列拥有同一个碱基的代价。 这些代价值可以固定为先验的代价表,以确保建树方法在计分时对每一种取代都使用确定的代价值。固定的代价方阵是典型的静态权重方阵,MP建树方法使用的就是这种方阵。如果使用这种权重,那么这个方法就会被称为“加权节约”。又如,ML建树方法,代价值是由即时的速率方阵得到的,这个方阵代表了各种取代可能会发生的概率的ML估计值。MP权重方阵只涉及简单的算术,而应用距离和ML速率方阵则可以引入复杂的代数。为了避免盲目使用不适当的方法,建议大家熟悉其内部的基本原理(见Li, 1997, and / or Swofford et al., 1996a)。 实际上,“前进”和“反向”取代速率被认为是相同的;这个取代模型被称为是“时间可逆”;这个模型拥有“静态”的性质,因为在所有的碱基频率中没有预知的变化。在系统发育的特殊历史中,不同序列中的碱基频率不同表明,前进和反向速率实际上可能会不同;而传统的取代权重或者速率方阵不能包容这个“非静态”环境;本节的结尾将讨论一个基于非静态取代模型(“log�det”)的建树方法,这个方法将会提供一个可供选择的计算方法。 通常,特征符状态的权重方阵都会或多或少地通过观察进行过估值,当然也可以从速率矩阵衍生得到。比如,如果假定两个转化的其中一个,发生的频率是每个颠换的两倍,那么据此就可以确定一个权重方阵,比如,A�G的转换代价为1,而A�T的颠换代价为2。(节约方法规定对角线元素值,或者说是不同序列中拥有相同碱基的代价值为零。这是节约方法的一个缺点��详见下文)在随后的建树步骤中,这套假定会把颠换的总数降至最低值,而力求把那些主要差异是转换的序列集中在一起。 任何一种“时间可逆”的核苷酸取代模型都可以用图9.5所示的方阵刻画,只是其中一个速率和其它速率的差异;在任意组合中,最多可以达到只有六个参数,其中每一个速率参数都是独立的(Swofford et al., 1996a;Li, 1997)。如果平衡的碱基频率不相等,则需要额外的参数;如果平衡的碱基频率不相等,但是却假定这些频率相等,那么系统发育进化树的最终结果将会出错(Li, 1997)。 侧线(paralinear)(Lake, 1994)和“log�det”(Lockhart et al., 1994)做了一些修正(见Swofford et al., 1996a)来满足非静态环境的需要;这个方法只适用于距离进化树的建立;在这个方法中,对于每一个序列匹配,各种类型和变化方向的原始取代的数目都会计算在一个4×4的方阵中(如图9.6)。每个方阵都会有一个代数行列式,这个行列式的log值是评估序列差异性的一个要素,因此被称为“log�det”。对那些拥有各种各样的碱基频率的序列进行双重比较,就会得到各种各样的方阵,也就会得到各种各样的行列式值;因此,在评估序列两两之间的距离的时候,就要受到序列两两之间的行列式值的影响,而且序列两两之间的比较允许适用不同的取代模型,因此沿着系统发育进化树的不同树枝,将会产生多元化。Log�det尤其对位点之间的速率差异(见下文)敏感,因此,碱基频率的偏向可能只存在于那些承受变化的位点。 位点内速率差异模型 除了取代模型的多元化以外,一个序列中各个不同位点之间取代速率的差异也会对建立进化树的结果产生很深远的影响(Swofford et al., 1996a);关于位点之间的速率差异(或者叫做位点异质性),有一个最明显的例子,就是在一个编码序列中,三联体编码的位点差异:在三联体编码中,第三个编码位点比另外两个位点更加容易发生变化;正是出于这个原因,许多系统发育分析方法在分析编码序列时,都会把第三个编码位点排除在外;但是在某些情况下,速率差异模型会更加敏锐(比如,对应于蛋白质或者rRNA的保守序列)。 对位点差异的取代速率进行估值的方法有非参数化模型(W.M. Yang et al., 1996),不变式模型和gamma分布模型(Swofford et al., 1996a)。非参数化方法源于特异位点的相对速率的范畴;这个方法可以在MP建树方法中使用,只要根据相对的变异频率对特异位点进行简单加权就可以了,当然进行加权时需要有关于真实进化树的预备知识;这个方法同样也可以适用于ML建树方法,但是在计算上被认为是不切实际的(W.M. Yang et al., 1996)。不变式模型对一定比例的位点进行估值,这些位点并不能自由变化;剩余的位点假定为等概率变化。至于gamma模型方法,它假定一个给定的序列变化的概率是遵守gamma分布规律的,据此指定位点的取代概率;gamma分布的形状(有形状参数α描述)描述了一个序列中各个位点的取代频率的分布(Swofford et al., 1996a, p. 444, Figure 13; cf. Li, 1997, p. 76, Figure 3.10; 注意尺度差异)。在一个混合方法中,可以假定一部分位点是不变的,而剩余的位点则是按照gamma分布变化的。 实际上,gamma修正可以是连续的,离散的或者自离散的(W.M. Yang et al., 1996)。连续gamma的意思是各个位点沿一条连续的概率曲线变化;目前,这个方法在绝大多数情况下无法计算。离散gamma逼近方法指定各个位点的概率,使得这些(大量的)概率值逼近gamma曲线。自离散模型假定相邻的位点的变化速率是相关联的;许多组位点被分为许多类,其中每一类中的位点的变化速率可能被假定为常量或者异类值。 进化树建立程序使用各种各样的位点速率差异修正方法。对于核酸数据,PAUP 4.0在单独或者混合使用时间可逆的距离建树方法和最大似然建树方法时,既使用不变式gamma模型,也使用离散的gamma模型;在使用log�det距离建树方法时,使用不变式模型(见下)。对于核酸,氨基酸和编码子数据,PAML使用连续的,离散的和自离散的gamma模型。对于核酸和氨基酸数据,PHYLIP使用一种离散的gamma模型。
现在已经有一些程序可以用来评估数据中的系统发育信号和进化树的健壮性(Swofford et al., 1997)。对于前者,最流行的方法是用数据信号和随机数据作对比实验(偏斜和排列实验);对于后者,可以对观察到的数据重新取样,进行进化树的支持实验(非参数自引导和对折方法)。似然比例实验可以对取代模型和进化树都进行评估。 随机进化树(偏斜实验) 模拟研究表明,通过随机的数据集所产生的随机的MP进化树的长度的分布是对称的,但是使用系统发育信号的数据集,其分布将是不对称的(图9.9; Hillis and Huelsenbeck, 1992)。在偏斜实验中,g1统计的临界值随着分类群数目的不同和序列中位点的不同而不同。这个实验并不评估一个特定的拓扑结构的可靠性,而且这个实验对其它的随机数据集中所呈现的信号都敏感,哪怕只是很少的一点。如果数据很明显地支持某些分组,而这些分组中的分类群被有选择地删除,那么这个实验可以用来决定系统发育信号是否还保留着,当然至少要为测试提供10种不同的特征符和5个分类群。PAUP中包含了这个程序。 随机的特征符数据(排列实验) 随机数据方法决定了一个从真实数据得到的MP进化树或者其中的一部分是否可以偶然得到。实际上,数据并非真正地随机化了,只是在每一个比对列中以不同次序排列,使得初始数据的共变性被消除了;结果产生了一个非随机序列的序列比对;正确地说,这些序列中的每一个位点都是从那些在整个比对中占据这个位点的碱基群体中随机得到的。排列结尾几率实验(PTP, the permutation tail probability test)对MP进化树的分值和那些通过对每一个位点都进行大量的排列组合多得到的数据所推算出的进化树的分值进行比较,从而决定在原始数据中是否存在着系统发育信号。一个依赖于拓扑结构的实验(T�PTP, topology-dependent test)对特殊的进化树的分值进行比较,从而决定这些差异是否可以产生偶然性;这个方法并不评估这个进化树或者其中的一部分是否正确(Faith and Trueman, 1996; Swofford et al., 1996b)。值得注意的是,T�PTP实验看来似乎是进一步地确认了进化树中那些同MP进化树很接近但是并不在其中的一些分组;这是因为这个方法探测集体的信号,这些信号可以把一个分类群放置在正确(如果不能说是准确,至少也是近似)的位置;结果可以通过附加的程序使用相关的数据子集进行调整(Faith and Trueman, 1996)。PAUP中包含了这个程序。 自引导方法 自引导方法是对进化树重新取样的评估方法,可以对距离建树方法、节约建树方法、似然进化方法以及衍生出的其它任何方法进行评估。这个方法是在1979年(Efron, 1979)提出的,并且由Joe Felsenstein将其引入(Felsenstein, 1985),作为系统发育分析中的进化树评估方法。典型的自引导分析结果是一个数字,这个数字同一个系统发育进化树的一个特定树枝相关,而这个系统发育进化树则给出了支持单源进化分支的自引导的重复比例。 那么在实际操作中应该怎么做呢?自引导方法的操作过程可以分为两个步骤,第一步先从原始数据集中产生(许多)新的数据集,然后经过计算得到一个数值,表征一个特定的数值(比方说,一个分类群)在进化树中出现的次数的比例;这个数值通常被称为自引导数值。从原始数据集中产生新的数据集的具体做法是重新取样,即从原始数据集中随机地“可以替换”地抽取各个列中的特征符作为新的样本。“可以替换”的意思是说每一个位点都可以重新取样,其抽取几率同其它任何位点的抽取几率都一样;结果是每一个新建的数据集同原始数据集的位点总数相同,但是某些位点重复了两次或者三次,而某些位点则丢失了;当然新建的数据集也有可能同原始数据集完全相同��或者走向另外一个极端,只有一个位点被反复抽取,总数达到500次,而原始数据集中其它499个位点都被丢失了。 虽然自引导评估方法已经成为整个系统发育分析中的一个普遍手段,但是对于这个方法究竟计算了什么仍然有一些争论。刚开始的时候,有人提议说自引导数值计算了重合性(Felsenstein, 1985)。在最近的阐述中,自引导程序被认为是计算了精确性��这个生物学相关的参数给出了得到真实的系统发生史的可能性(Felsenstein and Kishino, 1993)。模拟研究表明,在合适的条件下(各种替换速率基本相等,树枝基本对称),如果自引导数值大于70%,那么所得到的系统发育进化树能够反映真实的系统发生史的可能性要大于95%(Hillis and Bull, 1993)。同理,如果条件不是很合适,那么如果自引导数值大于50%,则精确性的评估就会过高(Hillis and Bull, 1993)。在某些条件下,如果自引导数值较高,可能会使系统发生史看起来很好,从而得出错误的结论。 带参数的自引导方法 带参数的自引导方法同不带参数的自引导方法不同,前者使用了模拟的但是仍然真实的复制品,而不是虚假的复制品。在进行系统发育的序列分析中,与原始数据集大小相同的复制数据集是通过一个特殊的序列进化模型得到的,数据集还包括根据这个模型得到的最佳进化树拓扑结构(Huelsenbeck et al., 1996a),然后就可以使用感兴趣的方法对每一个数据集进行分析。对实验进化树的树枝是否支持的判定方法同无参数的自引导方法大体相同。 作为一种还没有被其它方法(诸如进化树中所显示的分类群中的任何分组的单个谱系)(Huelsenbeck et al., 1996a, 1996b)检验过的测试假定,带参数的自引导方法还不能算是无参数的自引导方法之外的一种选择。在每一个复制品的分析中,“真实的”进化树(假定能够产生模拟数据的进化树)的分值可能会比每一个复制品的最好的进化树的分值都要大(或者小)。分值差异图则显示了样本偏差的一个真实的正规分布。任何期望的可选的进化树拓扑的分值差异的重要价值可以由这个正规图来决定。带参数的自引导方法可以同任何建树方法协同使用。目前,这个方法的局限因素在于程序生成模拟数据的可行性。有一个程序,能够在以下情况下模拟序列数据:模型包括两个取代类型(转化和颠换),碱基频率不相等,对于内部位点的速率差异设定或者没有设定gamma修正。这个程序可以在作者的Web站点上找到,这个站点由Berkeley维护(参见本章最后所附的Internet资源列表)。 似然比例实验 正如方法名称所暗示的那样,似然比例实验适用于ML分析。评估一个次优化的似然值对于最优化模型中的正规的误差分布极为重要。在理想情况下,误差曲线被假定为一个chi-平方分布,因此实验统计值应该是最优化数值和实验数值之差的两倍,而其自由度则是不同的参数的数目。 应用chi�平方实验来选择系统发育进化树存在不少问题,尤其是因为“参数空间的不规则性”(Z.Yang et al., 1995),但是如果取代模型之间的参数数目已知的话,这个方法可以用来评估取代模型的最优性。一旦我们用上述最大似然程序评估了一个取代模型和进化树,我们就可以用较少的参数对这个进化树进行评估(比方说,把位点内的速率设置为相同:图9.7)。 Kishino�Hasegawa实验 如果给定了伴随着决定进化树误差的不确定性,另外一个方法就可以决定与比对中每一个位点相联系的取样误差(Kishino and Hasegawa, 1989)。这个程序在PAUP中执行,可以用来测试一个特定的次优化的ML或者MP拓扑结构同最优化的拓扑结构相比是否显著不同,当然必须假定用来产生最优化的进化树的模型是正确的。这个方法不能被用来评估任意选择的拓扑结构:因为不同的拓扑结构可能会拥有不同的似然功能,在某个模型下,一个统计学意义较差的进化树在另一种模型下,其统计学意义可能会变得很好。可以把本方法同带参数的自引导方法结合起来(模型和进化树已经预先最优化了),以避免这个问题(见Sullivan et al., in press)。 约束进化树搜寻 评估进化树的一个最有效的方法是比较无约束搜索和有约束的搜索,约束条件是必须搜索同一个特殊拓扑结构相联系的最优化进化树。除了比较简单分值外,还可以把约束进化树同排列实验、似然比例实验、Kishino�Hasegawa实验以及带参数的自引导评估方法结合起来。