CN1717585A - 利用微分方程从时序基因表达数据推断基因调控网络 - Google Patents

Abstract

通过利用时间过程表达数据和一组线性微分方程，本方法的实施例能够被用来估计生物体基因之间的网络关系。通过确定哪些元素为0或者哪些元素在所研究的条件下不显著地改变，能够使用Aikaike信息准则和掩模工具来减少矩阵中元素的数目。使用最大似然估计和新的统计方法来估计所提出的网络关系的显著性。

Description

利用微分方程从时序基因表达数据推断基因调控网络

相关申请：

本申请要求在35 U.S.C§119(e)于2002年11月25日提交的美国临时专利申请No.60/428,827的优先权。其全部内容结合于此作为参考。

技术领域

本发明涉及确定生物体基因之间关系的方法。特别地，本发明包括利用微分方程线性***，从时间过程基因表达数据(time course gene expressiondata)推断基因调控网络的新方法。

背景技术

在生命科学、医学、药物发明和发展及制药行业中，当前研究和发展的最重要的方面中的一个为，对发展用于解释大量的第一手数据并基于这些数据引出结论的方法和器件的需求。生物信息学对***生物学的理解做出了重大贡献并保证产生对生命***成分之间的复杂关系更重大的理解。特别地，随着用于迅速探测被表达的基因和定量基因表达的新方法的出现，生物信息学甚至在不确定的情况下能够用来预测潜在的治疗目标，特定的基因在生物体生态学中可能起的确切作用。

遗传***的模拟是***生物学的中心主题。因为模拟能够基于生物学知识，通过预测或者推断以前未知的关系，网络评估方法能够支持生物模拟。

特别地，微阵列技术允许来自多种生物体的很多基因表达的研究。大量的第一手数据能够从来自生物体的许多基因得到，并且能够通过干扰或者通过突变、疾病或药物来研究基因表达。在特定的疾病中或者响应特定的干扰而特定的基因表达增加的发现，可以使人相信基因直接牵涉到疾病过程或者药物响应。然而，在生物有机体中，基因很少通过任何这样的干扰独立地调控，因为许多基因能够被一特定的干扰影响。因为很多不同的基因可以被如此影响，在这样的研究中了解基因之间的因果关系是非常困难的。因而正花费大量的努力来发展用于确定基因之间因果关系的方法，哪种基因对一种生物现象重要，哪种基因的表达对正在研究的生物过程不重要。尽管这样的不重要的基因表达作为生物或病理生理条件下的标记可能有用，如果这样的基因对于生理学或者病理生理学条件不重要，基于这样的基因发展药物可能不值得努力。相反，对于对一个过程识别为重要的基因，药物或者其它干扰的发展可能对发展用于与改变的基因表达关联的条件的治疗至关重要。

微阵列技术允许在同一时间测量很多基因的基因表达水平。使用互补DNA(cDNA)，微阵列分析能够容易地实现，但是RNA微阵列也能够用来研究基因表达。随着可利用的基因表达数据的数量迅速地增加，用来分析这种数据的技术仍在发展中。数学方法被日益地用来确定表达基因之间的关系。然而，从基因表达数据精确地导出基因调控网络会很困难。

在时序基因表达测量中，能够通过在少量的时间点测量基因表达水平来研究基因表达的时效模式。例如，已经在酿酒酵母菌(Saccharomycescerevisiae)的细胞周期测量了周期变化的基因表达水平(见参考文献1)。基因对缓慢改变的环境的响应已经在同种酵母菌的二次转变期间进行了测量(见参考文献2)。其它的实验测量响应生物体环境突然变化的时效基因表达模式。作为例子，在外部光强突然移位之后，测量了藻青菌集胞藻6803(Synechocystis sp.PCC 6803)的基因表达响应。

已经提出了几种方法来从表达数据推断基因相互关系(见参考文献2、5和6)。在簇分析中，基因基于他们的基因表达谱的相似点集群在一起。从测量的基因表达数据来推断布尔(Boolean)或贝氏(Bayesian)网络(见参考文献7、8、9、10、11和美国专利申请No：10/259.723及申请题目“NonlinearModeling of Gene Networks From Time Series Gene Expression Data”，2003年11月18日提交；律师案号No：GENN 1008 US1DBB，两个申请都全部结合于此作为参考)和使用任意的微分方程***模拟基因表达数据(见参考文献12)先前已经公开了。然而，为了可靠地推断这样一个任意的微分方程***，需要长序列的时序基因表达数据，这在当前常常还不可获得。

发明内容

为了克服现有技术的缺点，在本发明的某些方面，我们发展了使用微分方程线性***及从基因表达数据导出的信息来推断基因网络的方法。这种途径保持了微分方程的定量及内在因果关系的优点，而又足够简单易于计算处理。我们也发展了用于检验涉及基因调控网络的假定的新方法。

附图说明

参照其具体的例子，描述本发明的各个方面。本发明的其它特征能够通过参照附图理解，其中：

图1描绘来自枯草芽孢杆菌(Bacillus subtilis)的5条基因簇的基因表达随时间的图表。

图2描绘利用本发明的方法导出的在图1中描绘的5条基因簇的基因网络。

具体实施方式

