一种逐步回归基因调控网络的推断方法 【技术领域】
本发明涉及计算机领域，涉及一种逐步回归基因调控网络的推断方法。
背景技术
从大规模基因表达测量数据集中推断基因调控网络无论从计算还是实验的角度都是一项难题。主要原因在于，即使在不考虑具体的生化反应动力学特性(如基因调控网络的动态变化)的情况下，由一定数目的基因构成的网络结构也有相当多种可能性。因此基因网络预测或构建算法所面临的最大挑战在于：数据维数过大，实验样本有限而与之相关的网络结构却存在多种可能性。
从基因表达数据中成功地构建基因调控网络，往往需要复杂的计算，或者进行昂贵费时的前期实验设计，以产生足够多的高质量数据来避免过多的计算。因此，如何解决这两者间的矛盾，尽量达到一个平衡是大多数方法需要面临的难题。通常采用的解决方法是将问题进行简化，例如，假设整个系统处于稳定状态，使用稳态数据来推测网络，或者给每个基因与其它基因的连接个数加一个限制。
由于基因数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.对基因表达数据矩阵和基因扰动数据据矩阵分别进行数据归一化，使基因表达数据矩阵和基因扰动数据据矩阵均构成标准化数据，公式为：
xij=xij-x‾jLjj(x),i=1,2,...,N;j=1,2,...,M]]>
pij=pij-p‾jLjj(p),i=1,2,...,N;j=1,2,...,M]]>
其中，xij为基因表达数据矩阵中的元素；
为xj的离差平方和；
N为基因的数量，M为基因扰动的实验次数；
pij为基因扰动数据矩阵中的元素；
为pj的离差平方和。
D.分析标准化数据，用多重线性回归计算所有基因间相关系数矩阵，请参照图2，其具体步骤如下：
D1.从N个基因中的基因i选择第一个调控子，请参照图3，其具体步骤如下：
D11.对基因i建立N个只包含一个自变量的回归方程，回归方程为：
y^i=aij(1)xij(1),j=1,2,...,N]]>
其中，为yi的估计值，yi代表基因i的扰动值；
aij为调控子j对基因i的调控系数；
D12.对N个方程进行F检验，F检验值最大的基因j被暂时确定为基因i的第一个调控子；
D13.对第一个调控子j进行偏F检验，偏F检验值为Fj(1)，给定选入变量的显著性水平F进，判断是否成立，如果成立，则转步骤D14，如果不成立，则认为没有任何基因为基因i的调控子，回归结束，偏F检验的公式为：
Fj(1)=ΔSSRj(1)/1SSE(1)/(N-k)]]>
其中：ΔSSRj(1)为xj的偏回归平方和；
SSE(1)为残差平方差和；
k为基因i的调控子个数，自由度为N-k；
D14.基因j作为第一个自变量被选入回归模型。
D2.对基因i选择第二个调控子，即从N-1个基因中选择基因i的调控子，请参照图4，其具体步骤如下：
D21.对基因i建立N-1个包含两个自变量的回归方程，其中一个自变量在步骤D1中已被确定，另外一个自变量需要在以下步骤中确定。
D22.对N-1个回归方程进行F检验，F检验值最大的方程中的另一个自变量，即基因j被暂时确定为基因i的第二个调控子；
D23.对基因i的第二个调控子进行偏F检验，偏F检验值为Fj(2)，给定选入变量的显著性水平F进，判断是否如果是，转步骤D24，否则停止变量选择；
D24.基因i的第二个调控子，即变量，可被引入进回归方程中，
D25.对基因i的第一个调控子重新进行偏F检验，给定剔除变量的显著性水平F出，判断是否如果是，转步骤D26，否则结束该步骤。
D26.将基因i的第一个调控子从回归方程中剔除；
D3.终止引入新的调控子，完成基因i的回归模型的建立，得到基因i与其调控子之间的调控系数，请参照图5，其具体步骤如下：
D31.继续引入新的变量进入回归方程，将其作为基因i的调控基因；
D32.当既无显著的变量被选入回归方程，也无不显著的变量从回归方程中剔除时，得到的调控子集合为最优回归子集；
D33.得到包含了k个调控子的关于基因i回归方程，求解得到基因间调控系数，包含k个调控子的回归方程为：
-pi=a~iX]]>
其中，pi为1×M的向量，代表基因i在M次实验中的被扰动水平；
为1×k的向量，代表k个基因对基因i的调控系数；
X为k×M的矩阵，代表k个基因在M次实验中的表达水平。
E.将基因间相关系数矩阵可视化成网络，得到基因调控网络图。
参照图6，示出了本发明的一种逐步回归基因调控网络的推断方法与现有技术的代表性算法，即NIR算法，对不同规模的基因调控网络进行构建的敏感性对比图。本实验随机产生了节点数分别为10、20、50、100、200和500的基因调控网络各100个，平均度为1，分别使用一种逐步回归基因调控网络的推断方法和NIR算法进行基因调控网络的推断，统计两种算法的敏感性。从实验所得的对比图表明，由于NIR算法的时间复杂性较高，对于节点数为100及以上的网络无法计算，而计算节点数小于100的网络时，NIR算法的敏感性成下降趋势，本发明的一种逐步回归基因调控网络的推断方法却表现出了良好而稳定的敏感性，说明在敏感性方面本方法更优。
参照图7，示出了本发明的逐步回归基因调控网络的与现有技术的代表性算法，即NIR算法，对不同规模的基因调控网络进行构建的特异性对比图。本实验随机产生了节点数分别为10、20、50、100、200和500的基因调控网络各100个，平均度为1，分别使用一种逐步回归基因调控网络的推断方法和NIR算法进行基因调控网络的推断，统计两种算法的特异性。从实验所得的对比图表明，虽然计算节点数小于100的网络时，NIR算法的敏感性成上升趋势，但由于NIR算法的时间复杂性较高，因此对于节点数为100及以上的网络无法进行计算，而本发明的一种逐步回归基因调控网络的推断方法却表现出了良好而稳定的特异性，说明在特异性方面本方法更优。
参照图8，示出了本发明的一种逐步回归基因调控网络的推断方法与现有技术的代表性算法，即NIR算法，对不同规模的基因调控网络进行构建的覆盖率对比图。本实验随机产生了节点数分别为10、20、50、100、200和500的基因调控网络各100个，平均度为1，分别使用和NIR算法进行基因调控网络的推断，统计两种算法的覆盖率。从实验所得的对比图，虽然计算节点数小于100的网络时，NIR算法的覆盖率成上升趋势，但由于NIR算法的时间复杂性较高，因此对于节点数为100及以上的网络无法进行计算，而本发明的一种逐步回归基因调控网络的推断方法却表现出了良好而稳定的覆盖率，说明在覆盖率方面本方法更优。
参照图9，示出了本发明的一种逐步回归基因调控网络的推断方法与现有技术的代表性算法，即NIR算法，对不同规模的基因调控网络进行构建的精度对比图。本实验随机产生了节点数分别为10、20、50、100、200和500的基因调控网络各100个，平均度为1，分别使用一种逐步回归基因调控网络的推断方法和NIR算法进行基因调控网络的推断，统计两种算法的特异性。从实验所得的对比图表明，由于NIR算法的时间复杂性较高，因此对于节点数为100及以上的网络无法进行计算，且本发明的一种逐步回归基因调控网络的推断方法在特异性方面高于NIR算法。
参照图10，指出了本发明的一种逐步回归基因调控网络的推断方法与现有技术的代表性算法，即NIR算法，对不同规模的基因调控网络进行构建的计算时间对比图。本实验随机产生了9个测试用网络，节点数分别为5、10、20、30、40、50、60、80和100，且网络的平均度各不相同，使用一种逐步回归基因调控网络的推断方法和NIR算法并分别统计它们的计算时间。随着网络规模的逐渐增大，本方法相对于NIR方法的计算速度优势逐渐突出，尤其当网络节点数大于50时，本方法运行时间的增加明显小于NIR。因此，本方法在计算速度方面更优。
综上所述，从图6～9的对比图表明，说明本发明的一种逐步回归基因调控网络的推断方法的整体计算精度高于现有的代表算法NIR算法，从图10的对比图表明，该方法在计算速度方面也更优。
本发明所述的方法并不限于具体实施方式中所述的实施例，本领域技术人员根据本发明的技术方案得出的其它的实施方式，同样属于本发明的技术创新范围。