计算分子进化第四章
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
f ( X h ) x0 L0 ( x0 )
x0
0
(4.5)
x 是根的核苷酸为 x 的(先验)概率，由该模型下核苷酸平衡频率 0
给出。(具体可见书中 P105例) 修剪算法是一种主要的节省计算时间的方法。
8
4.2.3 时间可逆性、树根
分子系统发育分析中使用的大多数置换模型描述时间可逆马尔 可夫链。对于这类链，转换概率对任意i,j和t,满足 i pij (t ) j p ji (t ) 。 可逆性的一个重要结果是根可以沿着树任意移动而不影响似然率, Felsenstein(1981)称之为滑轮原理。例如，对公式(4.2)中的 x px x (t6 ) 置换 x px x (t6 ) ，注意
0 6 7 t7
t6
t8 t3
8
t1
1:T
t2
t4
4:C
t5
2:C 3:A
图4.1 一个用来演示似然函数计算的五物种树。顶端 示出所观测的一个位点上的核苷酸，枝长 t1 t8 的测度 5:C 是每位点核苷酸置换的期望数。
3
观测数据为TCACC,祖先节点用数字0,6,7,8标出，0为树根。连 到节点i的分支长度用 ti 表示，定义为平均每个位点核苷酸置换的期 望数目。 模型中的参数包括枝长和转换/颠换比率 ，用集合表示 t1, t2 , t3 , t4 , t5 , t6 , t7 , t8 , 由于假定位点间独立进化，整个数据集的概率是数据在单个位点的 概率之积，对数似然率则等价于对序列中的所有位点求和
r 1e r g (r; , ) ( )
(4.11)
具有均值α/β及方差 2 ，令α=β，使得均值为1。性状参数 与位点间置换率变异的大小有关。 如同在离散置换率模型中的一 样，我们不知道该位点的置换率，而不得不对置换率分布求平均
f (x h ) g (r ) f (x h r ; )dr
f (r x h ; ) f (r ) f (x h r ; ) f (x h )
(4.13)
参数θ可用其估计值代替(如最大似然估计)，已知这是一种经验贝 斯方法。在连续置换率模型下，人们可用后验均值作为该位点置 换率的估计值(Yang and Wang,1995)。在离散置换率模型下，公 式(4.13)中的置换率r取K可能值中的一个，有 f (r ) pk。人们可以 计算后验均值或者用具有最高后验概率的置换率作为最佳估计值 (Yang,1994a)。
Li ( xi ) pxi x j (t j ) L j ( x j ) pxi xk (tk ) Lk ( xk ) xj xk
(4.4)
这是对应于两个子节点j和k的两项的乘积。
7
该计算过程从树的顶端到根部，每个节点的所有后裔节点都被 访问过后再访问该节点，即后序树遍历(post-order tree traversal) 方法。在访问树上所有节点并计算根的概率向量之后，该位点的 数据概率为
6 6 0
0
0 6
p
我们有
x6 x0
(t6 ) px0 x8 (t8 ) px6 x8 (t6 t8 )
f (X h ) [ x6 px6 x7 (t7 ) px6 x8 (t6 t8 ) px7T (t1 )
x6 x7 x8
px7 C (t2 ) px6 A (t3 ) px8C (t4 ) p x8C (t5 )]
0
5
4.2.2 修剪算法
对祖先状态的所有组合求和很费力，在计算这类求和时，有一 种重要的技术很有用，它找出公因子且只对它们计算一次，这就 是嵌套法则或称Horner法则。对公式(4.2)应用嵌套法则意味着将 加号尽可能远地向右移，得
f (X h ) x0 px0 x6 (t6 ) px6 x7 (t7 ) px7T (t1 ) px7C (t2 )) px6 A (t3 ) x0 x7 x6 px0 x8 (t8 ) px8C (t4 ) px8C (t5 ) x8
p0 p1 f (x h r r1 ; ) f (x h ) p1 f (x h r r1; )
假如位点是恒定的 假如位点是可变的 (4.10)
14
4.3.1.2 伽马置换率模型
用一个连续分布去近似位点的可变置换率。最常用的分布是伽 马分布。伽马密度为
(4.6)
9
如果根在节点6，两个分支0-6和0-8合并到一个分支6-8，枝 长 t6 t8。结果树如图4.3所示 t6 t8 公式(4.6)也明示出图4.1中的模型 是过度参数化的，多了一个枝长。 6 t7 对 t6 和 t8 的任意组合，只要 t6 t8 7 8 相同，似然率就相同。数据 t4 t t1 t t3 2 5 不包含分别估计t6 和 t8 信息， 1:T 2:C 3:A 4:C 5:C 只有它们之和是可以估计的。 图4.3 滑轮原理也可用于简化对一些小树在理论研究中的似然率计 算。例如图4.4(a)中树的似然率计算设计对两个祖先节点的祖先状 态i和j求和。然而，将根移到物种1和物种2的共有祖先更为简便， 这样人们只要对一个位点的祖先状态求和。在仅仅一个位点上， 数据 x1 x2 x3 的概率则成为
由于我们并不知道每个位点属于哪一类置换率，因而在任一 位点观测数据的概率是对位点类别的加权平均
f (x h ) pk f (x h r rk ; )
k
(4.8)
似然率还是用位点间概率的乘积来计算。 由置换率r给出的观测数据 x h 的条件概率 f (x r ; ) 只是单置换 率模型下的概率，用所有枝长乘以r。而具有K类置换率的可变置 换率模型就是单置换率模型计算K次。 考虑图4.4(b)中的树，我们有
第四章 最大似然法
4.1 引言 4.2 树的似然计算 4.3 复杂模型下的似然计算 4.4 祖先状态重建 4.5 最大似然估计的数值算法 4.6 似然法的近似 4.7 模型选择与稳健性
1
4.1 引言
在本章中讨论一棵系统发育树上的多重序列似然计算。 系统发育分析中似然率的两种用途： 第一种，估计进化模型中的参数； 第二种，估计树拓扑结构，通过估计枝长和其他置换参数来实 现每棵树的对数似然率最大化，并将优化的对数似然率用做比较 不同树的树分值。
h
f ( x1 x2 x3 r ; ) j p jx1 (t1r ) p jx2 (t1r ) p jx3 [(2t0 t1 )r ]
j
(4.9)
但该模型被认为K不应超过3或4。Yang(1995a)实现了表4.1的普 适模型。
13
一个特例是不可变位点模型，它假设两类位点：一类为不可 变位点类，其置换率为 r0 0 ；另一类具有一个恒定置换率 r1 。由于 平均置换率为1，有r1 1 (1 p0 ) ，这里的不可变位点p0 比例是模型中 仅有的参数。注意，在对位排列中，一个可变位点不能有 r0 0 的 置换率，因而一个位点的数据概率为
2
4.2 树的似然计算
4.2.1 数据、模型、树及似然函数
已知s个已对位排列的同源序列，每个序列为n个核苷酸长，用 s×n阶矩阵X={x jh }表示，其中x jh 表示第j个序列中的第h个核苷酸. 令 X h 表示数据矩阵中的第h列，这里我们采用K80核苷酸置换 模型。下面我们用一个例子来解释似然率计算。如图
l log( L) log{ f (X h )}
h 1 n
(4.1)
ML方法通过对数似然率 l 最大化来估计 ,通常采用数值优化算法。 这里,我们考虑参数 给定时 l 的计算。
4
我们集中考虑一个位点。用 xi 表示祖先节点i的状态，标注下 标h,由于该位点的数据 (x h TCACC) 可由祖先核苷酸 x0 , x6 , x7 , x8 的任 意组合产生，计算 f (x h ) 就必须对已灭绝祖先的所有可能的核苷酸 组合求和
10
f ( x1 x2 x3 ) i pij (t0 t1 ) p jx1 (t1 ) p jx2 (t1 ) pix3 (t0 )
i j
j p jx1 (t1 ) p jx2 (t1 ) p jx3 (2t0 t1 )
j
(4.7)
t0
i
2t0 t1
j j
0
(4.12)
这里，参数 的集合包括 和枝长以及置换模型中的其他参数 (如 )。Yang(1993)、Gu等(1995)以及Kelly和Rice(1996)介绍了 连续置换率模型(如公式(4.12))下计算似然函数的一种算法。
15
位点置换率的贝斯估计
由许多序列组成的大型数据集使得我们可能估计每个位点的相 对置换率。在离散和连续置换率模型中，位点置换率都是从似然 函数中积分获得的随机变量。下面用位点数据给出的置换率条件 (后验)分布来估计置换率
f (X h ) [ x0 px0 x6 (t6 ) px6 x7 (t7 ) px7T (t1 )
x0 x6 x7 x8
(4.2)
px7C (t2 ) px6 A (t3 ) px0 x8 (t8 ) px8C (t 4 ) px8C (t5 )]
方括号中的各项为顶点数据TCACC和祖先节点数据 x0 , x6 , x7 , x8 的概 率，等于树根(节点)具有 x0 的概率，它由K80模型下 x 1 4 时树 上8个枝的8个转换概率相乘得到。
(4.3)
因而，我们在 x6 之前加 x7 ，在x0 之前对x6 和 x8 求和。换而言之， 我们对一个节点的祖先状态求和仅仅是在计算其所有后裔节点之 后。
6
用公式(4.3)计算f (x h ) 就形成了Felsenstein(1973b,1981)的修 剪法，它是动态规划算法的一种变体，其精髓是逐步计算许多子 树上的数据概率。下面介绍条件概率 Li ( xi ) ，给定节点i上的核苷酸 为 xi ，定义条件概率为节点i的子裔顶点上观测数据的概率。例如, 顶点1,2,3是节点6的子裔，设节点6具有状态x6 T ,则 L6 (T ) 为观测 到 x1, x2 , x3 =TCA的概率。 x 若节点i为一个顶点，其后裔点点仅包括节点i，i 为观测到的核 苷酸，有Li ( xi ) 1，否则为0。若节点i是一个具有子节点j和k的内节 点,我们有
16
4.4祖先状态重建
4.4.1 概述 4.4.2经验和等级贝斯重建 4.4.2.1边界重建(marginal reconstruction) 我们计算一个祖先节点的形状状态的后验概率。图4.2中考 虑根节点0，其核苷酸为 x0 ，给定位点数据，节点0的后验概率为
x0 L0 ( x0 ) f ( x0 ) f (x h x0 ; ) f ( x0 x h ; ) f (x h ) x0 L0 ( x0 )
表4.1 离散置换率模型 1
p1
位点类别 概率 置换率
2
p2
3
p3
…… …… ……
K
pk
r1
r2
r3
rk
两个约束，频率之和为1， pk 1 。其次，平均置换率固定为1， 即 pk rk 1 ，使得枝长可以用平均每个位点的核苷酸置换期望数 目来度量。后一个约束对于避免使用过多的参数是必要的。 12
t1
t1
t1
x1
x2
(a)
x3
x1
源自文库x2
(b)
x3
图4.4 (a)在一个三物种有根树上，似然计算必须对两个祖先节点上的祖先状态 i和j求和。(b)通过将树根移到物种1和物种2的祖先并对新树根上状态j求和可以 完成相同的计算。
11
4.3 复杂模型下的似然计算
4.3.1 位点间可变置换率模型 在实际序列中置换率在位点间常常是可变的，而忽略位点间 的置换率变异可能对系统发育分析产生重要影响。 4.3.1.1 离散置换率模型 在这个模型中，假定位点中具有不同置换率的若干类(如K类). (见表4.1).
x0
0
(4.5)
x 是根的核苷酸为 x 的(先验)概率，由该模型下核苷酸平衡频率 0
给出。(具体可见书中 P105例) 修剪算法是一种主要的节省计算时间的方法。
8
4.2.3 时间可逆性、树根
分子系统发育分析中使用的大多数置换模型描述时间可逆马尔 可夫链。对于这类链，转换概率对任意i,j和t,满足 i pij (t ) j p ji (t ) 。 可逆性的一个重要结果是根可以沿着树任意移动而不影响似然率, Felsenstein(1981)称之为滑轮原理。例如，对公式(4.2)中的 x px x (t6 ) 置换 x px x (t6 ) ，注意
0 6 7 t7
t6
t8 t3
8
t1
1:T
t2
t4
4:C
t5
2:C 3:A
图4.1 一个用来演示似然函数计算的五物种树。顶端 示出所观测的一个位点上的核苷酸，枝长 t1 t8 的测度 5:C 是每位点核苷酸置换的期望数。
3
观测数据为TCACC,祖先节点用数字0,6,7,8标出，0为树根。连 到节点i的分支长度用 ti 表示，定义为平均每个位点核苷酸置换的期 望数目。 模型中的参数包括枝长和转换/颠换比率 ，用集合表示 t1, t2 , t3 , t4 , t5 , t6 , t7 , t8 , 由于假定位点间独立进化，整个数据集的概率是数据在单个位点的 概率之积，对数似然率则等价于对序列中的所有位点求和
r 1e r g (r; , ) ( )
(4.11)
具有均值α/β及方差 2 ，令α=β，使得均值为1。性状参数 与位点间置换率变异的大小有关。 如同在离散置换率模型中的一 样，我们不知道该位点的置换率，而不得不对置换率分布求平均
f (x h ) g (r ) f (x h r ; )dr
f (r x h ; ) f (r ) f (x h r ; ) f (x h )
(4.13)
参数θ可用其估计值代替(如最大似然估计)，已知这是一种经验贝 斯方法。在连续置换率模型下，人们可用后验均值作为该位点置 换率的估计值(Yang and Wang,1995)。在离散置换率模型下，公 式(4.13)中的置换率r取K可能值中的一个，有 f (r ) pk。人们可以 计算后验均值或者用具有最高后验概率的置换率作为最佳估计值 (Yang,1994a)。
Li ( xi ) pxi x j (t j ) L j ( x j ) pxi xk (tk ) Lk ( xk ) xj xk
(4.4)
这是对应于两个子节点j和k的两项的乘积。
7
该计算过程从树的顶端到根部，每个节点的所有后裔节点都被 访问过后再访问该节点，即后序树遍历(post-order tree traversal) 方法。在访问树上所有节点并计算根的概率向量之后，该位点的 数据概率为
6 6 0
0
0 6
p
我们有
x6 x0
(t6 ) px0 x8 (t8 ) px6 x8 (t6 t8 )
f (X h ) [ x6 px6 x7 (t7 ) px6 x8 (t6 t8 ) px7T (t1 )
x6 x7 x8
px7 C (t2 ) px6 A (t3 ) px8C (t4 ) p x8C (t5 )]
0
5
4.2.2 修剪算法
对祖先状态的所有组合求和很费力，在计算这类求和时，有一 种重要的技术很有用，它找出公因子且只对它们计算一次，这就 是嵌套法则或称Horner法则。对公式(4.2)应用嵌套法则意味着将 加号尽可能远地向右移，得
f (X h ) x0 px0 x6 (t6 ) px6 x7 (t7 ) px7T (t1 ) px7C (t2 )) px6 A (t3 ) x0 x7 x6 px0 x8 (t8 ) px8C (t4 ) px8C (t5 ) x8
p0 p1 f (x h r r1 ; ) f (x h ) p1 f (x h r r1; )
假如位点是恒定的 假如位点是可变的 (4.10)
14
4.3.1.2 伽马置换率模型
用一个连续分布去近似位点的可变置换率。最常用的分布是伽 马分布。伽马密度为
(4.6)
9
如果根在节点6，两个分支0-6和0-8合并到一个分支6-8，枝 长 t6 t8。结果树如图4.3所示 t6 t8 公式(4.6)也明示出图4.1中的模型 是过度参数化的，多了一个枝长。 6 t7 对 t6 和 t8 的任意组合，只要 t6 t8 7 8 相同，似然率就相同。数据 t4 t t1 t t3 2 5 不包含分别估计t6 和 t8 信息， 1:T 2:C 3:A 4:C 5:C 只有它们之和是可以估计的。 图4.3 滑轮原理也可用于简化对一些小树在理论研究中的似然率计 算。例如图4.4(a)中树的似然率计算设计对两个祖先节点的祖先状 态i和j求和。然而，将根移到物种1和物种2的共有祖先更为简便， 这样人们只要对一个位点的祖先状态求和。在仅仅一个位点上， 数据 x1 x2 x3 的概率则成为
由于我们并不知道每个位点属于哪一类置换率，因而在任一 位点观测数据的概率是对位点类别的加权平均
f (x h ) pk f (x h r rk ; )
k
(4.8)
似然率还是用位点间概率的乘积来计算。 由置换率r给出的观测数据 x h 的条件概率 f (x r ; ) 只是单置换 率模型下的概率，用所有枝长乘以r。而具有K类置换率的可变置 换率模型就是单置换率模型计算K次。 考虑图4.4(b)中的树，我们有
第四章 最大似然法
4.1 引言 4.2 树的似然计算 4.3 复杂模型下的似然计算 4.4 祖先状态重建 4.5 最大似然估计的数值算法 4.6 似然法的近似 4.7 模型选择与稳健性
1
4.1 引言
在本章中讨论一棵系统发育树上的多重序列似然计算。 系统发育分析中似然率的两种用途： 第一种，估计进化模型中的参数； 第二种，估计树拓扑结构，通过估计枝长和其他置换参数来实 现每棵树的对数似然率最大化，并将优化的对数似然率用做比较 不同树的树分值。
h
f ( x1 x2 x3 r ; ) j p jx1 (t1r ) p jx2 (t1r ) p jx3 [(2t0 t1 )r ]
j
(4.9)
但该模型被认为K不应超过3或4。Yang(1995a)实现了表4.1的普 适模型。
13
一个特例是不可变位点模型，它假设两类位点：一类为不可 变位点类，其置换率为 r0 0 ；另一类具有一个恒定置换率 r1 。由于 平均置换率为1，有r1 1 (1 p0 ) ，这里的不可变位点p0 比例是模型中 仅有的参数。注意，在对位排列中，一个可变位点不能有 r0 0 的 置换率，因而一个位点的数据概率为
2
4.2 树的似然计算
4.2.1 数据、模型、树及似然函数
已知s个已对位排列的同源序列，每个序列为n个核苷酸长，用 s×n阶矩阵X={x jh }表示，其中x jh 表示第j个序列中的第h个核苷酸. 令 X h 表示数据矩阵中的第h列，这里我们采用K80核苷酸置换 模型。下面我们用一个例子来解释似然率计算。如图
l log( L) log{ f (X h )}
h 1 n
(4.1)
ML方法通过对数似然率 l 最大化来估计 ,通常采用数值优化算法。 这里,我们考虑参数 给定时 l 的计算。
4
我们集中考虑一个位点。用 xi 表示祖先节点i的状态，标注下 标h,由于该位点的数据 (x h TCACC) 可由祖先核苷酸 x0 , x6 , x7 , x8 的任 意组合产生，计算 f (x h ) 就必须对已灭绝祖先的所有可能的核苷酸 组合求和
10
f ( x1 x2 x3 ) i pij (t0 t1 ) p jx1 (t1 ) p jx2 (t1 ) pix3 (t0 )
i j
j p jx1 (t1 ) p jx2 (t1 ) p jx3 (2t0 t1 )
j
(4.7)
t0
i
2t0 t1
j j
0
(4.12)
这里，参数 的集合包括 和枝长以及置换模型中的其他参数 (如 )。Yang(1993)、Gu等(1995)以及Kelly和Rice(1996)介绍了 连续置换率模型(如公式(4.12))下计算似然函数的一种算法。
15
位点置换率的贝斯估计
由许多序列组成的大型数据集使得我们可能估计每个位点的相 对置换率。在离散和连续置换率模型中，位点置换率都是从似然 函数中积分获得的随机变量。下面用位点数据给出的置换率条件 (后验)分布来估计置换率
f (X h ) [ x0 px0 x6 (t6 ) px6 x7 (t7 ) px7T (t1 )
x0 x6 x7 x8
(4.2)
px7C (t2 ) px6 A (t3 ) px0 x8 (t8 ) px8C (t 4 ) px8C (t5 )]
方括号中的各项为顶点数据TCACC和祖先节点数据 x0 , x6 , x7 , x8 的概 率，等于树根(节点)具有 x0 的概率，它由K80模型下 x 1 4 时树 上8个枝的8个转换概率相乘得到。
(4.3)
因而，我们在 x6 之前加 x7 ，在x0 之前对x6 和 x8 求和。换而言之， 我们对一个节点的祖先状态求和仅仅是在计算其所有后裔节点之 后。
6
用公式(4.3)计算f (x h ) 就形成了Felsenstein(1973b,1981)的修 剪法，它是动态规划算法的一种变体，其精髓是逐步计算许多子 树上的数据概率。下面介绍条件概率 Li ( xi ) ，给定节点i上的核苷酸 为 xi ，定义条件概率为节点i的子裔顶点上观测数据的概率。例如, 顶点1,2,3是节点6的子裔，设节点6具有状态x6 T ,则 L6 (T ) 为观测 到 x1, x2 , x3 =TCA的概率。 x 若节点i为一个顶点，其后裔点点仅包括节点i，i 为观测到的核 苷酸，有Li ( xi ) 1，否则为0。若节点i是一个具有子节点j和k的内节 点,我们有
16
4.4祖先状态重建
4.4.1 概述 4.4.2经验和等级贝斯重建 4.4.2.1边界重建(marginal reconstruction) 我们计算一个祖先节点的形状状态的后验概率。图4.2中考 虑根节点0，其核苷酸为 x0 ，给定位点数据，节点0的后验概率为
x0 L0 ( x0 ) f ( x0 ) f (x h x0 ; ) f ( x0 x h ; ) f (x h ) x0 L0 ( x0 )
表4.1 离散置换率模型 1
p1
位点类别 概率 置换率
2
p2
3
p3
…… …… ……
K
pk
r1
r2
r3
rk
两个约束，频率之和为1， pk 1 。其次，平均置换率固定为1， 即 pk rk 1 ，使得枝长可以用平均每个位点的核苷酸置换期望数 目来度量。后一个约束对于避免使用过多的参数是必要的。 12
t1
t1
t1
x1
x2
(a)
x3
x1
源自文库x2
(b)
x3
图4.4 (a)在一个三物种有根树上，似然计算必须对两个祖先节点上的祖先状态 i和j求和。(b)通过将树根移到物种1和物种2的祖先并对新树根上状态j求和可以 完成相同的计算。
11
4.3 复杂模型下的似然计算
4.3.1 位点间可变置换率模型 在实际序列中置换率在位点间常常是可变的，而忽略位点间 的置换率变异可能对系统发育分析产生重要影响。 4.3.1.1 离散置换率模型 在这个模型中，假定位点中具有不同置换率的若干类(如K类). (见表4.1).