Chen从理论上考虑了使用线性微分方程模拟生物数据(见参考文献13)。在这个模型中，mRNA和蛋白质浓度都由线性微分方程***描述。这样的***能够描述为：

$\frac{d}{\mathrm{dt}}\underset{\‾}{x}\left(t\right)=\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}\·\underset{\‾}{x}\left(t\right),---\left(1\right)$

其中矢量 x(t)作为时间的函数包含mRNA和蛋白质浓度，并且矩阵

为常量，以[秒]^-1为单位。这个方程能够被看作布尔网络模型的一般形式，其中水平的数目无限而不是二元的。

在cDNA微阵列实验中，当蛋白质浓度未知时，通过测量相应的mRNA浓度通常仅只确定基因表达水平。因此我们集中在仅描述基因相互作用的微分方程***上。矩阵元Λ_ij于是代表基因j对基因i的作用，[Λ_ij]^-1为反应时间。

为了从所测量的数据推断微分方程***中的系数，先前已提出离散微分方程***(见参考文献13)，代替所测量的mRNA和蛋白质浓度，并求解因而产生的方程线性***以得到线性微分方程***中的系数Λ_ij。该方程***通常是不确定的。使用附加的必要条件，即基因调控网络应该是稀疏的，Chen表示模型能够以o(m^h+1)次构造，这里m是基因的数目以及h是***中每个微分方程所允许的非零系数的数目(见参考文献13)。

参数h特别选定，这会有两个意外的结果。因为矩阵 中每行将精确地具有h个非零元素，网络中每个基因或者蛋白质具有h个亲本基因或蛋白质，并且因此没有基因或者蛋白质能够存在于网络的顶点。其次，每个基因不可避免地为一个反馈环的成员。虽然反馈环可能存在于基因调控网络中，他们的存在应该从测量的数据中来确定而不是人为地创造。

另一方面，贝式网络，不允许环的存在。贝氏网络依赖于所估计的网络的联合概率分布，以便可分解为条件概率分布的积。仅在不存在环时这种分解是可能的。我们进一步注意到贝氏网络往往包含许多参数，并且因此需要大量的数据以进行可靠的估计。

因此，我们旨在找到允许网络中存在环的方法，但是不一定需要他们存在。利用方程1，通过限制可能出现在***中的非零系数的数目，我们构造了一个稀疏矩阵。代替特别选定这个数目，我们通过利用Akaike信息准则(Akaike’s Information Criterion)(AIC)从数据中估计在相互作用矩阵中哪些系数为0，允许基因调控路径的数目对于每个基因不相同。

我们的方法的多个方面能够用于寻找单独基因之间的网络，也可用于寻找基因簇之间的调控网络。作为例子，我们能够利用枯草芽孢杆菌的时间过程数据推断基因簇之间的基因调控网络。簇能够利用k平均聚类演算法(k-means clustering algorithm)创造。簇的生物功能能够从属于每个簇的基因的功能范畴来确定。

在一些实施例中，我们按照微分方程(方程1)线性***考虑m个基因之间的调控网络，这里矢量 x(t)包含在t时刻m个基因的表达率。这个微分方程***能够求解为：

$\underset{\‾}{x}\left(t\right)=\exp \left(\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}t\right)\·{\underset{\‾}{x}}_0,---\left(2\right)$

其中x₀包含在零时刻的基因表达率。在这个方程中，矩阵指数按照泰勒(Taylor)展开定义为：

$\exp \left(\underset{\‾}{\underset{\‾}{A}}\right)\≡\underset{i=0}{\overset{\∞}{\mathrm{\Σ}}}\frac{1}{i!}{\underset{\‾}{\underset{\‾}{A}}}^i,---\left(3\right)$

由于方程2非线性地依赖于

按照测量的数据 x(t)求解

将很困难。通过由差分方程：

$\frac{\mathrm{\Δ}\underset{\‾}{x}}{\mathrm{\Δt}}=\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}\·\underset{\‾}{x},---\left(4\right)$

或

$\underset{\‾}{x}(t+\mathrm{\Δt})-\underset{\‾}{x}\left(t\right)=\mathrm{\Δt}\·\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}\·\underset{\‾}{x}\left(t\right),---\left(5\right)$

替换微分方程(方程1)，可得到近似解。差分方程4或5为Chen所考虑的形式(见参考文献13)。为了统计地确定矩阵

的稀疏性，我们明确地添加一个误差 ε(t)，其将不变地出现在数据中：

$\underset{\‾}{x}(t+\mathrm{\Δt})-\underset{\‾}{x}\left(t\right)=\mathrm{\Δt}\·\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}\·\underset{\‾}{x}\left(t\right)+\underset{\‾}{\mathrm{\ϵ}}\left(t\right),---\left(6\right)$

通过利用这个方程，我们能够按照多维线性马尔可夫(Markov)模型有效地描绘基因调控网络。

能够假定误差有一个不依赖于时间的正态分布，如下所示：

$f(\underset{\‾}{\mathrm{\ϵ}}\left(t\right);{\mathrm{\σ}}^2)={\left(\frac{1}{\sqrt{2{\mathrm{\π\σ}}^2}}\right)}^m\exp \{-\frac{\underset{\‾}{\mathrm{\ϵ}}{\left(t\right)}^T\·\underset{\‾}{\mathrm{\ϵ}}\left(t\right)}{2{\mathrm{\σ}}^2}\},---\left(7\right)$

