系统发育树(Phylogenetic Tree)
建树方法
邻连法
在生物信息学中,邻连法(Neighbor-joining)是由 Naruya saiitou 和 Masatoshi Nei 于1987年提出的一种自下而上(聚集性)的聚类方法,用于创建系统发育树。该算法构建树通常用于基于DNA或蛋白质序列数据,需要知道每对类群(如物种或序列)之间的距离,以形成树。
邻接法以距离矩阵作为输入,表示每对类群之间的距离。该算法从一个完全未解析的树(A)开始,它的拓扑结构对应于星形网络,并迭代以下步骤,直到树完全解析并知道所有分支长度:
基于当前的距离矩阵计算矩阵Q (n为taxa个数);
在Q(i,j)种寻找到最小值,并将两个点(taxa)连接到一个新的中间节点(u);例如将f,g连接到新的u,u为中间节点;
计算每个类群到这个新节点的距离;例如计算a-e节点到新节点u之间的距离
计算taxa之外的一些节点到新节点的距离;例如计算w到v的距离
迭代上述过程,用新节点替换连接的邻居对,并使用前一步计算的距离。
n个物种使用邻接法重建系统发育树的步骤如下[2]:
(1)对每个物种i计算的净分歧度ri,实际上是以i为中心的星状树树长:
式中,n为终端节点数;dik为i和k之间的距离,从事先计算好的距离矩阵中读出。
(2)计算并确定最小速率校正距离(rate-corrected distance)M:
使得M最小的两个物种记为i和j。
(3)定义一个新节点u,u是节点i和j的父节点。节点u与节点i和j的距离为:
而节点u与其他节点的距离定义为:
(4)从距离矩阵中删除物种i和j,添加物种u,令n=n-1。
(5)如果尚余2个以上物种,返回到步骤(1)继续计算,直至系统树完全建成。
替换模型
两个序列间的替换数(K)是分子进化分析中常用的变量。如果序列比较相似,两个序列比对之间只有较少的替换,那么只要简单数一下替换个数就可以确定K值。然而,如果序列间的差异较大,某些位点可能会发生多次突变(图 7 6),除了我们可以观察到的变化之外,还可能发生了一些“隐藏”突变,因此直接用序列比对的替换个数就可能低估了两序列在最近的共同祖先之后发生的替换数。 图 7-6 两个序列间的多重突变
构建进化树需要估算DNA或蛋白质序列进化距离(distance),即观察到两条序列间替换数所需进化时间的定量测量。在经常发生替换的区域,有些位点可能会发生多重替换。计算进化距离需要考虑核苷酸或氨基酸之间替换的可能替换模型,从而估计实际替换数。蛋白质序列与DNA序列分别有不同的替换模型可选择,实际建树应用中,核酸序列的替换模型一般选择Kimura-2 parameter (Kimura-2参数),而蛋白质序列的替换模型一般选择Poisson correction(泊松修正)。
核苷酸替换模型
- (1) Jukes-Cantor模型 最早在1969年由Thomas Jukes和Charles Cantor提出一种估计每个位点核苷酸替换数K的方法。Jukes-Cantor模型(JC69)假设一组具有较少整体变异的序列比另一组具有大量变异的序列,在一个进化时间内,在任何特定位点经历多次替换的概率更少。如果α表示每单位时间的变化率,并且t表示一个小的单位时间,那么在给定的核苷酸处发生三种可能的取代之一(例如,对于A核苷酸,可能发生A→C,A→G,或A→T)的概率是3αΔt。这个模型假设所有变化都是同等可能的:A被T、G或C替换的可能性是一样的。在对这些变量进行一些操作之后,对可能发生在序列a和b之间的替换数估计(包括可观察的和隐藏的替换),我们可以将替换数K的计算公式表示为:
在该等式中,D是观察到的替换比例(总替换数/总核苷酸数)。然后我们可以直接检查序列比对,在这模型下计算K,并获得可以与其他一些对齐序列对的距离进行比较的进化距离估计。
- (2)Kimura双参数模型 在Jukes-Cantor模型中是以同样概率对待所有的核苷酸替换。但是,后来研究发现不同核苷酸类型的替换率相差较大,特别转换发生的频率至少是颠换的3倍。Motoo Kimura考虑到转换发生的频率高于颠换的情况,并于1980年提出Kimura双参数模型(K-2)。Kimura模型在估算K时同样考虑了“隐藏”替代,并引入另一个参数来解释这种差异,设定每年每个位点转换替代率(α)不同于颠换替代率(β)。 假设观察到的转换比例为S(转换/总替换)和颠换比例为V(颠换/总替换)。两个比对序列之间发生的替换数K的计算公式如下:
在该等式中,S (tranSitions)是序列中转换核苷酸的比例(转换/总替换),V (transVersions)是颠换核苷酸的比例(颠换/总替换)。如果不区分转换与颠换(Dab = S + V),这个方程就可简化成了Jukes-Cantor公式。
- (3)多参数模型 前面讨论了两个比较流行的模型,还存在许多其它模型,考虑更多的影响因素。由于不同物种的GC含量变化很大(如人类基因组约为40%,而一些细菌基因组超过60%),这会引起替换速率的差异。Timura三参数模型就在考虑转换/颠换偏差的基础上,还增加GC含量偏差的参数。Tamura-Nei模型考虑到并非所有转换都以相同的频率发生,如从A到G或G到A的转换与从C到T或T到C的转换的频率不同;Felsenstein模型允许不同位置的突变频率差异; 而Hasegawa-Kishino-Yano(HKY85)模型允许不同碱基都有不同的转换和颠换频率。 在具体应用中要选什么模型合适,还是与研究的序列有关,如GC含量高的序列推荐用Timura三参数模型。但一般认为使用单、双参数模型反而会比多参数模型得到更可靠结果。
氨基酸替换模型
两个蛋白质序列间的不同氨基酸的比例(p)可根据下面简单公式推算:
p = n / L
其中,n代表各种氨基酸在两序列间差别的数量,L是比对序列的长度。
与DNA序列一样,回复突变会导致严重低估替换数目。但是,因为有20个氨基酸,要准确计算两个蛋白质序列间发生的替换数,比估算DNA序列的替换数还要困难。某些氨基酸替换发生的频率会大于其它氨基酸,而且从一个氨基酸转变为另一个氨基酸的替换路径的长度也各不相同,如脯氨酸的CCC密码子转变成亮氨酸CUC密码子只要经历一个突变,而要转变成异亮氨酸AUC密码子,至少要经历两次变化。氨基酸替换对蛋白质功能的影响也不一样,会随氨基酸不同而变化,这令问题更加复杂。
泊松校正(Poisson correction)模型是精确估计氨基酸替代数的一个简单模型,假设给定位点氨基酸替代数的发生频率遵循泊松分布,从而把p转换成蛋白质序列间的氨基酸替换数。每个位点的氨基酸替换数(d)可以用下面公式估算:
References
- Dan Graur & Wen-Hsiung Li, Fundamentals of Molecular Evolution (2nd edn), 2000, Sinauer Associates Inc.
- 朱天琪、杨子恒, 系统发育树与溯祖分析(节选自马占山主编的生物信息学:计算技术和软件导论)