生物信息学第六章系统发生分析

第一节系统发１、系统发生分析系统发生（phylgny）是指生物形成或进化的历史，系统发生学(phylgnecs)研究物种之间的进化关系，其基本思想是比较物种的特征，并认为特征相似的物种在遗传学上接近。系统发生研究的结果往往以系统发生树（phylgneicree）表示，用它描述物种之间进化关系。系统发生分析早在时代就已经开始，经典系统发生学涉及的特征主要是物理或表型特征，如生物的大核酸序列或蛋白质分子。现今世界上存在的核酸和蛋白质分子都是从共同的祖先经过大量的细微改变进化而来，这些生物分子的序列中隐藏着大量关于物种进化的信息。现代系统发生分析通过比较生物分子序列，分析序列之间的关系，构造系统发生树，进而阐明各个物种的进化关系(FtchndMrolah,16)。在现代分子进化研究中，根据现有生物或物种多样性来重建生物的进化史是一个非常重要的问题。一个可靠的系统发生的推断，将揭示出有关生物进化过程的顺序，有助于了解生物进化的历史和进化机制。根据核酸和蛋白质的序列信息，可以推断物种之间的系统。其基本原理非分析，一般假设序列变化的速率相对恒定。关于地球上现代人的研究是一个典型而有趣的例子，科学家分析了取自世界不同地区许多人的线粒体DNA。线粒体DNA非常适合于系统发生分析，因为线粒体DNA从完全传到子代，不与父代DNA重组。分析结果表明，所有现代人都是一个非洲女性的后代。由于ＤＮＡ分子非常稳定，所以可以通过ＤＮＡ分析活着的生物，也可以分析死去的生物（Hermanandumel,194，甚至分析已经绝种的生物。当然用细胞核来研究系统发生关系时，遇到的一个严重的问题是，常常会被，导致在组中，一个可能有若干个拷贝。在进化过程中，这些拷贝各自演变，形成两个或的相似。在对不同物种的进行比较时，如果选择这类，其分析结果将存在问题。所有的生物都可以追溯到共同的祖先，生物的产生和分化就象树一样地生长、分叉，以树的形式来表示生物之间的进化关系是非常自然的事。可以用树中的各个分枝点代表一类生物的相对时间，两个分枝点靠得越近，则对应的两群生物进化关系越密切。对生物学家来说，重要的是知道某些物种或种群的历史及在每一次后趋异时间，当这些历史事件以系统发生树的形式表现时，此时的系统发生树称为物种（或种群）树。对任何一类生物，要想知道确切的物种（或种群）树是非常的，但可以检测这类生物所包含的一些的进化关系来推断物种（或种群树(常青等,1998)。如果系统发生树是基于一个的核酸或氨基酸序所立，此的统发树称为个树树时物种（或种群树并不一致，特别是在基因组中存在２个或的同样的拷贝时。但是，无论是物种（或种群）树还是树，都是用树一类的拓扑结构表示出来。２、系一般来说，系统发生树是一种二叉树。所谓树是一个无向非循环图。系统发生树由一系列节点和分支组成，其中每个节点代表一个分类单元（物种或序列，而节点之间的连线代表物种之间的进化关系。树的节又分为外部节点和节点。在一般情况下，外部点代表实际观察到的分类单元。节点又称为分支点，它代表了进化事件发生的位置，或代表分类单元进化历程中的祖先。分类单元亦称为运筹分类单位（opertonlaoomcnts，简称OTU，是一种由研究者选定的基本单位，在一项研究中分类单元一般应当一致。在下面的中，大多以序列（蛋白质、DNA）作为分类单位。树节点间的连线称为分支，其中一端与叶节点相连的为外支，不与叶节点相连的为内支。系统发生树有许多形式：可能是有根树，也可能是无根树；可能是一般的树，也可能是二叉树；可能是值的树（树中标明分支的长度，也可能是无权值树。在一棵有根树中，有一个唯一的根节点，代表所有其它节点的共同祖先，这样的树能够反应进化层次，而无根树则没有层次结构。二叉树是一种特殊的树，每个节点最多有两个子节点。在值的树中，分支的长度（或权值）一般与物种之间的变化成正比，它是关于生物进化时间或者遗传距离的一种度量形式。一般假设存在一个分子钟，进化的速率恒定。系统发生树具图6.1(a)所示的是一棵有根树，而图6.1(b)显示的是一棵无根树，图中的A、B、C、D为所研究对于给定的分类单元数，有很多棵可能的系统发生树，但是只有一棵树是正确的，分析的目标就是要寻找这棵正确的树。树112112 3ABC21B1321B13 6.1a)有根树；(b)３、系统发生分析步系统发生分析的目的是探讨物种之间的进化关系，其分析的对象往往是一组同源的序列。这些序列取自于不同生物组的共同位点。序列比对是进行同源分析的一种基本，是进行系统发生分析的基础，一般采用基于两两比对渐进的多重序列比对方法（FengandDoolittle，1996；Hein，，如无论是DNA序列，还是蛋白质序列，都是由特定字母表中的字符组成的，计算序列之间距离的一个前提条件是要有一个字符替换模型。替换模型影响序列多重对比排列的结果，影响系统发生树的构造结果。在具体的分析过程中，需要选择一个合理的字符替换模型，参见第三章的各种打分模型或代价、距离模型。（或者相似度是反映序列之间关系的一种度量，是建立系统发生树时所常用的一类数据。在计算距离之前，首先进行序列比对，然后累加每个比对位置的得分。可以应用第三章介绍的关于序列比较方法，直接计算序列之间的距离。如果在进行序列比较时使用的是相似性度量函数，则需要将相似度（或者得分）转换成距离。令S(i,j)是序列i和序列j比对位置得分的和，一种归一化的距离计算公式为：d(i,j)1S(i,j)Sr(i,Smax(i,j)Sr(i,

（6-其中Sr(i,j)是序列i和j随机化之后的比对得分的和，Smax(i,j)是可能的最大值（当两条序列相同时取最大值。两个序列归一化距离的值处于0和1之间，当两个序列完全一致时，距离为0，当两个序列差异很大时，距离接近于1。如果在上式中令Sr(i,j)=0，则计算公式变为：d(i,j)1 S(i,

（6-为了适合于处理相似性较小的序列，可以进一步修改距离计算公式：d(i,j) S(i,Smax(i,

（6-序列比对得分和可以根据常用的打分矩阵获得。如果待处理的序列是蛋白质，则用PAM矩阵、BLOSUM矩阵等，如果待处理的序列是DNA或者RNA，则用等价矩阵、核苷酸转换矩阵或者其它具有距离是系统发生分析时所使用的一类数据，另一类数据是所谓的离散特征数据。离散特征数据可分为二态特征与多态特征。二态的离散特征只有2种可能的状况，即具有与不具有某种特征，通常用“0”或“1”表示。例如，DNA序列上的某个位置如果是剪切位点，其特征值为1，否则为0。多态离散特征具有两种以上可能的状态，如核酸的序列信息，对序列中某一位置来说，其可能的核苷酸碱基有Ａ、Ｔ、Ｇ或Ｃ4种。系统发生树的构建方法很多种（Pennyetal192.;Flestin196；Aise994Swfordeta.，196；i197；Bociei2001。根据所处理数据的类型，可以将系统发生树的构建方法分为两大类。一类是基于距离的构建方法，利用所有物种或运筹分类单位间的进化距离，依据一定的原则及算法构建系统发生树。基本思路是列出所有可能的序列对，计算序列之间的遗传距离，选出相似程度比较大或非常相关的序列对，利用遗传距离进化关系。这类方法有非组平均法（unwigtdarrupmthdwiharthmeticmas、邻近归并法（neihorjinngmehd、Fitrglash法、最小进化方法（minmmvouion）等。另一类方法是基于离散特征的构建方法，利用的是具有离散特征状态的数据，如DNA序列中的特定位点的核苷酸。建树时着重分析运筹分类单位或序列间每个特征（如核苷酸位点）的进化关系等。属于这一类的方法有最大简约法（umarionymethod、最大似然法（umlklihoodmethod、进化简约法（evluinayarioymthod、相容性方法（comptbiiy）等。对相似性和距离数据，在重建系统发生树时只能利用距离法。离散特征数据通过适当的方法可转换成距离数据，因此，对于这类数据在重建系统发生树时既可以用距离法，亦可以采用离散特征法。根据建树算法在执行过程中采用的搜索方式，系统发生树的构建方法也可以分为以下3类。第一类是穷尽搜索方法：即产生所有可能的树，然后根据评价标准选择一棵最优的树。需要注意的是，系统发生树可能的个数随序列的个数急剧增加，假设要为n个物种建立系统发生树，则不同拓扑结构的树有(2n-3)!!棵。如果n为20，则大约有1021棵可能的树。假设序列的个数为7，则所有可能的有根树的个数为10,395，所有可能的无根树的个数为795。当序列的数目超过7时，由于树的个数太多，不可能采用穷尽搜索的方式来求取最优树。由于可能的树非常多，从计算量来看，这种方法只能处理很少的物种。第二类方法是分支约束方法：即根据一定的约束条件将搜索空间限制在一定范围内，产生可能的树，然后择优。这是人工智能技术中的一种空间搜索策略，这种搜索方式不需要搜索整个树空间，可大大提高搜索效率。第三类是启发式或经验性方法：根据先验知识或一定的指导性规则压缩搜索空间，提高计算速度。这种方法能够处理大量的分类单元，虽然不能保证所构建的树是最优的，但实际结果往往接近于最优解。当待分析的对象个数比较多时，必须采用分支约束方法或者启发式的方法。在构造系统发生树时需要考虑进化假设和进化模型(SwofordadOlen1990,Liandraur191)。系统发生树的类型可能是有根树，这意味着其中的一个序列代表其他所有序列共同的祖先，到达该序列的分支是树的最底层分支。另一方面，系统发生树可能是无根树，意味着没有共同的祖先。一般认为样本序列是随机进化的，序列中的所有位点的进化也是随机的而且是独立的。在进行具体的系统发生分析时，一般还要作一些假设（BaxeaisandFrancis，1998：序列必须是正确源（或平行进化）序列，序列之间的系统发育史是相同的，在序列比对中，不同序列的同一个位点都是同源的。另外，当两个物种在系统发生树上分化后，各自独立进化发展。对系统分析的样本也有一定的要求，要求样本足以反映感的问题，样本序列之间的差异包含了足以解决感的问题的系统发生信号。通过某种算法构造好一棵系统发生树之后，需要对树的合理性和可靠性进行分析。对于若干个序列，如果利用多种不同的分析方法进行系统发生分析，并且得到相似的进化关系，那么分析结果具有较高的度。系统发生分析一般是建立在“分子钟”基础上的。生物随着时间的推进而演化，进化的速率被视为进化研究中的基本问题之一（等，2000。进化速率就是在某一段时间内的遗传改变量。分子进化速率相关的分子钟的概念源于对蛋白质序列的研究。科学家们在比较几种动物的血红蛋白、细胞色素C的序列后注意到：这些蛋白质的氨基酸取代速率在不同的种系间大致相同，即分子水平的进化存在恒速现象（ZuckrnlandPaulig，162，95;rwntl.,1972。于是提出了进化在分子水平存在“时加。现在已知道不仅是氨基酸序列，ＤＮＡ序列也存在这种规律。Kimura进一步提出了具体的分子进化观点（Kimura983：对于各物种的每个蛋白质，如果用每个位点每年发生的氨基酸替换次数作为衡量分子进化的速率，则该速率是大致恒定的；功能上次要的分子（或者分子部分）的进化速率比功能重要的分子（或者分子部分）进化速率快；对现有分子结构或者功能破坏小的氨基酸替换比破坏力大的氨基酸替换发生得更加频繁。第二节基于距离的系统发生树构建立系统发生树的基本任务是：在给定的条件下（包括物种、物种的特征值或者序列，构造一棵最优的系统发生树。这里重点针对DNA序列或者蛋白质序列构建系统发生树。基于距离的系统发生树构建方法的基本思路是：给定一种序列之间距离的测度，在该距离测度下构建一棵系统发生树，使得该树能够最好地反映已知序列之间的距离。这种方法采用两两距离，建立一个距离矩阵，如表6.1所示。6.110(Saitou,1996)12345678923456789这里的距离代表两条序列之间各位点核苷酸替换数的（代价，也可DNADNA序列，可以采用等价矩阵（Li1997，也可以采用核苷酸转换矩阵（MichnerandSkl197或者采用其他具有非对称置换频率的矩阵（Felsesen，188；aoetl.，195）。对于蛋白质序列，一般采用PAM矩阵比较合适，因为该矩阵就是针对研究蛋白质序列之间的进化距离而设计的。为了便于分析，首先定义续加和距离函数，在该函数下，两个物种之间的距离与系统发生树中连接这两个物种的分支总长度成正比。这样，如果物种a和物种b由经过中间节点v的两条边相连，两条边的长度分别为dav和dbv，则它们之间的距离为dav+dv。这样，可以在系统发生树中确定a和b的相对位置。进一步，假设三个物种之间的距离分别为dabdacdc,可以通过求解线性方程计算出系统发生树的距离dav、dv和dcv。例，如果有三个物种，其两两距dab=dac=dbc=通过求解方程，得到如图6.2所示的u v ba图6.2根据方程求解结果构造系图6.2是一个简单的例子，在实际工作中，所处理的物种可能很多，因而需要求解的线性方程也很多，难以求解，或者方程组的求解过程存在着不确定性，因此需要采用数学近的方法。1、最小二乘法统计学上近的方法之一是最小二乘法。的目标是构造一棵树T，该树的叶节点代表物种，用该树物种之间的距离。通过优化，使下式最小化：nSSQ(T)Wij(Dijdij)i1

（6-这里，Dij为物ij的实际观察距离（或序列之间的计算距离）dij是物种i和j在系统发生树T中的距离，Wij是与物种i和j相关的权值。SSQ（T）是树T所有值与实际观察值偏差的累加和。权值Wij一般为1，或

=1/D （6-一般来说，寻找一棵最小方差树是一个NP-完全问题（DayWHE，1986），需要采用近似的算法。下面几种计算复杂度为多项式的启发式方法，即连锁聚类方inkageclustering、非组平均法（UPGMA）和邻近归并法（NeighborJoining2、连锁聚n个叶节点表示n个物种（序列，然后将进化过程中的祖先赋予树的节点。假设若干个序列是从一个共同的祖先进化而来，则系统S代表一个物种集合，定义在该集合上的距离测度d(x,y)（x,yS）满>0for（6-=0for（6-=（6- （6-6.3连锁聚类结果超距离测度d(x,y)满足下述条件： （6-给定序列，通过序列之间的两两对比排列，计算序列之间的进化距离，并假设得到的距离矩阵是一种满足超距离测度条件，则系统发生树的聚类构造过程如下：首先从距离矩阵中选择距离最小的一对序列，然后将这两个序列合二为一，形成一个新的对象（代表这两个序列的祖先，并重新计算这个新的对象与其它序列的距离。不同的实现方案采用不同的计算公式：单连锁（6-最大连（6-平均连（6-其中x代表y和z的合并，u代表任意其它对象。重复上述过程，直到所有的物种都被合并到一类为止。图6.3就是一个通过平均连锁聚类分析的结果。其一般过程如为：首先计算所有序列的两两比对，得到距离值，然后利用连锁聚类算法构造树，采多重比对问题转化为序列两两比对问题，3、非分组平均法（UnweightedPairGroupMethod,非分组平均法是一种较为常用的聚类分析方法(MichnerndSkal,15)，最早是用来（即若干个分类单元形成的集合）之间距离的计算公式不一样。当用来构建系统发生树时，其假定的前提条件是：在进化过程中，每一世系发生趋异的次数相同，即核苷酸或氨基酸的替换速率是均等且恒定的。通过非组平均法所产生的系统发生树可以说是物种树的简单体现，在每一次趋异发生后，从共同祖节点到两个分类单元间的分支长度一样。因此，这种方法较多地用于物种树的重建。非加权分组平均法在算法上较简单。在聚类之前，令所有的分类单元各自成为一类；然后反复将最邻近（距离最小）的两个类形成一个新的聚类，其分支点位于原来2个分类距离的1/2处，计算新的分类与分类间的平均距离；如此反复，直到所有的分类单元都聚为一类，最终得到一个完整的系统发生树。系统发生树的结构与分类过程相对应。设d是物种之间的距离函数，对于两个分类Ci和Cj，定义它们之间的距离DijD d(p,ni ni

（6-这里，ni=|Ci|,nj=|Cj|,分别为两个分类中分类单元的个为了给树中的每个分支赋予分支长度，使SSQ(T)最小，UPGMA通常能够获得满意的结果。该算初始化：使每个物种自成一类，如果有n个物种，则开始时共有n个类，每个类的大小为1，n个叶节点代表每个类；寻找具有最小距离Dij的两i、建立一个新的聚类（ij),该聚类i和类j的合ij(ij)，生长两个新的分支，将i和j连接到(ijDij/2；计算新的分类到其它类的距离，计算公式为下列的平均距离，

nin

nnin

（6-在距离矩阵中删除与类i和类j相应的行和列，为类(ij)加入新的行和列，其值按上述公运用UPGMA算法，得到一棵有根的系统发生树，从树根到任何叶节点的分支长度全都一样。也起点）出发，踏着同样的节律，沿不同的路径，演化成为当前的形式。而任意两个物种之间的距离是连接这两个物种所在叶节点路径上的距离和。在解决实际问题时，如果确实存在分子钟，则UPGMA能够保证找到最优的答案。该算法的示例见图6.4。图6.5是针对表6.1中的距离矩阵而构建系统发生树。 图6.4UPGMA和邻近归并法(NJ)分析结果(Shamir 图6.5UPGMA所构造的系统邻归并法是一种快速的类方法，方法Saitou和Nei于1987年首次(SitouandNei,197)。在构建系统发生树时，它取消了非分组平均法所作的假定，不需要关于分子钟的假设，在进化分支上，发生趋异的次数可以不同。最近的计算机模拟已表明它是最有效的基于距离数据重建系统发生树的方法之一。与非分组平均法相比，邻近归并法在算法上相对较复杂，它的是树上的节点而不是分类单元。这种方法的基本思想是：在进行类的合并时，不仅要求待合并的类是相近的，同时还要求待合并的类远离其它的类。在聚类过程中，根据原始距离矩阵，根据所有节点间的平均趋异程度对每两个节点间的距离进行调整，即将每个分类单元的趋异程度标准化，从而形成一个新的距离矩阵。重建时将距离最小的两个叶节点连接起来，合并这两个叶节点所代表的分类，形成一个新的分类。在树中增加一个父节点，并在距离矩阵中加入新的分类，同时删除原来的两个分类。随后，新增加的父节点被看成为叶节点，重复上一次循环。在每一次循环过程中，都有两个叶节点被一个新的父节点所取代,两个类被为一个新类。整个循环直到只剩一个类为止。从所得到的系统发生树来看，两个聚在一起的分类单元其所在的叶节点到父节点的距离并不一定相同。在每一次的循环中，在树中寻找两个物种的直接祖先。对于节点i，到其它节点的距离ui按下ui=ki(Dik/(n- （6-这里Dik是分i和分类k之间的距离，是动态更新的距离矩阵D中的元素。为了使所有分所示。新分支的长度dk,(ij)用求解线性方程的方法进行计算（6.2j图 寻找一对节点，使这两个节点靠近，但同时远离其它节算法如初始化（UPGMA一样循对于所有的分类单元iuiki(Dik/(n-选择一对分类单元i和j，使Dij-ui-uji和j归并为新的类(ij)，在树中添加一个新的节点，代表新生成的分类，计算从idi,(ij）=1/2Di,j+1/2(ui-uj),dj,(ij）=1/2Di,j+1/2(uj-D(ij）,k=1/2(Di,k+Dj,k-删除聚类ijDi,j的分支连接剩余的两个类。邻近归并法的处理示意如图6.4所示，针对表6.1中距离矩阵的实际计算结果见表6.2，构造的系统发生树如图6.7所示。详细的邻近归并算法见（等，1997。6.2邻近归并算法计算结果(OUT9=14.632(4.530E-10(1.407E-OUT2=-0.065(-2.011E-3(3.311E-Node12=4.463(1.382E-4(8.339E-Node13=4.281(1.325E-5(3.575E-OUT8= (1.680E-=(2.820E-OUT7= (5.932E-=(3.822E-Node14=5.102(1.580E-6(1.646E-(Last1.675(5.186E-0.907(2.808E-1.410(4.364E-Node或OUT后面的数字代表节点（或分类单元）与节点之间分支的长度，该长度对应于序列比较区域的核苷酸替换（括号里的数字对应每个位点的替换。2 45 9 图6.7利用邻近归并算法构造的系统在基于距离的建树方法中还有Fitargoliash法（FM,FitchandMrglash,97，该方将树中的序列两两距离的合适度极大化。另法是最小进化方法（ME,MinimumEvolution最小进化概念首先由Cavalli-Sforza等人用于邻近归并方法（Cavalli-SforzaandEdwards1967小进化方法先使用Fitch和Margoliash方法估计系统发生树上的分支长度，然后进行优化，求出所有分支长度和最小的树。利用最小二乘法进行优化，将观察到的距离相对于进化树距离的偏差的平方最小化（RzhetskyandNei,1992;Swoffordetal.,1996a;Felsenstein,1997FM方法不同，ME方法并不使用所有可能的序列两两距离和所有可能的相关路径长度，而是先根据树中的节点到叶节点之间的距离固定进化树节点的位置，然后根据这些观察点之间的最小计算误差，对的树枝长度进行优化。最小进化方法的具体求解过程是：首先利用邻近归并法构造出候选的树，再搜索与候选树具有6.1ME-REE进行计算，其结果如表6.3所示。所形成的系统发生树与图6-7相似，但是略有不同，优化以后的树的分支长度总和比图6.7中的系统发生树的分支长度总和短分分支长度置信度1—2—-3—4—5—6—7—8—9—10—11—12—13—14—15—16—17—第三节基于特征的系基于特征的系统发生分析的一般问题是：给定n个物种，m个用以描述物种的特征，以及每个物种所对应的特征值，构建一棵系统发生树，使得某个目标函数最大。输入一般为n×m的矩阵M，其中Mij代表第i物种之第j个特征的取值，如表6.4所示。在构建系统发生树假设特征是相互独立的，6.4分类单甲CAGGTA乙CAGACA丙CGGGTA丁TGCACT戊TGCGTA对于给定的条件，要在很多可能的树中找一个最佳的树。在实际应用中，不可能穷尽搜索所有可能的树，必须按照一定的方法、一定的策略在较短的时间内得到比较好的结果。1、最大简约法(um最大简约法最早是基于形态特征分类的需要而发展起来的，因使用的算法不同而有着许多版本，其中有些已被广泛地用于分子进化研究中根据离散特征数据构建系统发生树，如对DNA序列数据分析。最大简约法的目标是构造一棵反映分类单元之间最小变化的系统发生树。最大简约法利用的只是对简约分析能提供信息的特征，如在DNA序列数据中，利用的只是存在于核苷酸序列差异（至少有两种不同类型的核苷酸）的位点，这些位点称为简约信息位点（parsmonyifomtiesie对于系统发生树最直观的代价计算就是沿着各个分支累加特征变化的数目，而所谓简约就是使代价最小。利用最大简约方法构建系统发生树，实际上是一个对给定分类单元所有可能的树进行比较的过程，针对某一个可能的树，首先对每个位点祖先序列的核苷酸组成做出推断，然后统计每个位点用来阐明差异的核苷酸最小替换数目。在整个树中，所有简约信息位点最小核苷酸替换数的总和称为树的长度。通过比较所有可能树，选择其中长度最小的树作为最终的系统发生树，即最大简约树（umparsimonytree假设图6.8是根据表6.4所建立的一棵系统发生树。对于位点1，由于甲、乙、丙都取“C”值，边上有一个变化（左分支T变为C，或者右分支C变为T），就能解释各个分类单元在该位点的数据。对于位点2，甲、乙取“A”值，而丙、丁、戊取“G”值，通过“节点2”可以将它们分开，因此，“节点2”所连接的左分支上的一个变化（G变为A）就能解释各个分类单元在该位点的差异。位AGAG需要在“节点1”和“节点3”所连接的边上分别设置一个变换。因此位点4需要两个变的分支。总的来说，对于图6.8所示的整个系统发生树，共需要8个变节点节点节点 节点节点节点对于表6.4中所列出的5个分类单元，还可以形成其它拓扑结构的系统发生树，这些树为解释表6.4中各个位点的数据，可能需要或者更少的变化。而最大简约法就是要寻找一棵最小变化的系统发生树。最大简约法的针对待比较的物种，选择核酸或蛋白质序列。有些分子比其它分子变化慢，适合于进行距离分析，例如哺乳类的线粒体DNA、管家蛋白质等；根据每个序列比对的位置（即多重对比排列的每一列，确定相应的系统发生树，该树用最少的进化动作产生序列的差异，最终生成完整的树。对于一棵系统发生树T，假设树中的节点用V(T)表示，树的E(T)，以uj、vj分别表u和v的第j个特征，则树T的代价为S(T) {j|vjuj(u,v)E(T

给定一棵有根的系统发生树，所关心的是：如何计算该树的最小变化？如何标定树的节点？假设以vc代表v节点特征c的值，单特征（即每个物种所对应的序列仅含一个字符，作为单个特征的值）Fitch算法（itch，971v，其值vc取对应物种的特征值。然后执行下面两步：给每个节点v赋予一个集合Sv：如果v是叶节点，则Sv={vc}；如果v是节点，并且u、w是其子节点，如果SuSw，则Sv=SuSw；否则S(v)=SuSw。这个过程是从叶节点开始，直给定集合Sv，为每个节点v的特征c赋予值vc。如果v有一个父节点u满足ucSv,则将uc赋予vc，否则任取一个tSv赋予vc。这个过程的执行方向刚好与上一个过程相反，即从树根出发，直至叶节点为止，最后得到完全标定的树。应按前序遍历方式依次处理每个节点。最后树中变化的个数等于第一步计算得到空交集SuSw的个数。图6.9是对5个序列（每个序列仅含一个字符，作为物种特征）执行过程（1）的结果，最小的总代价为3，图点标注‘*’表明该节点的子节点交集SuSw的为空。 {CG6.9最大简约法处理结果示意（Shamir.上述算法对特征值之间变换的代价处理，不考虑其中的差别。其实各种特征值之间变换的代价是不一样的，在建立数学模型时应该区分这些差别。C

c代表特征c的特征值从i变化为j的代价。将上述算法改进为求最大简约树，这就比Fitch算法更一般化的Sakof（Snkf,17)算法。该算法可以处理多个特征，并考虑不同特征值变换代价的差异。计算过程如下：cv和特征t，计算St(v)，它代表树根为v的子vc=t的情况下最小代价值。Fitch算法中的第一步一样，后序遍历。对每个叶节点v：ccSc(v) v （6-c 对于每个节点v，设其两个子节点分别为 和w，则Sc(v)min{CcSc(u)}min{CcSc （6- 树根r的最小总代价等于S(T)minSct

（6-这里，m为特根据第一步的计算，确定每个特征在节点的最佳值对于根节点rrargminSc （6- 对于其他节点v，假设其父节点为u，vargmin(CcSc （6- 前面的工作实际上是针对已知拓扑结构的系统发生树，确定其中的节点标记，即确定树的节点特征取值。这仅仅是解决问题的一个步骤，或者是后面的一个步骤。关键的问题是系统发生树的拓扑结构又是如何确定的呢？实际上是要在众多可能的拓扑结构中选择一棵与给定数据相一致的树。这是一个典型的搜索问题，而且搜索空间非常大。分支约束（BranhndBundB&B）技术可用于解决这个问题(HenyandPeny,182)。分支约束是在一个复杂的空间中进行搜索的通用技术，搜索空间以从一个分层树的根节点至叶节点的一系列路径表示。该技术由Hendy和Peny首先引入到简约法之中(yd在的情况下，该算法的时间复杂度与穷尽搜索相近。最简单的算法形式是按照一定的顺序遍历搜索树，保存到目前为止到达叶节点路径的最小的代价，并作为代价的上限，记为B。在后续搜索过程中，如果到达搜索树的某个节点，并且到达该节点路径的代价大于B，则不再搜索以该节点为树根的子树，搜索树在此节点处不再继续扩展，因为即使搜索该节点的子树，也不会得到比B更好的结果。上述行为称为搜索树在某个节点剪枝。在开始搜索之前，可以根据期望值，或根据其他线索设定初始的上限值B。考虑这样一棵搜索树，以搜索树的节点代表相应的系统发生树。在搜索树的第k层上的节点代表具有k个叶节点，对应于前k个物种的所有可能的系统发生树，第k层上的节点的子节点将代表在上述系统发生树中加入第k+1个物种而形成的新系统发生树。搜索树的树根代表空系统发生树。由于在一棵系统发生树上增加一个新节点不会减小简约代价，所以搜索树满足单调条件，分支约束技术可以帮助对搜索树进行剪枝，从而压缩搜索空间，提高搜索效率。另一种优化方法技术是局部搜索，其基本思想是对于一个搜索空间定义邻近关系，在搜索空间遍历时，从一个点转移到它的邻近点。这样可以使用许多启发式算法，如贪婪算法、模拟退火算法等。总之，最终得到的是一个局部最优结果，并且接近于全局最优结果，或者就是全局最优结果。有许多定义邻近关系的方法,例如，在对于一棵系统发生树，分解它的任何子树，并重新拼装成与原树不同的树，那么这些新树都是原树的邻近树。2、最大似然法（um前面介绍的系统发生分析方法隐含地使用了各种概率模型，说明生物分子序列是如何进化的，通过系统发生树研究序列之间的进化关系。最大似然法（Felsesen981;Toreetl.19;etal.95;ulenekndanla97）明确地使用概率模型，其目标是寻找能够以较高概率产生观察数据的系统发生树。最大似然法的基本模型有两个部分：一是分支过程，它描述物种在进化历程上是如何的；二是采样模型，它描述生物学家是如何选择一组的物种进行研究的。采样模型比分支过程更难形式化，因此常常被忽略。一般的采样模型假设待分析的物种是被随机、独立地选中的。在生物学研究的许多领域，最大似然法常被用来检验一些假设，目前较多地用在遗传作图与临床检验上。因该法是建立在许多假定的基础上，并且获得最大似然估计的解较复杂，因此限制了它的推广应用。在DNA分子标记研究中，该技术已被用来重建由DNA序列数据及限制性位点数据产生的系统发生树。利用最大似然法来推断一组序列的系统发生树，需首先确定序列进化的模型。然后基于一定的模型考虑两个物种序列间的关系，计算分支的长度。这个过程需要寻找在某一进化距离上由第一种序列真正转换成第二种序列的可能性，并确定在最大可能下的进化距离。接着生成多个物种所构对于一棵给定的树，希望有一种评价的方法（KisinoandHaegaa，189。可以用可能性得分评估所作出的假设，即评价所得到的系统发生树T。对于给定的一组物种，假设它们的观察值为（M为向量，可以选择一棵树，使得P(M|T)最大，即最大似然法。在下面的中假定已经知道树的拓扑结构，目标是寻找最优的分支长度。在基于DNA或蛋白质序列的系统发生分析方面，与最大简约法相似，最大似然法首先依赖于一个合理可靠的多重序列的比对，然后检测每一列的变化。对于每一个可能的树，计算在每一列发现真实序列变化的可能性，将每个排列位置的概率相乘，其结果作为每棵树的可能性。具有最大似然值的树就是最可能的树。假设特征是两两独立的，分支是一个Markov过程，即节点拥有一个给定标识的概率仅仅是其父节即在时间tvu内，从状态x转换到状态y的概率，其中v、u代表系统发生树中的节点。假设在进化过程中特征的出现频率是恒定的，其值P(x)。假设有一个矩阵M，它是关于n个物种实际观察值，M描述每个物种m个特征的具体取值。同时假设存在一棵树T，其叶节点（如vu对应于这些物种，而树中的分支代表物种之间的距离tvu，求该树的似然值L=PM|T)。首先处理最简单的情况，即每个物种仅有一个特征。由于树种节点的标记未知，需要考虑所有可能的树，并加和对应的结果。例如，对于图6.10所示的树，可以写出公式：LP(M|T)P(r)Prs(trs)Prv(trv)Pvu(tvu)Pvw(tvw 其中r、v是节点可能的标识r v wu图6.10系统发生树似然值的计算将上述计算公式推广到多个特征，只需分别对每个特征进行反复的计算，然后将结果相乘（因为假定特征两两独立。推广的计算公式为：LP(M|T) P(Mj|Tcharacter P(Mj,R|Tcharacterjreconstruction P(root)Puv(tuv （6-characterjreconstruction 利用最大似然法得到树一般是有根树，然而如果特征替换是可逆的，即Pxy(tPyx(t)，则所对于一个特征j，定Cj(x,v)P(根为vCj(x,v)是关于子树节点v的条件概率，即节点v的第j个特征具有标识x成为子树的可能性。下首先进行初始化工作，对于所有的叶节点v和标识符号xC(x,v)

if(vj （6-u和w，x计算 Cj(x,v)Cj(y,u)Pxy(tvu)Cj(y,w)Pxy(tvw （6- 最终结m LCj(x,root)P( （6-j1 下面如何根据给定的树的拓扑结构寻找最优分支长度。首先假设除trv之外，所有的分支长度已知。如果r根节点，则 logLlogP(x)Cj(x,r)Pxy(trv)Cj(y, （6-j x, 可以采用许多不同的方法使logL最大，如Newton-Raphson方法，或者最大期望EM算法。如果r不是树根，也可以应用相似的方法。如果对于所有x、y、t，Pxy(t)=Pyx(t),则树根可以被赋予任何节点，而不影响L的值。换句话说，为了寻找节点r和v之间最优的分支长度，仅需要处理来自r接下来的一步是在没有任何关于分支长度先验信息的情况下，寻找最优的分支长度。主要的问题是，如果某个分支的长度发生变化，不能保证其它分支依然是最优的长度。在实际处理中，可以采用局部处理方式，在某个时刻仅优化一条边，通过遍历树，逐步处理所有的边，其最终结果与最优的系统发生树接近。3、相容性方法在解决系统发生问题时，定义目标函数考虑的另一方面是相容性(comatiilt)，即与一棵树相一致的特征个数。很显然，相容的特征数越多越好。相容性方法实际上是简约方法的一种简化，在所有的特征都是二值的情况下，这种方法非常有用。假设一个二值特征有两个可能的取值“0”和“1S代表一组分类单元，T是关于S的系统发生树。如果在解释叶节点中分类单元的特征数据时，只需要沿T的一条边变化，则称该特征与系统发生树T是相容的。图6.11是一棵具有4个分类单元的系统发生树，其中在图6.1(a)中特征是相容的，只需要一条边变化即可（标记星号的边；而在图6.11(b)中，特征是不相容的，因为需要两条边发生变化。*** *** 图6.11具有4个分类单元的系统发生树：(a)相容特征；(b)不相容定义：对于一个具有k个值的特征c，称c与树T相容，当且仅当T的节点存在一种标记方式，使得c变化总数为k-1。在简约方法中，力求找到最少的分支变化的系统发生树，使之与给定的各个分类单元特征值匹配。如果某个特征值在各个叶节点中出现n种状态，则至少需要有n-1次分支变化。相容性方法将树分成两种类型，即只需要最小变化的树和其它树，而最大简约法要区分各种不同的树。相比之下，相容性方法仅需要维护一小部分关于树相对性能的信息。给定一个DNA单碱基位点，如果所有物种在此位点或T‘C’,则该特征的取值数为2，而不是4。在下面的中，假设所有特征的取值都是二态的，即只有两种取值，分别用0和1可以证明，在二态的情况下，相容性是简约性的一种特例。只要在最大简约法中，对于每个变化至少两次的特征，规定其代价为2。对于一棵树Tni代表变化i次的特征个数，n0是叶节点中值恒定的特征个数，n1是需要发生一次变化的特征个数，则树T的总代价S(T)=n1+2(n-n0-n1)=2n-2n0- （6-其中n是总的特征个数，n和n0是固定不变的。于是代价最小化等价于使n1最大，而1征的个数。这样给出了一种计算最优相容树的思路，但这需要解决复杂的简约问题，而简约问题又是一个NP-完全问题，所以必须考虑其他的方法。处理相容问题的第一步与最大简约法相似。给定一棵叶节点已标记的树T，搜索最佳的相容代用Fitch算法M的情况下，寻找最佳的系统发生树。Sc1c2c1、c2相容的系统发生树？假设有两个特征c1c2,如果存在一棵树T，使得c1和c2各自与T相容，则称特征c1c2是两两相PC(c1,c2kc1、c2、…、ckT，使得c1、c2、…、ck都分别与Tc1、c2、…、ck是连带（jointlycompatible）6.1（两两相容性检Wilson,1965）M是分类单元特征取值矩阵，对于特征i,j，定义集合Sij, ={(x,y)|存在分类单元k使得Mki=x并且Mkj=y} 则当且仅当S{0,1}2，PC(c,c’)成立（PC(c,c’)代表c和c’两两相容。这里 014 证明：假设S{0,1}2，集合S最多有3元素。在第一种情况下，S仅有一个元素，则特征c和c’各自有一个单独的状态，但是这不可能，因为所的特征都是二值特征。在第二种情况下，Scc’有两个元素，则实际上可以将上述两个特征视为一个二值特征。 S有三个元素，可以假设{0,1}2-

={(x,y)}，即集合

中不包括特征值对(x,y)。图6.12与两个特征相容的树的基本结构，其特征组合值为（x,y）,(x,y),(x,y)。这里符号“”代表“取反”操作。图6.12中每个三角形表示一棵子树，每棵子树中两个特征的值不变。沿粗线段标记的y图6.12一棵与两个特征相容的树，有三个组合值举例说明定理6.1。假设有5个种属A、B、C、D、E，两个特征c1和2，各个种属特征取值情况见表6.(Flentin,198)。为c1和2建立一个22的特征值组合表（相当于定理6.1中的），表元素（i，j）的值反应两个特征取值的组合情况。如果某个种属的两个特征值分别为i和j，则特征值组合表元素（ij的值为“1”；如果找不到这样的组合，则表元素（ij的值为“0（如表6.6所示）。对于这样的特征值组合表，特征c1和c2是两两相容的当且仅当特征值组合表至少存在一个“0”元素。表6.5ABCDE1011000110表6.6c1=c1=c2=11c2=01假设特征值组合表中存在“0”元素，那么在所有4种特征取值组合中，有的组合没有出现。（i，j）T，c1=1-iT1下，而将所有特征值c1=i的种属放在其它子树中。令c1=i，c2=1-j为T的祖先节点。显然，c1T因为只需要在连接子树T1根节点的分支上一个变化即可。另外，当c1=1-i时，对于特征c2只T16.5T1A、C、D。Tc1=0，c2=1-1=0。很显然，Tc1的取值分为左右两个部分，c1TT1c2的取值也分为两个部分，因此c2与T也是相容的。 A(1,0) 6.136.2（两两相容性定理,Estabrook,1976）在一个二值特征集合C中，所有特征是连带相容的，c,c’C,PC(c,c’)成立。利用数学归纳证明。如果特征数为2，存在两个特征c1和c2，则定理6.1和上面的例子已经说明可以构造一棵树T，使得c1、c2分别与T相容。这表明归纳基础成立。假设定理6.2对于所有特征数m<n皆成立，则对于n个特征，任意选择一个特征c，根据c的取值将所有物种分为两个部分。根据条件，每个特征c’与c都是两两相容的，由定理6.1可知，特c’与c的某种取值组合一定不会在S中的物种出现。令这样的取值组合为（c=i，c’=j），于是对于所有特征c值等i的种属，它们的特征c’值只能为1-j。这样根据c将特征C分为两个集合，C0C1。C0包含所有满足下列条件的特征：如果c=0，这些特征将取多个值。C1包含所有满足下列条件的特征：如果c=1，这些特征将取多个值。显然，属于C0的特征在c=1的情况下只能取C1c=0的情况下也只能取一个值。S0S1c=0c=1v0v1C0的（或C1）的特征c’，v0（或者v1）的特征c’的取单个值，其值等于所有c=1（c=0）的种属的c’值。由于所有C0v0的状态与集合S1中的其它种属相符合，因此C0中的特征在种属集合S0{v0}中是两两相容的。同样，C1中的特征在种属集合S1{v1}中是两两相容的。由于C0和C1的特征个数小于n，根据归纳假设，C0、C1中特征是连带相容的。令T0（T1）是定义在S0{v0}（S1{v1}）的树，在这样的树中，属于C0（C1）的特征是连带相容的。为了构造一棵树T，使所有属于C的特征连带相容，在T0T1中寻找种属v0v1的位置，然后用一条公共的边eT0和T1e替代原来连接到v0和v1的边。通过上述操作形成树T，Tc相容，ceC0c’，T1中的所有种属与v0具有同样的值，因此如果v0T1替换，不需要额外的变化。这样，如果c’与T0相容，c’也与T相容。对于属于C1的特征也有同样的结论。由此可以断定，每个属于C的特征与T是相容，集合C中的特征是连带相容的。根据上述定理，解决“最大相容”的问题被简化为寻找最大的特征连带相容集合，归结为在两两相容图中寻找最大完全子图。两两相容图G定义如下:G=(V,E);V={v1,v2,…,vm}; （6-这里，V代表图中顶点的集合，而E代表分支或者边的集合。图论中有许多寻找最大完全子图的有效算法，比如基于分支约束的寻找最大完全子图算法。一般而言，对于相同的数据，相容性方法运算速度快于简约法。i*在相容图中找到最大完全子图，意味着找到最大的连带相容特征集合C*，接下来的任务是构造系统发生树。只要按照C*中的每个特征相继分割物种集合。在下面的算法中，给定矩阵M和特征集合C*，反复构造无根树。在每一步，利用C*中的一个特征扩展树T。集合C保存到目前为止尚未用到的特征。T的每个节点或者是没有标记，或者标记为Lv{1,2,…,n}，代表种属。对每个特征c和标记的节点v,令Lv()L{|k=i}，k代表特征c的取值。算法过程如下：i*初始初始T，使T仅有一个节点r，标记为循C中选择一个特征

r={1,2,…,n}，并C=Cn为种属的对于V(T)vLv元素的个数超过1，则 如果|Lvi(c)|0（i=0,1），则在T中增加两个新v0v1，分别标L0(c)、L1(c)，v0和v1到v 修改C,使得C=C重复执行上述循环过程，直到C为空算法执行的结果是产生无根节点的系统发生树，叶节点全部被标记，至于节点的标记工作相对简单，可以用Fitch算法。第四节１、各种方法比较在距离法中，连锁聚类方法比较简单，非分组平均法比较实用，当使用的距离数据是来源于多个的分析结果时，利用非分组平均法能得到可靠的系统发生树。对于离散特征分析方法，在不同世系间进化速率相差较大，并且进化速率恒定而树的分支很短的情况下，最大简约法并不能对一个真正的系统发生树作出始终一致的判断。即使有时最大简约法能得到一个始终一致的判断，但它获得一个正确树的效率通常要比邻近和最大似然法低(常青等,1998)。但在序列趋异程度较小、核苷酸替换的速率或多或少的恒定、没有很高的转换与颠换比及很强的Ｇ＋Ｃ含量偏差、所分析的核苷酸数量较多（大于几千）的情况下，最大简约法仍是一种较好的系统发生树构建法。另一方面，与距离法和最大似然法不同，最大简约法能利用序列中碱基的与缺失信息。在实际使用最大简约法和最大似然法构建系统发生树时，当给定的分类单元数量m（比如m<10可通过计算机对所有可能的树作穷尽搜索，确定最理想的树。这样做虽然非常耗时，但还是有可能的。当m大于10时，就不1)分支约束法，从一组具有潜在可能的树中确定优化的系统发生树。该方法能保证获得很好的结果，但当m大于或等于20时计算量仍然很大；(2)启发式搜索法，该法在分析中只对少部分的可能树进行比较，故m可取较大值，但此法不能保证发现最好的树。２、系统发生在系统发生推断中，统计分析的系统误差和随机误差均影响所建树的可靠性。对随机误差的影响，常采用一定的统计检验来分析获得的系统发生树的可靠性。一种是利用某一参量来对所获得树及其相近树进行结构差异检验。在最大似然法中常利用似然值。这种方法是一种保守检验，而且检验的程序非常复杂，需要很大的计算机内存。另一种类型是分析每个分支可靠性，其中常用的方法有(常青等,198)：（confidenceprobabilityCP，CP值越高，分支的长度也就越可靠。通常，当≥0.95或0.99时，可认为该分支的长度在统计上有效自举检验，这是一种重抽样技术，可用来估计在取样分布不知道或难以分析得到的情况下分支与统计有关的变异性。通过自举检验，可得到一个自举置信水平（boottapcnieneleel，简称BL。计算机模拟已表明当BCL＞09时，CP值与BCL值二者是非常相近的。与自举检验相近的另一种重抽样方法是弃半检验。有研究表明，在研究的核苷酸数量较少的情况下，即使CP或BCL值达到9%，所获得的结果仍然不十分。因此，在研究中应从不同的中尽可能分析较多数量的核苷酸，特别在研究不同生物间进化关系时，因为不同遭受的进化压力不同。对系统误差，需要降低它们对系统发生分析影响，增加所建系统发生树的可靠性。通常采用的（1）重新考虑分析时的假定，变换分析方法；（2）除去树中的长分支，因为一棵树中如果具有许多长分支，将会使分析中的误差复杂化；（3）（4）对某些特征或某一特征状态进行处理等。第五节本节通过实例说明系统发生分析的过程。这里使用莫斯科国立大学的系统发生树分析的对象是13条来自不同物种的同源蛋白质。对这些蛋白质序列进行多重序列比对，其结果见图6.14所示。根据序列比对结果计算序列之间的距离，生成距离矩阵。然后分别利用聚类方法和拓扑学方法建立系统发生树。拓扑学算法首先优化系统发生树的结构，确定树节点之间的分支连接，然后再计算分支的长度。该算法的基础是拓扑相似原理。图6.15是通过聚类分析得到的系统发生树，用Phylip格式表示如下：HumanHumanFLRTPKIVSGKDYNVTANSKLVIITAGARQChickenFLKTPKITSGKDYSVTAHSKLVIVTAGARQDogfishFLHTAKIVSGKDYSVSAGSKLVVITAGARQLampreyFLKTAKIVADKDYSVTAGSRLVVVTAGARQBarleyFLPRVRI-SGTDAAVTKNSDLVIVTAGARQMaizeyFLPRTRLVSGTDMSVTRGSDLVIVTAGARQLacto_caseiFTSPKKIYSA-EYSDAKDADLVVITAGAPQBacillus_steaFAPKPVDIWHGDYDDCRDADLVVICAGANQLacto_plantFTAPKKIYSG-EYSDCKDADLVVITAGAPQTherma_mariFTRRANIYAG-DYADLKGSDVVIVAAGVPQBifidoFYPTVSIDGSDDPEICRDADMVVITAGPRQThermus_aquaFAHPVWVRSGW-YEDLEGARVVIVAAGVAQMycoplasmaSLPFPISVSRYEYKDLKDADFIVITAGRPQ6.14|||||||||||||||||||||||||||||||| 6.156.14||

||||||||||||||||||||||| |||||||6.166.14156):0.049009):0.024465)。图6.16是利用拓扑方法构造的系统发生树，用Phylip格式表示如下：0.000500,Human:0.000500)参考文献常青,,1998.分子进化研究中系统发生树的重建.生物多样性,6(1):55-，，姜成林.1997.构建微生物分子分类系统进化树的快速运算法与数据结构.微生，张.2000.分子钟及其存在的问题.人类学学报，19（2）：151-AdachiJ,HasegamaM.1996.MOLPHYVersion2.3.ProgramsforMolecularphylogeneticsbasedonumlikelihood.Tokyo:InstituteofStatusticalMathematics.AviseJC.1994.MolecularMarkers,NaturalandEvolution.NewYork:Chapman&BaxevanisAD,FrancisBF.1998.Bioinformatics:Apracticalguidetotheysisofgnesandproteins.JohnWiley&Sons,NewYork.BrocchieriL.2001.Phylogeneticinferencesfrommolecularsequences:reviewandcritique.Theor.Popul.Biol.59(1):27-40.BrownRH,etal.1972.Theaminoacidsequenceofcytochromecfromhelixaspersa.Biochem.J.,128:971-974.Cavalli-SforzaLL,EdwardsAW.1967.Phylogeneticysis.Modelsandestimationprocedures.AmJHumGenet1967May;19(3):Suppl19:233.DayWHE.1986.Computationalcomplexityofinferringphylogeniesfromdissimilaritymatrices.BulletinofMathematicalBiology,49:461–467.FengDF,DoolittleRF.1996.Progressivealignmentofaminoacidsequencesandconstructionofphylogenetictreesfromthem.MethodsEnzymoi.266:368-382.FelsensteinJ.1981.EvolutionarytreesfromDNAsequences:aumlikelihoodapproach.J.Mol.Evol.,17:368-376.FelsensteinJ.1988.Phylogeniesfrommolecularsequences:inferenceandreliability.AnnualsRev.Genetics,22:521-565.FelsensteinJ.1996.Inferringphylogeniesfromproteinsequencesbyparsimony,distance,andlikelihoodmethods.MethodsEnzymol1996;266:418-27.FelsensteinJ.1997.ternativeleast-approachtoinferringphylogeniesfrompairwisedistances.Syst.Biol.46：101-111.FelsensteinJ.1998.InferringPhylogenies.ASUWPublishing,Seattle,FelsensteinJ.1998.Phylogeny,II:CompatibilityandumFitchWM,MargoliashE.1967.Constructionofphylogenetictrees.Science,155:279-FitchWM,MargoliashE.1967.Amethodforestimatingthenumberofinvariantaminoacidcodingpositionsinageneusingcytochromecasamodelcase.BiochemGenet.,1(1):65-71.FitchWM.1976.Towarddefiningthecourseofevolution:minimumchangeforaspecifiedtreetopology.SystematicZoology,20:406–416.HarverPH,PagelMD.1991.ThecomparativeMethodinEvolutionaryBiology.OxfordUniversitypress.HeinJ.1990.Unifiedapproachtoalignmentandphylogenies.MethodsEnzymol.183：trees.MathematicalBiosciences,60:133–142.HermannB,HummelS,editors.1994.AncientDNA.Springer-Verlag,NewHuelsenbeckJP,RannalaB.1997.Phylogeneticmethodscomeofage:testinghypothesesinanevolutionarycontext.Science1997Apr11;276(5310):227-32.KimuraM.(ed.)1983.Theneutraltheoryofmolecularevolution.CambridgeUniversityKishinoH,HasegawaM.1989.Evaluationofthe umlikelihoodestimateoftheevolutionarytreetopologies.Comput.Appl.Biosci.，10：189-191.LiW,GraurD.1991.Fundamentalsofmolecularevolution.SinauerAssociates,106-LiWH.1997.MolecularEvolution.Sunderland,MA:SinauerLockhartPJ,SteelMA,HendyMD,PennyD.1994.Recoveringevolutionarytreesunderamorerealisticmodelofsequenceevolution.Mol.Biol.Evol.，11：605-612.MaddisonDR.1991.Thediscoveryandimportanceofm