第４３卷第６期 ２０１３年１１月 东南大学学报（自然科学版） ＪＯＵＲＮＡＬ ＯＦ ＳＯＵＴＨＥＡＳＴ Ｉ ⅡＶＥＲＳＩＴＹ（Ｎａｔｕｒａｌ Ｓｃｉｅｎｃｅ Ｅｄｉｔｉｏｎ） Ｖｏ１．４３ ＮＯ．６ Ｎｏｖ．２０１３ ｄｏｉ：１０．３９６９／ｊ．ｉｓｓｎ．１００１—０５０５．２０１３．０６．００３ 基于基因扰动及变分逼近技术的基因网络推断 董自健 宋铁成 袁 创。 （ 东南大学信息科学与工程学院，南京２１００９６） （。淮海工学院电子工程学院，连云港２２２００５） （ 理工大学医疗科技及信息学系，） 摘要：为了有效提高基因网络推断的精度，基于基因表达数据和基因扰动数据，将基因 网络建模为结构方程模型，并进一步转化为验证性因子分析（ＣＦＡ）模型，然后使用贝叶斯方法 求解ＣＦＡ模型参数．在贝叶斯分析中，为减少计算量，不采用常用的马尔科夫一蒙特卡洛抽样方 法，而是采用变分逼近技术对参数的联合后验分布进行因式化，并获得参数的后验包含概率分布 及参数的后验分布．同时使用重要性抽样技术对ＣＦＡ模型的推断参数进行加权平均．仿真结果 表明，ＣＦＡ模型和变分逼近技术是有效和可靠的．根据实验数据，使用所提算法推导了具有３５ 个基因的酵母基因网络． 关键词：基因网络；验证性因子分析模型；变分逼近；重要性抽样 中图分类号：Ｑ７８ 文献标志码：Ａ 文章编号：１００１—０５０５（２０１３）０６—１１４７４３５ Ｖａｒｉａｔｉｏｎａｌ ａｐｐｒｏｘｉｍａｔｉｏｎ ｉｎｆｅｒｅｎｃｅ ｆｏｒ ｇｅｎｅ ｒｅｇｕｌａｔｏｒｙ ｎｅｔｗｏｒｋｓ ｆｒｏｍ ｇｅｎｅ ｐｅｒｔｕｒｂａｔｉｏｎｓ Ｄｏｎｇ Ｚ０ｉａｎ ‘ Ｓｏｎｇ Ｔｉｅｃｈｅｎｇ Ｙｕａｎ Ｃｈｕａｎｇ （ Ｓｃｈｏｏｌ ｏｆ Ｉｎｆｏｒｍａｔｉｏｎ Ｓｃｉｅｎｃｅ ａｎｄ Ｅｎｇｉｎｅｅｒｉｎｇ，Ｓｏｕｔｈｅａｓｔ Ｕｎｉｖｅｒｓｉｔｙ，Ｎａｎｊｉｎｇ ２１００９６，Ｃｈｉｎａ） （ Ｓｃｈｏｏｌ ｏｆＥｌｅｃｔｒｏｎｉｃ Ｅｎｇｉｎｅｅｒｉｎｇ，Ｈｕａｉｈａｉ Ｉｎｓｔｉｔｕｔｅ ｏｆＴｅｃｈｎｏｌｏｇｙ，Ｌｉｎｙｕｎｇａａｎｇ ２２２００５，Ｃｈｉｎａ） （ Ｄｅｐａｒｔｍｅｎｔ ｏｆ Ｈｅａｌｈ Ｔｅｃｈｎｏｌｔｏｇｙ ａｎｄ Ｉｎｆｏｒｍａｔｉｃｓ，ｔｈｅ Ｈｏｎｇ Ｋｏｎｇ Ｐｏｌｙｔｅｃｈｎｉｃ Ｕｎｉｖｅｒｓｉｔｙ，Ｈｏｎｇ Ｋｏｎｇ，Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ：Ｔｏ ｉｍｐｒｏｖｅ ｔｈｅ ｉｎｆｅｒｅｎｃｅ ａｃｃｕｒａｃｙ ｏｆ ｇｅｎｅ ｒｅｇｕｌａｔｏｒｙ ｎｅｔｗｏｒｋｓ（ＧＲＮ），ｕｓｉｎｇ ｂｏｔｈ ｇｅｎｅ ｐｅｒｔｕｒｂａｔｉｏｎｓ ａｎｄ ｇｅｎｅ ｅｘｐｒｅｓｓｉｏｎ ｄａｔａ，ＧＲＮ ｉｓ ｍｏｄｅｌｅｄ ａｓ ａ ｓｔｒｕｃｔｕｒａｌ ｅｑｕａｔｉｏｎ ｍｏｄｅｌ（ＳＥＭ），ａｎｄ ｆｕｒｔｈｅｒ ｔｒａｎＳｆｏｒｍｅｄ ｉｎｔｏ ｔｈｅ ｃｏｎｆｉｍａｔｒｏｒｙ ｆａｃｔｏｒ ａｎａｌｙｓｉｓ（ＣＦＡ）ｍｏｄｅ１．Ｔｈｅ Ｂａｙｅｓｉａｎ ａｐｐｒｏａｃｈ ｉｓ ｕｓｅｄ ｔｏ ｉｎｆｅｒ ｔｈｅ ｐａｒａｍｅｔｅｒｓ ｏｆ ｔｈｅ ｒｅｇｕｌａｔｏｒｙ ｎｅｔｗｏｒｋｓ．Ｉｎｓｔｅａｄ ｏｆ ｔｈｅ Ｍａｒｋｏｖ ｃｈａｉｎ Ｍｏｎｔｅ Ｃａｒｌｏ （ＭＣＭｃ）ｍｅｔｈｏｄ，ｔｈｅ ｖａｒｉａｔｉｏｎａｌ ａｐｐｒｏｘｉｍａｔｉｏｎ ｍｅｔｈｏｄ（ＶＡＭ）ｉｓ ａｐｐｌｉｅｄ ｆｏｒ ｉｔｓ ｌｏｗｅｒ ｃｏｍｐｕｔａ— ｔｉｏｎ ｃｏｓｔ，ｗｈｉｃｈ ｆａｃｔｏｒｉｚｅｓ ｔｈｅ ｉｏｉｎｔ ｐｏｓｔｅｒｉｏｒ ｄｉｓｔｒｉｂｕｔｉｏｎ ｏｆ ｐａｒａｍｅｔｅｒｓ，ａｎｄ ｏｂｔａｉａｓ ｔｈｅ ｐｏｓｔｅｒｉｏｒ ｉｎ— ｃｌｕｓｉｏｎ ｐｒｏｂａｂｉｌｉｔｙ ｄｉｓｔｒｉｂｕｔｉｏｎ ａｎｄ ｔｈｅ ｐｏｓｔｅｒｉｏｒ ｄｉｓｔｒｉｂｕｔｉｏｎ ｏｆ ｐａｒａｍｅｔｅｒｓ．Ａｎ ｉｍｐｏｒｔａｎｃｅ ｓａｍｐｌｉｎｇ ｔｅｃｈｎｉｑｕｅ ｉｓ ｔｈｅｎ ａｐｐｌｉｅｄ ｔｏ ｏｂｔａｉｎ ｔｈｅ ｗｅｉｇｈｔｅｄ ａｖｅｒａｇｅ ｏｆ ｔｈｅ ＣＦＡ ｉｎｆｅｒｒｅｄ ｐａｒａｍｅｔｅｒｓ．Ｓｉｍｕｌａｔｉｏｎｓ ａｒｅ ｃａｒｒｉｅｄ ｏｕｔ ｔｏ ｖｅｒｉｆｙ ｔｈｅ ｅｆｆｅｃｔｉｖｅｎｅｓｓ ａｎｄ ｒｅｌｉａｂｉｌｉｔｙ ｏｆ ｔｈｅ ＣＦＡ ｍｏｄｅｌ ａｎｄ ｔｈｅ ｖａｒｉａｔｉｏｎａｌ ａｐｐｒｏｘ— ｉｍａｔｉｏｎ．Ｂａｓｅｄ ｏｎ ｔｈｅ ｅｘｐｅｒｉｍｅｎｔａｌ ｄａｔａ，ｔｈｅ ｒｅｇｕｌａｔｏｒｙ ｉｎｔｅｒａｃｔｉｏｎｓ ａｍｏｎｇ ａ ３５ ｙｅａｓｔ ｇｅｎｅｓ ｎｅｔｗｏｒｋ ｒｅ ｉｄｅｎｔａｉｆｉｅｄ ｗｉｔｈ ｔｈｅ ｐｒｏｐｏｓｅｄ ＶＡＭ ａｌｇｏｒｉｔｈｍ． Ｋｅｙ ｗｏｒｄｓ：ｇｅｎｅ ｒｅｇｕｌａｔｏｒｙ ｎｅｔｗｏｒｋ（ＧＲＮ）；ｃｏｎｆｉｍａｔｒｏｒｙ ｆａｃｔｏｒ ａｎａｌｙｓｉｓ（ＣＦＡ）ｍｏｄｅｌ；ｖａｒｉａ— ｔｉｏｎａｌ ａｐｐｒｏｘｉｍａｔｉｏｎ；ｉｍｐｏｒｔａｎｃｅ ｓａｍｐｌｉｎｇ 基因网络（ＧＲＮ）结构的推导对理解基因 功能和复杂疾病的致病原因至关重要．基因网 量的实验．由于计算技术和计算能力的快速提高， 现在通常不是单纯使用实验手段，而是根据有限的 实验数据，使用统计方法推断基因网络． 络的推断需要识别出基因之间的相互作用，而仅使 用实验方法确定基因之间的相互作用需要进行大 由于引入了扰动信息，基因扰动技术可有效提 收稿日期：２０１３－０４－２５．作者简介：董自健（１９７３一），男，博士，副教授；宋铁成（联系人），男，博士，教授，博士生导师，ｓｏｎｇｔｃ＠ｓｅｕ．ｅｄｕ．ｃｎ． 基金项目：国家自然科学基金资助项目（６１２７１２０７）、江苏省自然科学基金资助项目（ＢＫ２０１１３９８）、江苏省高校优秀中青年骨干教师和校长 境外研修计划资助项目． 引文格式：董自健，宋铁成，袁创．基于基因扰动及变分逼近技术的基因网络推断［Ｊ］．东南大学学报：自然科学版，２０１３，４３（６）：１１４７— １１５１．［ｄｏｉ：１０．３９６９／ｊ．ｉｓｓｎ．１００１—０５０５．２０１３．０６．００３］ １１４８ 东南大学学报（自然科学版） 第４３卷 高ＧＲＮ推断的精度．基于基因表达数据，可以把 基因网络建模为图结构…，其中基因用节点 表示，节点之间的边表示基因之间的关系．另 外高斯图模型 、线性回归模型等 都可以用于 基因网络的推断．但这些方法的准确度都不 （２）作如下分解： ｙｌ ● ／Ｉ１ ● ● ： ： ： ｙ＝ ● ／ｌ ● 力＋ ● （３） ： ： ： 高．结构方程（ＳＥＭ）模型也可用于推断基因 网络 ，如Ｌｏｇｓｄｏｎ等 使用ＳＥＭ对基因网络 进行建模，并使用Ａｄａｐｔｉｖｅ Ｌａｓｓｏ（ＡＬ．ｂａｓｅｄ）算法 推导网络参数．通过仿真，作者认为ＡＬ．ｂａｓｅｄ算法 ｙ ｊ ＬＡＮ．ｊ 解／ｌ中参数，因此问题分解为 Ｌ 由于已知数据均在 和ｙ中，所以可逐行求 优于已有的其他各种算法． 基于基因表达和ｅＱＴＬｓ数据，本文把基因 网络建模为ＳＥＭ结构．由于这２类数据均已知，故 不存在潜变量，因此本文把ＳＥＭ模型转化为验证性 因子分析（ＣＦＡ）结构，并进而简化为线性回归模型． 为降低算法对数据的敏感性，本文使用贝叶斯推断 方法代替常用的最大似然以及各种Ｌａｓｓｏ方法求解 线性回归模型；为了降低计算量，不使用通常的马尔 科夫一蒙特卡洛（ＭＣＭＣ）方法，而是采用变分逼近 方法（ＹＡＭ）进行贝叶斯推断．仿真结果表明ＣＦＡ 模型和ＶＡＭ算法都是有效和可靠的．ＶＡＭ算法的 性能较好，在误检率和检测率上都远高于ＡＬ—ｂａｓｅｄ 算法；与最大似然等方法相比，变分逼近方法无需准 确估计初始参数，对数据适应性较好． １基因网络结构的ＣＦＡ模型 假设有 个基因，Ⅳ０个标记，Ｊ７ｖ个测量样 本．基因网络结构服从如下结构方程模型： Ｙ＝ＡＸ＋ＢＹ＋￡ （１） 式中，ｙ为Ｎ ｘＮ矩阵，其中ｙ 『为基因ｋ在第 次 测量中所获得的表达数据； 为Ｎｇ ＸＮｇ矩阵，包含 所有的网络参数； 为Ｎｑ ＸＮ ｅＱＴＬｓ数据矩阵，Ⅳ口 ≥ ；Ａ为 Ｘ Ｎｑ矩阵，表示每个ｅＱＴＬｓ对每个 基因的影响； 为Ｎｇ ｘＮ矩阵，其中 表示第 次 测量中基因ｋ的表达误差． 假设基因自身不存在相互作用，因此矩阵 的对角元素值为０．另外，假设所有的ｅＱＴＬｓ位点 已根据已有的方法确定出来，但是每个ｅＱＴＬｓ的 影响未知，因此Ａ中有Ⅳ０个位置已知的未知参 数，剩下的Ⅳ 一Ⅳｑ个元素都是０．不失一般性， 假设Ⅳ ＝Ｎｑ． 式（１）中所有未知参数都在 和Ａ中，因此模 型可变为 Ｙ＝ＡＯ＋ （２） 式中，Ａ＝［Ｂ Ａ］， ＝［ｙ ｒ．现在模型已变为 一类特殊形式的ＳＥＭ模型，即ＣＦＡ模型．将式 ｈｔｔｐ：／／ｊｏｕｒｎａ１．ｓｅｕ．ｅｄｕ．ｃｎ Ｙｋ＝Ａ ＋Ｙｋ ｋ＝１，２，…，Ｎｇ （４） ２ ＣＦＡ模型的变量逼近推断 假设 矩阵稀疏，即 中大部分数据为０，因 此可假设Ａ中所有参数服从均值为０的高斯分 布，且 服从Ｎ（Ｏ， Ｉ）． 将式（４）转化为 Ｄ ＝∑ ＋ ｋ＝１，２，…，Ｎｇ（５） 式中， ，为 的第 行参数，Ｄ＝２Ｎｇ．现在问题转 化为一个标准的线性回归问题．如前所述，在贝叶 斯推断中，如果采用ＭＣＭＣ方法，在变量较多的情 况下，计算量将非常大．因此本文使用变分逼近替代 ＭＣＭＣ方法用于网络参数的推断．变分逼近技 术已得到很多应用，如Ｌｏｇｓｄｏｎ等 和Ｃａｒｂｏｎｅｔｔｏ 等 把该方法用于多位点的基因组关联分析． 引入一个二进制向量Ｚ ＝｛Ｚ Ｚ也，…，Ｚ 。｝，用 于表示／ｌ 中每个元素的取值情况．如果Ｚ ，＝１，则 以目以概率１不为零；如果Ｚ村＝０，则Ａ 概率１为 零．定义参数组０＝（ｚ ，／ｌ ），假设该参数组的先验 概率Ｐ（Ｚ ，＝１）＝丌， ｆ～Ｎ（０， ＾２）． 求解的目标就是根据已知数据计算以 的后验 概率分布，以及每一个以 是否包含在模型中的概 率，即Ｐ（ｚ ＝ｌ Ｉ ｙ ，ｎ）．结合Ａ 后验分布以及 Ｐ（Ｚ ，＝１ Ｉ ， ），即可确定基因．『是否对基因ｋ有 作用，以及系数的大小． ２．１参数后验概率计算 Ｂａｙｅｓｉａｎ推断需要计算后验概率Ｐ（０ ｌ ，ｙ ， ‘ ），其中 是｛７ｒ， ， ２＾｝的一组先验分布参 数值．该计算通常需要多重积分，计算量非常大． 根据变分理论＿８ Ｊ，Ｐ（０ Ｉ ， ，ｔ， ’）一 ｑ（Ａ ，ｚ ），且ｑ（Ａ ，ｚ ）可分解为 Ｄ Ｄ ｑ（ａ ，Ｚｋ）＝Ｉ－Ｉｑ（Ａ灯，Ｚｋｊ）＝Ｉ－１ｑ（Ａ衄Ｉ Ｚｋｊ）ｇ（ｚ ） ｋ＝ｌ ＝１ （６） 其中，每个分项可表示为如下形式 Ｊ： 第６期 董自健，等：基于基因扰动及变分逼近技术的基因网络推断 １１４９ ／Ｌｗ ／ｚｋｊ，ｓｋ２ ｉ Ｉ ｚｋｊ）＝ａｋＮ（Ａｋｊ［村…组合．为简化起见，假设重要性抽样 （ｔ，‘ ）为均 ）ｍ 、 匀分布，所以ｐ（ｏ‘ ）＝１／６００，ｍ＝１：１：６００． ｑ（Ｚ村＝１）＝ 村 （７） 式（１３）中的Ｐ（Ｙｋ ｌ ｌｆ２，ｔ，）很难计算，Ｃａｒｂｏｎｅｔ— ｔｏ等 使用Ｊｅｓｓｅｎ不等式获得了其下界： Ｉ Ｉ根据变分逼近理论可得出迭代参数为 ，ｍ ２ 丽、（８） ｐ（Ｙ ｌｎ，ｕ）：一 Ｎｌ。ｇ（２耵 ）一 灯＝ （（ ）Ｊ一∑（ ） 蛐）（９ ％ 等 石 丽 式（８）一（１０）构成一个完整的迭代．在算法实现 中，可采用一种简单的迭代结束判断方式，即如果 连续２次迭代中产生的 差值很小，如∑（ 一 ｈ－ ） ＜ ，即可停止迭代．其中 是预设的一个 值，需要根据待求参数的精度要求及算法复杂度要 求而设定． ２．２基于重要性抽样技术的包含概率计算 为考察基因ｋ和基因．『之间的关系，还需判断 每个参数（ ）是否包含在模型中，即计算 Ｐ（Ｚ＾，＝１ ｌ ，Ｙｋ）＝ Ｊｐ（Ｚ ｆ＝１ ｌ ，Ｙｋ，ｖ）ｐ（ｔ，ｌ力， ）ｄｖ（１１） 该积分计算非常复杂，需要近似计算．已知 Ｐ（Ｚ ＝１ ｌ２／，ｙｔ，０）＝Ｏｔ材，为了计算Ａｔｆ是否包含在 模型中，需要根据式（１１）计算其后验概率．这里使 用重要性抽样，即 Ｍ ｐ（ｚ ＝１ Ｉ ， ） ｍ＝１ ∑ ｗ（ｖ ’）（１２） 式中，ｗ（ｖ ）为重要性抽样系数， （Ｙ Ｉ ，ｔ，‘ ）ｐ（ｖ‘ ） ｐ（ｖ‘ ’、 Ｐ（ Ｉ ， ‘ ）ｐ（ｖ‘ ） （１３） （ ‘ ’） 其中，１．１ ’为参数组（７ｒ， ， ２）的一组抽样点 （ ， ’， ’）；ｐ（ ’）为 ， ｔ，ｚ 的先验分 布参数组（仃， ， ２）在样本为 ’时所获得的抽 样值．这里假设Ｐ（Ｚ村＝１）：７ｒ，７ｒ～ｂｅｔａ（Ｓ ，ａ，ｂ）， 其中ａ，ｂ为超参数，需根据数据情况设定； 一Ⅳ （０，ｒｏ Ｉ）， 一ｇａｍｍａ（ ，Ｇ。，Ｇ ），其中Ｇ ，Ｇ６为 超参数，需预先设定；Ａ ｆ—Ｎ（Ｏ， ２ｄ）， 取自分布 ａ／（１＋ａｓ ） ，ａ需要预设．在参数推导过程中，需 要对每个超参数输入一个参数抽样序列，比如对 ｂｅｔａ分布，可取 ＝０．Ｏ１：０．Ｏ１：０．１；对ｇａｍｍａ分 布，取Ｓ ＝０．１：０．１：１；而Ｓ＾可取０．０５：０．０５：０．３． 因此，｛ｔ，‘¨，ｔ，（ ，…， （ ’｝共有１０×１０×６＝６００种 ２ 一 ２ ：耄Ｘ．－ａｌ（、 “ Ｄ叫。ｇ（等）一 ｃ ）ｌｏｇ（ ）＋ Ｄ警［１＋１。ｇ（轰）一 １（１４） 式中， ．ｘ／．ｔ 表示 和 两个行向量的点乘． 在计算时可使用这个下界代替Ｐ（Ｙｋ『力，１３）． 算法的整体思路是对每个参数组７．３ ’获得相应 的（７ｒ，ｏｒ ， ２），然后代人式（８）一（１０）进行迭代．迭 代收敛后，再根据式（１３）、（１４）计算相应的重要性 加权系数．在所有的ｕ ’计算结束后，使用式（１２）计 算后验包含概率．算法最后输出Ａ 估计，即 Ａ灯 ∑ ｗ（ ） （１５） 算法具体过程如下： 输入：ｘ，ｙ，抽样样本序列｛ ， ，ＳＡ｝，以及 ‘ ， ‘ ， ． １根据式（２）计算ｎ； ２从（ ， ， ）的联合分布中获得重要性抽样｛”（”，…， Ｖ（Ｍ ｝及其均匀分布值｛ （”（ ），…，ｐ（ｖ（ ）｝； ３ ｆｏｒ ｋ：１：Ⅳｇ（对参数矩阵逐行求解） 置ｔ，＝＂（ ；令ｆｌａｇ＝１；根据式（８）计算 ； ｗｈｉｌｅｆｌａｇ ｆｏｒ ＝１：Ｄ 根据式（９）、（１Ｏ）分别更新 和 ｅｎｄ Ｊ； 如果连续２次迭代中所产生的 差值∑（ 一 ） ＜ ，令ｆｌａｇ＝０，停止迭代． 根据（ ， ， ）的先验分布，令”＝ （ ，获得Ｐ （ｔ，（ｍ ） 根据式（１３）、（１４）获得重要性系数ｗ（”（ ）； 根据式（１２）获得后验概率ｐ（ｚ ＝１ ｌ ｎ， ）； 根据式（１５）计算网络参数Ａ ； ３ 仿真结果 令ⅣＥ表示 中边的数量，Ⅳｍ表示经算法推 导所获得的曰 中边的数量，Ⅳ；ａ。 表示存在于曰 但 不存在于 中边的数量，Ⅳ 。表示同时存在于 和 中边的数量，因此Ⅳ ＋Ⅳ ＝ ．检测率 （ｐｏｗｅｒ ｏｆ ｄｅｔｅｃｔｉｏｎ，ＰＤ）和误检率（ｆａｌｓｅ ｄｅｔｅｃｔｉｏｎ ｈｔｔｐ：／／ｊｏｕｍａ１．ｓｅｕ．ｅｄｕ．ｃａ １ｌ５０ ｒａｔｅ，ＦＤＲ）定义如下： ＰＤ： ’Ｅ 东南大学学报（自然科学版） 第４３卷 在检测率和误检率２个方面，ＶＡＭ算法均优于 ＡＬ－ｂａｓｅｄ算法． ’ＩＥ ，ＦＤＲ： 本文分别仿真了含有３Ｏ个基因的有向无环和 有向有环网络．每个基因节点平均有３条边，即ⅣＦ ＝４酵母基因网络的推断 本节使用ＶＡＭ算法对一酵母基因网络进行 推断．数据来自对酵母菌株ＢＹ４７１６和ＲＭ１１—１Ａ 进行交叉实验所获得的１ １２组基因表达数据，以及 基因标记数据 ．数据包含５ ７２７个基因．由 ３．如果从节点 到节点ｋ有边，则从（～１， ０．５）ｕ（０．５，１）中均匀抽样；否则 ，置为０． 一ｅＱＴＬ中每个 以概率（０．２５，０．５，０．２５）从集合 ｛１，２，３｝中抽取．测量噪声服从Ⅳ（０，０．１）分布． 不失一般性，设Ａ为单位阵．ｙ通过式（１）计算 获得．在仿真中，如果算法返回的Ｐ（Ｚ ，＝１ Ｉ力， Ｙ ）＞０．９５，且Ａ ，＞０．１０，则认为从节点．，到节点 ｋ之间有一条边，其他情况均判为２个节点之间 无关系． 本文主要与ＡＬ．ｂａｓｅｄ算法在检测率和误检率 上进行性能比较，图１给出了仿真结果．可以看出， １ ０ 金０ ０ ０ 抽样数量 （ａ）有向无环图的检测率 Ｏ．６ 重 ； Ｏ．３ １ 抽样数量 （ｂ）有向无环图的误检率 抽样数量 （Ｃ）有向有环图的检测率 ： 。 。 抽样数量 （ｄ）有向有环图的误检率 图１ ＶＡＭ算法在３０个基因网络中的性能 ｈｔｔｐ：／／ｊｏｕｒｎａ１．ｓｅｕ．ｅｄｕ．ｃｎ 于测量数据样本有限，且存在一定的误差，不可能 用这些数据较为准确地推断出全部５ ７２７个基因 之间的关系，因此需要使用预处理技术，选择 具有较强表达数据差异的基因进行分析＿５ｊ．本文 采用Ｌｏｇｓｄｏｎ等 提供的经过预处理的基因表达 数据以及相应的ｅＱＴＬｓ数据进行网络参数推导． 这些数据共包含３５个基因的相关测量数据，具体 酵母基因名称见图２． ＮＵＰ６０ 图２部分酵母基因网络结构图 所推导的网络结构如图２所示．其中８个 基因与其他基因之间没有相互关系，其他２７个基 因之间有４３条有向边．实线表示基因之间是正调 控，共１８条，其中ＯＲＣ５与ＮＵＰ６０之间有相互作 用的正；虚线表示负，共２５条边． 为了验证算法推导的结果，本文使用ＹＥＡＳ— ＴＲＡＣＴ网站的Ｇｅｎｅｒａｔｅ Ｒｅｇｕｌａｔｉｏｎ Ｍａｔｒｉｘ工 具¨ ，得到了包含上述３５个基因的网络图， 其中包含３个基因对，分别是ＨＡＬ９ Ｈ１Ｍ１［１３ ３，ＰＨＤ１ＹＥＬ０７３Ｃ¨ ，ＰＨＤ１ ＦＹＶ６ １ ，这些基因关系都经过已有实验验 证．ＶＡＭ算法推导出的４３条边中，包含了ＨＡＬ９ ＨＩＭ１和ＰＨＤ１ＹＥＩ￣７３Ｃ这２组关 系对． 第６期 董自健，等：基于基因扰动及变分逼近技术的基因网络推断 １１５１ ５ 结语 本文把基因网络建模为没有潜变量的验 证性因子分析模型，并把模型分解为线性回归模 型，逐行求解．在处理线性回归模型时，首先采用变 分逼近技术推导参数的后验概率分布，在此基础 上，使用重要性抽样技术计算参数包含在线性模型 中的后验概率，并得出网络参数的加权平均．将 ＶＡＭ算法与ＡＬ—ｂａｓｅｄ算法进行比较，验证了 ＶＡＭ算法的有效性和可靠性．最后使用实验数 据，针对一种酵母基因网络进行了推导，得出 了一个具有４３条边的参考基因网络结构． 参考文献（Ｒｅｆｅｒｅｎｃｅｓ） ［１］Ｆｒｉｅｄｍａｎ Ｎ．Ｉｎｆｅｒｉｒｎｇ ｃｅｌｌｕｌａｒ ｎｅｔｗｏｒｋｓ ｕｓｉｎｇ ｐｒｏｂａｂｉ— ｌｉｓｉｔｃ ｇｒａｐｈｉｃａｌ ｍｏｄｅｌｓ［Ｊ］．Ｓｃｉｅｎｃｅ，２００４，３０３ （５６５９）：７９９—８０５． ［２］Ｓｃｈａｆｅｒ Ｊ，Ｓｗｉｍｍｅｒ Ｋ．Ａｎ ｅｍｐｉｒｉｃａｌ Ｂａｙｅｓ ａｐｐｒｏａｃｈ ｔｏ ｉｎｆｅｒｒｉｎｇ ｌａｒｇｅ—ｓｃａｌｅ ｇｅｎｅ ａｓｓｏｃｉａｔｉｏｎ ｎｅｔｗｏｒｋｓ［Ｊ］． Ｂｉｏｉｎｆｏｒｍａｔｉｃｓ，２００５，２１（６）：７５４—７６４． ［３］Ｓｃｈａｆｅｒ Ｊ，ＳＷｉｍｍｅｒ Ｋ．Ａ ｓｈｒｉｎｋａｇｅ ａｐｐｒｏａｃｈ ｔｏ ｌａｒｇｅ— ｓｃａｌｅ ｃｏｖａｒｉａｎｃｅ ｍａｗｉｘ ｅｓｔｉｍａｔｉｏｎ ａｎｄ ｉｍｐｆｉｃａｔｉｏｎｓ ｆｏｒ ｆｕｎｃｔｉｏｎａｌ ｇｅｎｏｍｉｃｓ［Ｊ］．Ｓｔａｔｉｓｔｉｃａｌ Ａｐｐｌｉｃａｔｉｏｎｓ ｉｎ Ｇｅ— ｎｅｔｉｃｓ ａｎｄ Ｍｏｌｅｃｕｌａｒ Ｂｉｏｌｏｇｙ，２００５，４（１）：１１７５— １１８９． ［４］Ｘｉｏｎｇ Ｈ，Ｃｈｏｅ Ｙ．Ｓｔｒｕｃｔｕｒｌａ ｓｙｓｔｅｍｓ ｉｄｅｎｔｉｉｆｃａｔｉｏｎ ｏｆ ｇｅｎｅｔｉｃ ｒｅｇｕｌａｔｏｒｙ ｎｅｔｗｏｒｋｓ［Ｊ］．Ｂｉｏｉｎｆｏｒｍａｔｉｃｓ，２００８， ２４（４）：５５３—５６０． ［５］Ｌｏｇｓｄｏｎ Ｂ Ａ，Ｍｅｚｅｙ Ｊ．Ｇｅｎｅ ｅｘｐｒｅｓｓｉｏｎ ｎｅｔｗｏｒｋ ｒｅ— ｃｏｎｓｔｒｕｃｔｉｏｎ ｂｙ ｃｏｎｖｅｘ ｆｅａｔｕｒｅ ｓｅｌｅｃｔｉｏｎ ｗｈｅｎ ｉｎｃｏｒｐｏｒａ— ｔｉｎｇ ｇｅｎｅｔｉｃ ｐｅｒｔｕｒｂａｔｉｏｎｓ［Ｊ／ＯＬ］．ＰＬｏＳ Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｂｉｏｌｏｇｙ，２０１０，６（１２）．ｈｔｔｐ：／／ｗｗｗ．ｐｌｏｓｃｏｍｐｂｉｏ１． ｏｒｇ／ａｒｔｉｃｌｅ／ｉｎｆｏ％３Ａｄｏｉ％２Ｆ１０．１３７１％２Ｆｊｏｕｍａ１．ｐｃｂｉ． １００１０１４． ［６］Ｌｏｇｓｄｏｎ Ｂ Ａ，Ｈｏｆｆｍａｎ Ｇ Ｅ，Ｍｅｚｅｙ Ｊ Ｇ．Ａ ｖａｒｉａｉｔｏｎａｌ Ｂａｙｅｓ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｆａｓｔ ａｎｄ ａｃｃｕｒａｔｅ ｍｕｌｉｔｐｌｅ ｌｏｃｕｓ ｇｅ— ｎｏｍｅ—ｗｉｄｅ ａｓｓｏｃｉａｔｉｏｎ ａｎａｌｙｓｉｓ［Ｊ／ＯＬ］．ＢＭＣ Ｂｉｏｉｎ— ｆｏｒｍａｔｉｃｓ，２０１０，１１．ｈｔｔｐ：／／ｗｗｗ．ｂｉｏｍｅｄｃｅｎｔｒ１ａ．ｃｏｍ／ １４７１—２１０５／１ １／５８／． ［７］Ｃａｒｂｏｎｅｔｔｏ Ｐ，Ｓｔｅｐｈｅｎｓ Ｍ．Ｓｃａｌａｂｌｅ ｖａｉｒａｔｉｏｎａｌ ｉｎｆｅｒ－ ｅｎｃｅ ｆｏｒ Ｂａｙｅｓｉａｎ ｖａｒｉａｂｌｅ ｓｅｌｅｃｔｉｏｎ ｉｎ ｒｅｇｒｅｓｓｉｏｎ，ａｎｄ ｉｔｓ ａｃｃｕｒａｃｙ ｉｎ ｇｅｎｅｔｉｃ ａｓｓｏｃｉａｉｔｏｎ ｓｔｕｄｉｅｓ［Ｊ］．Ｂａｙｅｓｉａｎ Ａｎａｌｙｓｉｓ，２０１２，７（１）：７３—１０７． ［８］ＪｏｒｄａｎＭＩ，ＧｈａｈｒａｍａｎｉＺ，ＪａａｋｋｏｌａＴ Ｓ，ｅｔａ１．Ａｎｉｎ— ｒｔｏｄｕｃｔｉｏｎ ｔｏ ｖａｒｉａｔｉｏｎａｌ ｍｅｔｈｏｄｓ ｆｏｒ ｇｒａｐｈｉｃａｌ ｍｏｄｅｌｓ ［Ｊ］．Ｍａｃｈｉｎｅ Ｌｅａｒｎｉｎｇ，１９９９，３７（２）：１８３—２３３． ［９］ＡｔｉｔｓａＨ．Ｉｎｄｅｅｐｎｄｅｎｔｆａｃｔｏｒａｎａｌｙｓｉｓ［Ｊ］．ＮｅｕｒａｌＣｏｍ— ｐｕｔａｔｉｏｎ，１９９９，１１（４）：８０３—８５１． ［１０］Ｂｅａｌ Ｍ Ｊ．Ｖａｒｉａｔｉｏｎａｌ ａｌｇｏｒｉｈｔｍｓ ｆｏｒ ａｐｐｒｏｘｉｍａｔｅ Ｂａｙｅｓｉｎａ ｉｎｆｅｒｅｎｃｅ［Ｄ］．Ｌｏｎｄｏｎ：Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｌｏｎ— ｄｏｎ，２００３． ［１１］Ｓｔｒｎａｇｅｒ Ｂ Ｅ，Ｆｏｒｒｅｓｔ Ｍ Ｓ，Ｄｕｎｎｉｎｇ Ｍ，ｅｔ ａ１．Ｒｅｌａ— ｉｔｖｅ ｉｍｐａｃｔ ｏｆｎｕｃｌｅｏｔｉｄｅ ａｎｄ ｃｏｐｙ ｎｕｍｂｅｒ ｖａｒｉａｔｉｏｎ ｏｎ ｇｅｎｅ ｅｘｐｒｅｓｓｉｏｎ ｐｈｅｎｏｔｙｐｅｓ［Ｊ］．Ｓｃｉｅｎｃｅ，２００７，３１５ （５８１３）：８４８—８５３． ［１２］Ｋｎｏｗｌｅｄｇｅ Ｄｉｓｃｏｖｅｒｙ＆Ｂｉｏｉｎｆｏｒｍａｔｉｃｓ Ｒｅｓｅａｒｃｈ Ｇｒｏｕｐ．ＹＥＡＳＴＲＡＣＴ［ＥＢ／ＯＬ］．［２０１３－０３—１０］．ｈｔ． ｔｐ：／／ｗｗｗ．ｙｅａｓｔｒａｃｔ．ｃｏｍ／ｉｎｄｅｘ．ｐｈｐ． ［１３］Ｌｅｅ Ｔ Ｉ，Ｒｉｎａｌｉｄ Ｎ Ｊ，Ｒｏｂｅ￣Ｆ，ｅｔ ａ１．Ｔｒｎａｓｃｒｉｐｔｉｏｎａｌ ｒｅｇｕｌａｔｏｒｙ ｎｅｔｗｏｒｋｓ ｉｎ Ｓａｃｃｈａｒｏｍｙｃｅｓ ｃｅｒｅｖｉｓｉａｅ［Ｊ］． Ｓｃｉｅｎｃｅ，２００２，２９８（５５９４）：７９９—８０４． ［１４］Ｂｏｍｅｍａｎ Ａ Ｒ，Ｌｅｉｇｈ—Ｂｅｌｌ Ｊ Ａ，Ｙｕ Ｈ，ｅｔ ａ１．Ｔａｒｇｅｔ ｈｕｂ ｐｒｏｔｅｉｎｓ ｓｅｒｖｅ ａｓ ｍａｓｔｅｒ ｒｅｇｕｌａｔｏｒｓ ｏｆ ｄｅｖｅｌｏｐｍｅｎｔ ｉｎ ｙｅａｓｔ［Ｊ］．Ｇｅｎｅｓ ａｎｄ Ｄｅｖｅｌｏｐｍｅｎｔ，２００６，２０（４）： ４３５—４４８． ｈｔｔｐ：／／ｊｏｕｒｎａ１．ｓｅｕ．ｅｄｕ．ｃｎ