在所有的时刻对于所有的基因有一个相等的标准偏差σ。对于一系列在时刻t_i，i∈{1，...，n)，n个时间点的时序测量值 x _i的对数似然函数于是为：

$L(\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}},{\mathrm{\σ}}^2)=-\frac{\mathrm{nm}}{2}\ln \left[2\mathrm{\π}{\mathrm{\σ}}^2\right]-\frac{1}{2{\mathrm{\σ}}^2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\underset{\‾}{\overset{^}{\mathrm{\ϵ}}}}_i^T\·{\underset{\‾}{\overset{^}{\mathrm{\ϵ}}}}_i,---\left(8\right)$

其中

${\underset{\‾}{\overset{^}{\mathrm{\ϵ}}}}_i={\underset{\‾}{x}}_i-{\underset{\‾}{x}}_{i-1}-(t_i-t_{i-1})\·\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}\·{\underset{\‾}{x}}_{i-1},---\left(9\right)$

为在时刻t_i从测量数据估计的测量误差。

方差σ²的最大似然估计能够通过关于σ²取对数似然函数的最大值得到。这得出：

${\widehat{\mathrm{\σ}}}^2=\frac{1}{\mathrm{nm}}\underset{i=1}{\overset{1}{\mathrm{\Σ}}}{\underset{\‾}{\overset{^}{\mathrm{\ϵ}}}}_i^T\·{\underset{\‾}{\overset{^}{\mathrm{\ϵ}}}}_i\·---\left(10\right)$

将这个式子代入对数似然函数(方程8)得到：

$L(\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}},{\mathrm{\σ}}^2={\widehat{\mathrm{\σ}}}^2)=-\frac{\mathrm{nm}}{2}\ln \left[2\mathrm{\π}{\widehat{\mathrm{\σ}}}^2\right]-\frac{\mathrm{nm}}{2},---\left(11\right)$

为了得到矩阵 的最大似然估计 我们使用方程9将总平方误差

表示为：

${\widehat{\mathrm{\σ}}}^2=\frac{1}{\mathrm{nm}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}[({\underset{\‾}{x}}_i^T-{{\underset{\‾}{x}}_{i-1}}^T)\·({\underset{\‾}{x}}_i-{\underset{\‾}{x}}_{i-1})+{(t_i-t_{i-1})}^2{{\underset{\‾}{x}}_{i-1}}^T\·{\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}}^T\·{\underset{\‾}{x}}_{i-1}-2({\underset{\‾}{x}}_i^T-(t_i-t_{i-1}){{\underset{\‾}{x}}_{i-1}}^T)\·\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}\·{\underset{\‾}{x}}_{i-1}],$

(12)

并且关于 求导。我们得到关于

的线性方程：

$\underset{\‾}{\overset{^}{\underset{\‾}{\mathrm{\Λ}}}}=\underset{\‾}{\underset{\‾}{B}}\·{\underset{\‾}{\underset{\‾}{A}}}^{-1},---\left(13\right)$

其中矩阵 和 定义为：

$\underset{\‾}{\underset{\‾}{A}}\≡\underset{i=1}{\overset{n}{\mathrm{\Σ}}}[{(t_i-t_{i-1})}^2\·{\underset{\‾}{x}}_{j-1}\·{{\underset{\‾}{x}}_{i-1}}^T],---\left(14\right)$

$\underset{\‾}{\underset{\‾}{B}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}[(t_i-t_{i-1})\·({\underset{\‾}{x}}_i-{\underset{\‾}{x}}_{i-1})\·{{\underset{\‾}{x}}_{i-1}}^T],---\left(15\right)$

当不存在误差时，估计矩阵 等于真实矩阵 从生物学中我们得知基因调控网络并且因此 为稀疏的。然而，由于噪声的存在，在估计矩阵

中所有的元素都可能非零，即使在真实矩阵

中相应的元素为零。

在某些实施例中，如果总平方差因此增加很小，如方程12所给出的，能够设置该矩阵元等于零。正式地，我们将使用Akaike信息准则(见参考文献15、16)：

AIC＝2·[估计的模型的对数似然值]+2·[估计的参数的数目]，    (16)来决定哪个矩阵元素应该设置为等于零。通过比较模型中使用的参数的数目和所估计的模型中的总误差，AIC能够用来避免模型到数据的过拟合(overfitting)。具有最低AIC的模型被认为是最佳的。AIC基于信息论并且广泛地用于统计模型识别，尤其用于时序模型拟合(见参考文献17)。

我们于是能够使用掩模

来设置

的矩阵元等于零：

这里о表示哈达玛(Hadamard)(基于元素的(element-wise))乘积，并且掩模

是一个元素不为1就为0矩阵。相应的总平方误差

能够通过在方程12中用

替换 得到。给出掩模 能够通过求解一组方程18最小化总平方误差，

如果

