计算分子进化1.1-1.2.ppt

Computational Molecular Evolution 计算分子进化,杨子恒（Ziheng Yang）/著,第1部分 分子进化建模,第1章 核苷酸置换模型 第2章 氨基酸和密码子置换模型,2.5.1.4 rbcL基因应用实例,应用NG86方法来估计黄瓜和烟草中叶绿体蛋白1,5-二磷酸核 酮糖羧化酶/加氧酶大亚基(rbcL)基因间的 和 图2.5列举了数 据的一些基本统计值表2.5 黄瓜和烟草rbcL基因的基本统计量,第 1 章 核苷酸置换模型,1.1 引言 1.2 核苷酸置换和距离估计的马尔可夫模型 1.3 位点间可变的置换率 1.4最大似然估计 1.5马尔可夫链和广义模型下的距离估计,1.1引言,序列间的距离计算重要性 成对距离的计算是系统发育重建中距离矩阵方法的首要步骤，他采用聚类算法将距离矩阵转化为一个系统发育树； 用于距离计算的核苷酸置换马尔可夫过程模型构成了针对多重序列系统发育关系的似然分析与贝斯分析的基础两条序列间的距离被定义为平均每个位点核苷酸置换的期望数,简化的距离测度，即差异位点比例有时称为P距离 例 长度为100个核苷酸的两条序列间有10个位点差异，则p=10%=0.1 缺点：只适合关系非常近的序列，如p5%. 多重置换：一个可变位点可能是一次以上置换的结果，甚至一个不变位点也可能经历过回复（harbour）或平行（parallel）置换。

如下图,图1.1 当前两序列间只观察到两个差异位点，因而差异位点的比例值p=2/8=0.25， 而事实上，根据图中可得到真是的距离为平均每个位点10/8=1.25次置换核苷酸置换的马尔可夫模型,为了估计置换数目，需要建立一个概率模型来描述核苷酸间的变化而时间连续的马尔可夫链常被用于此目的 假定：核苷酸位点彼此独立进化， 任何特定位点上的置换都可以用一个马尔可夫链来描述，其中链状态为4个核苷酸 {T C A G}. 主要特性：无记忆性，给定一个当前状态，未来状态并不依赖于过去的状态换言之，链上核苷酸跳到其他核苷酸的概率依赖于当前状态，但与当前状态从何而来无关 对核苷酸间的置换率加以约束由此可以产生不同的核苷酸置换模型图1.2在核苷酸置换的3个马尔可夫链模型下核苷酸的相对置换率，这3个模型是JC69,K80和HKY85.线条粗细表示置换率，而圆圈大小表示稳定状态的分布1.2核苷酸置换和距离估计的马尔科夫模型,JC69模型（Juckes andCantor,1969） 假定每个核苷酸变成其他任何一个核苷酸都是相同速率 λ我们用 表示核苷酸i变成核苷酸j时的瞬时速率其 中（i,j=T,C,A,G）,则置换率矩阵为,为了描述马尔可夫链，对任意时间t0,我们需要一个相似的概率，这就是转换概率 即给定核苷酸i在时间t之后变成j的概率。

矩阵 则称为转换概率矩阵设想一条每个位点为核苷酸i的长序列，我们设一个时期t内每个位点都进化，则对j=T,C,A,G,核苷酸j在系列中的比例将为,我们将转换概率矩阵的两个数值不同的元素 和 标在图1.3上首先，P(t)在每行的总和为1； 其次，单位矩阵P(0)=Ι反映了无进化的情形（t=0）； 第三，速率λ和时间t只以一个乘积的形式λt在转换概率中出现； 最后，当t →∞时，对所有的i和j，有 图1.3,极限分布（limiting distribution）,t →∞时，链在状态j的概率用 表示，分布 是一个极限分布（limiting distribution） 若链的状态已经是极限分布，该链将继续保持该分布， 故极限分布也是一种“稳态分布”（steady-state distribution） 或“平稳分布”（stationary distribution）. 马尔可夫链模型是通过公式 计算转换概率得到 的，该公式考虑了可能经历的进化过程中的所有可能的路 径，尤其是一个马尔可夫链的转换概率满足查普曼-科尔莫哥 洛夫定理（C-K方程）,i ↓t1 k=(T,C,A,G) ↓t2 j,任意一个核苷酸i经过时间t1+t2变为核苷酸j的转换概率 是在一个中间时间点t1时所有可能状态k的总和。

