专利名称：一种逐步回归基因调控网络的推断方法
技术领域：
本发明涉及计算机领域，涉及一种逐步回归基因调控网络的推断方法。
背景技术：
从大规模基因表达测量数据集中推断基因调控网络无论从计算还是实验的角度 都是一项难题。主要原因在于，即使在不考虑具体的生化反应动力学特性(如基因调控网 络的动态变化)的情况下，由一定数目的基因构成的网络结构也有相当多种可能性。因此 基因网络预测或构建算法所面临的最大挑战在于数据维数过大，实验样本有限而与之相 关的网络结构却存在多种可能性。 从基因表达数据中成功地构建基因调控网络，往往需要复杂的计算，或者进行昂 贵费时的前期实验设计，以产生足够多的高质量数据来避免过多的计算。因此，如何解决这 两者间的矛盾，尽量达到一个平衡是大多数方法需要面临的难题。通常采用的解决方法是 将问题进行简化，例如，假设整个系统处于稳定状态，使用稳态数据来推测网络，或者给每 个基因与其它基因的连接个数加一个限制。 由于基因数N远大于单基因扰动实验次数M，要推断N个基因之间的调控关系，构 建基因调控网络非常困难。因此，以多重回归的网络识别算法，即NIR算法为代表的现有算 法提出假设每个基因最多有k个调控基因，也称为调控子。通过这样的假设，降低了数据 维度，将基因调控系数矩阵转换为稀疏矩阵，这确实符合生物网络为稀疏网络的理论，但同 时也对每个基因的调控子的选择从数量到质量都造成缺陷 1 、从生物学意义考虑，每个基因的调控子数量是不可能相同。如果对每个基因固 定最大非零连接的数量k，而事实上某个基因的调控子数量大于k，那么显然很多基因在入 选为调控子之前就失去了被选择为调控子的机会。 2、从统计意义考虑，在构建好基因调控网络的回归模型后，仅仅对整体模型进行 统计显著性检验是缺乏说服力的。因为整体的基因表达观测值对回归模型具有统计显著 性并不代表每个入选基因都分别具有统计显著性，要判断一个基因是否可以真正被选进模 型，还必须对选入进回归模型的基因作偏回归平方和的F检验。 3、从计算角度出发，若对所有可能的网络结构进行计算，这种采取遍历策略的串 行算法会造成过高的时空复杂度。
发明内容
本发明的目的在于为了解决以上诸多的技术问题提供一种逐步回归基因调控网 络的推断方法，该方法能推断出符合真实基因调控网络特性的稀疏网络，预测出网络中最 具统计显著性的边，同时避免过高的假阳性，克服高维小样本的实验数据问题，避免以往算 法中对目标基因强制设定最大连接数的不合理性。 为达到上述目的，本发明的构思是首先根据最优子集的回归原则，首先对目标基 因挑选最具统计显著性的调控子，然后对选入的调控子基因作偏回归平方和的F检验。每增加一个新的入选调控子，都要重新对先前入选的调控子进行F检验，若通过检验，则可以 继续被保留在回归模型中，否则将被剔除，重复这个过程直到未被选入的调控子中没有可 以再被选入的，并且已选入的调控子中没有可被剔除的，最后对所有选入的调控子用多重 线性回归计算基因间相互调控系数。 根据上述发明构思，本发明采用下述技术方案 —种逐步回归基因调控网络的推断方法，其特征在于具体操作步骤如下
A.读取基因表达数据矩阵和基因扰动数据矩阵； B.确定基因表达矩阵和基因扰动矩阵是否均为标准化数据，若基因表达数据矩阵 和基因扰动数据据矩阵均为标准化数据，则转步骤D，否则转步骤C ; C.对基因表达数据矩阵和基因扰动数据据矩阵分别进行数据归一化，使基因表达 数据矩阵和基因扰动数据据矩阵均构成标准化数据； D.分析标准化数据，用多重线性回归计算所有基因间相关系数矩阵； E.将基因间相关系数矩阵可视化成网络，得到基因调控网络图。 本发明的一种逐步回归基因调控网络的推断方法与已有技术相比较，具有如下面
显而易见的突出实质性特点和显著有点 1.该方法克服了高维小样本的实验数据问题。由于该方法采用的是每次只引入或 剔除一个调控子，因此在实验次数小于基因个数(M《N)的情况下，也能进行最优回归子集 的选择，因此解决了高维小样本实验数据的问题； 2.该方法避免了现有算法中对目标基因强制设定最大连接数的不合理性。本方法 不是事先统一规定每个基因的调控子个数，而是通过逐步回归的构思，为目标基因逐步选 择最具显著影响的调控子，因此每个基因的调控子个数是不同的，符合了基因调控网络的 真实情况； 3.随着基因调控网络规模的增大和网络稀疏度的增加，该方法在计算精度、计算 效率上均优于同类方法。
图1是本发明的一种逐步回归基因调控网络的推断方法的流程图。 图2是图1中步骤D所述的分析标准化数据，计算所有基因间相关系数矩阵的具
体流程图。 图3是图2中步骤D1所述的从N个基因中为基因i选择第一个调控子的具体流 程图。 图4是图2中步骤D2所述的对基因i选择第二个调控子的具体流程图。 图5是图2中步骤D3所述的终止引入新的调控子并完成回归模型的建立的具体
流程图。 图6是本发明的方法与现有的NIR算法分别对不同规模的基因调控网络进行构建 的敏感性对比图。 图7是本发明的方法与现有的NIR算法分别对不同规模的基因调控网络进行构建 的特异性对比图。 图8是本发明的方法与现有的NIR算法分别对不同规模的基因调控网络进行构建的覆盖率对比图。 图9是本发明的方法与现有的NIR算法分别对不同规模的基因调控网络进行构建 的精度对比图。 图10是本发明的方法与现有的NIR算法分别对不同规模的基因调控网络进行构 建的计算时间对比图。
具体实施例方式
以下结合附图对本发明的实施例作进一步详细说明。 本实施例中，本发明的逐步回归基因调控网络的推断方法的实验在上海大学系统 生物研究所的集群计算机上运行，该集群由14台IBM HS21刀片服务器和2台x3650服务 器组成计算和管理节点，网络连接采用千兆以太网和infiniband 2. 5G网。
本发明的一种逐步回归基因调控网络的推断方法，如图1所示，包括以下步骤
A.读取基因表达数据矩阵和基因扰动数据矩阵； B.确定基因表达矩阵和基因扰动矩阵是否均为标准化数据。若基因表达数据矩阵 和基因扰动数据据矩阵均为标准化数据，则转步骤D，否则转步骤C ; C.对基因表达数据矩阵和基因扰动数据据矩阵分别进行数据归一化，使基因表达 数据矩阵和基因扰动数据据矩阵均构成标准化数据，公式为 ~=^^, "l,2,…,W;y、l,2,…,M
v丄; &=^^， ^1，2"."A^ = 1,2,.."M 其中，Xij为基因表达数据矩阵中的元素； #为Xj的离差平方和； N为基因的数量，M为基因扰动的实验次数；
Pij为基因扰动数据矩阵中的元素； #为Pj的离差平方和。 D.分析标准化数据，用多重线性回归计算所有基因间相关系数矩阵，请参照图2， 其具体步骤如下 Dl.从N个基因中的基因i选择第一个调控子，请参照图3，其具体步骤如下
Dll.对基因i建立N个只包含一个自变量的回归方程，回归方程为
免.—')a:(/), j' = 1,2,...,W 其中J,为yi的估计值，yi代表基因i的扰动值；
aij为调控子j对基因i的调控系数； D12.对N个方程进行F检验，F检验值最大的基因j被暂时确定为基因i的第一 个调控子； D13.对第一个调控子j进行偏F检验，偏F检验值为F/1)，给定选入变量的显著性 水平Fa，判断《')> i^是否成立，如果《"> F^成立，则转步骤D14，如果《')> i^不成立，则认为没有任何基因为基因i的调控子，回归结束，偏F检验的公式为
<formula>formula see original document page 7</formula>
其中A SSR/1)为Xj的偏回归平方和； SSE(1)为残差平方差和； k为基因i的调控子个数，自由度为N_k ; D14.基因j作为第一个自变量被选入回归模型。 D2.对基因i选择第二个调控子，即从N-l个基因中选择基因i的调控子，请参照 图4，其具体步骤如下 D21.对基因i建立N-l个包含两个自变量的回归方程，其中一个自变量在步骤Dl 中已被确定，另外一个自变量需要在以下步骤中确定。 D22.对N-l个回归方程进行F检验，F检验值最大的方程中的另一个自变量，即基 因j被暂时确定为基因i的第二个调控子； D23.对基因i的第二个调控子进行偏F检验，偏F检验值为F/2)，给定选入变量的
显著性水平F^判断是否F/" 〉Fa，如果是，转步骤D24，否则停止变量选择； D24.基因i的第二个调控子，即变量，可被引入进回归方程中， D25.对基因i的第一个调控子重新进行偏F检验，给定剔除变量的显著性水平
F出，判断是否F/1) < & ，如果是，转步骤D26，否则结束该步骤。 D26.将基因i的第一个调控子从回归方程中剔除； D3.终止引入新的调控子，完成基因i的回归模型的建立，得到基因i与其调控子 之间的调控系数，请参照图5，其具体步骤如下 D31.继续引入新的变量进入回归方程，将其作为基因i的调控基因； D32.当既无显著的变量被选入回归方程，也无不显著的变量从回归方程中剔除
时，得到的调控子集合为最优回归子集； D33.得到包含了 k个调控子的关于基因i回归方程，求解得到基因间调控系数，包 含k个调控子的回归方程为
—/>, = 其中，Pi为1XM的向量，代表基因i在M次实验中的被扰动水平；
5,为1Xk的向量，代表k个基因对基因i的调控系数；
X为kXM的矩阵，代表k个基因在M次实验中的表达水平。
E.将基因间相关系数矩阵可视化成网络，得到基因调控网络图。
参照图6，示出了本发明的一种逐步回归基因调控网络的推断方法与现有技术的 代表性算法，即NIR算法，对不同规模的基因调控网络进行构建的敏感性对比图。本实验随 机产生了节点数分别为10、20、50、100、200和500的基因调控网络各IOO个，平均度为1， 分别使用一种逐步回归基因调控网络的推断方法和NIR算法进行基因调控网络的推断，统 计两种算法的敏感性。从实验所得的对比图表明，由于NIR算法的时间复杂性较高，对于节 点数为IOO及以上的网络无法计算，而计算节点数小于100的网络时，NIR算法的敏感性成 下降趋势，本发明的一种逐步回归基因调控网络的推断方法却表现出了良好而稳定的敏感 性，说明在敏感性方面本方法更优。
参照图7，示出了本发明的逐步回归基因调控网络的与现有技术的代表性算法，即 NIR算法，对不同规模的基因调控网络进行构建的特异性对比图。本实验随机产生了节点数 分别为10、20、50、100、200和500的基因调控网络各100个，平均度为l，分别使用一种逐 步回归基因调控网络的推断方法和NIR算法进行基因调控网络的推断，统计两种算法的特 异性。从实验所得的对比图表明，虽然计算节点数小于100的网络时，NIR算法的敏感性成 上升趋势，但由于NIR算法的时间复杂性较高，因此对于节点数为100及以上的网络无法进 行计算，而本发明的一种逐步回归基因调控网络的推断方法却表现出了良好而稳定的特异 性，说明在特异性方面本方法更优。 参照图8，示出了本发明的一种逐步回归基因调控网络的推断方法与现有技术的 代表性算法，即NIR算法，对不同规模的基因调控网络进行构建的覆盖率对比图。本实验随 机产生了节点数分别为10、20、50、100、200和500的基因调控网络各100个，平均度为l，分 别使用和NIR算法进行基因调控网络的推断，统计两种算法的覆盖率。从实验所得的对比 图，虽然计算节点数小于100的网络时，NIR算法的覆盖率成上升趋势，但由于NIR算法的 时间复杂性较高，因此对于节点数为IOO及以上的网络无法进行计算，而本发明的一种逐 步回归基因调控网络的推断方法却表现出了良好而稳定的覆盖率，说明在覆盖率方面本方 法更优。 参照图9，示出了本发明的一种逐步回归基因调控网络的推断方法与现有技术的 代表性算法，即NIR算法，对不同规模的基因调控网络进行构建的精度对比图。本实验随机 产生了节点数分别为10、 、50、100、200和500的基因调控网络各100个，平均度为l，分别 使用一种逐步回归基因调控网络的推断方法和NIR算法进行基因调控网络的推断，统计两 种算法的特异性。从实验所得的对比图表明，由于NIR算法的时间复杂性较高，因此对于节 点数为IOO及以上的网络无法进行计算，且本发明的一种逐步回归基因调控网络的推断方 法在特异性方面高于NIR算法。 参照图IO，指出了本发明的一种逐步回归基因调控网络的推断方法与现有技术的 代表性算法，即NIR算法，对不同规模的基因调控网络进行构建的计算时间对比图。本实验 随机产生了 9个测试用网络，节点数分别为5、 10、 20、 30、40、 50、60、80和100，且网络的平均 度各不相同，使用一种逐步回归基因调控网络的推断方法和NIR算法并分别统计它们的计 算时间。随着网络规模的逐渐增大，本方法相对于NIR方法的计算速度优势逐渐突出，尤其 当网络节点数大于50时，本方法运行时间的增加明显小于NIR。因此，本方法在计算速度方 面更优。 综上所述，从图6 9的对比图表明，说明本发明的一种逐步回归基因调控网络的 推断方法的整体计算精度高于现有的代表算法NIR算法，从图10的对比图表明，该方法在 计算速度方面也更优。 本发明所述的方法并不限于具体实施方式
中所述的实施例，本领域技术人员根据 本发明的技术方案得出的其它的实施方式，同样属于本发明的技术创新范围。
权利要求
一种逐步回归基因调控网络的推断方法，其特征在于，该方法具体步骤如下A.读取基因表达数据矩阵和基因扰动数据矩阵；B.确定基因表达矩阵和基因扰动矩阵是否均为标准化数据，若基因表达数据矩阵和基因扰动数据据矩阵均为标准化数据，则转步骤D，否则转步骤C；C.对基因表达数据矩阵和基因扰动数据据矩阵分别进行数据归一化，使基因表达数据矩阵和基因扰动数据据矩阵均构成标准化数据；D.分析标准化数据，用多重线性回归计算所有基因间相关系数矩阵；E.将基因间相关系数矩阵可视化成网络，得到基因调控网络图。
2. 根据权利要求1所述的一种逐步回归基因调控网络的推断方法，其特征在于，上述 步骤C中所述的将基因表达数据矩阵和基因扰动数据矩阵分别进行数据归一化，其归一化其中，Xij为基因表达数据矩阵中的元素^丄5为Xj的离差平方和；N为基因的数量，M为基因扰动的实验次数；Pij为基因扰动数据矩阵中的元素；V4f为Pj的离差平方和。
3. 根据权利要求2所述的一种逐步回归基因调控网络的推断方法，其特征在于，上述 步骤D中所述的分析标准化数据，用多重线性回归计算所有基因间相关系数矩阵，其具体 步骤如下Dl. N个基因中的基因i选择第一个调控子； D2.对基因i选择第二个调控子；D3.终止引入新的调控子，完成基因i的回归模型的建立，得到基因i与其调控子之间 的调控系数；
4. 根据权利要求3所述的一种逐步回归基因调控网络的推断方法，其特征在于，上述 步骤D1中所述的为基因i选择第一个调控子，其具体步骤如下Dll.对基因i建立N个只包含一个自变量的回归方程，自变量为基因i的调控子，回归其中J,为yi的估计值，yi代表基因i的扰动值；aij为调控子j对基因i的调控系数； D12.对N个方程进行F检验，F检验值最大的基因j被暂时确定为基因i的第一个调 控子；D13.对第一个调控子j进行偏F检验，偏F检验值为F/"，给定选入变量的显著性水平 Ftt，判断i^ > Fj4是否成立，如果《"> 4成立，则转步骤D14，如果《"> Fjj不成立，则认为没有任何基因为基因i的调控子，结束，偏F检验的公式为<formula>formula see original document page 2</formula>公式为方程为其中ASSR/1)为Xj的偏回归平方和；SSE("为残差平方差和；k为基因i的调控子个 数，自由度为N_k ;D14.基因j作为第一个自变量被选入回归模型。
5. 根据权利要求4所述的一种逐步回归基因调控网络的推断方法，其特征在于，上述步骤D2中所述的为基因i选择第二个调控子，其具体步骤如下D21.对基因i建立N-l个包含两个自变量的回归方程，其中一个自变量在步骤Dl中已 被确定，另外一个自变量需要在以下步骤中确定；D22.对N-1个回归方程进行F检验，F检验值最大的方程中的另一个自变量，即基因j 被暂时确定为基因i的第二个调控子；D23.对基因i的第二个调控子进行偏F检验，偏F检验值为F/2)，给定选入变量的显著性水平F进，判断是否F/2) 〉F进，如果是，转步骤D24，否则停止变量选择； D24.基因i的第二个调控子，即变量，可被引入进回归方程中；D25.对基因i的第一个调控子重新进行偏F检验，给定剔除变量的显著性水平F^判 断是否《(1) 〈&，如果是，转步骤D26，否则结束该步骤； D26.将基因i的第一个调控子从回归方程中剔除。
6. 根据权利要求5所述的一种逐步回归基因调控网络的推断方法，其特征在于，上述 步骤D3中所述的终止引入新的调控子，并完成回归模型的建立，其具体步骤如下D31.继续引入新的变量进入回归方程，将其作为基因i的调控基因；D32.当既无显著的变量被选入回归方程，也无不显著的变量从回归方程中剔除时，得到的调控子集合为最优回归子集；D33.得到包含了 k个调控子的关于基因i回归方程，求解可得到基因间调控系数，包含了k个调控子的回归方程为其中，Pi为lXM的向量，代表基因i在M次实验中的被扰动水平；4为1Xk的向量，代表k个基因对基因i的调控系数；X为kXM的矩阵，代表k个基因在M次实验中的表达水平。
全文摘要
本发明公开了一种逐步回归基因调控网络的推断方法，该方法步骤如下A.读取基因表达数据矩阵和基因扰动数据矩阵；B.确定基因表达矩阵和基因扰动矩阵是否均为标准化数据；C.对基因表达数据矩阵和基因扰动数据据矩阵分别进行数据归一化均构成标准化数据；D.分析标准化数据，计算所有基因间相关系数矩阵；E.将基因间相关系数矩阵可视化成网络，得到基因调控网络图。该方法能进行最优回归子集的选择，因此解决了高维小样本实验数据的问题；2.为目标基因逐步选择最具显著影响的调控子，符合了基因调控网络的真实情况；3.随着基因调控网络规模的增大和网络稀疏度的增加，在计算精度、计算效率上均优于同类方法。
文档编号G06F19/00GK101719195SQ200910199850
公开日2010年6月2日 申请日期2009年12月3日 优先权日2009年12月3日
发明者何冰, 宋安平, 张律文, 张武, 肖梅, 谢江 申请人:上海大学