Phylogenetic tree系统发育树(也称进化树)是一个分支图或一棵树,显示了各种生物物种或其他实体之间基于其物理或遗传特征的相似性和差异的进化关系。地球上所有的生命都是一个系统发育树的一部分,表明有共同的祖先。
在一棵有根的系统发育树中,每个有子代的节点代表推断出的这些子代的最近的共同祖先(MRCA),一些树枝长可以被解释为时间估计。每个节点被称为一个分类单位。内部节点一般被称为假设的分类单位,因为它们不能被直接观察到。树在生物信息学、系统学和系统发育学等生物学领域很有用。无根树只说明叶子节点的关联性,不需要知道或推断出祖先的根。
系统发育树是一种表示生物体之间进化关系的图示。系统发育树是假说,而不是确定的事实。系统发育树中的分支模式反映了物种或其他群体是如何从一系列共同的祖先演变而来的。在树上,如果两个物种有一个较近的共同祖先,那么它们的关系就比较近,如果它们有一个较远的共同祖先,那么它们的关系就比较差。系统发育树可以用各种方式绘制。将一棵树围绕其分支点旋转并不改变其承载的信息。Anatomy of a phylogenetic tree在系统发育树中,感兴趣的物种或群体位于被称为树枝顶端。例如,下面的系统发育树代表了五个物种之间的关系,A、B、C、D和E,它们位于树枝的末端。
树枝连接的模式代表了我们对树中物种如何从一系列共同祖先演变而来的理解。每个分支点(也称为内部节点)代表一个分化事件,或将一个群体分裂成两个后代群体。
**most recent common ancestor**在每个分支点上,都有从该分支点下来的所有群体的最近的共同祖先。例如,在产生物种A和B的分支点,我们会发现这两个物种的最近的共同祖先。在树根上方的分支点,我们会发现树上所有物种(A、B、C、D、E)的最近的共同祖先。
树上的每条横线代表一系列的祖先,一直到其末端的物种。例如,通向物种E的线代表该物种的祖先,因为它与树上的其他物种发生了分歧。同样,根代表一系列的祖先,直到树中所有物种的最近的共同祖先。
Which species are more related?在系统发育树中,两个物种的相关度(relatedness)有非常具体的含义。如果两个物种有一个较近的共同祖先,那么它们的亲缘关系就比较好,如果它们有一个较远的共同祖先,那么它们的亲缘关系就比较差。可以用一个相当直接的方法来寻找任何一对或一组物种的最近的共同祖先。在这个方法中,从承载两个感兴趣的物种的分支末端开始,在树上 '倒着走',直到找到物种线的交汇点。例如,假设我们想知道A和B或B和C的关系是否更密切。要做到这一点,我们将在树上沿着两对物种的线条向后走。由于在我们向后移动时,A和B首先在一个共同的祖先处汇合,而B在与A的交界点之后才与C汇合,因此我们可以说A和B比B和C关系更近。重要的是,有一些物种的亲缘关系我们不能用这种方法来比较。例如,我们不能说A和B是否比C和D关系更密切。这是因为,在默认情况下,树的横轴并不直接代表时间。因此,我们只能比较发生在同一枝系(从树根开始的同一枝系)的分支事件的时间,而不能比较发生在不同枝系的分支事件。
Some tips for reading phylogenetic trees你可能会看到许多不同形状的系统发育树。
上面的三棵树代表了物种A、B、C、D和E之间完全相同的关系。这些看起来不同的树中的相同信息提醒我们,在一棵典型的树中,有意义的是分支模式(而不是分支的长度)。
关于这些树的另一个关键点是,如果你旋转结构,使用其中一个分支点作为支点,你不会改变关系。因此,就像上面的两棵树一样,尽管它们的格式不同,但显示的关系是一样的,下面所有的树都显示了四个物种之间的相同关系。
如果你不能马上看出这是真的,就把注意力放在关系和分支点上,而不是放在图的顶部的物种排序(W、X、Y和Z)上。这种排序实际上并没有给我们提供有用的信息。相反,每张图的分支结构告诉我们理解树所需要的东西。到目前为止,我们所看的所有树都有漂亮、干净的分支模式,每个分支点只出现两个世系(血统)。然而,你可能会看到具有多分支(polytomy)的树,这意味着一个分支点有三个或更多的不同物种从它那里出来,一般来说,polytomy表明我们没有足够的信息来确定分支的顺序,可能需要在中间再添加一个物种来区分它们。
Where do these trees come from?为了生成系统发育树,我们经常比较和分析有关物种或其他群体的许多特征。这些特征可以包括外部形态学(形状/外观)、内部解剖学、行为、生化途径、DNA和蛋白质序列,甚至是化石的特征。
为了建立准确的、有意义的树,人们通常会使用许多不同的特征(减少任何一个不完美的数据导致错误的树的机会)。然而,系统发育树是假说,而不是确定的答案,它们取决于你制作用的的数据。随着时间的推移,树会被修改和更新,因为新的数据可以被添加到分析中。这在今天尤其如此,因为DNA测序提高了我们比较物种间基因的能力。
Building a phylogenetic tree系统发育树代表了关于一组生物体之间进化关系的假设。系统发育树可以利用物种或其他群体的形态学(体形)、生物化学、行为学或分子特征来建立。在建立树的过程中,我们根据共同的衍生性状(与群体祖先不同的性状)将物种组织成嵌套群体。基因或蛋白质的序列可以在物种之间进行比较,并用于建立系统发育树。关系密切的物种通常有很少的序列差异,而关系不大的物种往往有更多的差异。The idea behind tree construction我们如何建立一个系统发育树?其基本原则是达尔文的 '后代变异 '思想。基本上,通过观察现今生物体的修饰模式(新的性状),我们可以找出或者至少对它们从一个共同的祖先传下来的路径做出假设。
作为一个例子,让我们考虑下面的系统发育树(它显示了一组捏造的类似老鼠的物种的进化历史)。我们看到在这个群体的进化史上,有三个新的特征在不同的时间点出现:毛茸茸的尾巴、大耳朵和胡须。每种新性状都被产生该性状的祖先的所有后代物种所共享(用勾号表示),但在该性状出现之前分裂出来的物种则没有。
当我们建立系统发育树时,在一个群体的进化过程中出现的、与该群体的祖先的性状不同的性状被称为衍生性状。在我们的例子中,毛茸茸的尾巴、大耳朵和胡须是衍生性状,而瘦小的尾巴、小耳朵和没有胡须是祖先的性状。重要的一点是,一个衍生性状可以通过一个特征的损失或增益出现。例如,如果E谱系上有另一种变化导致尾巴的丧失,那么无尾就会被认为是一种衍生性状。
数据集中的物种或其他群体之间共享的衍生性状是帮助我们建树的关键。如上所示,共有的衍生性状往往形成嵌套模式,提供了物种进化过程中何时发生分支事件的信息。
当我们从一个数据集中构建系统发育树时,我们的目标是利用现今物种的共有衍生性状来推断其进化史的分支模式。然而,诀窍在于我们无法观察我们感兴趣的物种的演变,也无法看到每个谱系中何时出现了新的性状。
相反,我们必须向后找。也就是说,必须观察我们感兴趣的物种--如A、B、C、D和E--并找出哪些性状是祖先的,哪些是衍生的。然后,我们可以利用共同的衍生性状将物种组织成如上图所示的嵌套组合。以这种方式制成的树是关于物种进化史的假设--通常是一个最简单的分支模式,可以解释它们的性状。
Example: Building a phylogenetic tree作为我们研究的一部分建立一棵系统发育树,我们将不得不选择哪一组生物体来排列成一棵树。我们还必须选择这些生物的哪些特征作为树的基础(从它们许多不同的物理、行为和生物化学特征中)。
如果我们要为一个类别建立系统发育树,那么我们有可能会得到一组特征,通常是以表格的形式,我们需要将其转换为一棵树。例如,这个表格显示了各种特征的存在( )或不存在(0)。
7LXl3g接下来,我们需要知道每个特征的哪种形式是祖传的,哪种是衍生的。例如,肺的存在是一个祖先的特征,还是一个衍生的特征?作为提醒,祖先特征是我们认为存在于相关物种的共同祖先中的特征。衍生性状是我们认为在该祖先后代的某处出现的一种形式。如果没有能力查看过去(这很方便,但可惜不可能),我们如何知道哪些性状是祖先的,哪些是衍生的?
在家庭作业或考试的背景下,你要解决的问题可能会告诉你哪些特征是衍生的,哪些是祖先的。
如果你在做自己的研究,你可能有一些知识可以让你识别祖先和衍生性状(例如,基于化石)。
你可能会得到一个外群的信息,一个与感兴趣的物种的关系比它们彼此之间更远的物种,如果我们得到了一个外群,外群可以作为祖先物种的代表。也就是说,我们可以假设它的特征代表了每个特征的祖先形式。
外群,就像我们试图研究的物种一样,随着时间的推移继续进化。因此,外群物种有可能(独立地)获得一些在我们感兴趣的物种中发现的相同的衍生性状。在这种情况下,我们可能会意外地将一个衍生性状确定为祖先性状。
例如,在我们的例子中,鳗鱼是一种没有骨架的无颌鱼,是我们的外群。如表中所示鳗鱼缺乏所有列出的特征:它没有肺、颌、羽毛、砂囊或毛。根据这些信息,我们将假设没有这些特征是祖先的,而每个特征的存在是一种衍生性状。
现在,我们可以开始建立我们的树,根据它们共同的衍生特征对生物体进行分组。一个好的开始是寻找最多的生物体之间共享的衍生特征。在本例中,这就是颌的存在:除了外群物种(鳗鱼)之外,所有的生物都有颌。因此,我们可以在树上画出从其他物种中分离出来的鳗鱼枝系,我们可以把颌的出现放在非鳗鱼物种的分支上。
接下来,我们可以寻找下一个最大的生物群体所共有的衍生性状。那就是是肺,羚羊、秃鹰和鳄鱼都有,但海鲈鱼没有。基于这种模式,我们可以画出海鲈鱼的分支谱系,我们可以把肺的出现放在通往羚羊、秃鹰和鳄鱼的谱系上。
按照同样的模式,我们现在可以寻找下一个最大数量的生物体所共享的衍生特征。那就是鳄鱼和秃鹰共有的砂囊(羚羊没有)。根据这些数据,我们可以画出羚羊系从短吻鳄和秃鹰系中分离出来的分支,并将砂囊的出现放在后者身上。
那我们剩下的毛皮和羽毛的特征呢?这些性状是衍生的,但它们不是共享的,因为每一个都只在单一物种中发现。不共享的衍生性状不能帮助我们建立一棵树,但我们仍然可以把它们放在树上最可能的位置。对于羽毛,这是在通往秃鹰的枝系上(在与鳄鱼分化之后)。对于毛皮来说,这是在羚羊枝系上,在它与鳄鱼和秃鹰分化之后。
Parsimony and pitfalls in tree construction当我们在构建上述树时,我们使用了一种叫做简约法的方法。简约法本质上意味着我们要选择能够解释我们观察结果的最简单的解释。在做树的背景下,它意味着我们选择需要最少的独立遗传事件(性状的出现或消失)发生的树。例如,我们也可以用下面的树来解释我们看到的性状模式。
这一系列的事件也为我们在五个物种中看到的性状提供了一个进化的解释。然而,它不太合理,因为它需要更多独立的性状变化发生。因为我们把海鲈鱼放在哪里,我们必须假设颌独立出现了两次(一次在海鲈鱼,一次在羚羊、秃鹰和短吻鳄)。这使得该树共有6个标记,或性状变化事件,而在上述更简明的树中只有5个标记。
在这个例子中,似乎相当明显的是有一棵最好的树,而计算标记似乎不是很必要。然而,当研究人员将系统发育作为他们工作的一部分时,他们经常使用大量的特征,而这些特征的模式很少能100/100%,百分百地相互一致。相反,会有一些冲突,即一棵树更适合于一个特征的模式,而另一棵树更适合于另一个特征的模式。在这些情况下,研究人员可以使用简约法来选择最适合数据的一棵树(假设)。
你可能想知道。为什么这些树不都是彼此一致的,不管它们是建立在什么特征上?毕竟,一群物种的进化在过去确实以一种特定的方式发生。问题是,当我们建立一棵树时,我们是在用不完整的、有时是不完美的数据来重建这一进化历史。比如说。
我们可能并不总是能够区分反映共同祖先的特征(同源特征)和类似但独立产生的特征(通过趋同进化产生的类似特征)。
趋同进化:并非所有看起来相似的身体特征都是共同祖先的标志。相反,一些身体上的相似之处是类似的:它们在不同的生物体中独立进化,因为这些生物体生活在类似的环境中或经历类似的选择压力。这个过程被称为趋同进化。
例如,胡须出现了两个独立进化。如果我们不知道这个群体的真实历史并试图重建它,我们可能会把胡须解释为产生于一个单一事件。这样,胡须的数据就会与其他性状的数据相冲突。
在一个物种的进化史上,性状可以多次获得和丧失。一个物种可能有一个衍生性状,但在进化过程中又失去了这个性状(恢复到祖先的形式)。
在这棵树上,物种E经历了一个基因变化,导致它失去了茂密的尾巴,而获得了该群体祖先中的瘦小尾巴。如果我们不知道这个群体的真实历史,并试图重建它,我们可能会假设物种E是一个没有蓬松尾巴的祖先的后代。在这种假设下,尾巴数据将与其他性状的数据相冲突。
生物学家经常使用许多不同的特征来建立系统发育树,因为有像这样的错误来源。即使所有的特征都被仔细选择和分析,仍有可能使其中一些特征得出错误的结论(因为我们没有关于过去发生的事件的完整信息)。
Using molecular data to build trees一个已经彻底改变并将继续改变系统发育分析的工具是DNA测序。通过DNA测序,我们不是用生物体的物理或行为特征来建树,而是可以比较它们的直系同源基因(进化上相关)基因或蛋白质的序列。
这种比较的基本原则与我们上面所讲的相似:有一个祖先形式的DNA或蛋白质序列,在进化过程中可能发生了变化。然而,一个基因或蛋白质并不只是对应于一个存在于两种状态的单一特征。
相反,一个基因的每个核苷酸或一个蛋白质的氨基酸可以被看作是一个独立的特征,一个可以通过突变翻转到多种状态(例如,一个核苷酸的A、T、C或G)。因此,一个有300个核苷酸的基因可以代表300个不同的特征,存在于4个状态中!我们从序列中得到的信息量很大。因此,我们可以期望在系统发育树中得到的分辨率,要比我们使用物理性状时高得多。
为了分析序列数据并确定最可能的系统发育树,生物学家通常使用计算机程序和统计算法。
不过一般来说,当我们在物种之间比较一个基因或蛋白质的序列时。差异数越大,说明物种之间的关系越差较少的差异数对应于更相关的物种例如,假设我们比较人类和其他各种物种之间的血红蛋白(血液中的携氧蛋白)的β链。如果我们比较人类和大猩猩的蛋白质,我们会发现只有1个氨基酸的差异。如果我们比较人类和狗的蛋白质,我们会发现15个差异。人类与鸡相比,氨基酸差异高达45个,而人类与鳗鱼(一种无颌鱼)相比,我们看到了127个差异。这些数字反映出,在所考虑的物种中,人类与大猩猩的关系最大,与鳗鱼的关系最远。
Wikipedia:Phylogenetic treeA highly resolved, automatically generated tree of life, based on completely sequenced genomesHistory查尔斯-达尔文(1859年)制作了最早的进化树插图之一,并在其开创性著作《物种起源》中关键性地普及了进化 '树 '的概念。一个多世纪后,进化论生物学家仍然使用树状图来描述进化,因为这种图有效地传达了这样一个概念,即通过适应性的和半随机的血统分裂来实现物种进化。随着时间的推移,物种分类已不再是静态的,而是更加动态的。
Tree of life基于rRNA基因的系统发育树,显示了三大界:细菌界、古细菌界和真核生物界。系统发育树底部的黑色分支将生物体的三个分支连接到最后的普遍共同祖先。在没有外群的情况下,根是推测的。
PropertiesRooted tree有根的系统发育树是一棵有指向性的树,它有一个独特的节点--根,对应于树上所有实体的(通常是假定的)最近的共同祖先。根节点没有父节点,但是作为树中所有其他节点的父节点。因此,根是一个度数为2的节点,而其他内部节点的最小度数为3(这里的 '度数 '是指传入和传出的边的总数)。
最常见的定根方法是使用一个没有争议的外群--足够接近于允许从性状数据或分子测序中进行推断,但又足够远,成为一个明确的外群。另一种方法是中点生根(midpoint rooting),或者通过使用非稳态替换模型也可以使树生根 。
Rooted treeUnrooted tree无根树说明了叶子节点的关联性,而没有对祖先进行假设。它们不需要知道或推断出祖先的根。无根树总是可以通过简单地省略根而从有根树中生成。相比之下,推断无根树的根需要某种手段来识别祖先。这通常是通过在输入数据中加入一个外群,使树根必然在外群和树中其他类群之间,或者通过引入关于每个分支的相对进化率的额外假设,如应用分子钟假说来实现。
Special tree typesDendrogram树状图是树的总称,无论是否是系统发育树。
Dendrogram of the phylogeny of some dog breedsCladogramA cladogram only represents a branching pattern; i.e., its branch lengths do not represent time or relative amount of character change, and its internal nodes do not represent ancestors.[13]
Phylogram进化分枝图只代表一个分支模式;也就是说,它的分支长度不代表时间或特征变化的相对数量,它的内部节点不代表祖先。
鳞翅目动物的年代图。在这种系统发育树类型中,分支长度与地质时间成正比。
Spindle diagram纺锤图或气泡图,因为它是由美国古生物学家Alfred Romer推广的。它表示分类学的多样性(水平宽度)与地质时间(垂直轴),以反映各种分类群的丰度随时间的变化。然而,纺锤图并不是进化树:分类学上的纺锤形掩盖了亲代分类群与子代分类群的实际关系,并有涉及亲代分类群的副系的缺点。
9Ctq7SCoral of life达尔文提到的珊瑚可能是一个比树更合适的比喻。系统发育的珊瑚图对于描绘过去和现在的生命很有用,而且它们比树有一些优势(允许吻合等) 。
The Coral of LifeConstruction使用计算系统发育学方法构建由非微不足道的输入序列组成的系统发育树。距离矩阵方法,如邻接法或UPGMA,从多个序列排列中计算遗传距离,是最简单的实现方法,但没有调用进化模型。许多序列比对方法,如ClustalW,也是通过使用较简单的建树算法(即基于距离的算法)来创建树木。最大简约法是另一种估计系统发育树的简单方法,但意味着隐含的进化模型(即简约法)。更高级的方法使用最大似然的最优性标准,通常是在贝叶斯框架内,并将明确的进化模型应用于系统发育树的估计。[3] 使用许多这些技术来识别最优树是NP困难的,[3] 因此启发式搜索和优化方法与树的评分函数结合使用,以确定符合数据的合理的好树。
树形构建方法可以根据以下几个标准来评估:[21]
效率(计算答案需要多长时间,它需要多少内存?)
功率(它是否很好地利用了数据,或者信息被浪费了?)
一致性(如果每次给定同一模型问题的不同数据,它是否会重复收敛于同一答案?)
稳健性(它是否能很好地应对违反基础模型假设的情况?)
可证伪性(当它不好用时,即当假设被违反时,它是否会提醒我们?)
File formatsNexus file一个NEXUS文件是由一个固定的标题#NEXUS和多个块组成。每个区块以**BEGIN block_name;开始,以END;**结束。关键词是不区分大小写的。Comment被括在方括号[...]内
对于常见的数据类型,有一些预定义的块名。例子包括:
TAXA块
TAXA块包含关于分类群的信息。
数据块
DATA块包含了数据矩阵(例如:序列比对)。
TREES块
TREES块包含使用Newick格式描述的系统发育树,例如:((A,B),C);。
下面的例子使用了上述三种块的类型:
#NEXUSBegin TAXA; Dimensions ntax=4; TaxLabels SpaceDog SpaceCat SpaceOrc SpaceElf;End;Begin data; Dimensions nchar=15; Format datatype=dna missing=? gap=- matchchar=.; Matrix [ When a position is a 'matchchar', it means that it is the same as the first entry at the same position. ] SpaceDog atgctagctagctcg SpaceCat ......??...-.a. SpaceOrc ...t.......-.g. [ same as atgttagctag-tgg ] SpaceElf ...t.......-.a. ;End;BEGIN TREES; Tree tree1 = (((SpaceDog,SpaceCat),SpaceOrc,SpaceElf));END;Newick formatExamples(,,(,)); no nodes are named(A,B,(C,D)); leaf nodes are named(A,B,(C,D)E)F; all nodes are named(:0.1,:0.2,(:0.3,:0.4):0.5); all but root node have a distance to parent(:0.1,:0.2,(:0.3,:0.4):0.5):0.0; all have a distance to parent(A:0.1,B:0.2,(C:0.3,D:0.4):0.5); distances and leaf names (popular)(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5)F; distances and all names((B:0.2,(C:0.3,D:0.4)E:0.5)F:0.1)A; a tree rooted on a leaf node (rare)当一棵无根的树用Newick符号表示时,一个任意的节点被选为它的根。无论是有根的还是无根的,通常一棵树的表示都是以内部节点为根。
一个有根二叉树是以内部节点为根的,每个内部节点正好有两个直系子节点。
一棵在任意内部节点上生根的无根二叉树,其根节点正好有三个直系子节点,而其他每个内部节点正好有两个直系子节点。
一棵从叶子上生根的二叉树,根节点最多只有一个直系子节点,而每个内部节点正好有两个直系子节点。
Tree: The full input Newick Format for a single treeSubtree: an internal node (and its descendants) or a leaf nodeLeaf: a node with no descendantsInternal: a node and its one or more descendantsBranchSet: a set of one or more BranchesBranch: a tree edge and its descendant subtree.Name: the name of a nodeLength: the length of a tree edge.Bayesian Estimation of Species Divergence Times Integrating Fossil and Molecular Information由于基因随着时间的推移以恒定的速度积累变化,两个物种之间的遗传距离,由积累的变化数量与物种分化的时间成正比。因此,分子可以作为一个时钟,通过积累的变化保持物种分化的时间。如果化石记录或地质事件可以用来给系统发育树上的物种分化事件指定一个绝对的地质时间,那么就可以把所有计算出来的遗传距离转换成绝对的地质时间。这种分子钟测年法的原理可以扩展到处理进化速率的局部变异中。对分子钟测年法的关键是使用化石信息来校准时钟。
地质时间 * 碱基替换速率 = 遗传距离
Approximate likelihood calculation for Bayesian estimation of species divergence times从分子数据中对分歧时间进行贝叶斯推断在计算上是很昂贵的。分歧时间的后验分布不能用分析法计算,所以必须依靠马尔科夫链蒙特卡洛(MCMC)方法来获得后验分布的近似值。在MCMC过程中计算分子序列排列的log-likelihood是后验计算中最耗费时间的部分。Thorne及其同事(1998, MBE 15:1647)提出了一个对数似然函数的正态近似,以加快MCMC期间的计算。然而,这种近似的准确性从未被充分测试过。在这个项目中,我们对该近似进行了详细的数学分析,并建议使用变量变换来提高准确性。使用真实的分子数据进行的广泛测试表明,该近似值是高度可靠的。我们新的基于ARCSINE的变换现在是MCMCTREE程序中默认的近似似然方法。你可以在这里阅读2011年2月发表在MBE上的关于近似似然法的论文。用传统的精确方法对几个物种的分化时间进行贝叶斯估计可能需要几天到几个月的时间,而用近似方法则可能在一到两天内完成分析。对于基因组规模的数据,精确的方法变得过于昂贵。MCMCTREE中的近似方法使我们能够使用超过2100万个位点的排列来估计哺乳动物物种的分化时间,并解决了分子进化中的一个长期争议(见上文和这里)。
The impact of fossil calibrations on divergence time estimation贝叶斯统计中多维先验的构建是一个复杂的问题。物种分化时间的贝叶斯估计涉及到分化时间的联合先验的构建。对于一组s个物种,有s-1个时间需要估计,必须构建一个s-1个多维联合时间先验。不同的贝叶斯分歧时间估计程序(如MULTIDIVTIME、BEAST或MCMCTREE)使用不同的策略将化石信息整合到系统发育中,以构建时间先验。这些不同的策略可能导致看似相似的化石校准密度产生不同的先验。反过来,尽管使用相同的分子数据和化石信息,不同的先验可能导致不同程序的后验时间估计值出现惊人的差异。
MCMCTree tutorialsMCMCTree在各种分子钟模型下,利用soft fossil constraints对物种分歧时间进行贝叶斯估计。该程序使用一个多序列比对(核苷酸或蛋白质)作为输入,一个带有化石校准的系统发育树,以及一个包含程序指令的控制文件(通常称为mcmctree.ctl)。MCMCTree是PAML软件包的一部分。
Tutorial 1: Divergence time of apes在本教程中,将分析7种猿类的线粒体蛋白编码基因的数据集。包含在PAML软件包的示例文件中(examples/DatingSoftBound)。文件 'mtCDNApri123.txt '包含了核苷酸序列比对信息。该比对方式被分为3个部分,分别对应于第1、第2和第3个密码子位点。文件 'mtCDNApri.trees '包含与化石标定的7个物种相关的系统发育树。
7 1((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);//end of file((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);((((human, (chimpanzee, bonobo)) 'B(.06, .08)', gorilla), (orangutan, sumatran)) 'B(.12, .16)', gibbon);第一行是物种的数量(7)和树的数量(1)。然后给出Newick格式的树。该树不能有分支长度。该树有两个化石定标,一个是人类-黑猩猩的最近共同祖先:'> .06 < .08',另一个是大猩猩的最近共同祖先:'> .12 < .16'。时间单位是100个百万年(Myr)。因此,第一项校准将人类和黑猩猩的共同祖先限制在6-8百万年前。MCMCTree中的校正界限是soft的,有一个小概率(默认为0.025)的界限可能被违反。请注意,该树在根部没有化石校准。MCMCTree总是需要对树的根部进行标定,当这种标定不存在于树的文件中时,必须在控制文件中进行标定(变量RootAge)。
控制文件 'mcmctree.ctl '包含运行MCMCTree程序的所有必要指令。你可以用文本编辑器打开这个文件:
seed = -1 seqfile = mtCDNApri123.txt treefile = mtCDNApri.trees mcmcfile = mcmc.txt outfile = out.txt ndata = 3 seqtype = 0 * 0: nucleotides; 1:codons; 2:AAs usedata = 1 * 0: no data; 1:seq like; 2:normal approximation; 3:out.BV (in.BV) clock = 2 * 1: global clock; 2: independent rates; 3: correlated rates RootAge = '<1.0' * safe constraint on root age, used if no fossil for root. model = 0 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85 alpha = 0 * alpha for gamma rates at sites ncatG = 5 * No. categories in discrete gamma cleandata = 0 * remove sites with ambiguity data (1:yes, 0:no)? BDparas = 1 1 0.1 * birth, death, sampling kappa_gamma = 6 2 * gamma prior for kappa alpha_gamma = 1 1 * gamma prior for alpha rgene_gamma = 2 20 1 * gammaDir prior for rate for genes sigma2_gamma = 1 10 1 * gammaDir prior for sigma^2 (for clock=2 or 3) finetune = 1: .1 .1 .1 .1 .1 .1 * auto (0 or 1): times, musigma2, rates, mixing, paras, FossilErr print = 1 * 0: no mcmc sample; 1: everything except branch rates 2: everything burnin = 2000 sampfreq = 10 nsample = 20000*** Note: Make your window wider (100 columns) before running the program.seed:这设置了程序所使用的随机种子。当设置为一个负整数时(比如-1,就像上面给出的例子),程序将使用计算机的当前时间来设置种子,所以每次你运行MCMCTree时,它将从一个不同的种子开始,MCMC的结果看起来也会不同。如果你需要可重复的结果,你可以将种子设置为一个奇数或一个偶数。你应该使用这个选项,并至少运行两次程序,以确定运行的结果非常相似。seqfile和treefile:这分别是序列比对文件和树文件的名称。mcmcfile:关于分析的参数的MCMC运行报告。outfile:一旦程序完成,它将把结果的摘要写到这个文件中。ndata:比对文件中分区的数量。在我们的例子中,有三个分区,分别对应于线粒体蛋白质中的三个密码子位置。usedata:当设置为1时,似然函数以正常方式计算,MCMC分析正常进行。当usedata=0时,不计算似然函数(它被设置为1),所以只计算先验。当usedata=2和=3时,将进行近似似然计算和分支长度的ML估计。它们可用于分析核苷酸、氨基酸和密码子序列;分别使用核苷酸、氨基酸和密码子替换模型。它规定了近似值计算,输入(梯度和Hessian矩阵等)在file中。 * 0: no data不使用,不会进行likelihood估算,会快速得到mcmc树,但分歧时间不可用; * 1:seq like,使用多序列比对数据进行likelihood估算,正常进行mcmc; usedata=1时model无法选择; * 2:normal 进行正常的approximation likelihood分析,不读取多序列比对数据,直接读取当前目录的in.BV文件,in.BV是由usedata = 3时生成的out.BV重命名得来; * 3:程序利用多序列比对数据调用baseml/codeml命令对数据进行分析,生成out.BV文件。seqtype:变量,为核苷酸序列设置0,为密码子序列设置1,为氨基酸序列设置2。clock:要使用的时钟模型。这里我们将使用独立速率模型(clock=2),其中速率遵循对数正态分布(也就是说,速率的对数是正态分布)。RootAge:如果树文件中没有提供,则使用根的校准。这里我们使用'<1.0'即所有类人猿最近的共同祖先最大限制为100Myr。model, alpha and ncatG:要使用的替换模型。在这个例子中,我们将使用JC69(它的计算非常快)。因为α=0,所以我们在本例中不使用速率变化的gamma模型。cleandata=0:意味着在似然计算中,比对gap和模糊字符将被视为缺失数据。= 1意味着在分析之前,任何至少有一个序列有比对gap或模糊字符的位点都将被删除。这个变量用于usedata=1和3,如果usedata=2,则没有影响。BDparas:控制出生-死亡过程的参数。出生-死亡过程用于构建树中没有化石标定的节点的时间先验。这里我们使用默认的1 1 0.1,它产生了统一的节点年龄先验。kappa_gamma和alpha_gamma:设置kappa(转换/颠换比率)的GAMMA分布参数;设置GAMMA形状参数alpha的GAMMA分布参数。rgene_gamma: 设置序列中所有位点平均[碱基/密码子/氨基酸]替换率的Dirichlet-GAMMA分布参数:alpha=2、beta=20、初始平均替换率为每100个million年(取决于输入有根树文件中的时间单位)1个替换。若时间单位由100Myr变换为1Myr,则要设置成'2 2000 1'。总体上的平均进化速率为:2 / 20 = 0.1 个替换 / 每100Myr,即每个位点每年的替换数为 1e-9。sigma2_gamma:设置所有位点进化速率取对数后方差(sigma的平方)的Dirichlet-GAMMA分布参数:alpha=1、beta=10、初始方差值为1。当clock参数值为1时,表示使用全局的进化速率,各分枝的进化速率没有差异,即方差为0,该参数无效;当clock参数值为2时,若修改了时间单位,该参数值不需要改变;当clock参数值为3时,若修改了时间单位,该参数值需要改变。print:如果设置为1,MCMC的输出和结果的摘要将被写入硬盘(MCMC被写入mcmc.out文件,摘要被写入上述设置的outfile文件)。如果设置为2,每个基因座(分区)的分支率将被附加到输出中。请注意,对于宽松的时钟模型(时钟=2或3),输出的速率相当大,可能需要更多的时间将它们全部写入输出文件(打印=2)。因此,根据你使用的模型,你可以决定 'print '是设置为1还是2。如果设置为0,结果将只被打印到屏幕上。你不希望这样。burnin、sampfreq和nsample:在我们的例子中,程序将放弃前2000次迭代作为burn-in,然后它将每10次迭代采样一次,直到收集到20000个样本。总的来说,MCMC将运行2,000 10×20,000=202000次迭代。通常情况下,你应该收集10,000到20,000个样本以获得良好的统计总结。大样本量(例如100,000)往往会浪费大量的硬盘空间,提供很少的统计改进。程序也可能需要很长的时间来总结结果。如果你需要增加MCMC的长度(以提高收敛性),增加sampfreq,但保持nsample在一个合理的值。运行程序
mcmctree mcmctree.ctl一旦MCMC完成(达到100%),程序将对结果进行总结,并将总结打印到屏幕上。该程序还将生成几个输出文件。'out.txt'、'SeedUsed'、'mcmc.txt '和 'FigTree.tre'。
第一棵树只包含节点标签。第二棵树包含以时间为单位的分支长度。第三棵树包含以时间为单位的分支长度,以及节点年龄的可信区间。最后三棵树(只有当你把 '打印 '设置为2时才会写在输出文件中)有碱基替换率而不是分支长度,每棵树代表每个基因座(在我们的例子中,三个密码子位置中的每一个)。在树之后,后验平均值、相应的95%可信区间和最高后验概率密度95%可信区间(HPD,即在后验下具有给定概率的最短可能区间,本例中为0.95)。
out.txtSpecies tree for FigTree. Branch lengths = posterior mean times; 95% CIs = labels((((1_human, (2_chimpanzee, 3_bonobo) 12 ) 11 , 4_gorilla) 10 , (5_orangutan, 6_sumatran) 13 ) 9 , 7_gibbon) 8 ;((((human: 0.066223, (chimpanzee: 0.028385, bonobo: 0.028385): 0.037838): 0.025932, gorilla: 0.092155): 0.053388, (orangutan: 0.052619, sumatran: 0.052619): 0.092923): 0.009781, gibbon: 0.155323);((((human: 0.066223, (chimpanzee: 0.028385, bonobo: 0.028385) [&95%={0.0196606, 0.0356899}]: 0.037838) [&95%={0.0467423, 0.0799248}]: 0.025932, gorilla: 0.092155) [&95%={0.0655926, 0.109516}]: 0.053388, (orangutan: 0.052619, sumatran: 0.052619) [&95%={0.0371283, 0.0650519}]: 0.092923) [&95%={0.104711, 0.168268}]: 0.009781, gibbon: 0.155323) [&95%={0.110733, 0.185412}];FigTree.tre#NEXUSBEGIN TREES; UTREE 1 = ((((human: 0.066223, (chimpanzee: 0.028385, bonobo: 0.028385) [&95%={0.0196606, 0.0356899}]: 0.037838) [&95%={0.0467423, 0.0799248}]: 0.025932, gorilla: 0.092155) [&95%={0.0655926, 0.109516}]: 0.053388, (orangutan: 0.052619, sumatran: 0.052619) [&95%={0.0371283, 0.0650519}]: 0.092923) [&95%={0.104711, 0.168268}]: 0.009781, gibbon: 0.155323) [&95%={0.110733, 0.185412}];END;Tutorial 2: Divergence time of apes with approximate likelihood calculation使用近似似然方法估计分歧时间有两个步骤。在第一步中,通过最大似然法估计分支长度,以及最大似然函数的gradient和Hessian。在第二步中,使用MCMC进行分歧时间的估计。
创建一个名为 'Hessian '的新文件夹,并将树、比对和控制文件复制到这个文件夹中。打开控制文件 'mcmctree.ctl'。将usedata = 3。
运行程序
mcmctree mcmctree.ctlMCMCTree将为比对中的每个分区创建四个临时文件。'tmp000X.txt '包含第一个分区的比对,'tmp000X.trees '包含第一个比对的树,'tmp000X.ctl '是BASEML程序分析 'tmp000X.trees '的控制文件,'tmp000X.out '是分析的总结;X为分析的编号(如1,2,3,...)。然后MCMCTree将调用BASEML三次,对三个分区的分支长度、gradient和Hessian进行最大似然估计。
现在,打开名为 'out.BV '的文件。这个文件包含每个分区的三组分支长度、gradient和Hessian:
7(((human: 0.025136, (chimpanzee: 0.013241, bonobo: 0.010461): 0.014365): 0.013406, gorilla: 0.029460): 0.025520, (orangutan: 0.022252, sumatran: 0.020639): 0.041605, gibbon: 0.062501); 0.025520 0.013406 0.025136 0.014365 0.013241 0.010461 0.029460 0.041605 0.022252 0.020639 0.062501 0.000000 0.006237 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000Hessian -9.652e 04 -3713 -1.037e 04 -1.176e 04 -1.373e 04 -2.065e 04 -8997 -4832 -2037 -2380 -4792 -3713 -1.766e 05 -4315 -143.4 -1.372e 04 -9111 -2.129e 04 -7018 2411 -5781 -9587第一行是物种的数量(7),然后是带有分支长度的(无根)树,然后是2×7-3=11个分支长度的矢量,然后是gradient通常都是零值)和Hessian矩阵(11×11=121值)。如果你向下滚动文件,你会看到另一个区块,其中有另一棵树、一组分支长度等,对应于第二个分区;再往下,最后一个区块对应于第三个分区。每当BASEML完成对一个分区的分析,它就会把树、gradient和Hessian写入一个叫做 'rst2 '的文件中。MCMCTree会收集 'rst2 '文件,并将它们合并到更大的 'out.BV '文件中。
创建一个名为' approx01 '的新文件夹。在这个文件夹中复制树、比对、控制和 'out.BV '文件。进入新文件夹,将 'out.BV '文件重命名为 'in.BV'。打开 'mcmctree.ctl '文件,修改usedata = 2。
mcmctree mcmctree.ctlMCMCTree现在将进行分歧时间估计,但这次是使用gradient和Hessian来计算近似似然函数。与任何MCMC一样,你需要再次运行分析以检查收敛性(convergence)。创建一个新的文件夹,将其称为 'actrox02',并重复分析的MCMC步骤(没有必要重复第一步对分支长度、gradient和Hessian的估计)。
Tutorial 3: Approximate likelihood with protein data如果输入的比对文件是氨基酸序列,如果我们想使用似然的方法,还需要一些额外的步骤。我们将使用 'examples '目录下的 'abglobin.aa '比对文件。这个文件包含了五种哺乳动物的球蛋白序列。创建一个名为 'mcmctreeglobin '的新目录,并将 'abglobin.aa '文件复制到其中。创建树文件,并将其称为 'abglobin.trees'。
修改配置文件
seed = -1 seqfile = abglobin.aa treefile = abglobin.trees mcmcfile = mcmc.txt outfile = out.txt ndata = 1 seqtype = 2 * 0: nucleotides; 1:codons; 2:AAs usedata = 3 * 0: no data; 1:seq like; 2:normal approximation; 3:out.BV (in.BV) clock = 2 * 1: global clock; 2: independent rates; 3: correlated rates RootAge = '<1.0' * safe constraint on root age, used if no fossil for root.model = 0 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85 alpha = 0 * alpha for gamma rates at sites ncatG = 5 * No. categories in discrete gamma cleandata = 0 * remove sites with ambiguity data (1:yes, 0:no)?BDparas = 1 1 0.1 * birth, death, sampling kappa_gamma = 6 2 * gamma prior for kappa alpha_gamma = 1 1 * gamma prior for alpha rgene_gamma = 2 20 1 * gamma prior for rate for genes sigma2_gamma = 1 10 1 * gamma prior for sigma^2 (for clock=2 or 3) finetune = 1: .1 .1 .1 .1 .1 .1 * auto (0 or 1): times, rates, mixing, paras, RateParas, FossilErr print = 1 * 0: no mcmc sample; 1: everything except branch rates 2: everything burnin = 2000 sampfreq = 10nsample = 20000mcmctree mcmctree.ctlMCMCTree将生成 'tmp0001.ctl'、'tmp0001.trees'、'tmp0001.out '和 'tmp0001.txt '文件,并调用CODEML来生成蛋白质数据的Hessian矩阵。然而,MCMCTree将使用最简单的蛋白质模型(Poisson and no gamma rates),这对实际数据分析不是很有用。删除生成的 'out.BV '和 'rst '文件。将 'wag.dat '文件从 'dat '目录中复制到 'mcmctree-globin'。打开 'tmp0001.ctl',并编辑它
seqfile = tmp0001.txt treefile = tmp0001.trees outfile = tmp0001.out noisy = 3 seqtype = 2 model = 2 * 2: Empirical aaRatefile = wag.dat fix_alpha = 0 alpha = .5 ncatG = 4 Small_Diff = 0.1e-6 getSE = 2 method = 1这个新的控制文件将以empirical rate matrix(WAG)和位点间的gamma rates运行。现在你可以调用CODEML
codeml tmp0001.ctl来生成合适的Hessian矩阵,使用WAG Gamma。将文件 'rst2 '重命名为 'in.BV',现在你有了一个用WAG Gamma计算的Hessian矩阵。编辑 'mcmctree.ctl '并设置 usedata = 2
mcmctree mcmctree.ctlFigTree.tre#NEXUSBEGIN TREES; UTREE 1 = ((((rabbit: 0.813250, rat: 0.813250) [&95%={0.406218, 1.27432}]: 0.268244, human: 1.081494) [&95%={0.62304, 1.55951}]: 0.251728, goat-cow: 1.333223) [&95%={0.819077, 1.78561}]: 0.467555, marsupial: 1.800778) [&95%={1.70009, 1.90029}];END;Figtree查看结果
当然,想要做成像文章里的图一样,可能需要iTol,ggtree等进行绘图。
下篇文章预告:
CAFE5分析基因家族扩张和收缩
ggtree进化树绘制
Referencehttps://en./wiki/Phylogenetic_treeBuilding a phylogenetic tree (article) | Khan AcademyPhylogenetic trees | Evolutionary tree (article) | Khan Academyhttps://yanzhongsino./2021/03/25/bioinfo_caculate.divergence.time/http://abacus.gene./software/MCMCtree.Tutorials.pdf