${\underset{\‾}{\underset{\‾}{M}}}_{\mathrm{ij}}=1:{[{\underset{\underset{\‾}{\‾}}{\overset{^}{\mathrm{\Λ}}}}^{\′}\·\underset{\‾}{\underset{\‾}{A}}]}_{\mathrm{ij}}=B_{\mathrm{ij}};$

如果 $M_{\mathrm{ij}}=0:{\widehat{\mathrm{\Λ}}}_{\mathrm{ij}}^{\′}=0;---\left(18\right)$

得出最大似然估计

在这个方程中，

和

利用测量的基因表达水平 x _i从方程14和15确定。我们于是通过将根据方程11的替换估计的对数似然函数代入到方程16中来计算相应于

的AIC：

估计的参数为

和我们所允许的非零的矩阵

元素。从这个方程，能够看到当平方误差减少时，随着非零元素的数目增加，AIC可能增加。现在可以通过找到对于AIC来说产生最低值的掩模

根据基因表达数据推断基因调控网络。

对于任何但是最平凡的情形，可能的掩模 的数目非常大，使得进行彻底地搜寻以找到最佳掩模不可行。替代地，我们能够使用贪婪(greedy)搜寻方法。开始，以对于每个掩模元素为1和0相等的概率，我们能够随机地选择掩模。通过改变每个掩模元素M_ij，能够减小AIC。这个过程能够继续直到找到最后的掩模，对于这个掩模不能实现进一步AIC的减小。这种算法能够从不同的(例如，随机的)初始掩模重复开始，并且能够用来确定具有最小的相应AIC的最终掩模

如果在几十次试验中都找到这个最佳的掩模，能够合理地断定没有更好的掩模存在。

我们已经描述并证明了从测量基因表达数据以微分方程线性***的形式来推断基因调控网络的方法。由于进行典型测量的时间点的数目有限，得到基因调控网络通常是一个不确定的难题。因为生物学上作为结果的基因调控网络预期为稀疏的，我们设置某些矩阵项等于零，并且仅利用非零项推断网络。非零项的数目，并且因此，网络的稀疏性，利用Akaike信息准则从数据确定，而不利用任何特别参数。

按照微分方程描述基因网络至少具有三个优点。第一，微分方程组描述基因之间的因果关系：系数矩阵的系数Λ_ij确定基因j对基因i的作用。第二，它以清楚的数字形式描述基因相互作用。第三，由于微分方程***中呈现的大量信息，其它网络形式能够容易地从它导出。另外，我们能够将推断的网络连接到其它分析或者可视化工具，如GON(Genomic ObjectNet)。

在先前描述的方法中，或者不能够找到任何环(如在贝氏网络模型中)或者该方法人为地在网络中产生回路。尽管这里描述的方法允许网络中出现环，但并不需要他们的存在。仅只在有数据保证时可找到环。例如，当利用MMGE培养基中枯草芽孢杆菌的时间过程数据推断基因簇之间的调控网络时，我们发现一些簇是环的部分，而其它的不是(参看下面的例子和图2)。

如果基因的数目m等于或者大于试验的数目n，方程18中的矩阵

是奇异的。这个问题于是就不确定，并且能够找到具有零总误差

和为-∞的AIC的相互作用矩阵 通过对足够小数量的基因或者基因簇应用这种方法，或者通过限制网络中亲本的数目，能够避免我们方法的这种不足。

用于评估网络关系的统计显著性的方法

在本发明的其它实施例中，提供了用于确定网络关系分析的统计显著性的方法。在无效假设(null hypothesis)下，能够假定基因不被实验操作影响。在不同时间点的测量对数比(log-ratio)于是相等。我们能够进一步假设对数比具有零平均的正态分布。在一些情形，统计测试，如学生的t-检验(Student’s t-test)，会在每个时间点进行以确定哪个对数比显著地不同于零。然而，学生的t-检验会不可靠，因为仅只少许测量的数据组。因此，在包括每个时间点仅两次测量的数据组的一些实施例中，我们设计了一种新的统计测试，将多个时间点的测量合并在一起。特别地，如例2中所示，我们将这种方法应用到来自所有的8个时间点的数据上。令人满意的是该方法能够用于其它类型的实验，并且将在下面描述。

实现本方法的步骤在下面描述。

步骤1：在每个时间点，计算平均对数比为

${\stackrel{\‾}{x}}_{\mathrm{ji}}=\frac{1}{2}\underset{k=\mathrm{1,2}}{\mathrm{\Σ}}{x}_{\mathrm{ji}}\left[k\right].---\left(21\right)$

在无效假设下， x_j·(在一个时间点两个基因表达对数比的平均)是具有零平均正态分布并且估计的标准偏差为 的随机变量。

步骤2：接着从所有的测量估计标准偏差(如，对于如例1中包括的数据组8×2＝16)：

${\widehat{\mathrm{\σ}}}_{j\left|H_0\right|}\sqrt{\frac{1}{2n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=\mathrm{1,2}}{\mathrm{\Σ}}{\left(x_{\mathrm{ji}}\right[k\left]\right)}^2},---\left(20\right)$

其中x_ji[k]表示对于基因j在时间点i测量k的数据值。

步骤3：就绝对值而言比测量值 x_ji大的 x_j·的联合概率为：

