第6章《系统发生分析》.docx

第6章 系统发生分析系统发生（或种系发生、系统发育，phylogeny）是指生物形成或进化的历史，系统发生学(phylogenetics)研究物种之间的进化关系，其基本思想是比较物种的特征，并认为特征相似的物种在遗传学上接近。系统发生研究的结果往往以系统发生树（phylogenetic tree）表示，用它描述物种之间的进化关系。通过对生物学数据的建模提取特征，进而比较这些特征，研究生物形成或进化的历史。在分子水平上进行系统发生分析具有许多优势，所得到的结果更加科学、可靠。分子系统发生分析主要分成三个步骤：分子序列或特征数据的分析、系统发生树的构造以及结果的检验。6.1 分子系统发生与系统发生树6.1.1 分子系统发生分析系统发生学是进化生物学的一个重要研究领域，系统发生分析早在达尔文时代就已经开始。从那时起，科学家们就开始寻找物种的源头，分析物种之间的进化关系，给各个物种分门别类。经典系统发生学研究所涉及的特征主要是生物表型(phenotype)特征，所谓的表型特征主要指形态学的（结构的）特征，如生物体的大小、颜色、触角个数，也包括某些生理的、生化的以及行为习性的特征。通过表型比较来推断生物体的基因型(genotype)，研究物种之间的进化关系。但是，利用表型特征是有局限性的。有时候关系很远的物种也能进化出相似的表型，这是由称为趋同进化(convergent evolution)的过程造成的。例如，如果一个生物学家按照生物体是否有眼睛来构建进化树，那么他可能将人类、两翼昆虫和软体动物放在同一个进化组中，因为它们都有光探测器官。在这个例子中，很明显这三种生物体并不具有密切的关系，在其它特征上有天壤之别，就是它们眼睛的构造也大相径庭。这说明表型有时候会误导我们，表型相似并不总是反映基因相似。用表型来判定进化关系的另一个问题是，对于许多生物体很难检测到可用来进行比较的表型特征。例如，即使用显微镜检查，细菌也几乎没有明显的特性。当我们试图比较关系较远的生物体的时候，第三个问题又出现了，即什么样的表型特征能用来比较呢？例如分析细菌、蠕虫和哺乳动物，它们之间的共同特征实在是少之又少。随着人们对生物的认识从宏观发展到微观，科学家对物种分类的依据也从宏观上的形态发展到了微观上的分子，并且有了突破性的进展，系统发生分析进入分子层次。科学家认为现今世界上存在的核酸和蛋白质分子都是从共同的祖先经过不断的进化而形成的，作为生物遗传物质的核酸和作为生命机器的蛋白质分子中存在着关于生物进化的信息，可用于系统发生关系的研究。在分子水平上进行分析具有许多表型分析所没有的优势，所得到的结果更加科学、可靠。分子系统发生分析直接从核酸序列或蛋白质分子提取的信息，作为物种的特征，通过比较生物分子序列，分析序列之间的关系，构造系统发生树，进而阐明各个物种的进化关系。当然，这些分子不仅在序列上保留进化的痕迹，它们的结构也保留着进化的痕迹。在分子水平上研究生物之间的关系早在20世纪初就开始了。科学家发现，当把一个生物体的血液注射到另一个测试生物体中时，产生免疫反应的程度与这两个生物体的进化关系直接相关，由此最早正确地推断出人和猿比它们和其它灵长类动物有更近的共同祖先。直到20世纪中期，分子数据才开始被广泛应用于系统发生研究。蛋白质电泳使得我们可以在一些浅层特征上，如分子大小和电荷，来分离和比较相关的蛋白质。20世纪60年代，蛋白质测序成为可能；20世纪70年代，研究者开始能够获得基因组信息，特别是DNA序列。蛋白质序列和DNA序列为分子系统发生分析提供了可靠的数据。在现代分子进化研究中，根据现有生物基因或物种多样性来重建生物的进化史是一个非常重要的问题。一个可靠的系统发生的推断，将揭示出有关生物进化过程的顺序，有助于我们了解生物进化的历史和进化机制。根据核酸和蛋白质的序列信息，可以推断物种之间的系统发生关系。其基本原理非常简单，从一条序列转变为另一条序列所需要的变换越多，那么，这两条序列的相关性就越小，从共同祖先分歧的时间就越早，进化距离就越大；相反，两个序列越相似，那么它们之间的进化距离就可能越小。为了便于分析，一般假设序列变化的速率相对恒定。关于地球上现代人起源的研究是一个典型而有趣的例子，科学家分析了取自世界不同地区许多人的线粒体DNA，分析结果表明，所有现代人都是一个非洲女性的后代。线粒体DNA非常适合于系统发生分析，因为线粒体DNA从母体完全传到子代，不与父代DNA重组。由于DNA 分子非常稳定，所以我们既可以通过DNA 分析活着的生物，也可以分析死去的生物，甚至分析已经绝种的生物。当然，用细胞核基因来研究系统发生关系时，遇到的一个严重的问题是，基因常常会被复制，导致在个体基因组中，一个基因可能有若干个拷贝。在进化过程中，这些拷贝各自演变，形成两个或更多的相似基因。在对不同物种的基因进行比较时，如果选择这类基因，其分析结果的可靠性将存在问题。所有的生物都可以追溯到共同的祖先，生物的产生和分化就象树一样地生长、分叉，以树的形式来表示生物之间的进化关系是非常自然的事。可以用树中的各个分枝点代表一类生物起源的相对时间，两个分枝点靠得越近，则对应的两群生物进化关系越密切。系统发生分析一般是建立在分子钟(molecular clock)基础上的。生物随着时间的推进而演化，进化的速率被视为进化研究中的基本问题之一。进化速率就是在某一段时间内的遗传改变量。分子进化速率相关的分子钟的概念源于对蛋白质序列的研究。在长期的进化过程中，有着相似的功能约束的位点的分子进化速率则几乎完全一致。20世纪60年代最早由Emile Zuckerkandl 和Linus Pauling所做的蛋白质序列的比较研究表明，蛋白质同系物的替换率就算过了千百万年也能保持恒定，因此他们将氨基酸的变异积累比做分子钟。科学家们在比较几种动物的血红蛋白、细胞色素C的序列后注意到：这些蛋白质的氨基酸取代速率在不同的种系间大致相同，即分子水平的进化存在恒速现象。分子时钟在不同的蛋白质中运行的速率是不同的，但是两个蛋白质同系物的差异始终和它们独立分化的时间成正比。两序列间稳定的变异速率，不仅有助于确定物种间系统发生关系，而且能够像利用放射性衰变考察地质年代那样，准确测定序列分化发展的时间。不同物种间的蛋白质氨基酸序列差异随着分歧时间的加大而增加，DNA 序列也存在这种规律。Kimura进一步提出了具体的分子进化观点：对于各物种的每个蛋白质，如果用每个位点每年发生的氨基酸替换次数作为衡量分子进化的速率，则该速率是大致恒定的；功能上次要的分子（或者分子部分）的进化速率比功能重要的分子（或者分子部分）进化速率快；对现有分子结构或者功能破坏小的氨基酸替换比破坏力大的氨基酸替换发生得更加频繁。尽管以上的分析给我们带来了很多希望，但是，Zuckerkandl 和Pauling的分子时钟假说还是有争议的。经典进化学家们认为形态的进化不够稳定，这与分子以稳定的速度变异不一致。关于分化时间也有不同意见，这些意见对这个假说的核心即进化率是稳定的表示质疑。6.1.2 系统发生树一般来说，系统发生树是一种二叉树。所谓树，实际上是一个无向非循环图。系统发生树由一系列节点（nodes）和分支（branches ）组成，其中每个节点代表一个分类单元（物种或序列），而节点之间的连线代表物种之间的进化关系。树的节点又分为外部节点（terminal node）和内部节点（internal node）。在一般情况下，外部节点代表实际观察到的分类单元，而内部节点又称为分支点，它代表了进化事件发生的位置，或代表分类单元进化历程中的祖先。分类单元是一种由研究者选定的基本单位，在同一项研究中，分类单元一般应当一致。在下面的讨论中，我们基本上以序列（DNA序列或蛋白质蛋白质）作为分类单位。树节点间的连线称为分支，其中一端与叶节点相连的为外支，不与叶节点相连的为内支。系统发生树有许多形式：可能是有根树（rooted tree），也可能是无根树（unrooted tree）；可能是一般的树，也可能是二叉树；可能是有权值的树（或标度树，scaled tree，树中标明分支的长度），也可能是无权值树(或非标度树，unscaled tree)。在一棵有根树中，有一个唯一的根节点，代表所有其它节点的共同祖先，这样的树能够反映进化层次，从根节点历经进化到任何其它节点只有唯一的路径。系统发生分析中一个重要的差别是，有的能由系统发生树推断出共同祖先和进化方向，而有的却不能。无根树没有层次结构，无根树只说明了节点之间的关系，没有关于进化发生方向的信息。但是，通过使用外部参考物种（那些明确地最早从被研究物种中分化出来的物种），可以在无根树中指派根节点。例如，在研究人类和大猩猩时，可用狒狒作为外部参考物种，树的根节点可以放在连接狒狒与人和大猩猩共同祖先的分支上。二叉树是一种特殊的树，每个节点最多有两个子节点。在有权值的树中，分支的长度（或权值）一般与分类单元之间的变化成正比，它是关于生物进化时间或者遗传距离的一种度量形式。一般假设存在一个分子钟，进化的速率恒定。系统发生树具有以下性质：（1） 如果是一棵有根树，则树根代表在进化历史上是最早的、并且与其它所有分类单元都有联系的分类单元；（2） 如果找不到可以作为树根的单元，则系统发生树是无根树；（3） 从根节点出发，到任何一个节点的路径均指明进化时间或者进化距离。图6.1(a)所示的是一棵有根树，而图6.1(b)显示的是一棵无根树，图中的A、B、C、D为所研究的分类单元。 树根 A C 1 1 2 1 1 2 1 1 1 3 B 3 A B C D D (a) (b)图6.1 系统发生树。 (a)有根树；(b)无根树。对于给定的分类单元数，有很多棵可能的系统发生树，但是只有一棵树是正确的，分析的目标就是要寻找这棵正确的树。基于单个同源基因差异构建的系统发生树称为基因树（gene tree），这比称作物种树（species tree）更为合理。因为这种树代表的仅仅是单个基因的进化历史，而不是它所在物种的进化历史。物种树一般最好是从多个基因数据的分析中得到。基因树和物种树之间的差异是很重要的。例如，假设只用HLA的等位基因来构建物种树，许多人将与大猩猩分在一起，而不是和其他人分在一起。6.1.3 距离和特征用于构建系统发生树的分子数据分成两类：（1）距离（distances）数据，常用距离矩阵描述，表示两个数据集之间所有两两差异；（2）特征(characters)数据，表示分子所具有的特征。分子系统发生分析的目的是探讨物种之间的进化关系，其分析的对象往往是一组同源的序列。这些序列取自于不同生物基因组的共同位点。序列比对是进行同源分析的一种基本手段，是进行系统发生分析的基础，一般采用基于两两比对渐进的多重序列比对方法，如ClustalW程序。通过序列的比对，可以分析序列之间的差异，计算序列之间的距离。无论是DNA序列，还是蛋白质序列，都是由特定字母表中的字符组成的。计算序列之间距离的一个前提条件是要有一个字符替换模型，替换模型影响序列多重对比排列的结果，影响系统发生树的构造结果。在具体的分析过程中，需要选择一个合理的字符替换模型，参见第3章的各种打分模型或代价、距离模型。距离（或者相似度）是反映序列之间关系的一种度量，是建立系统发生树时所常用的一类数据。在计算距离之前，首先进行序列比对，然后累加每个比对位置的得分。可以应用第3章介绍的关于序列比较方法，直接计算序列之间的距离。如果在进行序列比较时使用的是相似性度量函数，则需要将相似度（或者得分）转换成距离。令S(i,j)是序列i和序列j比对位置得分的加权和，一种归一化的距离计算公式为：（6-1）其中，Sr(i,j)是序列i和j随机化之后的比对得分的加权和，Smax(i,j)是可能的最大值（当两条序列相同时，取最大值）。两个序列归一化距离的值处于0和1之间，当两个序列完全一致时，距离为0；当两个序列差异很大时，距离接近于1。如果在上式中令Sr(i,j)=0，则计算公式变为：（6-2）为了适合于处理相似性较小的序列，可以进一步修改距离计算公式：（6-3）序列比对得分的加权和可以根据常用的打分矩阵获得，如果待处理的序列是蛋白质，则用PAM矩阵、BLOSUM矩阵等；如果待处理的序列是DNA或者RNA，则用等价矩阵、核苷酸转换-颠换矩阵或者其它具有非对称置换频率的矩阵。距离是系统发生分析时所使用的一类数据，另一类数据就是所谓的离散特征数据。离散特征数据可分为二态特征与多态特征。二态的离散特征只有2种可能的状况，即具有与不具有某种特征，通常用“0”或“1”表示。例如，DNA序列上的某个位置如果是剪切位点，其特征值为1，否则为0。多态离散特征具有两种以上可能的状态，如核酸的序列信息，对序列中某一位置来说，其可能的核苷酸碱基有、共4种。可以将特征数据转换为距离数据。如果建立所有可能状态之间相似性的度量，特征数据就很容易被转换成距离数据。6.1.4 分子系统发生分析过程分子系统发生分析主要分成三个步骤：（1）分子序列或特征数据的分析；（2）系统发生树的构造；（3）结果的检验。其中，第一步的作用是通过分析，产生距离或特征数据，为建立系统发生数提供依据。系统发生树的构建方法很多种。根据所处理数据的类型，可以将系统发生树的构建方法大体上分为两大类。一类是基于距离的构建方法，利用所有物种或分类单元间的进化距离，依据一定的原则及算法构建系统发生树。基本思路是列出所有可能的序列对，计算序列之间的遗传距离，选出相似程度比较大或非常相关的序列对，利用遗传距离预测进化关系。这类方法有非加权组平均法（unweighted pair group method with arithmetic means）、邻近归并法（neighbor joining method）、Fitch-Margoliash法、最小进化方法（minimum evolution）等。另一类方法是基于离散特征的构建方法，利用的是具有离散特征状态的数据，如DNA序列中的特定位点的核苷酸。建树时，着重分析分类单位或序列间每个特征（如核苷酸位点）的进化关系等。属于这一类的方法有最大简约法（maximum parsimony method）、最大似然法（maximum likelihood method）、进化简约法（evolutionary parsimony method）、相容性方法（compatibility）等。对相似性和距离数据，在重建系统发生树时只能利用距离法。离散特征数据通过适当的方法可转换成距离数据，因此，对于这类数据在重建系统发生树时，既可以用距离法，亦可以采用离散特征法。根据建树算法在执行过程中采用的搜索方式，系统发生树的构建方法也可以分为三类。第一类是穷尽搜索方法，即产生所有可能的树，然后根据评价标准选择一棵最优的树。需要注意的是，系统发生树可能的个数随序列的个数急剧增加。假设要为n个分类单元建立系统发生树，则可能的有根树个数（NR）和无根系统发生树个数（NU）可用下面的算式计算得到：（6-4）（6-5）可以看到，随着n的增加，可能的有根系统发生树和无根系统发生树的数目迅速增加。表6.1中列出了一些n值，以及对应的有根树和无根树的数目。当n大于等于15时，可能的系统发生树数目变得很惊人，但是只有其中的一棵树代表了待分析基因或者生物体之间的真实进化关系，我们的目的就是找出这棵反映真实进化关系的树。表6.1 对不同的n，可能的有根树和无根树数目数据数目有根树数目无根树数目21133141535105151034,459,4252,207,02515213,458,046,767,8757,905,853,580,625208,200,794,532,637,891,559,375221,643,095,476,699,771,875从计算量来看，这种方法只能处理很少的分类单元。当分类单元个数n大于一定值（如15），几乎不可能采用穷尽搜索的方式来求取最优树。第二类方法是分支约束方法，即根据一定的约束条件将搜索空间限制在一定范围内，产生可能的树，然后择优。这是人工智能技术中的一种空间搜索策略，这种搜索方式不需要搜索整个树空间，可大大提高搜索效率。第三类是启发式或经验性方法，即根据先验知识或一定的指导性规则压缩搜索空间，提高计算速度。这种方法能够处理大量的分类单元，虽然不能保证所构建的树是最优的，但实际结果往往接近于最优解。当待分析的对象个数比较多时，必须采用分支约束方法或者启发式的方法。在构造系统发生树时需要考虑进化假设和进化模型。系统发生树的类型可能是有根树，这意味着其中的一个序列代表其他所有序列共同的祖先。另一方面，系统发生树可能是无根树，意味着没有共同的祖先。一般认为序列是随机进化的，序列中的所有位点的进化也是随机的而且是独立的。在进行具体的系统发生分析时，一般还要作一些假设：序列必须是正确无误的，待分析的序列是同源的，所有的序列都起源于同一个祖先序列，并且它们不是共生同源（或平行进化）序列，在序列比对中，不同序列的同一个位点都是同源的。另外，当两个物种在系统发生树上分化后，各自独立进化发展。对系统分析的样本也有一定的要求，要求样本足以反映感兴趣的问题，样本序列之间的差异包含了足以解决感兴趣的问题的系统发生信号。通过某种算法构造好一棵系统发生树之后，需要对树的合理性和可靠性进行分析。对于若干个序列，如果利用多种不同的分析方法进行系统发生分析，并且得到相似的进化关系，那么分析结果具有较高的可信度。6.2 基于距离的系统发生树构建方法建立系统发生树的基本任务是：在给定的条件下（包括分类单元、分类单元的特征值或者序列），构造一棵最优的系统发生树。这里重点讨论针对DNA序列或者蛋白质序列构建系统发生树。基于距离的系统发生树构建方法的基本思路是：给定一种序列之间距离的测度，在该距离测度下构建一棵系统发生树，使得该树能够最好地反映已知序列之间的距离。这种方法采用两两距离，建立一个距离矩阵，如表6.2所示，根据距离矩阵构造系统发生树。表6.2 10条核酸序列的距离矩阵。12345678920.051630.05500.003140.04830.02210.025350.05820.06510.06850.054960.00940.04160.04500.03840.054970.01250.05840.06190.05510.06510.015780.02840.06870.07220.06540.07540.03170.028590.09250.12210.12590.11850.13700.08200.07860.0927100.19210.21830.22280.20540.23090.17980.17950.18330.1860这里的距离代表两条序列之间各位点核苷酸替换数的估计值。6.2.1 最小二乘法在很多情况下，这种距离矩阵传达了大部分进化信息。这里的距离可以是序列比对的距离（代价），也可以是序列相似性得分。运用前面介绍的序列比较程序，可以计算出序列之间的距离。但是，在进行序列比较时，应根据比较的对象是DNA序列还是蛋白质序列，正确选用不同的打分矩阵。为了便于分析，首先定义一种连续加和距离函数，在该函数下，两个分类单元之间的距离与系统发生树中连接这两个分类单元的分支总长度成正比。这样，如果分类单元a和分类单元b由经过中间节点v的两条边相连，两条边的长度分别为dav和dbv，则它们之间的距离为dav+ dbv。这样，可以在系统发生树中确定a和b的相对位置。进一步，假设三个分类单元之间的距离分别为dab、dac、dbc,如果分类单元a和分类单元b由经过中间节点v的两条边相连，再经过节点u与分类单元c相连,则可以通过求解线性方程计算出系统发生树的各种内部距离。例，如果有三个分类单元，其两两距离如下：dab = 0.5;dac = 0.9; dbc = 0.9假设分类单元a和分类单元b的分歧起始时间是相同的，根据分子时钟假说，和的值应该是相等的，进一步假设节点u到其它节点的距离相同，则通过求解方程，得到如图6.2所示的一棵树。 u0.450.2 v0.250.25 a b c 图6.2 根据方程求解结果构造系统发生树 图6.2是一个简单的例子，但是，在实际工作中，所要处理的分类单元可能很多，因而，需要求解的线性方程也很多，难以求解，或者方程组的求解过程存在着不确定性。因此，需要采用数学逼近的方法。统计学上逼近的方法之一是最小二乘法。我们的目标是构造一棵树T，该树的叶节点代表分类单元，用该树预测分类单元之间的距离。通过优化，使下式最小化：（6-6）这里，Dij为分类单元i和j的实际观察距离（或序列之间的计算距离），dij是分类单元i和j在系统发生树T 中的距离，Wij是与分类单元i和j相关的权值。SSQ（T）是树T所有预测值与实际观察值偏差的累加和。权值Wij一般为1，或（6-7）一般来说，寻找一棵最小方差树是一个NP-完全问题，需要采用近似的算法。下面讨论几种计算复杂度为多项式的启发式方法，即连锁聚类方法（linkage clustering）、非加权组平均法（UPGMA）和邻近归并法（Neighbor Joining）。6.2.2 连锁聚类方法及非加权分组平均法连锁聚类属于一般的聚类分析方法，当用来构建系统发生树时，其假定的前提条件是：在进化过程中，核苷酸或氨基酸的替换速率是均等且恒定的，在每一次分歧发生后，从共同祖节点到两个分类单元间的分支长度一样。在构建系统发生树时，首先用n个叶节点表示n个分类单元（序列）， 每个分类单元自成一类，然后将进化过程中的祖先赋予树的内部节点，通过反复的聚类使所有的分类单元都聚为一类，最终得到一个完整的系统发生树。假设若干条序列是从一个共同的祖先进化而来，则系统发生树将是一个有根树，并且从根节点出发到所有叶节点路径的长度相同。对于给定的序列，通过序列之间的两两比对，计算序列之间的进化距离，然后根据距离矩阵构造系统发生树。算法的基本思路是首先从距离矩阵中选择距离最小的一对分类单元（序列），令它们分别为x和y，然后将这两个分类单元合二为一，形成一个新的对象（代表这两个分类单元的祖先，记为z），并重新计算这个新的对象与其它分类单元(或对象，以u表示)之间的距离d(z,u)。不同的实现方案采用不同的计算公式：单连锁聚类: d(z,u)=min(d(x,u),d(y,u)（6-8）最大连锁聚类： d(z,u)=max(d(x,u),d(y,u)（6-9）平均连锁聚类： d(z,u)=1/2(d(x,u)+d(y,u)（6-10）其中z代表x和y的合并，u代表任意其它对象。每次合并所形成的新对象实际上是一个聚类,以一个内部节点表示，该节点到x、y 所在节点的距离相同，其值等于d（x，y）的一半，而到其它节点的距离按照上述公式计算。每次合并后，修改距离矩阵。重复上述过程，直到所有的分类单元都被合并到一类为止。连锁聚类算法的执行过程如下：（1） 初始化：使每个分类单元自成一类，如果有n个分类单元，则开始时共有n个类，每个类的大小为1，分别用n个叶节点代表每个类；（2） 执行下列循环：l 寻找具有最小距离Dxy 的两个类x、y；l 建立一个新的聚类z, 该聚类是类x和类y的合并；l 在树中建立一个新的内部节点z，生长两个新的分支，将x 和y 连接到z，并使这两个分支到x、y 中各个叶节点的长度为d（x，y）/2；l 按照公式（6-8）、（6-9）或（6-10）计算新的分类到其它类的距离l 在距离矩阵中删除与类x和类y 相应的行和列，为类z 加入新的行和列，其值按上述公式计算；重复循环，直到仅剩一个类为止。系统发生树的结构与分类过程相对应。图6.3就是一个通过平均连锁聚类分析的实例。距离矩阵如图6.3（a）所示。初始化以后，共有5个分类，分别对应于5个分类单元A、B、C、D及E。第一次聚类时，选择距离最小的两个分类A和 C进行合并，形成一个新的分类（AC）。然后修改距离矩阵，（AC）到B的距离等于（8+8）/2=8，其它距离类推。最终结果见图6.3 (b)。在系统发生树中两个分类单元之间的距离为连接这两个分类单元所在叶节点对应路径的长度。ABCDB8C48D686E 8488（a） 11 2 3 2 2 2 2 B E D A C 图6.3 平均连锁聚类：（a）距离 矩阵；（b）聚类结果 （b）运用连锁聚类算法，得到一棵有根的系统发生树，从树根到任何叶节点的分支长度全都一样，也就是说，所有物种的突变速率相同，存在一个固定节律的“分子钟”，各个物种从树根（代表某个历史起点）出发，踏着同样的节律，沿不同的路径，演化成为当前的形式。而任意两个物种之间的距离是连接这两个物种所在叶节点路径上的距离之和。如果上述过程处理的分类单元是序列的话，则聚类的过程实际上也是一个进行序列多重比对的过程，将序列多重比对问题转化为序列两两比对问题。非加权分组平均法（Unweighted Pair Group Method with Arithmetic mean, UPGMA）是一种较为常用的聚类分析方法，最早是用来解决分类问题的。使用该方法的前提条件也是稳定的进化速率。从分类过程来看，该方法与前面介绍的平均连锁聚类方法相似，但是关于各个分类（即若干个分类单元形成的集合）之间距离的计算公式不一样。在平均连锁聚类过程中，一个新类到其它类之间的距离就是简单的原距离平均值。这样的计算非常简单，但是如果类中分类单元个数不一样，原距离矩阵中各个距离值对新距离计算的贡献就不一样，或者说是经过“加权”的，称这样的聚类为加权分组平均（Weighted pair group method with Arithmetic mean, WPGMA）。在非加权分组平均法中，在计算新分类到其它分类之间的平均距离时按照各分类中分类单元的数目进行加权处理。令类x和类j所形成的新分类为（xy），新分类到其它类u的距离按照下述公式计算： （6-11）其中ni、nj、(ni+nj)分别为i类、j类、(ij)类的元素个数。这样，每个距离对最终结果的贡献一样，即是“非加权”的。图6.4是针对表6.2中的距离矩阵利用UPGMA 而构建系统发生树。0.1 0.5 1.0图6.4 UPGMA所构造的系统发生树167823459106.2.3 距离变换法连锁聚类和UPGMA算法的一个缺陷是假定所有家系的进化速率是相同的，但是实际情况并不总是这样。进化速率的变化对极易导致连锁聚类和UPGMA算法产生错误拓扑结构的树。假设有4个分类单元A、B、C和D ，其系统发生关系及各个分类单元之间的距离如图6.5所示，距离矩阵见表6.3。如果利用UPGMA进行分析，则首先合并A和C，（AC）到B的距离等于9/2+11/2=10，（AC）到D 的距离等于12/2+10/2=11；进一步合并（AC）和B ，（AC）B）到D 的距离等于211/3 + 114/3 = 12；最终得到图6.6所示的系统发生树，但是，该树显然与真实树有出入。 1 2 6 3 3 6 A B C D 图6.5 4个分类单元A、B、C和D之间的系统发生关系表6.3 四个分类单元的距离矩阵ABCB9C811D121410 1 1 6 5 4 4 A C B D图6.6 UPGMA根据表4.3构造的系统发生树有一些基于距离矩阵的方法考虑了不同的家系有不同的进化速率，其中最简单最早的算法是距离变换法（Transformed Distance Method）。这种方法充分利用了外群或外部参考物种(outgroup)，即先于其它所有被考虑的物种（称为内群或内部物种，ingroup）从它们的共同祖先中分化出来的那些物种。假设有4个物种A、B、C和D ，其中D是物种A、B和C的外部参考物种，并且已知关于这四个物种的距离矩阵。D可作为变换其它物种之间的距离的外部参考，变换式如下： （i，j=A，B，C）（6-12）其中dij是物种i和j之间的变换后距离，是外部参考物种与全体内部物种之间的平均距离，此例中= (dAD + dBD + dCD)/3）。引入只是为了确保所有变换后的距离为正值，因为从进化的角度看，负值是不可能的。用值和原来的距离矩阵可以为物种A、B和C产生一个新的距离变换矩阵。可以利用基本的UPGMA方法针对这个新的矩阵构造系统发生树，这样可以得到一棵具有正确拓扑结构的树。距离变换法的优势体现在那些很简单却容易被忽略的方面：内部物种只是在分化发生后进化分离出来的，所以它们积累的替换数目一定是从那以后才有了差异。此时，外部参考物种为比较它们替换速率提供了客观参考的框架。6.2.5 邻近归并法邻近归并法（Neighbor Joining）是另一种快速的聚类方法，该方法是Saitou和Nei于1987年首次提出的。在构建系统发生树时，该方法取消了非加权分组平均法所作的假定，不需要关于分子钟的假设，在进化分支上，发生趋异的次数可以不同。与非加权分组平均法相比，邻近归并法在算法上相对较复杂，它跟踪的是树上的节点而不是分类单元。这种方法的基本思想是：在进行类的合并时，不仅要求待合并的类是相近的，同时，还要求待合并的类远离其它的类。在聚类过程中，根据原始距离矩阵，根据所有节点间的平均趋异程度，对每两个节点间的距离进行调整，即将每个分类单元的趋异程度标准化，从而形成一个新的距离矩阵。重建时，将距离最小的两个叶节点连接起来，合并这两个叶节点所代表的分类，形成一个新的分类。在树中增加一个父节点，并在距离矩阵中加入新的分类，同时删除原来的两个分类。随后，新增加的父节点被看成为叶节点，重复上一次循环。在每一次循环过程中，都有两个叶节点被一个新的父节点所取代,两个类被合成为一个新类。整个循环直到只剩一个类为止。从所得到的系统发生树来看，对于两个聚在一起的分类单元，其所在的叶节点到父节点的距离并不一定相同。在每一次循环中，都要在树中寻找两个分类单元的直接祖先。对于节点x，到其它节点的距离dx 按下式进行估算： （6-13）这里Dxy是分类x和分类y之间的距离，是动态更新的距离矩阵D中的元素。 为了使所有分支长度的和最小（或称为最小进化原则），选择Dxy-dx-dy最小的一对节点x 和节点y 进行归并（如图6.7所示）。 x 树的其它节点 y 图6.7 寻找一对节点x和y，使这两个节点靠近，但同时远离其它节点。算法如下：（1） 初始化（与UPGMA一样）（2） 循环 l 对于所有的分类单元x，按公式（6-13）计算dx；l 选择一对分类单元x 和y ，使Dxy-dx-dy最小；l 将x和y归并为新的类 (xy)，在树中添加一个新的节点，将它与节点x和y连接 ，新节点代表新生成的分类，计算从x和y到新节点（xy）的分支长度； dx,(xy） = 1/2Dx,y +1/2(dx-dy), dy,(xy） = 1/2Dx,y +1/2(dy-dx) l 计算新类与其它类的距离；D(xy）,u = 1/2(Dx,u + Dy,u - Dx,y) l 删除聚类x和y，添加新类 (xy),更新距离矩阵；l 如果有两个以上的分类存在，则继续执行循环，否则用长度为Di,j的分支连接剩余的两个类 。例如，假设6个分类单元的距离矩阵如表6.4所示。首先分别用6个叶节点代表分类单元，计算每个节点到其它节点的距离：dA = 5+4+7+6+8=30/4=7.5dB = 10.5dC = 8dD = 9.5dE = 8.5dF = 11表6.4 6个分类单元的距离矩阵ABCDEB5C47D7107E6965F811898第一次归并时，选择使Dxy-dx-dy最小的一对分类单元x=A 和y=B，将这两个邻近的单元归并，以一个新节点（xy）代表，同时计算（xy）到x和y所在节点的距离：dx,(xy）= Dx,(xy）/2+(dx-dy)/2=5/2+(7.5-10.5)/2=1dx,(xy）= Dx,(xy）/2+(dy-dx)/2=5/2+(10.5-7.5)/2=4进一步计算新类与其它类的距离，更新距离矩阵，重复循环。最终结果见图6.8。图6.8 利用邻近归并算法构造的系统发生树 6.3 基于特征的系统发生树构建方法基于特征的系统发生分析要解决的一般问题是：给定n个分类单元，m个用以描述分类单元的特征，以及每个分类单元所对应的特征值，构建一棵系统发生树，使得某个目标函数最大。输入一般为一个nm的矩阵M，其中，Mij代表第i个分类单元之第j个特征的取值，如表6.5所示。在构建系统发生树时，假设特征是相互独立的，即一个特征的变化不影响另一个特征。另外，还假设在进化过程中，两个物种分叉后独立进化，互不影响。表6.5 分类单元特征矩阵分类单元位点1位点2位点3位点4位点5位点6甲CAGGTA乙CAGACA丙CGGGTA丁TGCACT戊TGCGTA对于给定的条件，要在很多可能的树中找一棵最佳的树。在实际应用中，不可能穷尽搜索所有可能的树，必须按照一定的方法、一定的策略在较短的时间内得到比较好的结果。6.3.1 最大简约法最大简约法(Maximum Parsimony)最早是基于形态特征分类的需要而发展起来的，因使用的算法不同而有着许多版本，其中有些已被广泛地用于分子进化研究中根据离散特征数据构建系统发生树，如对DNA序列数据分析。最大简约法的目标是构造一棵反映分类单元之间最小变化的系统发生树。最大简约法利用的只是对简约分析能提供信息的特征，如在DNA序列数据中，利用的只是存在于核苷酸序列差异（至少有两种不同类型的核苷酸）的位点，这些位点称为简约信息位点（parsimony informative site）。具体来说，信息位点就是指能由位点产生的突变数目把一棵树与其它树区分开来的位点。如果对于某个位点，所有序列都有同样的字符，则这个位点称为不变位点(invariant)。显然不变位点是非信息位点（uninformative site）。如果一个位点是信息位点，那么它至少有两种不同的核苷酸，并且这些核苷酸至少出现两次。所有的简约法程序在开始时都将这条简单的规则应用于输入数据集。显然，表6.5中位点6 是非信息位点，该位点将被舍弃，在简约法分析中不再被考虑。但是非信息位点对基于距离的方法中两两相似度的得分都有贡献，仅这一点差别就可能使这两类方法产生的结果有很大的不同。对于系统发生树最直观的代价计算就是沿着各个分支累加特征变化的数目，而所谓简约就是使代价最小。利用最大简约方法构建系统发生树，实际上是一个对给定分类单元所有可能的树进行比较的过程，针对某一个可能的树，首先对每个位点祖先序列的核苷酸组成做出推断，然后统计每个位点用来阐明差异的核苷酸最小替换数目。在整个树中，所有简约信息位点最小核苷酸替换数的总和称为树的长度或树的代价。通过比较所有可能的树，选择其中长度最小、代价最小的树作为最终的系统发生树，即最大简约树（maximum parsimony tree）。假设图6.9是根据表6.5所建立的一棵系统发生树。对于位点1，由于甲、乙、丙都取“C”值，而丁、戊取“T”值，它们分别处于树的两边，因此，只需要在系统发生树“根节点”所连接的两条边上有一个变化（左分支T变为C，或者右分支C变为T），就能解释各个分类单元在该位点的数据。对于位点2，甲、乙取“A”值，而丙、丁、戊取“G”值，通过“节点2”可以将它们分开，因此，“节点2” 所连接的左分支上的一个变化（G变为A），就能解释各个分类单元在该位点的差异。位点3的情况与位点1相同，位点4的情况比前几个位点复杂，甲、乙、丙、丁、戊的取值方式为“G、 A、G、A、G”，需要在“节点1”和“节点3”所连接的边上分别设置一个变换。因此位点4需要两个变化。同样位点5需要两个变化。对于位点6，显然只需要一个变化，其变化点处于“节点3”所连接的分支。位点6属于非信息位点，在简约分析时可以不再考虑它。总的来说，对于图6.9所示的整个系统发生树，共需要8个变化。 甲 乙 丙 丁 戊 节点1 节点3 节点2 根节点图6.9 对应于表6.4的系统发生树对于表6.5中所列出的5个分类单元，还可以形成其它拓扑结构的系统发生树，这些树为解释表6.5中各个位点的数据，可能需要更多或者更少的变化。而最大简约法的目的就是要寻找一棵最小变化的系统发生树。最大简约法的处理过程如下：（1） 针对待比较的物种，选择核酸或蛋白质序列。有些分子比其它分子变化慢，适合于进行进化分析，例如哺乳类的线粒体DNA、管家蛋白质等；（2） 比较各个序列，产生序列的多