内蒙科大生物信息学课件第6章 分子系统发生分析

第六章分子系统发生分析2022/9/10BIOINFORMATICS1本章提要：本章首先介绍了分子系统发生分析的基本内容，几种分子进化模型。重点学习分子系统树的构建方法，最后介绍分子系统发生树的分析软件。2022/9/10BIOINFORMATICS2

系统发生（或种系发生、系统发育，phylogeny)是指生物形成或进化的历史。系统发生学(phylogenetics)研究物种之间的进化关系，其基本思想是比较物种的特征，并认为特征相似的物种在遗传学上接近。2022/9/10BIOINFORMATICS3

系统发生研究的结果往往以系统发生树（phylogenetictree）表示，用它描述物种之间的进化关系。通过对生物学数据的建模提取特征，进而比较这些特征，研究生物形成或进化的历史。在分子水平上进行系统发生分析具有许多优势，所得到的结果更加科学、可靠。分子系统发生分析主要分成3个步骤：分子序列或特征数据的分析，系统发生树的构造以及结果的检验。2022/9/10BIOINFORMATICS4§6.1分子系统发生与系统发生树6.1.1分子系统发生分析系统发生学是进化生物学的一个重要研究领域，系统发生分析早在达尔文时期就已经开始。从那时起，科学家们就开始寻找物种的源头，分析物种之间的进化关系，给各个物种分门别类。经典系统发生学研究所涉及的特征主要是生物表型(phenotype)特征，所谓的表型特征主要指形态学的（结构的）特征，如生物体的大小、颜色、触角个数，也包括某些生理的、生化的以及行为习性的特征。2022/9/10BIOINFORMATICS5通过表型比较来推断生物体的基因型(genotype)，研究物种之间的进化关系。但是，利用表型特征是有局限性的。有时候关系很远的物种也能进化出相似的表型，这是由称为趋同进化(convergentevolution)的过程造成的。2022/9/10BIOINFORMATICS6例如，如果一个生物学家按照生物体是否有眼睛来构建进化树，那么他可能将人类、两翼昆虫和软体动物放在同一个进化组中，因为它们都有光探测器官。在这个例子中，很明显这3种生物体并不具有密切的关系，在其他特征上有天壤之别，就是它们眼睛的构造也大相径庭。这说明表型有时候会误导我们，表型相似并不总是反映基因相似。2022/9/10BIOINFORMATICS7用表型来判定进化关系的另一个问题是，对于许多生物体很难检测到可用来进行比较的表型特征。例如，即使用显微镜检查，也难以发现细菌的明显特性。当我们试图比较关系较远的生物体的时候，第三个问题又出现了，即什么样的表型特征能用来比较呢？例如，分析细菌、蠕虫和哺乳动物，它们之间的共同特征实在是少之又少。2022/9/10BIOINFORMATICS8随着人们对生物的认识从宏观发展到微观，科学家对物种分类的依据也从宏观上的形态发展到了微观上的分子，并且有了突破性的进展，系统发生分析进入分子层次。科学家认为，现今世界上存在的核酸和蛋白质分子都是从共同的祖先经过不断的进化而形成的，作为生物遗传物质的核酸和作为生命机器的蛋白质分子中存在着关于生物进化的信息，可用于系统发生关系的研究。2022/9/10BIOINFORMATICS9在分子水平上进行分析具有许多表型分析所没有的优势，所得到的结果更加科学、可靠。分子系统发生分析直接利用从核酸序列或蛋白质分子提取的信息，作为物种的特征，通过比较生物分子序列，分析序列之间的关系，构造系统发生树，进而阐明各个物种的进化关系。当然，这些分子不仅在序列上保留着进化的痕迹，它们的结构也保留着进化的痕迹。2022/9/10BIOINFORMATICS1020世纪60年代，蛋白质测序成为可能；20世纪70年代，研究者开始能够获得基因组信息，特别是DNA序列。蛋白质序列和DNA序列为分子系统发生分析提供了可靠的数据。2022/9/10BIOINFORMATICS11在现代分子进化研究中，根据现有生物基因或物种多样性来重建生物的进化史是一个非常重要的问题。一个可靠的系统发生的推断，将揭示出有关生物进化过程的顺序，有助于我们了解生物进化的历史和进化机制。根据核酸和蛋白质的序列信息，可以推断物种之间的系统发生关系。其基本原理非常简单，从一条序列转变为另一条序列所需要的变换越多，那么，这两条序列的相关性就越小，从共同祖先分歧的时间就越早，进化距离就越大；相反，两个序列越相似，那么它们之间的进化距离就可能越小。2022/9/10BIOINFORMATICS12为了便于分析，一般假设序列变化的速率相对恒定。关于地球上现代人起源的研究是一个典型而有趣的例子，科学家分析了取自世界不同地区许多人的线粒体DNA，分析结果表明：所有现代人都是一个非洲女性的后代。线粒体DNA非常适合于系统发生分析，因为线粒体DNA从母体完全传到子代，不与父代DNA重组。2022/9/10BIOINFORMATICS13由于DNA分子非常稳定，所以我们既可以通过DNA分析活着的生物，也可以分析死去的生物，甚至分析已经绝种的生物。当然，用细胞核基因来研究系统发生关系时，遇到的一个严重的问题是，基因常常会被复制，导致在个体基因组中，一个基因可能有若干个拷贝。在进化过程中，这些拷贝各自演变，形成两个或更多的相似基因。在对不同物种的基因进行比较时，如果选择这类基因，其分析结果的可靠性将存在问题。2022/9/10BIOINFORMATICS14所有的生物都可以追溯到共同的祖先，生物的产生和分化就像树一样地生长、分叉，以树的形式来表示生物之间的进化关系是非常自然的事。可以用树中的各个分支点代表一类生物起源的相对时间，两个分支点靠得越近，则对应的两群生物进化关系越密切。2022/9/10BIOINFORMATICS15系统发生分析一般是建立在分子钟(molecularclock)基础上的。生物随着时间的推进而演化，进化的速率被视为进化研究中的基本问题之一。进化速率就是在某一段时间内的遗传改变量。分子进化速率相关的分子钟的概念源于对蛋白质序列的研究。在长期的进化过程中，有着相似功能约束位点的分子进化速率则几乎完全一致。2022/9/10BIOINFORMATICS1620世纪60年代最早由EmileZuckerkandl和LinusPaining所作的蛋白质序列的比较研究表明，蛋白质同系物的替换率就算过了千百万年也能保持恒定，因此他们将氨基酸的变异积累比作分子钟。科学家们在比较几种动物的血红蛋白、细胞色素c的序列后注意到：这些蛋白质的氨基酸取代速率在不同的种系间大致相同，即分子水平的进化存在恒速现象。2022/9/10BIOINFORMATICS17分子时钟在不同的蛋白质中运行的速率是不同的，但是，两个蛋白质同系物的差异始终和它们独立分化的时间成正比。两序列间稳定的变异速率，不仅有助于确定物种间的系统发生关系，而且能够像利用放射性衰变考察地质年代那样，准确测定序列分化发展的时间。不同物种间的蛋白质氨基酸序列差异随着分歧时间的加大而增加，而DNA序列也存在这种规律。2022/9/10BIOINFORMATICS18Kimura进一步提出了具体的分子进化观点：对于各物种的每个蛋白质，如果用每个位点每年发生的氨基酸替换次数作为衡量分子进化的速率，则该速率是大致恒定的；功能上次要的分子的进化速率比功能重要的分子进化速率快；对现有分子结构或者功能破坏小的氨基酸替换比破坏力大的氨基酸替换发生得更加频繁。2022/9/10BIOINFORMATICS19尽管以上的分析给我们带来了很多希望，但是，Zuckerkandl和Paining的分子时钟假说还是有争议的。经典进化学家们认为形态的进化不够稳定，这与分子以稳定的速度变异不一致。关于分化时间也有不同意见，这些意见对这个假说的核心即进化率是稳定的表示质疑。2022/9/10BIOINFORMATICS206.1.2系统发生树一般来说，系统发生树是一种二叉树。所谓树，实际上是一个无向非循环图。系统发生树由一系列节点（nodes）和分支（branches）组成，其中每个节点代表一个分类单元（物种或序列），而节点之间的连线代表物种之间的进化关系。2022/9/10BIOINFORMATICS21系统发生树有许多形式：可能是有根树(rootedtree)，也可能是无根树（unrootedtree)；可能是一般的树，也可能是二叉树；可能是有权值的树（或标度树，scaledtree，树中标明分支的长度），也可能是无权值树（或非标度树，unscaledtree）。在有权值的树中，分支的长度（或权值）一般与分类单元之间的变化成正比，它是关于生物进化时间或者遗传距离的一种度量形式。一般假设存在一个分子钟，进化的速率恒定。2022/9/10BIOINFORMATICS22