$P=\underset{i=1}{\overset{n}{\mathrm{\Π}}}{P}_i=\underset{i=1}{\overset{n}{\mathrm{\Π}}}p\left(\right|{\stackrel{\‾}{x}}_j\·|>|{\stackrel{\‾}{x}}_{\mathrm{ji}}\left|\right)=\underset{i=1}{\overset{n}{\mathrm{\Π}}}[1-\mathrm{erf}\left(\frac{\left|{\stackrel{\‾}{x}}_{\mathrm{ji}}\right|}{{\widehat{\mathrm{\σ}}}_{j\left|H_0\right|/\sqrt{2}}}\right)],---\left(22\right)$

其中，erf是误差函数。对于这个乘积中的单个因素P_i，我们通常会挑选一个显著性水平α，并且如果P_i＜α则舍弃无效假设。

步骤4：采用准则：P＜αⁿ用于舍弃无效假设。这允许我们通过利用关于那个基因所有可利用的数据来确定在实验期间该基因的表达水平是否显著地改变。

步骤5：确定一个基因改变的表达水平是否显著。

用于确定基因之间网络关系的方法和新的统计方法能够用在研究，生物医学科学，包括诊断学中，以便发展新的诊断及用于在制药工业中选择先导化合物。

例子

下面的例子旨在阐明本发明的实施例，而不限制其范围。能够发展其它的实施例而不背离本发明的范围，并且本发明的方法及其变体能够在不存在不适当的实验下用来推断枯草芽孢杆菌或其它生物体中不同基因的调控网络。所有这样的实施例被认为是本发明的部分。

例1：枯草芽孢杆菌中的基因网络

最近在枯草芽孢杆菌的MMGE基因表达实验中测量了用于利用基因表达数据寻找基因调控网络的本发明实施例。MMGE是包含葡萄糖和谷氨酰胺(如碳和氮源)的合成的基本培养基。在这种介质中，诱导了小分子生物合成所需要的基因的表达，如氨基酸。在这个实验中，在一个小时的时间间隔中于八个时间点测量了4320 ORF的表达水平，在每个时间点进行两次测量。

数据准备和分析

为减少出现在数据中的测量噪声的影响，每个基因的表达水平与测量背景水平比较。无论在红色或是绿色通道中具有低于平均背景水平的平均基因表达水平的基因从分析中除去。

对3823个剩下的基因应用全局归一化，并且计算了基因表达率以2为底的对数。我们对测量的对数比应用统计测试来确定它们是否显著地不同于零。

以上描述的方法的流程图重现在下面的总结中。

步骤1：计算在每个时间点，每个基因表达的平均对数比；

步骤2：从所有的测量计算标准偏差；

步骤3：计算联合概率；

步骤4：采用用于统计显著性的准则；及

步骤5：确定一个基因改变的表达水平是否显著。

在这个例子中，我们选择一个显著性水平α＝0.00025以便预期的假阳性数(0.00025×3823＝1)可接受。通过应用这个准则到这3823个基因，我们发现684个基因显著地受到影响。

例2：枯草芽孢杆菌基因的聚类

使用k平均聚类，枯草芽孢杆菌的这684个基因顺序地群集成5组。使用欧几里得(Euclidean)距离测量基因之间的距离，而簇的质心(centroid)定义为簇中所有基因的中部(median)。挑选簇的数目以避免显著的重叠。k平均算法从不同的随机初始聚类开始而重复1000000次。找到最佳解81次。

完整的聚类结果可在下列网站得到：

http：//bonsai.ims.u-tokyo.ac.jp/-mdehoon/publications/Subtilis/ clusters.html.

为了确定所创造的簇的生物学功能，我们考虑每个簇中所有基因在枯草芽孢杆菌数据库中的功能范畴。表1列出所形成的5个簇的主要功能范畴。

图1表示对于每个簇，基因表达的对数比作为时间的函数。尽管在该时间过程期间，簇I、II和V的表达水平相当大地改变，簇II和III具有相当稳定的表达水平。特别地，簇IV能够看作为一个包罗万象的簇，被分配到其上的基因不能适合其它簇。

                                表1

利用k平均聚类创造的5个簇的主要功能范畴。

簇	基因的数目	主要功能范畴
			I	42	2.2：11基因；1.1：9基因
II	62	1.2：15基因；2.2：12基因
			III	187	5.1：30基因；6.0：23基因；1.2：22基因
IV	343	5.1：40基因；5.2：39基因；1.2：33基因
			V	50	1.2：15基因；2.1.1：15基因

                        基因的功能范畴

参考法国巴斯德研究所(Institute Pasteur)枯草芽孢杆菌数据库的功能范畴

1.1：      细胞壁1.2：      运输/结合蛋白质和脂蛋白2.1.1：    碳水化合物代谢和相关分子----特定通路2.2：      氨基酸代谢和相关分子5.1：      与来自枯草芽孢杆菌的未知蛋白质类似的5.2：      与来自其它生物体的未知蛋白质类似的6.0：      无相似性

图1表示对于每个簇，作为时间的函数的基因表达的对数比，其从测量的基因表达数据确定。

分段网络构造

从那12个基因的测量对数比中，我们构造矩阵

和 并计算矩阵 从一个随机的初始掩模开始，计算掩模