图1.4 查普曼-科尔莫哥洛夫定理的图示1.4）,估计距离,基于公式（1.1）得知，任意核苷酸的总置换率为3λ，若两条序列在时间t分开，则这两条序列间的距离为d=3λt 假设两条序列间n个位点中有x个差异，则差异位点的比例为 基于公式（1.3）子裔序列不同于祖先序列的核苷酸概率为 将其换算为观测比例 ，我们得距离估计值,（1.5）,（1.6）,这里，对数底为常数为e如果 ，则该距离公式不可用，两条随机序列应该有约75%的差异位点；如果 ，估计距离将为无穷 大是变量 的一个二项式比例值为了推断的方差，我们将 看作是 的函数并采用Δ技术（见附录），得,（1.7）,1.2.2 K80模型,转换（transition）：两个嘧啶（T↔C）或者两个嘌呤（A↔G）间的置换； 颠换（transvertion）:一个嘧啶和一个嘌呤（T,C↔A,G）间的置换 设转换率为α，颠换率为β，该模型称为K80模型（亦称 Kimura双参数模型）速率矩阵如下,（1.8）,任意一个核苷酸的总置换率为α+2β，两条序列经过时间t的距离为d= (α+2β). αt是平均每个位点的期望转换数目，2βt是平均每个位点的期望颠换数目。

人们可以将αt和βt看作模型的两个参数，但一般直接用距离d及转换、颠换比率κ=α/β更为方便转换概率矩阵为,（1.9）,这里矩阵的3个不同的元素为,,注意：,序列数据能被概括为转换（S）与颠换（V）差异位点的比例值根据模型的对称性（公式（1.9）），一个核苷酸位点上出现一个转换差异的概率为E(S)= , 和E(V)= 对观测比例值S和V导出两个容易求解的公式：,（1.11）,相应地，转换距离αt和颠换距离2βt的估计值为,（1.12）,该距离公式当且仅当1-2S-V0和1-2V0时适用S和V为多项式比例值，有Var(S)=S(1-S)/n,Var(V)=V(1-V)/n以及 Cov(S,V)=-SV/n,我们可以采用Δ技术来推导 和 变量-协 方差矩阵特别地， 的方差为,（1.13）,（1.14）,1.2.3HKY85、F84、TN93等模型 首先介绍TN93模型,Jukes和Cantor与Kimura的模型都有对称的置换率，即对 所有的i和j, ，这样的马尔可夫链对所有的i都有 的固定分布，即当置换过程打到平衡时，序列中4个核苷酸 的比例将趋于一致这一假设对于实际的每一套真是数据都 不成立。

这里介绍TN93模型下的置换率矩阵为,当置换率用参数 ， ， ， 来确定时，它们也给出稳态分布，其嘧啶和嘌呤的频率分布为,Q能写成 其中，U为非奇异矩阵，且有 λ为Q的特征值，有 得,（1.16）,（1.18）,（1.19）,将 代入公式(1.17)中，得,这里，,当t从0→∞时随着 ，对角线元素 从1降到 ，而非对角线元素 从0增加到 ，与起始核苷酸i无关现在考虑基于该模型的序列距离估计,首先，定义距离核苷酸i的置换率为 ，且在4个核苷酸状态不同当置换过程打到平衡时，马尔可夫链在4个状态T,C,A,G中所化的时间总和分别与平衡频率 成比例因此，平均置换率为 设S1为两个不同嘧啶出现的比例，S2为两个不同嘌呤出现的比例V为位点出现颠换差异的比例由 则有,(1.21),(1.22),(1.23),可得到如下估计值,,,这里,。