类群1、祖先类群(ancestralgroup)

如果一个类群(物种)至少有一个子裔类群，这个原始的类群就称为祖先类群。2、单系类群(monophyleticgroup)

包含一个祖先类群所有子裔的群组称为单系类群，其成员间存在共同祖先关系。2022/9/10BIOINFORMATICS233、并系类群(paraphyleticgroup)和复系类群(polyphyleticgroup)

不满足单系类群要求，各成员间又具有共同祖先特征的群组称为并系类群；各成员既不具有共同衍生特征也不具有共同祖先特征，只具有同型特征的分类群组称为复系类群。2022/9/10BIOINFORMATICS244、内类群(ingroup)和外类群(outgroup)

一项研究所涉及的某一特定类群可称为内类群；不包括在内类群中又与之有一定关系的类群可称为外类群。外类群的建立是为了比较内类群中各成员的特征差异。2022/9/10BIOINFORMATICS255、姊妹群(sistergroup)

与某一类群在谱系关系上最为密切的类群称为姊妹群。姊妹群是单系类群的一种常见类型。2022/9/10BIOINFORMATICS26

图6-1示出6个分类群(A-F)进行不同划分所产生的单系、并系和复系类群的例子。图6-1(a)中单系类群为：{A,B}，{E，D，F｝、{C,D,E,F｝、{A,B,C,D,E,F｝图6-1(b)中并系类群为:{C,D,E}、{B,C,D,E,F}等图6-1(c)中复系类群为：{B,D}、{B,D,E,F}图6-1单系类群(a)、并系类群(b)、复系类群(c)

图6-2示出了内类群、外类群和姊妹群的例子。图中G1={A，B，C}和G2={S，G1}可视为内类群；O1，O2为G2的外类群；S为G1的姊妹群图6-2内类群、外类群和姊妹群2022/9/10BIOINFORMATICS29

有根树和无根树

系统树是表达类群(或序列)间系统发育关系的一种树状图，可划分为以下几种类型。以外类群作为树根的系统树称为有根树，没有外类群作为树根的系统树称为无根树。有根树中有一个惟一的根节点，代表所有其他节点的共同祖先，这样的树能够反映进化层次，从根节点历经进化到任何其他节点只有惟一的路径。无根树没有层次结构，无根树只说明了节点之间的关系，没有关于进化发生方向的信息。2022/9/10BIOINFORMATICS30当分类单位数目为n(n>2)时，全部可能的有根树(二歧分支)的数目为：

当n≥3时，无根树(二歧分支)的数目为：(6-2)(6-1)2022/9/10BIOINFORMATICS31

随着n的增加，可能的有根树和无根树的数目迅速增加。表6-1中列出了一些n值下对应的有根树和无根树的数目。例如，当n＝10时，NR＝34459425，NU＝2027025。当n≥15时，可能的系统树的数目变得非常巨大。然而，这其中只有一棵树代表了真实的进化关系，我们的任务就是找出这棵反映真实进化关系的树。数据数目n有根树数目NR无根树数目NU211331415351051510344594252207025152134580467678757905853580625表6-1对不同的n，可能的有根树和无根树的数目AABBCCDDEE时间(a)(b)R

图6-3示出了有根树和无根树的例子，其中R为有根树(a)的树根。图6-3有根树(a)和无根树(b)2022/9/10BIOINFORMATICS34

标度树枝(scaledbranch)和非标度树枝(unscaledbranch)在图6-4示出的两个系统树中，(a)中树枝的长度代表了性状状态变异的数量，称为标度树枝系统树；而(b)中各树枝的长度并不表示性状状态变异的数量，但所有节点的位置仍与分化时间（divergencetime）相对应，可称为非标度树枝系统树。图6-4标度树枝(a)和非标度树枝(b)系统树2022/9/10BIOINFORMATICS36

基因树(genetree)和物种树(speciestree)基因树是由一个基因所构建的系统树，这种树代表的仅仅是单个基因的进化历史，而不是它所在物种的进化历史。物种树则表达了某一特定类群的进化路径(evolutionarypathway)。虽然基因树和物种树同属系统树，但两者至少在两个方面存在差异：2022/9/10BIOINFORMATICS37

①从两个不同的物种中获取的两个基因，其分化可能早于物种的分化；②基因树的拓扑结构可能与物种树不一致，因为两个或多个基因树之间很可能存在不一致之处。如何将由多个基因或基因组建立的基因树综合成为一个物种树，是分子系统学目前所面临的一个难题。2022/9/10BIOINFORMATICS38

图6-5示出了物种树和基因树间3种可能的关系。图中圆角矩形(代表物种形成)和实心圆点(代表基因分化)分别表示物种树和基因树。在(a)和(b)中，物种树的拓扑结构与基因树一致。其中(a)中基因间的分化时间大致等于物种间的分化时间；(b)中基因X和Y的分化时间大大早于物种间的分化时间；(c)中基因树与物种树的拓扑结构不同。对于中性位点而言，上述3种情况的概率均可通过物种分化时间(T=t1-t2)和有效种群规模(effectivepopulationsize)N来计算，分别标注在图6-5的树下。图6-5物种树（矩形）和基因树（实心圆点）2022/9/10BIOINFORMATICS406.1.3距离和特征用于构建系统发生树的分子数据分成两类：

(1)距离(distances)数据，常用距离矩阵描述，表示两个数据集之间所有的两两差异；

(2)特征(characters)数据，表示分子所具有的特征。2022/9/10BIOINFORMATICS41分子系统发生分析的目的是探讨物种之间的进化关系，其分析对象往往是一组同源的序列。这些序列取自于不同生物基因组的共同位点。序列比对是进行同源分析的一种基本手段，是进行系统发生分析的基础，一般采用基于两两比对渐进的多重序列比对方法，如ClustalW程序。通过序列的比对，可以分析序列之间的差异，计算序列之间的距离。2022/9/10BIOINFORMATICS42无论是DNA序列，还是蛋白质序列，都是由特定字母表中的字符组成的。计算序列之间距离的一个前提条件是要有一个字符替换模型，替换模型影响序列多重比对的结果，影响系统发生树的构建结果。在具体的分析过程中，需要选择一个合理的字符替换模型，参见第5章的各种打分模型或代价、距离模型。2022/9/10BIOINFORMATICS43距离（或者相似度）是反映序列之间关系的一种度量，是建立系统发生树时所常用的一类数据。在计算距离之前，首先进行序列比对，然后累加每个比对位置的得分。可以应用第5章介绍的关于序列比较的方法，直接计算序列之间的距离。2022/9/10BIOINFORMATICS44如果在进行序列比较时使用的是打分函数或相似性度量函数，则需要将相似度（或者得分）转换成距离。令S(i,j)是序列i和序列j各个比对位置得分的加权和，一种归一化的距离计算公式为(6-3)2022/9/10BIOINFORMATICS45其中，Sr(i,j)是序列i和j随机化之后的比对得分的加权和，Smax(i,j)是两条序列所有可能的比对的最大值（当两条序列相同时，取最大值）。两个序列归一化距离的值处于0和1之间，当两个序列完全一致时，距离为0；当两个序列差异很大时，距离接近于1。如果在上式中令Sr(i,j)=0，则：

(6-4)2022/9/10BIOINFORMATICS46为了适合于处理相似性较小的序列，可以进一步修改距离计算公式：

(6-5)2022/9/10BIOINFORMATICS47序列比对得分的加权和可以根据常用的得分矩阵获得，如果待处理的序列是蛋白质，则用PAM矩阵BLOSUM矩阵等；如果待处理的序列是DNA或者RNA，则用等价矩阵、核苷酸转换一颠换矩阵或者其他具有非对称置换频率的矩阵。2022/9/10BIOINFORMATICS48距离是系统发生分析时所使用的一类数据，另一类数据就是所谓的离散特征数据。离散特征数据可分为二态特征与多态特征。二态的离散特征只有两种可能的状况，即具有与不具有某种特征，通常用“0"或“1”表示。例如，DNA序列上的某个位置如果是剪切位点，其特征值为1，否则为0。2022/9/10BIOINFORMATICS49多态离散特征具有两种以上可能的状态，如核酸的序列信息，对序列中某一位置来说，其可能的碱基有A、T、G、C共4种。可以将特征数据转换为距离数据。如果建立所有可能状态之间相似性的度量，特征数据就很容易被转换成距离数据。2022/9/10BIOINFORMATICS50§6.2分子进化模型与序列分歧度计算6.2.1核苷酸序列进化Jukes-Cantor单参数模型

DNA序列进化的一个基本过程就是核苷酸随时间而变化(置换)。核苷酸置换模型可以用矩阵表示。图6-6示出了Jukes和Cantor(1969)单参数模型的基本假定，即核苷酸在每个方向上的置换率均为α。图6-6Jukes-Cantor单参数模型2022/9/10BIOINFORMATICS52对应的核苷酸置换率矩阵M为：

(6-6)2022/9/10BIOINFORMATICS53

利用该模型，可以计算序列间核苷酸置换概率。设某一位点上的碱基在起始时间(t=0)为A，在时间t=1时，同一位点上的碱基仍为A的概率为：

PA(1)=1-3α(6-7)

t=2时，该位点碱基为A的概率：

PA(2)=(1-3α)PA(1)+α[1-PA(1)]

(6-8)2022/9/10BIOINFORMATICS54

时间为t+1时，概率为：

PA(t+1)=(1-3α)PA(t)+α[1-PA(t)](6-9)定义时间t+1和t的概率差为ΔPA(t)有：

ΔPA(t)=PA(t+1)-PA(t)=-3αPA(t)+α[1-PA(t)]=-4αPA(t)+α

(6-10)2022/9/10BIOINFORMATICS55

上述公式适用于离散时间过程。对连续时间过程，可转换为：令PA(0)=1，有：PA(t)=1/4+3/4e-4αt(6-12)(6-11)Kimura的两参数模型

Kimura的两参数模型将转换(transition)和颠换(transversion)置换率分别设为α和β，如图6-7所示。图6-7Kimura两参数模型2022/9/10BIOINFORMATICS57对应的矩阵M为：(6-13)2022/9/10BIOINFORMATICS58

Jukes-Cantor单参数模型(JC)和Kimura两参数模型(KZP)在分子进化研究中应用较为广泛。其他模型有Kimura三参数模型(K3ST)、Felsenstein“等输入”模型(F81)、Lanave等“广义时间可反转”模型(GTR)、Hasegawa等五参数模型(HKY85)、Kishino和Hasegawa广义模型(SYM)以及Tamura和Nei六参数模型(TrN)等。近来，Yang与Adachi和Hasegawa还提出了广义可反转马尔可夫模型，等等。2022/9/10BIOINFORMATICS596.2.2蛋白质编码序列进化

Miyata和Yasunaga(1980)最早提出了估算两个蛋白质序列间置换数目的方法。在这些方法中，必须将同义置换(synonymous)和非同义置换(nonsynonymous)分开考虑。起始和终止密码子应排除在外，因为它们几乎不随时间变化。2022/9/10BIOINFORMATICS60

在研究蛋白质序列进化的工作中，最基本的参数计算步骤包括以下几步:首先，将核苷酸序列的位点分为非简并的(nondegenerate,如果在该位点所有可能的变化都是非同义的)、两重简并的(twofolddegenerate,如果3个可能变化中的一个是同义的）和四重简并的(fourfolddegenerate,如果所有可能的变化都是同义的)。计算这3种情况的数目，记为Li(i=0,2,4)。2022/9/10BIOINFORMATICS61

其次，比较两个密码子之间的同义与非同义变化。对于只有一个核苷酸差异的两个密码子来说，差异值容易计算；对于不止一个核苷酸差异的两个密码子比较，则必须考虑所有可能的进化途径。例如，AAT(Asn)和ACG(Thr)间的两种可能的途径是：

途径1：AAT(Asn)→ACT(Thr)→ACG(Thr)

途径2：AAT(Asn)→AAG(Lys)→ACG(Thr)2022/9/10BIOINFORMATICS62

途径1中同义与非同义变化各1次，而途径2中有2次非同义变化。假定两种途径的可能性相同(非加权)，同义差异值(MA)为(1+0)/2=0.5，非同义差异值(MA)为(1+2)/2=1.5。不少学者已对加权方法进行了讨论。2022/9/10BIOINFORMATICS63

最后，可以将核苷酸置换分为转换和颠换两种，差异值分别记为Si和Vi(i=0,2,4)。在两重简并的位点中，转换是同义的而颠换是非同义的，这在哺乳动物线粒体编码中是没有例外的。在通用编码中有两个例外：精氨酸密码子(CGA,CGG,AGA和AGG）的第一个位置，异亮氨酸(AUU,AUC和AUA）的最后一个位置。在这两种例外情况中，所有的同义变化包括在S2中，而所有的非同义变化包括在V2中。2022/9/10BIOINFORMATICS64

在实际应用中，常见的氨基酸置换有Dayhoff模型、Jones-Taylor-Thomton模型、mtREV模型等。近来，有关氨基酸置换模型的建立及应用是分子进化研究中最为活跃的领域之一。2022/9/10BIOINFORMATICS656.2.3核苷酸序列分歧度

DNA序列间的分歧度(sequencedivergence)是一种相异性指数，可通过序列成对比较获得碱基差异值，然后应用序列进化模型来校正（Li，1997）。例如，图6-8给出了两个同源序列间核苷酸置换的例子，其中存在12个置换但仅有3个（位点2,位点5和位点7)能被检测出来，用作计算碱基差异值。图6-8同源序列间的核苷酸置换2022/9/10BIOINFORMATICS67

设序列长度为L，序列间的差异值为N，P＝N/L。应用Jukes-Cantor单参数模型，计算K值的公式为：

K＝-3/4ln(1-4P/3)(6-14)当L足够大时，K值的取样误差为：(6-15)2022/9/10BIOINFORMATICS68

设序列间转换置换的差异值为N1，P＝N1/L；颠换置换的差异值为N2，Q＝N2/L。应用Kimura两参数模型，计算K值的公式为：K＝1/2ln(a)+1/4ln(b)(6-16)这里，a=1/(1-2P-Q)，b=1/(1-2Q)。当L足够大时，K值的取样误差为：

V(K)=[a2P+c2Q-(aP＋cQ)2]/L(6-17)这里c=(a＋b)/2。2022/9/10BIOINFORMATICS696.2.4蛋白质编码序列分歧度

设KS为两个序列间同义变化的分歧度，KA为非同义变化的分歧度，应用Jukes-Cantor单参数模型，可以计算：(6-18)(6-19)2022/9/10BIOINFORMATICS70

这里，Ms为两个序列间密码子比较获得的同义差异值，Ns为同义位点的数目；MA为非同义差异值，NA为非同义位点的数目。

令Pi=Si/Li，Qi=Vi/Li，应用Kimura两参数模型，可以计算不同简并情形下的分歧度：Ki＝Ai+Bi

(6-20)

Ai=ln(ai)/2-ln(bi)/4(6-21)Bi＝ln(bi)/2(6-22)2022/9/10BIOINFORMATICS71其中：ai=1/(1-2Pi-Qi)，bi=1/(1-2Qi)。方差估计为：(6-23)其中：ci＝(ai-bi)/2

(6-24)2022/9/10BIOINFORMATICS72

两个序列间同义和非同义变化的分歧度，可以分别计算:(6-25)(6-26)2022/9/10BIOINFORMATICS73其中个参数的定义同上。方差估计为：(6-27)(6-28)2022/9/10BIOINFORMATICS74

上述公式中未考虑转换与颠换比率(即加权问题)。如果用加权平均(L2A2+L4A4)/(L2+L4)作为两重简并和四重简并情形中转换与颠换比率的估计值，可以计算：(6-29)(6-30)2022/9/10BIOINFORMATICS75方差估计为：(6-31)(6-32)2022/9/10BIOINFORMATICS76§6.3分子系统树的构建6.3.1分子系统发生分析过程分子系统发生分析主要分成3个步骤：(1)分子序列或特征数据的分析；(2)系统发生树的构造；(3)结果的检验。2022/9/10BIOINFORMATICS77其中，第(1)步的作用是通过分析，产生距离或特征数据，为建立系统发生树提供依据。系统发生树的构建方法有很多种。根据所处理数据的类型，可以将系统发生树的构建方法大体上分为两大类。一类是基于距离的构建方法，利用所有物种或分类单元间的进化距离，依据一定的原则及算法构建系统发生树。其基本思路是列出所有可能的序列对，计算序列之间的遗传距离，选出相似程度比较大或非常相关的序列对，利用遗传距离预测进化关系。2022/9/10BIOINFORMATICS78这类方法有：非加权分组平均法(unweightedpairgroupmethodwitharithmeticmeans，UPGMA)邻接法(neighborjoiningmethod，NJ)Fitch-Margoliash（FM）法最小进化方法（minimumevolution，ME)等。2022/9/10BIOINFORMATICS79另一类方法是基于离散特征的构建方法，利用的是具有离散特征状态的数据，如DNA序列中的特定位点的核苷酸。建树时，着重分析分类单位或序列间每个特征（如核苷酸位点）的进化关系等。属于这一类的方法有最大简约法(maximumparsimonymethod，MP)、最大似然法(maximumlikelihoodmethod，ML)、进化简约法（evolutionaryparsimonymethod，EP）、相容性方法（compatibility）等。2022/9/10BIOINFORMATICS806.3.2距离矩阵法(distancematrixmethod)这里介绍的距离矩阵法包括两种方法：UPGMA法和邻接法(neighbor-joiningmethod)。在使用这两种方法前都必须获得一个对称距离矩阵(m阶方阵)D={dij}m×n，其中m为OTU数目。距离系数的公式很多。例如，Nei(1972)的遗传距离系数适用于限制性内切酶和同功酶数据，Jukes-Cantor单参数距离系数和Kimura两参数模型距离系数则广泛用于各种序列数据。2022/9/10BIOINFORMATICS81UPGMA法

UPGMA法是使用算术平均的不加权对群法的缩写，也称为类平均法，是目前聚类分析中使用得最多的一种聚合策略。设两个已聚合过的类群OTUp和OTUq中分别包含了np和nq个原始类群，再将OTUp和OTUq聚合后组成的OTUr与其他OTUi间的距离dr,i可以用以下公式计算：

(6-33)2022/9/10BIOINFORMATICS82ABCDEA0B80C480D6860E84880

例6-1对下面距离矩阵，用UPGMA法构建系统发生树2022/9/10BIOINFORMATICS83解：首先，根据最小距离值4(对角线元素0除外)，将A和C聚合为AC。AC与B、D、E间的距离分别为：

2022/9/10BIOINFORMATICS84ACBDEAC0B80D680E8480新的距离矩阵为：2022/9/10BIOINFORMATICS85根据最小距离值4(对角线元素0除外)，将B和E聚合为BE。BE与AC、D间的距离分别为：

2022/9/10BIOINFORMATICS86ACBEDAC0BE80D680新的距离矩阵为：2022/9/10BIOINFORMATICS87根据最小距离值6(对角线元素0除外)，将AC和D聚合为ACD。ACD与BE间的距离为：

根据以上结果，可以构建一个系统树

图6-9

用UPGMA法构建系统树2022/9/10BIOINFORMATICS89邻接法邻接法的原始思路由Saitou和Nei(1987)提出，后由Studier和Keppler(1988)修正。该方法的关键步骤是计算发散系数ri，二是生成一个速率校正距离阵(ratecorrecteddistancematrix）M＝{Mij}m×n。公式分别如下：2022/9/10BIOINFORMATICS90

Mij=dij-(ri+rj)/(N-2)(6-35)(6-34)这里，N是终端节点的数目。对所有的i和j，设j>i，然后找出Mij最小值所对应的i和j。最后，根据每一步骤的结果绘制系统树。例6-2下面是一个5个分类群5SrRNA的例子。分类群及缩写为：①Bacillussubtilis(Bsu)②B.stearothermophilus(Bst)③Lactobacillusviridescens(Lvi)④Acholeplasmamodicum(Amo)⑤Micrococcusluteus(Mlu))。距离系数半对角矩阵为：BsuBstLviAmoMluriri/(N-2)Bsu—0.17150.21470.30910.23260.92790.3093Bst-0.4766—0.29910.33990.20581.01630.3388Lvi-0.4905-0.4356—0.27950.39431.18760.3959Amo-0.4527-0.4514-0.5689—0.42891.35740.4525Mlu-0.4972-0.5535-0.4221-0.4441—1.26160.4205解：计算ri和ri/(N-2)值，以及Mij，组成表6-2。*N=5表6-2邻接法中间结果2022/9/10BIOINFORMATICS93

表6-2上半角矩阵为距离系数矩阵，由此计算ri和ri/(N-2)的值；下半角为速率校正距离矩阵M，其中MLvi,Amo＝-0.5689为最小值。设新节点(邻接节点)为u1，则：

SLvi,u1=SLvi,Amo/2-(rLvi-rAmo)/(2×(N-2))=0.2795/2+(0.3959-0.4525)/2=0.1114SAmo,u1=0.2795-0.1114=0.1681

Sbsu,u1=(Sbsu,Lvi+Sbsu,Amo-SLvi,Amo)/2=0.1222Sbst,u1=(Sbst,Lvi+Sbst,Amo-SLvi,Amo)/2=0.1798SMlu,u1=(SMlu,Lvi+SMlu,Amo-SLvi,Amo)/2=0.2719BsuBstMluu1riri/(N-2)Bsu—0.17150.23160.12220.52630.2631Bst-0.3701—0.20580.17980.55710.2785Mlu-0.3856-0.4278—0.27190.71030.3551u1-0.4278-0.3856-0.3701—0.57390.2869删去已连接的类群Amo和Lvi后，计算结果见表6-3

表6-3邻接法中间结果(N=4)2022/9/10BIOINFORMATICS95从表6-3中可以发现MBsu,u1=-0.4278为最小值。对节点u2,有：

SBsu,u2=0.1222/2+(0.2631-0.2869)/2=0.0492Su1,u2=0.1222-0.0492=0.0730SBst,u2=(SBst,u1+SBst,Bsu-SBsu,u1)/2=0.1146SMlu,u2=(SMlu,u1+SMlu,Bsu-SBsu,u1)/2=0.1912BstMluu2riri/(N-2)Bst—0.20580.11460.32040.3204Mlu-0.5116—0.19120.39700.3970u2-0.5116-0.5116—0.30580.3058删去Bsu和u2后的计算结果见表6-4表6-4邻接法中间结果(N=3)MBst,u2=-0.5116为3个最小值之一。对节点u3有

SBst,u3=0.1146/2+(0.3204-0.3058)/2＝0.0646Su2,u3=0.1146-0.0646＝0.0500最后，SMlu,u3=0.1412根据以上结果，绘制无根系统树。树枝长度用上述计算结果标注。图6-10用邻接法构建系统2022/9/10BIOINFORMATICS98

ABCDEFA054768B0710911C0768D059E08F0例6-3对下面距离矩阵，用邻近归并法构建系统发生树。ABCDEFriri/（N-2）A054768307.5B-1307109114210.5C-11.5-11.50768328D-10-10-10.5059389.5E-10-10-10.5-1308348.5F-10.5-10.5-11-11.5-11.504411解：发散系数ri，速率校正距离阵M＝{Mij}m×n。公式分别如下：

*N=62022/9/10BIOINFORMATICS100

其中MA,B＝-13为最小值。设新节点(邻接节点)为AB,SA,AB=SA,B/2+(rA-rB)/(2×(N-2))=5/2+(7.5-10.5)/2=1SB,AB=5-1=4SC,AB=(SC,A+SC,B-SA,B)/2=3SD,AB=(SD,A+SD,B-SA,B)/2=6SE,AB=(SE,A+SE,B-SA,B)/2=5SF,AB=(SF,A+SF,B-SA,B)/2=7ABCDEFriri/（N-2）AB03657217C-120768248D-10-10059279E-10-10-1208248F-10.67-10.67-10.67-10.6703210.67删去已连接的类群A和B后，计算结果见下表（N=5）2022/9/10BIOINFORMATICS102

其中MAB,C＝-12为最小值。设新节点(邻接节点)为ABCSAB,ABC=SAB,C/2+(rAB-rC)/(2×(N-2))

=3/2+(7-8)/2=1

SC,ABC=3-1=2

SD,ABC=(SD,AB+SD,C-SAB,C)/2=5

SE,ABC=(SE,AB+SE,C-SAB,C)/2=4

SF,ABC=(SF,AB+SF,C-SAB,C)/2=6ABCDEFriri/（N-2）ABC0546157.5D-12059199.5E-12-1308178.5F-13-12-1202311.5删去已连接的类群AB和C后，计算结果见下表（N=4）2022/9/10BIOINFORMATICS104

其中MD,E＝-13为最小值。设新节点(邻接节点)为DESD,DE=SD,E/2+(rD-rE)/(2×(N-2))=5/2+(9.5-8.5)/2=3SE,DE=5-3=2SABC,DE=(SABC,D+SABC,E-SD,E)/2=2SF,DE=(SF,D+SF,E-SD,E)/2=62022/9/10BIOINFORMATICS105ABCDEFriri/（N-2）ABC02688DE-140688F-14-1401212删去已连接的类群D和E后，计算结果见下表（N=3）2022/9/10BIOINFORMATICS106

设新节点(邻接节点)为ABCDE,SABC,ABCDE=SABC,DE/2+(rABC-rDE)/(2×(N-2))

=2/2+(8-8)/2=1

SE,DE=2-1=1SF,ABCDE=(SF,ABC+SF,DE-SABC,DE)/2=51根据以上结果，获得一个系统树 图6-11邻近归并法构建系统树2022/9/10BIOINFORMATICS108Fitch-Margoliash法Fitch和Margoliash(1976)提出了一种距离构树法，应用较广。结合例6-4，说明该方法的基本步骤如下。例6-4设A~D4个类群间的距离矩阵为：2022/9/10BIOINFORMATICS1091、从距离矩阵中挑选距离最近的两个类群(A,C)聚合，聚合节点为1(见下图)；2、在图中加入第三个类群(暂定为B)。标注A与1之间的距离为a，B与1之间的距离为b，C与1之间的距离为c；定义A与B之间的距离为A到除C之外的所有其他类群的平均距离，标注为dAB。图6-12Fitch-Margoliash法构建系统树（1）2022/9/10BIOINFORMATICS111

显然：dAB=a+b=(8+12)/2=10同理：dCB=c+b=(9+11)/2=10加上：dAC=a+c=7计算获得：

A=3.5c=3.5b=6.5其中a和c为最终结果，b为中间结果。2022/9/10BIOINFORMATICS112用UPGMA法中公式计算节点1到B和D的距离：d1,B=8.5d1,D=11.52022/9/10BIOINFORMATICS113构建新的距离矩阵：2022/9/10BIOINFORMATICS1143.继续加入第四个类群，重复步骤2。计算获得：

d＝8.5b＝5.5l'＝3

4.最后，获得一个无根树。注意，中间节点1和2之间的距离调整为d12=|l’-3.5|=0.5图6-13Fitch-Margoliash法构建系统树（2）2022/9/10BIOINFORMATICS1166.3.3基于特征的构树法基于特征的系统发生分析要解决的一般问题是：给定n个分类单元，m个用以描述分类单元的特征，以及每个分类单元所对应的特征值，构建一棵系统发生树，使得某个目标函数最大。一般情况下输入的是一个nm的矩阵M，其中，Mij代表第i个分类单元的第j个特征的取值。在构建系统发生树时，假设特征是相互独立的，即一个特征的变化不影响另一个特征。另外，还假设在进化过程中，两个物种分叉后独立进化，互不影响。2022/9/10BIOINFORMATICS117对于给定的条件，要在很多可能的树中找一棵最佳的树。在实际应用中，不可能穷尽搜索所有可能的树，必须按照一定的方法、一定的策略在较短的时间内得到比较好的结果。2022/9/10BIOINFORMATICS118

最大简约法(maximumparsimony，MP)最大简约法最早是基于形态特征分类的需要而发展起来的，具体的算法有许多版本，其中有些已被广泛地用于分子进化研究中，根据离散特征数据构建系统发生树。最大简约法的目标是构建一棵反映分类单元之间最小变化的系统发生树。最大简约法利用的只是对简约分析能提供信息的特征，对核苷酸序列，理论上每个位点均可构建3种可能的系统树。2022/9/10BIOINFORMATICS119然而，并非所有位点都能用于重建系统发育关系。能够用于构建最简约系统树的位点称为信息位点（informativesite）。例如，考虑如下4个序列：位点123456789序列1AAGAGTGCA序列2AGCCGTGCA序列3AGATATCCA序列4AGAGATCCG2022/9/10BIOINFORMATICS120位点1显然不是一个信息位点，因为所有序列在该位点上均为A，性状状态变化的数目为零。位点6和位点8的情形也相同。在位点2上，只有序列1为A，其余均为G。由于3种可能的树中都只有1个变化（从G到A），因而该位点也不能作为信息位点。位点3和位点4同样导致3种可能的树中性状状态变化相同（分别为2和3）（图6-14），不能视为信息位点。

在位点5上，树1中是1个变化，而树2和树3均为2个变化，该位点为一个信息位点。逐一讨论所有位点，可以确定位点5、位点7和位点9为信息位点。以位点5为例，最简约系统树为树1。2022/9/10BIOINFORMATICS121对于系统发生树最直观的计算代价方法就是沿着各个分支累加特征变化的数目，而所谓简约就是使代价最小。利用最大简约方法构建系统发生树，实际上是一个对给定分类单元所有可能的树进行比较的过程。针对某一个可能的树，首先对每个位点祖先序列利用最大简约方法构建系统发生树，实际上是一个对给定分类单元所有可能的树进行比较的过程。2022/9/10BIOINFORMATICS122针对某一个可能的树，首先对每个位点祖先序列的核苷酸组成做出推断，然后统计每个位点用来阐明差异的核苷酸最小替换数目。在整个树中，所有简约信息位点最小核苷酸替换数的总和称为树的长度或树的代价。通过比较所有可能的树，选择其中长度最小、代价最小的树作为最终的系统发生树，即最大简约树。2022/9/10BIOINFORMATICS123图6-14确定信息位点2022/9/10BIOINFORMATICS124随着序列数目和信息位点数目的增加，最简约系统树的数目将迅速增加，通常只能使用计算机程序来发现所有的最简约系统树。其算法步骤包括：

(1)确定所有的信息位点；

(2)对每个信息位点计算核苷酸变化最小数目，并对所有信息位点的最小数目求和；

(3)选取核苷酸变化总和最小的树作为最简约系统树。注意，对一个数据集可能产生多个最简约系统树。2022/9/10BIOINFORMATICS125

进化简约法

进化简约法也称为无变度法。该方法的具体步骤为：（1）、选择每4个序列为一组（2）、寻找两个具有嘌呤和两个具有嘧啶的位点（3）、考虑3种可能的树型组合（图6-15），分别称为树型X、Y和Z2022/9/10BIOINFORMATICS126（4）、找出序列A、B均为嘌呤或均为嘧啶的位点（C和D则相反），并计算支持位点的数目和相反树枝X顺序的总数，分别记为X+和X-；同样，用序列A、C和B、D计算Y+和Y-的总数，并用序列A、D和B、C计算Z+和Z-的总数。图6-154个序列的3种树型组合2022/9/10BIOINFORMATICS127（5）计算对树型X,Y和Z的净支持率；

X=X＋-X－

Y=Y＋-Y－

Z=Z＋-Z－（6）

用自由度f=1的χ2检验来估计统计显著性；

χx2＝X2/(X＋+X－)χy2＝Y2/(Y＋+Y－)χz2＝Z2/(Z＋+Z－)2022/9/10BIOINFORMATICS128例6-5对下列4个序列SA、SB、SC和SDSAAUCAGGCUUGCACUAACUGGSBAGGAGAAGUAAGGCCACUUCSCAGGUGUAAUCAGGGCAGAACSDAGGUAGCUUUUGCAACGAUA计算各参数：X＋＝3＋3＝6，X一＝0Y＋=0，Y一＝1Z十＝4，Z一＝2＋1＝32022/9/10BIOINFORMATICS129故：

X=X+-X-=6Y=Y+-Y-=-1Z=Z+-Z-=4-3=1

统计显著性检验：

χx2＝62/(6＋0)＝6χy2＝(-1)2/(1+0)=1χz2＝12/(4+3)=0.14在自由度为1的情况下，只有对χx2才具有显著性。因此，获得图6-16所示进化简约树。2022/9/10BIOINFORMATICS130图6-16进化简约法构建系统树例子2022/9/10BIOINFORMATICS1316.3.4其它方法

最大似然法（MaximumLikelihood，ML）最大似然法是一种重要的统计方法。早在1967年该方法就被用于基因频率数据分析。然而，由于计算复杂性等原因，该方法并未得到广泛的应用。直到20世纪80年代后，最大似然法才逐渐引起分子进化学家的重视。目前，已成为分子系统树构建的核心方法之一。最大似然法明确地使用概率模型，其目标是寻找能够以较高概率产生观察数据的系统发生树。最大似然法是一类完全基于统计的系统发生树重建方法的代表。2022/9/10BIOINFORMATICS132该方法在每组序列比对中考虑了每个核苷酸替换的概率。例如，转换出现的概率大约是颠换的三倍。在一个三条序列的比对中，如果发现其中有一列为一个C，一个T和一个G，我们有理由认为，C和T所在的序列之间的关系很有可能更接近。2022/9/10BIOINFORMATICS133由于被研究序列的共同祖先序列是未知的，概率的计算变得复杂；又由于可能在一个位点或多个位点发生多次替换，并且不是所有的位点都是相互独立，概率计算的复杂度进一步加大。尽管如此，还是能用客观标准来计算每个位点的概率，计算表示序列关系的每棵可能的树的概率。然后，根据定义，概率总和最大的那棵树最有可能是反映真实情况的系统发生树。2022/9/10BIOINFORMATICS134设物种数为N，对位排列后DNA或氨基酸序列的长度为n，用这些序列组成的矩阵为：

假定不同位点的进化是独立事件，根据该数据矩阵可以进行多种不同的似然估计。2022/9/10BIOINFORMATICS135计算构树数据的似然率对于一棵给定的树，我们希望有一种评价的方法。可以用可能性得分评估我们所作出的假设，即评价所得到的系统发生树T。对于给定的一组分类单元，假设它们的观察值为M（M为向量），可以选择一棵树，使得P(M|T)最大，即最大似然法。2022/9/10BIOINFORMATICS136令tvu代表节点v和u之间的分支长度，反映的是遗传距离或者进化时间。以概率Pxy(tvu)表示在时间tvu内，从状态x转换到状态y的概率。假设有一个矩阵M，它是关于n个分类单元的实际观察值，M描述每个分类单元m个特征的具体取值。同时假设存在一棵树T，其叶节点（如v、u）对应于这些分类单元，而树中的分支代表分类单元之间的距离tvu，求该树的似然值L=P(M|T)。2022/9/10BIOINFORMATICS1372.计算树与子树的似然率

设长度为n的部分似然率矩阵为q，定义qi＝Pix(t)，这里t为树枝长度。有：Qj是部分似然率的乘积。2022/9/10BIOINFORMATICS138计算树枝长度的似然率4.θ的最大似然估计

θ定义为maximizelogL(θ|X,T)θΘ,θ满足此外，用最大似然法还可以估计距离和距离矩阵等。2022/9/10BIOINFORMATICS139

相容法

在系统发育中，与简约法同时发展起来的还有相容法。该方法通过比较性状的相容性来构造若干性状状态树，然后将这些性状状态树综合构建成一个系统树。由于该方法需要发现相容性集团，因而又称为集团分析法。

一些学者讨论了如何确定两个核苷酸位点的相容性问题，同时指出并非所有的成对相容的位点联合起来一定相容，在这种情况下能否发现一个相容系统树是值得怀疑的。迄今为止，由于相容法在理论上还有一些疑问，加之计算上的复杂性，在分子系统发育分析的实际工作中应用极少。2022/9/10BIOINFORMATICS140

系统发育网络系统发育网络是基于系统发育的网络式样，而不是分支(树状)式样的。图6-17示出了一个系统发育网络，它是根据黑猩猩(Ch)，倭猩猩(Bo)，大猩猩(Go)和人类(Hu)的ABO血型基因的7号外显子序列构建的。图中数字表示对应边上的核苷酸位置；边长与核苷酸差异成比例；星号(HCG节点)表示人—黑猩猩—大猩猩共用的祖先节点。2022/9/10BIOINFORMATICS141图6-17基于ABO血型基因7号外显子序列的系统发育网络2022/9/10BIOINFORMATICS142

需要说明的是，造成系统发育网状式样的因素很多。人们正在发展新的统计模型和计算机模拟方法来检测杂交事件，以获得正确的系统发育式样。此外，多种数学工具已开始应用于构建系统发育网络。例如，Moulton等人尝试用直接循环图(directedacyclicgraph，DAG)和最大似然法构建人类T-淋巴细胞病毒(HTLV)的系统发育网络。2022/9/10BIOINFORMATICS143

神经网络方法(neuralnetworkmethod）

人工神经网络是近年来逐渐发展起来的一种建模与分析工具，已广泛应用于生物学研究。该方法在系统发育分析中具有一定的潜力，因为它可以通过“训练”过程来确定相关数据的关系，特别是二级和更高级的相关性。同时，适用于非线性生物系统的进化与分类研究。2022/9/10BIOINFORMATICS1446.3.5方法比较

构建分子系统树的方法很多，然而，目前还没有一种方法可以适合于各种数据或各种条件。从如下几个方面对各种方法进行比较，可能对实际工作中选择合适的方法具有参考价值。2022/9/10BIOINFORMATICS145.假设和计算时间

不同的构树方法一般是建立在不同假设下的。了解这些假设将有助于理解各种方法（特别是当假设条件不满足时应慎用的方法）的应用范围。例如，在距离矩阵方法中，UPGMA法隐含了所有分支中速率相等的假设。然而，很多实际情况违背这一假设，特别是序列长度较短时，该方法构建的系统树更容易造成错误。2022/9/10BIOINFORMATICS146

邻接法没有速率一致的假设，而是采用“校正（速率）”矩阵，它依赖于距离系数的准确性。当序列较短时，对距离的估测可能仍有较大统计误差。最大简约法并没有明晰的假设。当序列间的分歧度很小时，同型情形就少，简约法效果较好。反之，分歧度大时，同型情形很普遍，简约法效果较差。最大似然法对于进化速率和核苷酸置换型式的假设十分明确，同时对违背假定的情形也不敏感。但该方法在计算上仍然较为复杂。2022/9/10BIOINFORMATICS147

计算时间

距离矩阵方法最为省时。

而最大简约法需要比较大量可能的简约树，当OTU的数目大于11时，可能难以穷尽最简约的系统树。分枝界限法是一种节省运算次数的算法。当OTU数目很大时，可以使用“启发式方法”最大似然法最为费时。进化简约法在计算时间上要少于最大似然法和最大简约法，但多