的过程重复1000次。找到最佳的解55次。因此，不太可能存在其它具有较低AIC的掩模。注意到可能的掩模的总数目为2²⁵＝33,554,432。

图2表示所找到的网络。在网络中簇的亲本的数目在0和5之间变化。簇III和IV作为网络的顶点出现，而簇I、II和V连接到一个环中。注意到这种网络既不能由先前提出的方法(见参考文件13)产生，也不能由贝氏网络模型产生。

网络中两个最强的相互作用分别是簇IV对簇V和簇II的正效应和负效应。簇II和V的基因表达水平的相反行为极大可能地由簇IV引起，而不是簇II和V之间的直接相互作用。

图2表示如根据MMGE时间过程数据和本发明的方法所确定的5个基因簇之间的网络。如由相互作用矩阵

中相应元素所给出的一样，数值表示一个基因簇对另一个簇的影响有多强。实际上，这个矩阵代表基因表达水平有多迅速地彼此响应。作为例子，如果簇II、III和IV的表达水平不改变，则簇I的基因表达水平的变化将引起簇V的表达水平在1/(5.0小时^-1)＝12分钟之内相当大地改变。

参考文献

1. P.T.Spellman，G.Sherlock，M.Q.Zhang，V.R.Iyer，K.Anders，M.B.Eisen，P.O.Brown，D.Botstein，and B.Futcher，″Comprehensive identification of cellcycle-regulated genes of the yeast Saccharomyces cerevisiae by microarrayhybridization″Mol.Biol.Cell 9(1998)3273-3297.

2. J.L.DeRisi，V.R.Iyer，and P.O.Brown，″Exploring the metabolic and geneticcontrol of gene expression on a genomic scale″Science 278(1997)680-686.

3. Y.Hihara，A.Kamei，M.Kanehisa，A.Kaplan，and M.Ikeuchi，″DNAmicroarray analysis of cyanobacterial gene expression during acclimation to high light″The Plant Cell 13(2001)793-806.

4. M.J.L.de Hoon，S.Imoto，and S.Miyano，″Statistical analysis of a small set oftime-ordered gene expression data using linear splines″Bioinformatics，in press.

5. M.B.Eisen，P.T.Spellman，P.O.Brown，and D.Botstein，″Cluster analysis anddisplay of genome-wide expression patterns″Proc.Natl.Acad.Sci.USA 95(1998)14863-14868.

6. P.Tamayo，D.Slonim，J.Mesirov，Q.Zhu，S.Kitareewan，E.Dmitrovsky，E.S.Lander，and T.R.Golub，″Interpreting patterns of gene expression with self-organizingmaps：Methods and application to hematopoietic differentiation″Proc.Natl.Acad.Sci.USA 96(1999)2907-02912.

7. S.Liang，S.Fuhrman，and R.Somogyi，″REVEAL，a general reverseengineering algorithm for inference of genetic network architectures″Proc.Pac.Symp.on Biocomputing 3(1998)18-29.

8. T.Akutsu，S.Miyano，and S.Kuhara，″Inferring qualitative relations in geneticnetworks and metabolic pathways″Bioinformatics 16(2000)727-734.

9. N.Friedman，M.Linial，I.Nachman，and D.Pe′er，″Using Bayesian networks toanalyze expression data″J.Comp.Biol.7(2000)601-620.

10. S.Imoto，T.Goro，and S.Miyano，″Estimation of genetic networks andfunctional structures between genes by using Bayesian networks and nonparametricregression″Proc.Pac.Symp.on Biocomputing 7(2002)175-186.

11. S.Imoto，S.-Y.Kim，T.Goto，S.Aburatani，K.Tashiro，S.Kuhara，and S.Miyano，″Bayesian network and nonparametric heteroscedastic regression for nonlinearmodeling of genetic network″Proc.IEEE Computer Society Bioinformatics Conference(2002)219-227.

12. E.Sakamoto and H.Iba，″Evolutionary inference of a biological network asdifferential equations by genetic programming″Genome Informatics 12(2001)276-277.

13. T.Chen，H.L.He，and G.M.Church，″Modeling gene expression withdifferential equations″Proc.Pac.Symp.on Biocomputing 4(1999)29-40.

14. R.A.Hom and C.R.Johnson，Matrix Analysis.Cambridge University Press，Cambridge，UK(1999).

15. H.Akaike，″Information theory and an extension of the maximum likelihoodprinciple″Research Memorandum No.46，Institute of Statistical Mathematics，Tokyo(1971).In B.N.Petrov and F.Csaki(editors)，2nd Int.Symp.on Inf.Theory.AkadémiaiKiadó，Budapest (1973)267-281.

16. H.Akaike，″A new look at the statistical model identification″IEEE Trans.Automat.Contr.AC-19(1974)716-723.

17. M.B.Priestley，Spectral Analysis and Time Series，Academic Press，London(1994).

18. Microbial Advanced Database Organization (Micado).http：//www-mig.versailles.inra.fr/bdsi/Micado/.

19. I.Moszer，P.Glaser，and A.Danchin，″SubtiList：a relational database for theBacillus subtilis genome″Microbiology 141(1995)261-268.

20. I.Moszer，″The complete genome of Bacillus subtilis：From sequence annotationto data management and analysis″FEBS Letters 430(1998)28-36

