（概率论与数理统计专业论文）基于随机微分方程和结构em算法的系统发生树的构建.pdf

基于随机微分方程和结构e m 算法的系统发生树的构建中文摘要 中文摘要 在生命科学中，人们不断地探索物种之间的进化关系，而系统发生树成为描述 这种进化关系最好的手段之一随着许多物种的基因测序工程的完成和生物技术 的发展，人类拥有了大量的生物数据并且发现一些同源基因能够表征物种间的特 征所以人们又将系统发生树的构建着眼于物种的基因序列2 0 世纪8 0 年代以来， 许多分子水平上的系统发生树的构建方法相继提出，这些方法大多数是基于合理 基因序列之间距离的定义和在距离之下作出一定的假设来构建系统发生树 本文第三章考虑进化除了主要受进化自身规律的影响，还受各种内、外部因 素的影响，这种影响使得进化表现出某种波动特性，我们近似地认为这种波动是一 个b r o w n 运动所以我们借用随机微分方程的手段构建了序列问进化的随机微分方 程的模型在对扩散函数，漂移函数作出限定后，我们用g i b b s 抽样求解该随机微分 方程 第四章根据中性进化理论，假设序列间的进化是相互独立的，在此基础上我们 考虑这样的问题：如果对n 条d n a 序列，我们只知道部分序列之间的距离，如何寻 求这n 条d n a 序列的最优系统发生树我们借用结构e m 算法的理论，给出构建系 统发生树的结构尉算法，并且我们给出算法中q 函数的处理方法 关键词：系统发生树，随机微分方程，g i b b s 抽样，结构e i 算法 作者：陈爱忠 指导老师：汪四水( 副教授) 基于随机微分方程和结构e m 算法的系统发生树的构建 英文摘要 p h y l o g e n e t i ct r e e s c o n s t r u c t i o nb a s e do ns t o c h a s t i c d i f f e r e n t i a le q u a t i o na n ds t r u c t u r a le ma l g o r i t h m a b s t r a c t i nl i f e s c i e n c e s p e o p l ec o n s t a n t l ye x p l o r ee v o l u t i o n a r yr e l a t i o n s h i p sa m o n g s p e c i e s t h ep h y l o g e n e t i ct r e ei so n eo f t h eb e s tm e t h o d si na l lm e t h o d sw h i c hd e s e r i b e t h i se v o l u t i o n a r yr e l a t i o n s h i p s w i t ht h ec o m p l e t i o no fg e n o m e ss e q u e n c i n go fm a n y s p e c i e sa n dt h ea d v a n c eo fm i c r o a r r a yt e c h n o l o g y , w eb e g i nt op o s s e s saw e a l t ho f v a l u a b l eb i o l o g i c a ld a t a ，a n dd i s c o v e rt h a ts o m eh o m o m o r p h i cg e n e sc a nr e p r e s e n tt h e s p e c i e s f e a t u r e s oa tt h es a m et i m ep e o p l ec o n s t r u c tt h ep h y l o g e n e t i ct r e e sb a s e do n t h eg e n es e q u e n c e so fs p e c i e s f r o m19 8 0 s ，m a n ym e t h o d so fp h y l o g e n e t i ct r e e s c o n s t r u c t i o ni nm o l e c u l el e v e lh a v eb e e nr a i s e d m o s to fm e t h o d sa r eb a s e do nt h e r e a s o n a b l ed e f i n i t i o no fd i s t a n c e sa m o n gg e n es e q u e n c e sa n dc e r t a i nh y p o t h e s i su n d e r t h ed i s t a n c e i nc h a p t e r3 t h ee v o l u t i o ni sm a i n l ya f f e c t e dn o to n l yb yt h et h ee v o l u t i o nl a w i t s e l f , b u ta l s ob yav a r i e t yo fi n t e m a la n de x t e r n a lf a c t o r s t h i sa f f e c tm a k e st h a t e v o l u t i o ne x p r e s s e ss o m ef l u c t u a t i o nc h a r a c t e r i s t i c s w ea p p r o x i m a t i v e l yc o n s i d e rt h a t t h i sk i n do ff l u c t u a t i o ni sb r o w n i a nm o t i o n t h e r e f o r e ，w eu s et h em e a n so fs t o c h a s t i c d i f f e r e n t i a le q u a t i o n ( s d e ) t ob u i l das d em o d e lb e t w e e nt h ee v o l u t i o no fs e q u e n c e s a f t e rs e t t i n gal i m i tt ot h ed r i f ta n dd i f f u s i o nf u n c t i o n s ，w eg i v et h es o l u t i o no ft h i s e q u a t i o nb yg i b b ss a m p l i n g i nc h a p t e r4 ，b a s e do nt h eh y p o t h e s i so fn e u t r a le v o l u t i o n a r y t h e o r ya n dt h e e v o l u t i o n sb e t w e e ns e q u e n c e sb e i n gi n d e p e n d e n to fe a c ho t h e r , w ec o n s i d e rt h e f o l l o w i n gq u e s t i o n ：i fw eo n l yk n o wap a r to ft h ed i s t a n c e sa m o n gs e q u e n c e sf o r 刀d n as e q u e n c e s h o wd ow ec o n s t r u c tt h em o s tv a l u a b l ep h y l o g e n e t i ct r e e w eg i v e t h es t r u c t u r a le ma l g o r i t h mo ft h ec o n s t r u c t i o no fp h y l o g e n e t i ct r e e sa f t e rm a k i n gu s e o ft i l et h e o r yo fs t r u c t u r a le ma l g o r i t h ma n di n t r o d u c et h em e t h o dt od e a lw i t hq f u n c t i o ni nt h i sa l g o r i t h m k e y w o r d s ：p h y l o g e n e t i ct r e e ，s t o c h a s t i cd i f f e r e n t i a le q u a t i o n ，g i b b ss a m p l i n g ，e m a l g o r i t h m w r i t t e nb yc h e na i z h o n g s u p e r v i s e db ya s s o c i a t ep r o f w a n gs i s h u i 苏州大学学位论文独创性声明及使用授权声明 学位论文独创性声明 本人郑重声明：所提交的学位论文是本人在导师的指导下，独立 进行研究工作所取得的成果。除文中已经注明引用的内容外，本论文。 不含其他个人或集体已经发表或撰写过的研究成果，也不含为获得苏 州大学或其它教育机构的学位证书而使用过的材料。对本文的研究作 出重要贡献的个人和集体，均已在文中以明确方式标明。本人承担本 声母目的法律责任。 研究生签名：蛰：盈惠。，e l期：翌2 ：生皇 学位论文使用授权声明 苏州大学、中国科学技术信息研究所、国家图书馆、清华大学论 文合作部、中国社科院文献信息情报中心有权保留本人所送交学位论 文的复印件和电子文档，可以采用影印、缩印或其他复制手段保存论 文。本人电子文档的内容和纸质论文的内容相一致。除在保密期内的 保密论文外，允许论文被查阅和借阅，可以公布( 包括刊登) 论文的 全部或部分内容。论文的公布( 包括刊登) 授权苏州大学学位办办理。 研究生签名：，螫：蛊如e t 期：迎2 ：垒：垒 导师签名：趁空生e t期：丝翌：! ：芝 基于随机微分方程和结构e m 算法的系统发生树的构建第一章序言 第一章序言 1 1 生物背景知识的介绍 人们对生物的认识从宏观发展到微观，科学家对物种分类的依据也从宏观的 形态发展到微观上的分子近几年，随着不同的分子测序技术的飞速发展，系统发 生分析也进入分子层次，并且大量的d n a 分子数据的不断涌现，这给生物学家提供 了大量的数据，使其能够实现重构地球上所有生命的系统发生树的梦想 系统发生树是一个二叉树所谓树，实际上是一个无向非循环图系统发生树由 一系列的节点和分支构成，其中每个节点代表一个分类单元，而节点之间的连线( 分 支) 代表物种或序列之间的进化关系树的节点又分为外部节点和内部节点系统发 生树一般分为有根树和无根树 如图l ： n y ( a )( b ) 图l ：( a ) 为有根树，( b ) 为无根树 图l ( a ) 是一有根树，( b ) 是一无根树，a ，b ，c ，d ，e 称为外部节点，代表物种，或 者d n a 序列，我们也可称为叶；n 1 ，n 2 ，n 3 称为内部节点，节点之间的连线称为树的 枝，有时枝上也可以赋值( 我们称为距离) 我们有如下的结论： 若树的外部节点有n ( n 3 ) 个，则建立的有根的系统发生树的数目为： 基于随机微分方程和结构e m 算法的系统发生树的构建第一章序言 删( 刀) = 丽( 2 n - 3 ) ! ( 1 1 1 ) 而无根树的数目为： = 罴焉 n 坨， 下面我们用数学归纳法对这两个式子进行证明： 首先我们先明白下面的事实，对于有根树而言，当叶的数目为3 时，则建立的树 的节点个数为5 ，分支的个数为4 ；当叶的数目增加一个时，树的节点将增加2 个，同 样分支的个数增加2 枝，所以我们可以归纳出物种与节点和分支的关系： 当物种数为 时，有根树的节点个数为2 n 一1 ，分支个数为2 n 一2 同理 当物种数为刀时，无根树的节点个数为2 n 一2 ，分支个数为2 n 一3 下面我们对( 1 1 1 ) 式进行证明： 当刀= 3 时，有根树的形式有图2 中的3 种：设物种为a ，b ，c c a b 图2 ：三个物种的所有有根树的形式 将玎= 3 代入( 1 1 1 ) 式的右边，得3 所以刀= 3 时，( 1 1 1 ) 式成立 假设珂= 七时也成立，即 m 炉筹盖 这时每棵树的节点个数为2 七一l ，分支的个数为2 七一2 当刀= 七+ l 时，即物种增加一个，则就是在原有的树上的每个节点上都可以再增加 一个分支，即 基于随机微分方程和结构e m 算法的系统发生树的构建第一章序言 ”( 七+ 1 ) = 端( 2 七一1 ) 将力= 七+ 1 代入( 1 1 1 ) 式的右边有： ( 2 ( 七+ 1 ) 一3 ) 1 2 ( k + 1 - 2 ) ( 七+ 卜2 ) ! ( 2 k 一1 ) ! ( 2 k 一3 ) ! ( 2 七一2 ) ( 2 k 1 ) = ：一：= 一 2 k - x ( 七一1 ) 12 x 2 一2 ( 七一2 ) ! ( 七一1 ) ( 1 1 4 ) 式与( 1 1 3 ) 式相等! 所以( 1 1 1 ) 得证 接着我们对( 1 1 2 ) 式证明： 当玎= 3 时，无根树的形式有图3 中的1 种：设物种为a ，b ，c c a b ( 1 1 3 ) 图3 ：三个物种无根树的形式 将刀= 3 代入( 1 1 2 ) 式的右边，得1 所以r = 3 时，( 1 1 2 ) 式成立 假设刀= 七时也成立，即 一= 羲焉 这时每棵树的节点个数为2 k 一2 ，分支的个数为2 k 一3 当胛= k + 1 时，即物种增加一个，则就是在原有的树上的每个分支上都可以再增加 一个分支，即 俨砸+ 1 ) = 筹焉泓- 3 ) ( 1 1 5 ) 将即= 七+ 1 代入( 1 1 2 ) 式的右边有： ! 兰! 垒! ! 二三21 ：! 圣生二兰! ! ：( 2 k - 5 ) ! ( 2 k - 4 ) ( 2 k - 3 ) ：! 兰垄二皇! ! ( 2 k 一3 、l 一= 一= 一= 一 一l 2 * + 1 - 3 ( 尼+ 1 3 ) 12 * - 2 ( 七一2 ) 12 2 k - 3 ( 七一3 ) ! ( 忌一2 )2 k - 3 ( 七一3 ) ! 。7 ( 1 1 6 ) v 一23 一一二七 堕气 一2 基于随机微分方程和结构e m 算法的系统发生树的构建第一章序言 ( 1 1 5 ) 式与( 1 1 6 ) 式相等，所以( 1 。1 2 ) 式得证 1 2 序列间距离的确定 树的分支长度我们称为距离，系统发生树中距离的确定一般有下列两种方法： 1 2 1 基于序列比对 在计算距离之前，首先进行序列比对，然后累加每个比对位置的得分可以运用 序列比较的方法，直接计算序列之问的距离如果在进行序列比较时，使用的是打分 函数或相似性度量函数，则需要将相似度( 或者是得分) 转换成距离 令s ( f ，) 是序列i 和序列各个比对位置得分的加权和，则归一化的距离 d ( i ，_ ，) 的计算公式为： 斗眷耥 其中母瓴- ，) 是序列i 和序列随机化之后的比对得分加权和，( f ，) 是序列 f 和序列歹所有可能的比对的最大值( 当两个序列相同时，取最大值) 两个序列归一 化距离的值处于0 和l 之间，当两个序列完全一致时，距离为0 ，当两个序列差异很大 时，距离接近1 1 2 2 基于统计的方法【1 】 首先我们先介绍几个概念： 字母表：字母表是指一个集合，它的每个元素称为字母例如：我们所研究的 d n a 序列，s = a ，t ，g ，c ) ，我们就称s 为4 个字母的字母表 字：字是由字母表中的字母构成的给定长度的字符串，若字母表有t 个元，给定 长度为z ，则长为，的字的种类数为f 。 窗口和窗口字：给定长度为z ，我们就可以得到长度为，个单位的窗口其中每 4 基于随机微分方程和结构e m 算法的系统发生树的构建第一章序言 个单位允许填入1 个字符，当所有单位都填满后，我们就得到1 个字，我们称为窗1 ：3 字 窗口滑动：用给定长度为，的窗口，在长度为刀的符号串t 上滑动，可以得到 以一，+ 1 个窗口字 特征矩阵：统计不同种类的窗口字可以得到1 个2 x t 7 的矩阵形的第一行 记录这个字符序列中在窗1 ：3 滑动过程中出现的不同种窗口字；第二行记录对应窗 口字的数目 设两个d n a 序列a ，b ，长度分别为n i 刀：，窗口长度为z ，则不同类型的窗口字数 总和为4 7 根据两个序列的特征矩阵，我们可以定义两个序列a ，b 之间的距离d ( a ，b ) ： 州固= 粉 n ( a ，曰) ：刀。+ 刀：一2 1 + 2 ，d ( a ，占) ：4 ii 疗j 一刀刘 耐表示序列a 的特征矩阵的第i 种窗口字的数目，砰表示序列b 的特征矩阵 的同种窗v i 字的数目于是： 州瑚= 石士 1 2 2丽l 2 鲰i = l 一砰 ，l l + 一+ 一。 1 3 基于距离构建系统发生树的方法 基于距离的系统发生树构建的方法的基本思路是：给定一种序列之间距离的 度量，在该距离度量下构建一棵系统发生树，使得该树能够最好地反映已知序列之 间的距离2 0 世纪8 0 年代以来，基于分子水平上距离的系统发生树的构建方法相 继提出，l iw h ( 1 9 8 1 ) 【2 】 s a t i t o un ，n e im ( 1 9 8 7 ) t 3 1 ，s o u r d i sj ，k r i m b a sc ( 1 9 8 7 ) t 4 1 ， s t u d i e r , j a ，k e p p l e rk ( 1 9 8 8 ) 口】，使得系统发生树的构建成为一个热门课题 下面我们介绍两种在不同假设条件下构建系统发生树的方法： 基于随机微分方程和结构e m 算法的系统发生树的构建第一章序言 1 3 1 非加权分组平均法 当用来构建系统发生树时，其假定的前提条件是：在进化过程中，核苷酸的替换 速率是均等的且恒定的，则在每一次分歧发生后，从共同祖节点到两个分类单元间 的分支长度是相等的 其基本思路如下： ( 1 ) 使所研究的每一个核苷酸序列( 节点) 自成一类开始如果有刀个核苷 酸序列( 节点) ，则有栉个类，每个类中包含一个元素； ( 2 ) 寻找距离最小的两个类将其合并，产生新的类如：若类c f ，c ，之问的 距离最小，则定义c 。= c ，c j ； ( 3 ) 定义c 辨与其他类c ，之间的距离，以及g ，c ，之间的距离 烈c ，c j ) 2 赢荟善烈墨力 崛绷= 塑焉铲 其中撑c f ，# q 表示类c ，c j 中序列的数目 经过这样一步运算，就从刀个类降为万一1 类； ( 4 ) 回到第一步，反复迭代，直到降为只剩下一个类，而这个类包含了所有 的序列 1 3 2 邻近加入法 邻近加入法( n e i g h b o rj o i n i n g ) 是另一种快速的聚类方法，该方法是s a i t o u 和 n e i 于1 9 8 7 年首次提出的，并且由s t u d i e r , j a ，k e p p l e rk 发展在构建系统发生树 时，该方法不需要对进化速率做出假设，并且从共同的祖节点到两个分裂单元间的 分支长度也不一样与t f - 力n 权分组平均法相比，邻近加入法在算法上相对较复杂， 它跟踪的是树上的节点而不是分类单元 6 基于随机微分方程和结构e m 算法的系统发生树的构建第一章序言 狂聚类过程中，根据原始距离矩阵，以及所有节点问的平均趋异程度，对每两个 节点间的距离进行调整，即将每个分类单元的趋异程度标准化，从而形成一个新的 距离矩阵重建时，将距离最小的两个叶节点连接起来，合并这两个叶节点所代表的 分类，形成一个新的分类在树中增加一个父节点，并在距离矩阵中加入新的分类， 同时删除原来的两个分类随后新增加的父节点被看成为叶节点，重复上一次循 环在每一次，循环过程中，都有两个叶节点被一个新的父节点所取代，两个类被合成 为一个新类整个循环直到只剩一个类为止从所得到的系统发生树来看，对于两 个聚在起的分类单元，其所在的叶节点到父节点的距离并不一定相同 该方法步骤如下： ( 1 ) 使所研究的每一个核苷酸序列( 节点) 自成一类开始如果有刀个核苷酸 序列( 节点) ，则有，个类，每个类中包含一个元素； ( 2 ) 对于所有分类单元，计算以，以2 i 至d ( x ，少) ，x ，j ，为分类单元中的 类： ( 3 ) 选择一对分类单元x 和y ，使d ( x ，y ) - 以一d ，最小； ( 4 ) 将x 和y 归并成新的类( 砂) ，在树中添加一个新的节点，将它与节点x 和 y 连接，新的节点代表新生成的分类，计算从x 和y 到新的节点( 砂) 的分 支长度： d ( z ，( 砂) ) = 三d ( x ，j ，) + 三( 以一办) m ，( 砂) ) = 兰m ，y ) + 圭( d y 一以) ； ( 5 ) 计算新类( 训与其他类u 的距离： d ( ( 】吵) ，“) = ( d ( x ，“) + d ( y ，”) 一d ( x ，y ) ) ； ( 6 ) 删除聚类x ，y ，添加新类( x y ) ，更新距离矩阵； ( 7 ) 如果有两个以上的分类单元存在，则继续执行循环，否则合并剩余的两 个类，并且连接这两个类 基于随机微分方程和结构e m 算法的系统发生树的构建第二章 随机微分方程和e m 算法 第二章随机微分方程和e m 算法 2 1 随机微分方程 定义l ：任一族随机变量 _ ：f t ) 称为以丁为参数值的随机过程； 定义2 ：设 形：f r + ) 是聊维随机过程，它满足以下条件： ( 1 ) w o = 0 ，i ：1 8 ： ( 2 ) 正态性：若0 s f o o ，则 形一彬 n ( 0 ，( f s ) ) ，r ”的 正定矩阵： ( 3 ) 增量独立性：若0 j t ，则 形一形) 是相互独立的， 则称彬为一个m 维b r o w n 运动，当= ，( 单位矩阵) 时，称形为m 维标准 b r o w n 运动 众所周知：常微分方程( o d e ) 一般形式是： d x ( t ) = f ( t ，x ( t ) ) d t ，t j ( 2 1 1 ) 其中厂：，x r dj r d 凡确定地依赖于初始状态的时间系统，通常可表示为某个形如( 2 1 1 ) 的常微分 方程模型 然而，在现实生活中，由于随机干扰的作用，时问系统通常表现出不确定性，致使 形如( 2 1 1 ) 的方程不再是正确描述这类系统的有效工具这需要对方程( 2 1 1 ) 进行 适当的改造，使之能反映随机干扰的作用 我们通常在( 2 1 1 ) 右端加入一个b r o w n 运动，使其转化为随机微分方程( s d e ) ： d x ( t ) = f ( t ，x ( t ) ) d t + g ( t ，z o ) ) d f y ( f ) ，t j ( 2 1 2 ) 其中f ：j xr d 啼r d 与g ：j r d 专r d ”是给定的b o r e l 可测函数，通常称 为漂移函数，g 称为扩散函数，形= ( w l ，w 2 ，1 6 。) r 是给定的m 维b r o w n 运动 基于随机微分方程和结构e m 算法的系统发牛树的构建第二章随机微分方程和e m 算法 下面自然我们要问：方程( 2 1 2 ) l 基j 解是如何定义? 以及在什么条件- f ( 2 1 2 ) 的解存 在且唯一? 定义3 ：若随机过程“：f t ) 满足以下的条件： ( 1 ) x ( f ) 是连续适应过程； ( 2 ) f ( t ，x o ) ) z 1 ( j ，r d ) ，g ( t ，x ( f ) ) z 1 ( j ，r d ”) ； ( 3 ) x o ) 满足如下的随机积分： x ( f ) = 而+ 厂( j ，x ( s ) ) 凼+ g ( j ，x ( s ) ) d ( j ) 则称x ( f ) 为方程( 2 1 2 ) 的具有初始值的解 引理2 1 - 设函数f ( t ，x ( f ) ) 与函数g ( f ，x ( f ) ) 在d x r d 上满足以下条件： ( 1 ) 一致l i p s c h i t z 条件：存在常数c ，使得： l i f ( t ，功一f ( t ，y ) l l + l l g ( t ，x ) - g ( t ，y ) 0 c l l x - y l l ； ( 2 ) 线性增长条件：存在正常数m ，使得： i 矿( f ，x ) 1 1 2 + l l g ( t ，x ) 1 1 2m o + i l x l l 2 ) 其中j j - | j 欧几里德范数 则对任给的初值，方程( 2 1 2 ) 存在唯一解x o ) 引理2 1 的证明见l i a n g j i a nh u 等( 2 0 0 7 ) 6 】 注1 - 若f ( t ，0 ) 与g ( t ，0 ) 对任意t r 均有界，则条件( 1 ) 可推出条件( 2 ) 因为：妙( f ，x ) 1 1 2 - l v ( t ，o ) 1 1 2 + i v ( ，x ) - f ( t ，o ) 1 1 2 r + c l l x 0 2 s ，( 1 + i l x l l 2 ) s l = m a x t ，c ) 同理：i i g ( t ，x ) 1 1 2 s ：( 1 + t l x l l 2 ) 贝t jl f ( t ，x ) 1 1 2 + i g ( t ，x ) 1 2 ( s 。+ s ：) ( 1 + 0 x i l 2 ) = m ( 1 + x l l 2 ) 令m = s 。+ s ： 9 基于随机微分方程和结构e m 算法的系统发牛树的构建第二章随机微分方程和e m 算法 2 2 1 酬算法的步骤 2 2酬算法 e m 算法最初是由d e m p s t e r a p , l a i r d n ，r u b i nd8 ( 1 9 7 7 ) 1 7 1 提出的e m 算法是 在缺失数据情况下求解参数模型中参数的最大似然估计或后验众数的一种迭代算 法e m 算法对处理一些带有不完全数据的问题变得相当容易，使得存在缺失数据 时我们能有效的进行参数估计e m 算法还可以用来处理一些困难的统计问题，像 混合模型和分层模型等，这些统计问题可能并不涉及缺失数据，但是如果我们恰当 地引入缺失数据，这些问题就很容易处理了 令y = ( ，) 是完全数据，其中为观察数据，圪插为缺失数据，下面我们引 进一些符号：( p 为参数) p ( oi ) 表示基于观察数据的参数口的后验分布密度，称为观察后验分布； p ( od 表示基于完全数据的参数口的后验分布密度，称为添加后验分布； p ( i ，口) 表示在给定口和观察数据r o 如下潜在数据匕括的条件分布密度 e m 算法就是利用p ( oi 功，p ( i ，d 两个密度函数，将缺失数据和参数秒 紧密联系起来 完全数据】，的密度函数可以写为： p ( ol 即= p ( 】么i 口) p ( 匕插ij 幺，p ) 有了完全数据】，的密度，我们就可以写出它的对数似然： l o g p ( 0l ，) = l o g p ( y o b s1 秒) + l o g p ( y m 缸il j = 廊，p ) e m 算法是一种迭代算法，我们目标是计算观察后验分布p ( oly ) 的参数日，它 的每一次迭代分两步完成：e 步和m 步假定经过s 步，参数目的估计值为臼( ”，则 s + 1 步迭代过程为： e 步：将l o g p ( 8 i ，匕。) 关于缺失数据匕捃的条件分布求期望 1 0 基于随机微分方程和结构e m 算法的系统发生树的构建第二章 随机微分方程和e m 算法 q ( o0 j , k ) = l o g p ( 口i ，) p ( k 渺，) 吧= 研l o g p ( od 渺，】 m 步：将q ( 0 1 秒”，r o b , ) 极大化，即求0 件1 这样就形成了一次迭代口s 哼秒什n 将上述两步迭代直到0 n 一秒，l 或 l i q o 1 ip ，r o b , ) 一q ( o j i 伊p n ，k ) 0 充分小为止 2 2 2 酬算法的收敛性 e m 算法的主要目的是提供了一个简单的迭代算法来计算后验众数或者似 然函数的最大似然估计，则就会提出这样的问题：由e m 算法迭代得到的序列 秒5 ) 是否收敛，如果收敛，能否是l o g p ( oi ) 的最大值或者是局部最大值，下面的两个 引理解决了这两个问题： 引理l ：e m 算法的每一次迭代后均可提高观察后验密度函数值，即： p ( o 1 i ) p ( o 5 i ) 引理2 ：如果q 徊l 口”，) 关于p 和0 0 ) 均是连续的，则e m 算法得到的估计序列 秒5 ) 收敛到l o g p ( oi ) 的稳定点 引理1 和2 的证明见、u ，c f j ( 1 9 8 3 ) t 引 引理2 的条件在大多数的情况下是满足的，所以e m 算法的结果只能保证收敛 到局部极大值点，并不能保证收敛到最大值点我们在运用e m 算法时，通过选择不 同的初始值进行迭代，可以解决这样的问题 基于随机微分方程和结构e m 算法的系统发生树的构建第三章基于随机微分方程的系统发生树的构建 第三章基于随机微分方程的系统发生树的构建 在构建系统发生树时，如果我们能知道两个序列或者两个物种之间的距离( 进 化时间) ，我们就可以通过非加权分组平均法，邻近加入法等方法构建系统发生，下 面我们利用随机微分方程来求解序列之间的距离 3 1 模型的建立 设所研究的，l 条d n a 序列分别是墨，s ：，s 。，我们目标是建立这刀条d n a 序 列的系统发生树，寻找他们的进化关系( 如果我们要研究甩个物种之间的进化关系， 可以选取疗个物种之间具有物种表征的同源d n a 序y o ) 下面我们主要是求出这行条d n a 序列之间的两两距离，序列s ，和s ，在t 时刻 之间的距离记为y q ( t ) ( i 歹，f = 1 , 2 ，拧，j = 1 , 2 ，刀) ，我们主要从两个方面考虑序列 之间的进化： ( 1 ) 进化主要受自身规律的影响，则在很短的时间址内进化距离我们记为 y q ( t + a t ) 一y q ( t ) 所以我们可以建立进化距离和时间之间的微分关系： y q ( t + a t ) 一l ( f ) = f ( y q ( t ) ；8 ) a t ( 口o 为参数) 即：d y o ( t ) = f ( y q ( t ) ；o ) d t ， 其中( 巧( f ) ；口) 是反映序列s ，和s ，之间进化的漂移函数； ( 2 ) 当然，进化还受到内，外部因素的影响，如自然环境的突变，人为的对物种进 化的干预和物种之间的相互干预等，我们近似地认为内，外部因素对进化过程的影 响是一种标准的b r o w n 运动扩散系数我们记为g ( l ( f ) ；秒) ， 贝u 有d 巧( f ) = g ( 巧( f ) ；伊) d ( f ) 1 2 基于随机微分方程和结构e m 算法的系统发生树的构建 第三章基于随机微分方程的系统发生树的构建 n x ( n - i ， - 、 g ：r 0 专r _ 广， 形( f ) = ( ( f ) ，( f ) 。既x ( 川) ( f ) ) 丁是竺拿坐维标准 2 布朗运动 综合上述两个因素，我们可以建立如下随机微分方程： d y ，( t ) = f ( y o ( t ) ；o ) d t + g ( ( t ) ；o ) d w ( t ) 推广到我们所研究n 物种时，我们可以得到如下形式的随机微分方程组： d y l 2 ( f ) d r , 3 ( f ) d r , 。( f ) d 艺3 ( f ) 哝纠) 。( f ) 我们简记为 g ( y 1 2 ( f ) ；口) g ( 1 3 ( f ) ；d g ( 鬈。( ，) ；力 g ( 艺3 0 ) ；d g 化删。( f ) ；口) d y ( t ) = f ( y ( t ) ；o ) d t + g ( 】，( f ) ；秒) d 形( f ) 其中： d y ( t ) = d r l 2 0 ) d r , 3 ( f ) 媚。( f ) 妈( f ) d ( 一) 。( f ) d w ( t ) n x ( n - 1 ) n x ( n 1 ) n x ( n 1 ) 月( p 1 ) 一n x ( n - i ) f ：r 2 o 专r 2 g ：r 2 xo 专r 22 为了以后表达的方便，我t f - i r r ( 力为z ，g = 旦丛笋尘 d r , = f ( y , ；o ) d t + g ( y , ；8 ) d w ( t ) ( 3 1 1 ) 如果( 3 i 1 ) 式满足( b d me r a k e r ( 2 0 0 1 ) 9 1 ) ： ( 1 ) 漂移函数和扩散函数应满足l i p s c h i t z 条件和线性增长条件 即：| i f ( x ；o ) 一( y ；秒) 0 + 0 9 ( z ；秒) 一g ( y ；口) l i c l x - y l i +出 一、，j 秒 、，、，、j、， e 泐泐 一泐泖 似 )、， ) ) 月 枷 2 3 疗 3 加 瓯 基于随机微分方程和结构e m 算法的系统发生树的构建第三章基于随机微分方程的系统发生树的构建 l b e ( x ；o ) 1 1 2 + l l g ( = ；o ) 1 1 2 c 2 ( i + l l x l l 2 ) ( c 为正常数) ( 2 ) g ( x ；口) g ( z ；口) 是正定矩阵，v x r 9 贝l j ( 3 1 1 ) 式有解以后我们总是假定满足这两个条件 3 2 对随机微分方程的求解 利用随机微分方程来解决实际问题时，会涉及到缺失数据，这就对解随机微分 方程带来困难所以我们利用g i b b s 抽样来解该随机微分方程，可以避免这样的 问题 首先我们对( 3 1 1 ) 式进行离散，有： r = f ( y , ；0 ) k t 4 - g ( r ；口) 形( f )( 3 2 1 ) i jd 其中a w n ( o ，址)，! q 是z q 维的单位矩阵 注2 ：当f 足够小的时候，离散的近似解可以近似认为是连续解 下面我们引进缺失数据的问题： 我们将时间间隔【0 ，r 】分割成s = m t 个等距点，0 = f 0 t l l o g p ( o ，m ( 。) ) ，因此，在迭代过程中我们选择一个结构m ( 州) 使得 期望最大，利用该定理，我们可以改造m 步，即选择结构m ( 槲) ，使得 q ( m ( 州) im 。) ) q ( m ( 。) im 。) ) 利用定理4 1 1 在加上一些条件，我们就可以是结构 m ( 。) ) 的收敛性 定理4 1 2 ：假设m ( o ) ，m ( 1 ) ，m ( 2 ) 是结构e m 算法产生的一组结构序列如果 硝的结构总数有限，则 m m ) 收敛 4 2 似然函数的确定 在上一节中，我们知道在e 步中，涉及到完全数据的似然，这就需要首先对给定 一棵系统发生树的拓扑结构( d n a 序列的进化关系) ，如何定义关于缺失数据的似 然函数，即如何建立关于树的概率模型而在基于d n a 序列的系统发生分析方面， 概率模型( f e l s e n s t e i n ，j ( 1 9 8 1 ) 旧) 首先依赖于一个合理的可靠的多重序列比对，然 基于随机微分方程和结构e m 算法的系统发生树的构建第p l l 章基于结构e m 算法的系统发生树的构建 后检测每一列的变化 4 2 1 结构已知时的系统发生树的似然函数的确定 我们以图4 为例，建立模型： 4 5 图4 ：以节点0 为根的系统发生树 图4 中的树有五个外部节点，3 个内部节点，根为0 ，节点和根均为d n a 序列，所以 一共九个序列，我们记为 j 。，j ：，j ，) ，序列的长度为k ，分支长度为f 。f 。 ( 1 ) 假设两个序列不同的位点之间的进化是独立的，则在给定序列岛，s ，在第七个位 点，在时间t 内发生进化的概率为p 与( k ) s s ( 七) ( f ) ，所以序列s ，s ，在时间f 内发生进化 的概率为： ( 2 ) 过过该树的所有内部节点和外部节点的序列是已知的，则似然函数是： l = 7 1 s op ( t o ) p 丑( f 1 ) p 屯( f 2 ) p 蚋( f 8 ) p 曲曲( f 3 ) p 曲。，( f 7 ) p 即( t 4 ) p 。，( t 5 ) 1 9 、，l 似p似 焉 p k 兀树 = 、，- 0 而 p 基于随机微分方程和结构e m 算法的系统发生树的构建第四章基于结构e m 算法的系统发生树的构建 【4 2 1 ) 为根的先验分布 ( 3 ) 若所有的外部节点是已知的，而内部节点是未知的，则似然函数为： 三= l r , o p , o , 。( t 。) p j 6 ( ，。) p 铀( ，2 ) p 饷( f 8 ) p 即，( f 3 ) p 球，( f ，) 见，( f 。) 如屯( ，5 ) s o 5 1s l ( 4 2 2 ) ( 4 ) 下面我们通过迭代的思想简化，假设磋是第七个节点序列的似然函数，如果 两个子代序列分别为s ，s ，则 雀) _ ( p h ( f ，) 崔) ( p ( f ) 崔) s i s i 所以对于该树，我们可以先计算雀和磋，然后计算磋和线= 三 所以似然函数为： 三= ( 万j o l ( o ) 三。毛 s o = ( p j o ( 气) p 钠o 。) p 讷( f ：) ) ( p ( t s ) p 码( f ，) p 屯却( f ，) 嚣 s oj 6吨 = 万 q ( m ( 。) lm ( 。) ) 将上述两步进行直到 q ( m 川) im 。) ) = q ( m 。) im ( 。) ) 为止 我们得到的关于这刀条d n a 序列的系统发生树的结构 m ( 。) ) 是否收敛呢? 答 案是肯定的，因为根据我们在序言中证明的结论，n 条d n a 序列，可以建立的有根树 基于随机微分方程和结构e m 算法的系统发生树的构建第四章基于结构e m 算法的系统发生树的构建 的数目为脚，( 加东禹无根树的数目为蚺脚，= 栽禹所以不管 是有根树还是无根树，树的数目总数是有限的，即朋是有限集，根据定理4 1 2 ，我们 知 m 。) ) 是收敛的，- 记l i m m ( 。) - - m ，则m 就是我们所要求的这刀条d n a 序列的最 优的系统发生树 基于随机微分方程和结构e m 算法的系统发生树的构建 第五章结论与研究 第五章结论与研究 随着生物技术的进步和人类基因组计划的实施，人类拥有的生物数据越来越 多，使得从分子水平构建系统发生树成为可能，但数据的海量增长使得科学家面临 很大的挑战本文运用随机微分方程的方法刻画生物序列之间的进化距离，并运用 了g i b b s 抽样求解该方程，并且针对如果只给定所研究的d n a 序列之间部分的序 列之间的距离，运用结构e m 算法去寻求最优的系统发生树。 但由于生物实际问题本身的复杂性，以及统计方法在处理实际问题的局限性 在系统发生树构建方面关于还有许多需要进一步研究的问题： ( 1 ) 我们在基于随机微分方程的系统发生树的构建，模型建立一节中，假定外部 因素对进化过程的影响是一种标准的布朗运动，但在实际中的外部因素对 进化的影响可能是及其复杂的，我们如何更好的量化这些外部因素的影响 ( 2 ) 在基于结构e m 算法构建系统发生树，似然函数建立的这一节中，我们假定 任意的两个序列不同的位点之间的进化是独立的，这在实际问题中不一定 成立，如果把这个假设去掉，我们如何处理 ( 3 ) 本文提出运用随机微分方程和结构e m 算法求解d n a 序列的系统发生树 的理论，下面是如何运用实际数据来检验这两个理论的有效性 ( 4 ) 我们已经提出两种方法构建系统发生树，在加上先前的方法，我们如何来比 较这些方法在研究生物序列进化方面的准确性和效率 ( 5 ) 随着生物技术的进步，我们可以得到更多的基因表达序列，我们如何挑选这 些信息去构建系统发生树 基于随机微分方程和结构e m 算法的系统发生树的构建参考文献 参考文献 【l 】 汪浩统计方法求系统树云南大学学报( 自