21. T.W.Anderson and J.D.Finn，The New Statistical Analysis of Data.SpringerVerlag，New York(1996).

22. H.Matsuno，A.Doi，Y.Hirata，and S.Miyano，″XML documentation ofbiopathways and their simulation in Genomic Object Net″Genome Informatics 12(2001)54-62.Genomic Object Net is available at http：//www.GenomicObject.net.

Claims (28)

1.一种用于推断基因之间网络关系的方法，包括：

(a)提供一个用于生物体的一组基因的定量时间过程数据库，所述库包括基于在所述基因组中的每个基因表达的时间过程的表达结果，量化所述基因彼此间的平均作用和每个时间点的变异性的测量；

(b)从所述库创造一个稀疏矩阵，所述矩阵具有从其上除去的零系数；

(c)从所述矩阵产生一组线性微分方程；及

(d)求解所述方程组以产生所述网络关系。

2.如权利要求1所述的方法，其中利用Akaike信息准则(AIC)来识别所述零系数。

3.如权利要求1和2中任何一个所述的方法，其中所述微分方程为：

$\frac{d}{\mathrm{dt}}\underset{\‾}{x}\left(t\right)=\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}\·\underset{\‾}{x}\left(t\right),$

其中矢量 x(t)作为时间的函数包含所表达的cDNA的总量，并且矩阵

为常量，以[秒]^-1为单位。

4.如权利要求1-3中任何一个所述的方法，其中所述矩阵包含元素Λ_ij，Λ_ij代表基因j对基因i的作用，并且其中[Λ_ij]^-1代表所述基因j对基因i的作用的反应时间。

5.如权利要求1-4中任何一个所述的方法，其中所述求解的微分方程为：

$\underset{\‾}{x}\left(t\right)=\exp \left(\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}t\right)\·{\underset{\‾}{x}}_0,$

6.如权利要求1-5中任何一个所述的方法，其中所述指数Λt(exp(Λ))根据以下表达求解：

$\exp \left(\underset{\‾}{\underset{\‾}{A}}\right)\≡\underset{i=0}{\overset{\∞}{\mathrm{\Σ}}}\frac{1}{i!}{\underset{\‾}{\underset{\‾}{A}}}^i.$

7.如权利要求1-6中任何一个所述的方法，其中所述微分方程通过求解以下差分方程估计：

$\frac{\mathrm{\Δ}\underset{\‾}{x}}{\mathrm{\Δt}}=\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}\·\underset{\‾}{x}.$

8.如权利要求1-7中任何一个所述的方法，其中，根据以下公式，所述稀疏矩阵进一步包括一个估计的误差：

$\underset{\‾}{x}(t+\mathrm{\Δt})-\underset{\‾}{x}\left(t\right)=\mathrm{\Δt}\·\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}\·\underset{\‾}{x}\left(t\right)+\underset{\‾}{\mathrm{\ϵ}}\left(t\right).$

9.如权利要求1-8中任何一个所述的方法，其中，根据以下公式，所述误差具有不依赖于时间的正态分布：

$f(\underset{\‾}{\mathrm{\ϵ}}\left(t\right);{\mathrm{\σ}}^2)={\left(\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}^2}}\right)}^m\exp \{-\frac{\underset{\‾}{\mathrm{\ϵ}}{\left(t\right)}^T\·\underset{\‾}{\mathrm{\ϵ}}\left(t\right)}{2{\mathrm{\σ}}^2}\},$

其中标准偏差σ在所有时刻对于所述基因的每一个相等。

10.如权利要求1-9中任何一个所述的方法，其中，根据以下公式，方差σ²的最大似然估计通过关于σ²取对数似然函数的最大值确定：

${\widehat{\mathrm{\σ}}}^2=\frac{1}{\mathrm{nm}}\underset{i=1}{\overset{1}{\mathrm{\Σ}}}{{\widehat{\underset{\‾}{\mathrm{\ϵ}}}}_i}^T\·{\widehat{\underset{\‾}{\mathrm{\ϵ}}}}_i\·$

11.如权利要求1-10中任何一个所述的方法，其中所述方差σ²根据以下公式确定：

${\widehat{\mathrm{\σ}}}^2=\frac{1}{\mathrm{nm}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}[({{\underset{\‾}{x}}_i}^T-{{\underset{\‾}{x}}_{i-1}}^T)\·({\underset{\‾}{x}}_i-{\underset{\‾}{x}}_{i-1})+{(t_i-t_{i-1})}^2{{\underset{\‾}{x}}_{i-1}}^T\·{\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}}^T\·\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}\·{\underset{\‾}{x}}_{i-1}-2({{\underset{\‾}{x}}_i}^T-(t_i-t_{i-1}){{\underset{\‾}{x}}_{i-1}}^T)\·\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}\·{\underset{\‾}{x}}_{i-1}],$

12.如权利要求2-11中任何一个所述的方法，其中所述AIC根据以下公式取最小值：

AIC＝2·[估计的模型的对数似然值]+2·[估计的参数的数目].

13.如权利要求1-12中任何一个所述的方法，其中根据以下公式，掩模

用来设置

的矩阵元素等于0：

这里°表示基于元素的乘积，并且掩模

是一个其元素不为1就为0矩阵。

14.如权利要求13所述的方法，其中通过应用由最小化以下公式所产生的掩模 来设置矩阵元素为0：

如果 ${\underset{\‾}{\underset{\‾}{M}}}_{\mathrm{ij}}=1:{[{\hat{\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}}}^{\′}\·\underset{\‾}{\underset{\‾}{A}}]}_{\mathrm{ij}}=B_{\mathrm{ij}};$

如果 ${\underset{\‾}{\underset{\‾}{M}}}_{\mathrm{ij}}=0:{\hat{\underset{\‾}{\underset{\‾}{\mathrm{\Λ}}}}}_{\mathrm{ij}}^{\′}=0,$ 从而得出最大似然估计

15.如权利要求2所述的方法，其中所述AIC根据以下公式最小化：

16.如权利要求13所述的方法，其中选择所述掩模

以最小化AIC。

17.一种介质，其含有一个或多个利用上述权利要求所述的方法获得的基因之间网络关系的结果，所述结果存储在所述介质上。

18.一种用于确定网络关系的统计显著性的方法，包括：

(a)在每个时间点对于每个基因计算表达的平均对数比；

(b)从所有的测量中计算标准偏差；

(c)计算联合概率；及

(d)采用用于统计显著性的准则。

19.如权利要求18所述的方法，其中所述步骤(a)利用以下公式确定：

${\stackrel{\‾}{x}}_{\mathrm{ji}}=\frac{1}{2}\underset{k=\mathrm{1,2}}{\mathrm{\Σ}}{x}_{\mathrm{ji}}\left[k\right].$

20.如权利要求18或者19所述的方法，其中所述步骤(b)利用以下公式确定：

${\widehat{\mathrm{\σ}}}_{j\left|H_0\right|}=\sqrt{\frac{1}{2n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=\mathrm{1,2}}{\mathrm{\Σ}}{\left(x_{\mathrm{ji}}\right[k\left]\right)}^2},$

其中x_ji[k]为对于基因j在时间点i测量k的数据值。

21.如权利要求18-20中任何一个所述的方法，其中就绝对值而言比测量值 x_ji大的 x_j·的联合概率利用以下公式计算：

$P=\underset{i=1}{\overset{n}{\mathrm{\Π}}}{P}_i=\underset{i=1}{\overset{n}{\mathrm{\Π}}}p\left(\right|{\stackrel{\‾}{x}}_{j\·}|>|{\stackrel{\‾}{x}}_{\mathrm{ji}}\left|\right)=\underset{i=1}{\overset{n}{\mathrm{\Π}}}[1-\mathrm{ref}\left(\frac{\left|{\stackrel{\‾}{x}}_{\mathrm{ji}}\right|}{{\widehat{\mathrm{\σ}}}_{j\left|H_0\right|/\sqrt{2}}}\right)],$

其中erf为误差函数。

22.如权利要求18-21中任何一个所述的方法，其中选择一个显著性水平α。

23.如权利要求18-22中任何一个所述的方法，其中如果P_i＜α则舍弃无效假设。

24.如权利要求18-23中任何一个所述的方法，其中如果P＜αⁿ则舍弃无效假设，其中n为进行基因表达计算的时间点的数目。

25.一种用于确定网络关系的统计显著性的方法，包括：

(a)利用以下公式，在每个时间点对于每个基因计算表达测量的平均对数比：

${\stackrel{\‾}{x}}_{\mathrm{ji}}=\frac{1}{2}\underset{k=\mathrm{1,2}}{\mathrm{\Σ}}{x}_{\mathrm{ji}}\left[k\right].$

(b)利用以下公式，计算所述测量的标准偏差：

${\widehat{\mathrm{\σ}}}_{j\left|H_0\right|}=\sqrt{\frac{1}{2n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=\mathrm{1,2}}{\mathrm{\Σ}}{\left(x_{\mathrm{ji}}\right[k\left]\right)}^2}.$

其中 x_ji[k]为对于基因j在时间点i测量k的数据值。

(c)利用以下公式，计算就绝对值而言比测量值 x_ji大的 x_j·的联合概率：

$P=\underset{i=1}{\overset{n}{\mathrm{\Π}}}{P}_i=\underset{i=1}{\overset{n}{\mathrm{\Π}}}p\left(\right|{\stackrel{\‾}{x}}_{j\·}|>|{\stackrel{\‾}{x}}_{\mathrm{ji}}\left|\right)=\underset{i=1}{\overset{n}{\mathrm{\Π}}}[1-\mathrm{ref}\left(\frac{\left|{\stackrel{\‾}{x}}_{\mathrm{ji}}\right|}{{\widehat{\mathrm{\σ}}}_{j\left|H_0\right|/\sqrt{2}}}\right)],$

其中erf为误差函数；及

(d)应用用于统计显著性的准则以确定是否舍弃无效假设。

26.如权利要求25所述的方法，其中如果P＜αⁿ则舍弃无效假设，其中n为进行基因表达计算的时间点的数目。

27.一种用于推断如在这里充分描述的基因网络的方法。

28.一种用于确定如在这里充分描述的网络关系的统计显著性的方法。