应用推导方法和生物学约束条件估计基因网络的制作方法

专利名称：应用推导方法和生物学约束条件估计基因网络的制作方法
技术领域：
本发明涉及使用应用生物学约束条件的数学推导方法估计生物系统中的基因网络。具体的，本发明涉及使用贝叶斯估计(Bayesian estimation)估计基因网络，其搜索范围通过对遗传体系生物学约束条件的认识进行限定。
背景技术：
由基因表达测定值推导基因网络是系统生物学中的主要挑战。如果基因网络能够被正确推导，将能够导致更好的理解细胞过程，并由此应用于药物开发、疾病研究和其它领域。由表达水平的测定值估计基因网络是生物信息学研究的一个焦点。
常用的模拟基因网络的方法是贝叶斯网络(Bayesian networks)。在贝叶斯网络中，网络的行为通过生物体相关基因群中所有基因的表达水平的联合概率分布进行描述。联合概率分布可以通过使用被我们称为网络的定向非循环图(directed acyclic graph)分解到条件概率分布中。给出基因表达测量数据，许多评分函数(score function)已经被提议用于选择网络。致力于使用已有知识的新评分函数正在开发。给定数据，则应用基于网络可能性的评分函数评价网络。选择了评分函数，估计基因网络的任务则是找到具有最小分值的网络。
然而，由于NP-hard的基因网络估计的问题和超指数大小的搜索范围，一些人使用启发式算法，例如贪心算法(greedy algorithms)或模拟退火(simulated annealing)以估计基因网络。由于启发式算法不提供任何准确性的保证，因而仍然不清楚所用的计算方法的网络估计的偏差程度。
对于大约30个或更少基因的基因网络的一些数据，能够推导出最优基因网络模型。这个结果适合任何评估函数sG×2G，其指定基因g和g的一组亲本的分值。然而，虽然最优的基因网络模型是最可能的模型，仍存在非常不同的对于一些数据具有几乎同样可能性的模型，尤其是由于基因表达测量值通常仅提供基因网络的部分信息。另外，即使对于仅10个基因的基因网络，存在大约4.17*1018种可能的网络模型。因此，即使最优的基因网络模型通常也不匹配目标基因网络。
发明概述在本发明的某些方面，我们给出有效估计基因网络的方法，虽然其在有限的意义上提供所估计的网络的最优性。第一，创建可能的基因网络的推导模型。第二，选择搜索范围(space)的生物学相关子范围(subspace)。然后，我们通过反复应用计算最优小基因网络的算法找到所选子范围中最优的解决方案。子范围的选择可以通过使用生物学知识或者，如果没有已知的现有知识，例如通过使用聚类方法。
在某些实施方式中，该方法允许通过以生物学上意义微小但计算有效的方式限定搜索范围来利用生物学知识。该方法改善关于限定的搜索范围的最优性，并因此提供了通过启发/经验加以实施的那部分估计和最优实施的那部分估计之间的清楚的区别。
在其它实施方式中，使用贝叶斯网络，基因网络的行为被模拟为所有基因的联合概率分布。这提供了基因相互作用的非常通用的模拟。联合概率分布可以被分解为条件概率的产物，其表示基因被其它基因调节。该分解可以被表示为定向非循环图。贝叶斯网络模型可以被用于寻找生物学上看似可信的基因网络。
基因表达信息可以使用本领域公知的方法产生，包括例如微阵列分析。由于这些方法是本领域所公知的，因此我们不需要进一步描述它们。但是，可以理解的是网络关系的可靠性依赖于所获得的基因表达信息的准确性和/或可靠性。因此，在一些实施方式中，需要进行重复实验和对照研究以增加基因表达数据的可靠性。
在另一方面，我们提供用于以可能性顺序列举最优和次优网络的理论基础和算法。在另一实施方式中，我们比较了最优基因网络模型和有关现有基因网络的已知知识。在另一实施方式中，我们提供了获取最优基因网络模型的可靠部分的方法并应用该方法于枯草芽孢杆菌和大肠杆菌的已知数据从而比较相关物种的基因网络模型。
我们的结果显示我们的方法识别的部分网络明显比最优网络模型本身更好的符合现有知识。由于我们将我们的方法以最好的n基因网络模型为基础，因而可以从这些估计中比来自单独基于启发式方法的方法更加可靠地获得结论。
本发明还包括执行基因网络推导分析的系统，其至少包括适合接受来自大量不同基因的基因表达数据的输入设备，和适合运行推导基因网络关系的程序的处理器。在另外的实施方式中，系统包括用于适合显示、打印或发送基因网络分析结果给远程位置的输出设备。
附图简述参考发明实施方式的特别描述和附图描述本发明，其中

图1描述通过实施例1的算法3预测的基因网络。
图2a描述由枯草芽孢杆菌获得的时程数据。
图2b描述由枯草芽孢杆菌获得的破坏体(disruptant)数据。
图3描述从大肠杆菌获得的数据。
图4描述由大肠杆菌获取的基因网络关系基序(motif)。
图5描述由枯草芽孢杆菌获得的基因网络关系基序。
发明详述细胞中基因相互作用的阐明是近年来研究的重要焦点(20)。在某些实施方式中，可以使用贝叶斯网络模拟基因网络(4，5，6，11，12，14，16，19)。在贝叶斯网络中，网络的行为通过所有基因的表达水平的联合概率分布进行描述。联合概率分布必须使用定向非循环图(directed acyclic graph)或“网络”被分解到条件概率分布中。给出基因表达测量数据，许多评分函数(scorefunction)已经被提议用于选择网络(3，4，11)。致力于使用以前知识的新评分函数正在开发(13，23)。选择了评分函数，估计基因网络的任务则是找到具有最小分值的网络。
然而，由于NP-hard的基因网络估计问题(1)和超指数大小的搜索范围(17)，研究者使用启发式算法如贪心算法(greedy algorithms)或模拟退火(simulated annealing)以估计基因网络。由于启发式算法不提供任何准确性的保证，因而仍然不清楚计算方法对网络估计的偏差程度。
在本发明的某些方面，我们给出有效估计基因网络的方法，虽然其在有限的意义上提供网络的最优性。第一，创建可能的基因网络的推导模型。第二，选择搜索范围的生物学相关子范围(subspace)。然后，我们通过反复应用最优计算小基因网络的算法找到所选子范围中最优的解决方案(15)。子范围的选择可以通过使用生物学知识或者，如果没有已知的现有知识，例如通过使用聚类方法。
在某些实施方式中，该方法允许通过以生物学上意义微小但计算有效的方式限定搜索范围来利用生物学知识。该方法改善关于限定的搜索范围的最优性，并因此提供了通过启发/经验加以实施的那部分估计和最优实施的那部分估计之间的清楚的区别。
实施例下述实施例被用于说明本发明的各方面，而不被用于进行限制。此外，对于每个实施例，参考文献被引用在括号中和显示在每个实施例末尾的参考文献列表中。本发明的其它实施方式可以由本领域技术人员在不需要不适当实验的情况下实施。所有的那些实施方式都被包括在本发明的范围内。
实施例1使用生物学约束条件寻找最优基因网络在选择的子范围中寻找最优解决方案所需的运行时间为O(n)，其中n为基因的数量，因此允许任意大的基因网络的估计。我们注意到使用功能性sG×2G→IR的任意评分函数s可以使用该方法，其适用于贝叶斯网络框架中所有评分函数以及大部分其它评分函数的情况。
我们执行我们的方法并将它应用于酵母数据和枯草芽孢杆菌数据。枯草芽孢杆菌的估计网络和生物学知识的比较显示了明显的一致。
必要的符号在第2部分介绍，方法在第3部分详细解释，网络估计的结果在第4部分讨论。
1.1.预备整个实施例中，我们将使用G表示基因的集合，为其预测基因网络，并假定我们被提供了评分函数sG×2G→IR，其赋予分值给基因g∈G和一组亲本基因AG。对于网络N，N的分值被定义为分值(N)＝def∑g∈Gs(g，PN(g))，其中PN(g)表示N中g的亲本的集合。我们从(15)引入一些符号。
定义1.1我们为所有的g∈G和AG定义FG×2G→IR为F(g，A)＝defminBAs(g，B)。F(g，A)是当选择限定为A时g的亲本的最好的选择。假设我们使用π表示一组AG的数列π{1，...|A|}→A，∏A表示所有A的数列的集合。
定义1.2假设AG。我们对于所有的π∈∏A，定义QA∏A→IR为QA(π)=defΣg∈AF(g,{h∈A|π-1(h)<π-1(g)}).---(1)]]>定义1.3对于所有AG，我们定义M2G→AG∏A为M(A)=defminπ∈ΠAQA(π)---(2)]]>在(15)中，显示QA(π)为所有网络的范围中最优的网络，其中M(A)为拓扑排序，1为数列，其是基于A的至少一个最优网络的拓扑排序。
我们从(15)规定下列算法，我们将使用其为本工作的子程序。
算法1.1步骤1为所有g∈G计算F(g，φ)＝s(g，φ)。
步骤2对所有AG，A≠φ和所有g∈G计算F(g，A)为min{s(g，A)，minα∈AF(g，A-{a})}。
步骤3设定M(φ)＝φ。
步骤4对所有AG，A≠φ，进行下列两个步骤步骤计算4ag*＝arg mingA(F(g，A-{g}+QA-{g}(M(A-{g})))。
步骤对所有1≤i＜|A|，设定M(A)(i)＝M(A-{g*})(i)，和4bM(A)(|A|)＝g*。
步骤5返回QG(M(G))。
算法1.1在步骤1和步骤2计算函数F，在步骤3和步骤4利用F计算函数M。在步骤5，M被用于计算最优网络的分值。本算法被用于为多至约30个基因的基因网络寻找最优基因网络模型，但是对于大基因网络变得可靠性降低。
定理1.1(15)算法1.1使用O(n·2n)动态设计步骤寻找最优网络。
1.2方法我们首先证明定理，然后我们从它推导出两个用于估计大基因网络的算法。
1.2.1理论基础下面的推导为本实施例中我们的方法预备了基础。
定理1.2假设m，c∈IN。对每个g∈G，假设CgG-{g}为一组候选亲本。假定下列两种情况成立1.对所有g∈G，|Cg|≤m。
2.由Cg，g∈G诱导的图形的强连接组分不大于c。
然后关于所选候选亲本的最优网络能够在时间O(|G|)内被找到。
验证假设L＝(G，E)表示基于被Cg，g∈G诱导的基因集合G的图形。那表示E完全含有边缘(h，g)，对其h∈Gg.成立。假设S1，...SnG表示L的强连接组分。现在假设L’＝({S1，...，Sn}，E’)表示强连接组分的图形，其中E.定义为E’＝def{(Si，Sj)||i≠j，gεSj(g，h)εE}。
然后，通过强连接组分的定义，在L’.W.l.o.g.中没有循环。让我们假设S1，...，Sn为对L’的拓扑分类，其表示在L’中没有边缘(Si，Sj)，其中j＜i。
稍微修改的算法1.1可以被用于根据所选候选亲本，为G计算最优网络。为了将算法1.1应用于强连接组分Si，我们以下列方式修改算法1.在步骤1和步骤2的F的计算中，我们仅仅为所有的g∈Si和所有的ACg计算F(g，A)。
2.我们替换步骤4a中的F(g，A-{g})项(见算法1的定义)为F(g，(Cg＝Si)∪(Cg∩A))。
这两个改变将候选亲本的限定引入算法，并允许g∈Si以具有Si以外的亲本。这不影响算法1.1的正确性，因为可以以对算法1一样的方式对修改的算法1.1进行定理1.1(15)的验证。我们以随机顺序将这种方式修改的算法1.1用于每个强连接组分Si，对每个i∈1，...，n产生部分网络Ni＝(G，Ei)。然后我们返回网络N=(G,Uni=1Ei).]]>通过定理1.1，每个部分网络Ni最优化。从L’的非循环性(acyclicity)得出N是非循环的。因此，使用定义分值(N)＝3g∈Gs(g，PN(g))，N是最优的网络。
由于我们有|Si|≤c和|Cg|≤m，进行一次算法1的计算时间由一些常数限定。因此，将算法1.1应用于所有的强连接组分能够在O(|G|)内被完成，其完成验证。
我们注意到算法1.1和算法1.1作为验证中定义的子程序的应用可以以获得所有最优网络以及次优网络以等级顺序列表的方式执行。计算的时间将增加，其取决于用于计算的网络的数量，但是仍然是可行的。这对于在分值较小改变的情况下评估估计网络的稳定性颇有价值。
1.2.2没有已有知识的情况下预测大基因网络通过定理1.2，为了获得有效的网络预测，可以为每个基因选择候选亲本，从而导致图中含有所有候选相互作用的强连接组分可以被常数限制。由于基因网络被假定由稀疏地相互连接的高度连接的块(blocks)组成(10，18)，该限制看上去满足许多研究的设定。因此，定理1.2提供了基因网络估计的一般方法的基础。
在本工作中，我们提供了两个被设计用于使用定理1.2的算法的实例。第一个(算法1.2)假定没有在先知识，并通过聚类检测高度连接的块。第二个算法(算法1.3)应用生物学知识以满足定理1.2的限制。
算法1.2步骤1在G中聚类基因，使得没有聚类大于c个基因。
步骤2根据大小递减排序聚类C1，...，Cn。
步骤3对每个i∈{1，...，n}和对每个g∈Ci，从C1U...UCn中选择多至m个候选亲本。
步骤4使用定理1.2计算最优的基因网络模型。
由于所有基因g的候选亲本都选自聚类，g被包含在聚类之内或之前，图中没有被Cg，g∈G诱导的循环能横跨两个聚类。因此，强连接组分不大于c，因为对于所有聚类Ci而言|Ci|≤c，其显示了算法1.2的正确性。
属于同一个聚类的基因相互关联，并因此组成基因网络的紧密连接的部分，而属于不同聚类的基因具有较低的相关性，并因此不太可能在基因网络中连接。我们使用k均值聚类，以Pearson相关性作为距离度量。
用于根据大小递减排序聚类的理由如下所示。在排序开头部分的聚类具有对于候选基因的选择最强的限制，因为存在很少的在前聚类。为最大化从其中选出候选亲本的基因的数量，大聚类将被放在开头部分。
我们注意到步骤1和步骤2强加的子范围的限制仍然允许任意两个基因被连接，仅限制一些基因对边缘的方向。
对步骤3中基因g的候选亲本h的选择存在许多合理的标准。例如基因的相关性或s(g，{h})。在我们的执行中，我们使用了后者的标准。
1.2.3使用在先知识的大基因网络的预测如果存在关于基因网络的在先知识，所述在先知识涉及紧密连接的组分和不同组分中基因相互作用的方向，下列算法能够被应用。
算法1.3步骤1聚合G中的基因于群Ci中，其中|Ci|≤c，并根据生物学知识排列它们C1，...，Cn步骤2对每个i∈{1，...，n}和对每个基因g∈Ci，从C1U...UCi中选出多至m个的候选亲本。
步骤3使用定理2计算最优的基因网络模型。
如对算法1.2，算法1.3的正确性直接取决于步骤2中候选亲本的选择方式。不仅紧密连接组分，而且它们的顺序都已知的例子是在细胞周期的特殊阶段被激活的基因。在关于紧密连接组分和/或它们的顺序的知识是部分知识的情况下，可使用算法1.2和算法1.3的结合。
1.3结果我们已经实施了算法1.2和算法1.3，其使用现有的算法1.1的执行方式(15)，和公众可获得的聚类软件(8)。我们选择BNRC分值(11)作为评分函数s，因为它能够模拟非线性的基因相互作用并能够不需要离散化而处理基因表达数据。
1.3.1将算法1.2用于枯草芽孢杆菌数据我们将算法1.2用于枯草芽孢杆菌数据(9)。我们选择6个转录起始因子(sigma factors)和79个已知被所选转录起始因子调控的操纵子(21)。操纵子的表达率被定义为它们相应基因的平均表达率。因此，我们具有该85个基因相应操纵子的网络中的79个已知调节相互作用，并且如果算法1.2具有预测能力，将识别这些中的重要部分。我们设定参数m(所选候选亲本的最大数量)为15，以及参数c(聚类的最大大小)为25。在步骤1中计算的聚类的实际大小范围在8到24。对于该计算，我们使用具有96个CPUs，每个900MHz的Sun Fire 15K超级计算机运行约一小时。
为了评价所估计的网络的生物学意义，我们计算由算法1.2预测的网络中每个操纵子和每个转录起始因子之间的距离，其中我们定义该距离为连接操纵子和转录起始因子的无引导途径的边缘(edges)的最小数量。当正确的转录起始因子比其它转录起始因子更加接近操纵子时，我们将调节性相互关系评估为被检测到的。这使得能够为我们的估计结果计算p值，因为对于无意义的预测方法(predictor)，检测到调节性相互关系的概率至多为1/6。实际概率低于《方程式GIW20146》，因为我们将中断的连接(breaking ties)计算为未被检测到的。
表1.1将算法1.2用于枯草芽孢杆菌数据.

表1.1显示该评估方案的结果。我们观察到估计的基因网络与两个转录起始因子，20148和20149的生物学知识明显一致。这表明估计的结果含有有价值的生物学信息。我们注意到主要未被检测到的调节性关系是中断连接(breaking ties)，因此操纵子和它们正确的转录起始因子之间的距离仍然是最小的，但是一个或多个其它的转录起始因子同样接近。我们观察到对20150的特别强的显著性与来自[CITE:dehoon03b]对于微分方程网络(differentialequation network)的结果一致。因此，可能数据集对于20151附近的基因网络区域比对其它区域具有更加丰富的信息。
我们注意到很少有工作来解决以原则性方式评估基因网络估计质量的问题。根据我们的知识，仅有一个其它的出版物给出了良好定义的p值以证明估计[CITE:dehoon03b]的显著性。进行基因网络估计的原则性评估时会遇到几个问题。例如关于基因网络的知识不完整，以及不确定在何种试验条件下基因网络的哪个部分工作。另外，如果发现了真实基因网络和估计网络之间结构上的差异，将区分其是实质性的错误差异还是可接受的差异，例如过渡边缘(transitive edges)。上述在估计的网络中分析最短途径的方法是解决若干这些问题的方法之一，但是关于这些问题的工作还需要进一步继续。
1.3.2将算法1.3应用于细胞周期数据算法1.3是基于在先知识的应用以找到生物学上有意义的子范围。在此，我们将算法1.3应用于酿酒酵母的RNA微阵列数据，该数据是从使用同步细胞培养的时程试验(time-course experiments)中获得的(22)。该数据集被期望含有关于在酵母细胞周期期间工作的基因网络的一些信息。
已知三个基因对cln1/cln2，clb5/clb6，和clb1/clb2在细胞周期的各G1/S期、S期、M期被特异性激活(2)。通过基于该事实的聚类分析，几个基因被预测在这些期之一被明显激活(13)。我们选择一组43个基因，将它们分为三群，每群对应细胞周期的一部分。我们将这些群根据细胞周期阶段的时间顺序排序(算法1.3的步骤1)，设置候选亲本的最大数量(参数)为20，并将算法1.3用于估计基因网络。使用1.9GHz的单CPU进行不到一个小时的计算。图1显示估计的结果。
我们观察到总共87个预测的基因相互关系。表1.2显示了从第一基因到第二基因的直接相互关系的数量，所述基因通过基因所属的群(行)和(列)分开。我们看到87个相互关系中52个在群中，这很好的反应了预期的大部分相互作用发生于在细胞周期的相同阶段具有显著活性的基因之间。对于来自不同群的基因对估计了35个相互关系。对于每个时间上相邻的各G1/S和S/G2，S/G2和M，在两个群之间的相互关系的数量为13，但对于从G1/S期到M期的相互关系仅有9个。因此虽然基因网络估计是基于将基因的集合分解为三群，但群之间的相互关系也被很好的说明了。
表1.2基因群的预测的相互关系的数量

我们为每个基因计算相互关系的数量，其是20158的亲本的数量加上该基因子代的数量。在估计的基因网络(图1)中具有最高数量的预测相互关系的基因列于表1.3中，其还给出了酵母基因组数据库中注释的这些基因的分子功能。我们观察到所有具有至少7个预测的相互关系的基因在转录水平上具有调控活性或者具有至今仍未知的功能。在具有6个预测相互作用的三个基因中，一个基因(hcml)具有转录水平的调控功能，但其它两个是细胞周期依赖的信号蛋白。我们推断所有具有更高数量预测相互作用的基因具有已知的活跃调节作用或是具有未知功能。这显示估计的基因网络含有至少一定程度的生物学信息。
重要基因cln2和clb5在表1.3的基因群中具有更低等级的观察结果可以通过它们的信号功能进行解释，该功能可从mRNA测量值中仅仅以间接的方式观察到。有趣的是，swi5被预测为被fkh2所调控，这是已知的调控关系[CITEzhu00]。因此，其它基因在转录水平的活性被认为是更容易观察到的。
表1.3具有至少6个预测相互关系的基因的分子功能


正相反，表1.4显示具有最低数量的预测相互关系的基因。我们发现了几个具有代谢或信号功能的基因，一个具有未知功能的基因，一个mRNA结合基因(sto1)和一个转录因子(swi4)。由于sto1被认为涉及核剪切，较小可能具有调控功能。因此，在11个具有最低数量的预测相互关系的基因中只有一个具有已知的转录水平的调控功能，其在生物学上是看似合理的，因为许多信号活动从mRNA测量中是无法观察到的。我们认为与单纯的聚类结果相比基因网络估计可以提供更多的基因功能的信息。
表1.4具有最多2个预测相互作用的基因的分子功能


我们在此提供用于有效估计大基因网络的一般方法的理论基础。根据该方法的算法由两部分组成。第一部分是用于在搜索范围中识别有生物学意义的子范围的启发式或经验的方法。第二部分是在所选子范围中基因网络的最优估计。
我们已经表明该方法能够估计有意义的基因网络，适当和有效的应用生物学知识，然而计算的时间不会对基因数量产生实际限制。
由于方法不依赖于某些分值或某种基因表达数据，它是普遍可适用的。结合已有知识例如序列信息(23)或结构信息(13)建立新的分值是进一步增加预测能力的有效方法。其它进行启发式子范围选择的方法包括在启发式估计的基因网络周围构建子范围，或者为聚类平均基因表达值并应用算法1为聚类寻找最优的基因网络，其随后能够被用于寻找聚类的顺序(而并非如算法1.2的步骤2中的根据大小的简单排序)。
实施例1的参考文献1.Chickering，D.M.，Learning Bayesian networks is NP-complete，in D.Fisher and H.-J.Lenz，editors，Learning from DataArtificial Intelligence andStatistics V，Springer-Verlag，1996.
2.Cho，R.J.，Campbell，M.J.，Winzeler，E.A.，Steinmetz，L.，Conway，A.，Wodicka，L.，Wolfsberg，T.G.，Gabrielian，A.E.，Landsman，D.，Lockhart，D.J.，Davis，R.W.，A genome-wide transcriptional analysis of the mitotic cell cycle，Molecular Cell，265---73，1998.
3.Cooper，G.F.，Herskovits，E.，A Bayesian method for the induction ofprobabilistic networks from data.Machine Learning，9309--347，1992.
4.Friedman，N.，Goldszmidt，M.，Learning Bayesian networks with localstructure，Jordan，M.I.(ed.)，Kluwer Academic Publishers，421--459，1998.
5.Friedman，N.，Linial，M.，Nachman，I.，Pe′er，D.，Using Bayesiannetworks to analyze expression data，Journal of Computattional Biology，7601--620，2000.
6.Hartemink，A.J.，Gifford，D.K.，Jaakkola，T.S.，Young，R.A.，Usinggraphical models and genomic expression data to statistically validate models ofgenetic regulatory networks，Pacific Symposium on Biocomputing，WorldScientific，422--433，2001.
7.Hartemink，A.J.，Gifford，D.K.，Jaakkola，T.S.，Young，R.A.，Combininglocation and expression data for principled discovery of genetic regulatorynetwork models，Pacific Symposium on Biocomputing，World Scientific，7437--449，2002.
8.de Hoon，M.J.L.，Imoto，S.，Nolan，J.，and Miyano，S.，Open sourceclustering software，Bioinformatics，2003，in press.
http://bonsai.ims.u-tokyo.acjp/mdehoon/software/cluster/9.de Hoon，M.J.L.，Ott，S.，Imoto，S.，Miyano，S.，Validation of noisydynamical system models of gene regulation inferred from time-course geneexpression data at arbitrary time intervals，European Conference onComputational Biology，2003.
10.Hartwell，L.H.，Hopfield，J.J.，Leibler，S.，Murray，A.W.，Frommolecular to modular cell biology，Nature，402C47--C52，1999.
11.Imoto，S.，Goto，T.，Miyano，S.，Estimation of genetic networks andfunctional structures between genes by using Bayesian networks andnonparametric regression，Pacific Symposium on Biocomputing，World Scientific，175--186，2002.
12.Imoto，S.，Kim，S.，Goto，T.，Aburatani，S.，Tashiro，K.，Kuhara，S.，Miyano，S.，Bayesian network and nonparametric heteroscedastic regression fornonlinear modeling of genetic network，Journal of Bioinformatics andComputational Biology，1231--252，2003.
13.Nariai，N.，Kim，S.，Imoto，S.，Miyano，S.，Using protein-proteininteractions for refining gene networks estimated from microarray data byBayesian networks，Pacific Symposium on Biocomputing，World Scientific，inpress，2004.
14.Ong，I.M.，Glasner，J.D.，Page，D.，Modelling regulatory pathways inE.coli from time series expression profiles，Bioinformatics，18241--248，2002.
15.Example 2.
16.Pe′er，D.，Regev，A.，Elidan，G.，Friedman，N.，Inferring subnetworksfrom perturbed expression profiles，Bioinformatics，17215--224，2001.
17.Robinson，R.W.，Counting labeled acyclic digraphs，New Directionsin the Theory of Graphs，239--273，1973.
18.http://www.yeastgenome.org/19.Shen-Orr，S.S.，Milo，R.，Mangan，S.，Alon，U.，Network motifs in thetranscriptional regulation network of Escherichia coli，Nature Genetics，3164--68，2002.
20.Smith，V.A.，Jarvis，E.D.，Hartemink，A.J.，Evaluating functionalnetwork inference using simulations of complex biological systems，Bioinformatics，18216--224，2002.
21.van Someren，E.P.，Wessels，L.F.A.，Backer，E.，Reinders，M.J.T.，Genetic network modeling，Pharmacogenomics，3(4)507--525，2002.
22.Sonenshein，A.L.，Hoch，J.A.，Losick，R.，Bacillus subtilis and itsclosest relativesfrom genes to cells，ASM Press，Washington，D.C.，2001.
23.Spellman，P.，Sherlock，G.，Zhang，M.，Iyer，V.，Anders，K.，Eisen，M.，Brown，P.，Botstein，D.，Futcher，B.，Comprehensive identification of cellcycle-regulated genes of the yeast Saccharomyces cerevisiae by microarrayhybridization，Molecular Biology of the Cell，93273--3297，1998.
24.Tamada，Y.，Kim，S.，Bannai，H.，Imoto，S.，Tashiro，K.，Kuhara，S.，Miyano，S.，Estimating gene networks from gene expression data by combiningBayesian network model with promoter element detection，Bioinformatics，19227--236，2003.
25.Zhu，G.，Spellman，P.T.，Volpe，T.，Brown，P.O.，Botstein，D.，Davis，T.N.，Futcher，B.，Two yeast forkhead genes regulate the cell cycle andpseudohyphal growth，Nature，40690--94，2000.
实施例2为小基因网络寻找最优模型介绍从基因表达测定值推导基因网络是系统生物学的主要挑战。如果基因网络能够正确推导，将能够导致更好的理解细胞过程，并由此应用于药物开发、疾病研究和其它领域。
贝叶斯网络是模拟基因网络的最常用的方法(见参考文献3，4，7，9，11，13，14，17)。在贝叶斯网络中，基因网络的行为被模拟为所有基因的联合概率分布。这导致基因相互关系的非常普遍的模拟。联合概率分布可以分解为条件概率P(Xg|X1，...，Xn)的产物，其表示一些基因g1，...，gn对基因g的调控。该分解可以表示为定向非循环图。贝叶斯网络模型被显示可以寻找生物学上看似可信的基因网络(见参考文献4，9)。
然而，学习贝叶斯网络的困难在于它的大搜索范围。n个基因的基因网络的搜索范围是含有n个顶点的定向非循环图的范围。带有n个顶点(cn)的定向非循环图的数量的递归公式和渐近式来自Robinson(见参考文献15)。我们在此指定渐进式为
cn=n!·2n2·(n-1)r·zn;r~0.57436;z~1.4881]]>例如，存在大约2.34·1072种可能的含有20个基因的网络，和对于含有30个基因的基因网络有大约2.71·10158种可能的解决方案。甚至对于9个基因的基因网络(搜索范围大小大约1.21·1015)，强力方法也需要在超级计算机上运行几年的计算时间。另外，已知寻找最优网络的问题是NP-hard(见参考文献1)，甚至对于离散分值BDe(见参考文献2，3)和MDL(见参考文献3)也如此。因此，研究者至今仍使用启发式方法，例如模拟退火(见参考文献8)或贪心算法(见参考文献9)来估计贝叶斯网络(见参考文献18)。
然而，由于启发式的准确性是不确定的，难以根据启发式估计的网络得出结论。为了克服该问题，我们分析了超指数搜索范围的结构，并开发了在指数时间内于超指数搜索范围中寻找最优解决方案的算法。该方法对于20个或更多基因的基因网络是切实可行的，其取决于所使用的具体概率分布。此外，加入生物学合理的假设，可以为多至40个基因的基因网络推导最优网络。
克服启发式的不确定致使有可能比较统计学模型在推断生物学正确的基因网络方面的能力。而且，该方法是改进基因的已知功能群的基因网络的有效工具。
方法在第2.2部分中被提供。在第2.3部分中我们提供了应用该方法的结果，其显示它能够从生物学上准确地估计基因网络。
2.2方法2.2.1预备在该部分中，我们假定我们被提供了一组基因G和如由几个群所使用的网络评估函数(见参考文献4，9，17)，即函数sG×2G→R，其赋予基因g∈G和一组亲本基因AG分值。对于网络N，N的分值被定义为分值(N)＝3g∈Gs(g，PN(g))，其中PN(g)表示N中g的亲本的集合。
实施例
1.BDe分值(见参考文献2，3)给定数据，则分值与网络的后验概率成比例。当使用Bde分值时，需要离散微阵列数据。
2.MDL分值(见参考文献3)MDL分值使用最小描述长度法则(minimaldescription length principle)，并且也使用离散数据。
3.BNRC分值(见参考文献9)BNRC分值使用无参数回归以获取非线性基因相互关系。由于数据不需要离散，因此没有信息会丢失。推导网络的任务是为每个基因寻找一组亲本基因，从而使得所获得的网络是非循环的，并且网络的分值最小。
下面，我们介绍一些描述算法所需的符号。
定义2.1F我们为所有的g∈G和AG定义于G×2G→R为F(g，A)＝defminBAs(g，B)。
通过定义，F(g，A)的意义为当亲本必须选自子集A时基因g的亲本的最优选择。对于每个非循环图，存在顶点的排序，从而使得所有的边缘沿着排序的方向定位。相反的，当提供G的固定顺序，我们能够考虑遵循所给顺序的所有图的集合，就如我们在接下来的定义中所做的。
集合AG的排序可以被描述为数列π{1，...，|A|}→A。让我们使用∏A表示A的所有数列的集合。
定义2.2π-线性假设AG和 假设NA×A是网络。我们假定当且仅当所有(g，h)∈Nπ-1(g)＜π-1(h)成立时N是π-线性。
现在我们使用上面的定义并定义函数QA，其将使我们得以为所提供的π计算最佳π-线性网络的分值，就如我们下面所示的。
定义2.3QA假设AG。对于所有的 我们定义QA∏A→R为QA(π)=defΣg∈AF(g,{h∈A|π-1(h)<π-1(g)}).---(2)]]>
如果我们能够使用函数F和Q为特定的数列π计算最好的π-线性网络，那么我们为寻找最优网络所需要做的是寻找最优的数列π，其产生全局最小值。形式上，我们为该步骤定义函数M。
定义2.4M对于所有的AG，我们定义M2G→AG∏A为M(A)=defminπ∈ΠAQA(π)---(3)]]>算法2.1使用上述符号，算法2.1可以被如下定义。
步骤1对所有g∈G计算F(g，φ)＝s(g，φ)。
步骤2对所有AG，A≠φ和所有g∈G计算F(g，A)为min{s(g,A),mina~∈AF(g,A-{a})}.]]>步骤3设定M(φ)＝φ。
步骤4对所有AG，A≠φ，进行下列两个步骤步骤4a计算g*＝arg mingA(F(g，A-{g}+QA-{g}(M(A-{g})))。
步骤4b对所有1≤i＜|A|，设定M(A)(i)＝M(A-{g*})(i)，M(A)(|A|)＝g*。
步骤5返回QG(M(G))。
在步骤2和步骤4中提供的递归公式中，我们希望为基数m＝|A|的子集AG计算函数F各(resp.)M，并且需要为基数m-1的子集计算函数F各(resp.)M的函数值。因此，我们能够应用步骤2和步骤4的动态程序为增加基数的子集A计算函数F各(resp.)M。在步骤4的递归公式中，首先在步骤4a中计算数列M(A)的最后的元素g，然后在步骤4b中设定M(A)。
正确性和时间复杂度算法2.1的步骤2中的递归公式的正确性直接取决于F的定义。因此，在步骤1和步骤2运行后，所有基因g和所有子集A□G的函数F的值被存储到存储器中。在进行步骤3和步骤4之前，我们规定关于函数QA的意义的引理。
引理2.1假设AG和π∏A。假设N*A×A是具有最小分值的π-线性网络。那么，QA(π)＝score(N*)成立。
证明.在π-线性图中，基因g只能具有亲本h，其在由π编码的顺序中位于上游，也就是π-1(h)＜π-1(g)。因此，当为g选择亲本，我们限定为B＝{h∈A|π-1(h)＜π-1(g)}，且F(g，B)是在这种情况下的最优选择。由于在π-线性图中，所有的边缘遵循由π编码的顺序，我们能够通过该方法独立地为所有基因选择亲本，其证明了引理2.1。
使用引理2.1，我们证明函数M能够通过步骤4中提供的公式被计算。
引理2.2假设AG。假设g*＝arg ming∈A(F(g，A＝{g})+QA-{g}(M(A-{g}))).
通过π(i)＝M(A-{g*})(i)定义π∈∏A，且π(|A|)＝g。
那么，π＝M(A)。
证明.假设π-‘∈∏A。通过M的定义，我们不得不显示QA(π)QA(π’)。假设N*是最优的π’-线性网络，M*是最优的π-线性网络。那么，通过引理2.1，QA(π)QA(π’)等价于分值N*分值M*。让我们表示π’的最后项为h＝π′(|A|)。我们注意到对于所有BG，QB(M(B))是通过上述定义基于B的全局最优网络的分值。因此，我们有分值(M*)=s(h,PM*(h))+Σg∈A-{h}s(g,PM*(g))]]>≥s(h,PM*(h))+QA-{h}(M(A-{h}))]]>≥F(h，A-{h})+QA-{h}(M(A-{h}))≥minh∈AF(h，A-{h})+QA-{h}(M(A-{h}))=F(g*,A-{g*})+QA-{g*}(M(A-{g*}))]]>分值(N*)，其证明了引理2.2。
由于Q能够用F直接计算，算法能够在步骤5中计算QG(M(G))。最后，QG(M(G))通过定义是最优贝叶斯网络的分值，其显示了正确性。
如果每个基因和每个子集AG的最佳亲本的信息与F(g，A)共同存储，最优的网络可以在QG(M(G))的计算期间被构建。
定理2.1没有约束的最优网络算法使用O(nA2n)动态规划步骤寻找最优网络。
证明.在步骤1和步骤2中的动态规划需要O(nA2n)(n＝|G|)步骤，且在每步计算一个分值。在步骤3和步骤4的动态规划中需要O(2n)步，其中每步包括寻找一些先前存储的分值。注意到步骤4a中函数QA不需要被实际计算，因为在先前的步骤中QA-{g}能够与M(A-{g})共同存储。因此，总的时间复杂度是O(nA2n)。
在生物学现实中，虽然调控基因的子代的数量可以非常高，亲本的数量可以被假定是被限制的。当我们限制亲本的数量，分值计算的数量充分降低，其使得可以进行更大网络的计算。
我们规定下列推论，其在实践上非常有意义(见第2.3部分)。
推论2.1假设∈IN是常数。其中没有基因具有超过m个亲代的最优网络能够在O(n·2n)步动态规划步骤中被找到。[引证]如果我们不希望通过常数限制亲代的数量，但是改为能够为每个基因选择确定数量的候选亲本，复杂度改变如下。
推论2.2假设∈IN是常数。对于每个g∈G，假设CgG是Cg≤m的集合。其中每个基因g只在Cg中具有亲本的最优网络可以在O(2n)步动态规划步骤中找到。
证明.由于每个基因的亲本选自恒常大小的集合，在步骤1和步骤2中的动态规划的复杂度成为常数。因此，总的复杂度成为O(2n)。
我们注意到，动态规划在我们算法中的两个应用可以被执行为动态规划的单个应用，因为当我们为大小为m的集合计算函数M时，我们仅仅需要大小为m-1的集合的函数F的函数值。因此，只有函数F和M对大小为m-1和m的集合的函数值需要在同一时间被存储在存储器中。这在实践上对于减少所需的存储器数量很有意义。
我们也注意到算法能够被修改以便同样计算次优解决方案。计算第二最佳或第三最佳网络可能对于在分值少量改变的情况下评估推断的网络的稳定性是有价值的。
结果上述算法2.1作为C++程序被执行。作为评估函数，现有的BNRC分值、BDe分值和MDL分值的执行被使用。执行所有的三种方法(定理2.1，推论2.1和2.2)。
我们将程序应用于173个微阵列的数据集，测定酿酒酵母对各种应激状态(stress conditions)的反应(见参考文献5)。
对热休克数据的应用从数据集中我们选择从25℃到37℃热休克实验的15个微阵列和从不同温度到37℃的热休克实验的5个微阵列。然后我们选择涉及或推定涉及热休克反应的9个基因的集合。图2显示关于BNRC分值的最优网络。
我们观察到转录因子MCM1被估计调节三个其他基因，但它不被该组中一个基因调节，这是看似合理的。在我们的基因集合中的第二个转录因子HSF1被估计调节三个其他热休克基因。它也被估计被HSP70-蛋白(SSA1)调节，这在以前被报导过(见参考文献16)。这些基因中的另一个侣伴蛋白SSA3也好象在热休克反应中具有活跃作用，并与SSA1和HSP104相互作用，其与Glover和Lindquist的报导一致(见参考文献98)。
表2.1记录了该分析的结果。总体而言，该结果是生物学看似合理的，并且指出侣伴蛋白SSA1和SSA3在热休克反应期间的活跃作用。我们认为最优推导的基因网络是有意义的，并且对于基因调节的阐明是有用的。
表2.1涉及热休克反应的基因

计算的可能性和局限性虽然甚至像第2.3.1部分中所推导的网络那样的小范围网络也不能用强力方法(等式1)推导，它们能够通过我们的程序使用单个的奔腾以CPU1.9GHz运行约10分钟被最优推导。为了评估本方法的实际可能性，我们从数据集中选择20个已知在基因调节中具有活跃作用的基因(见参考文献12)，并使用全部173个微阵列估计具有最优BNRC分值的网络。通过使用具有90个CPUs，每个900MHz的Sun Fire 15K超级计算机，计算在约50个小时内完成。作为该计算实验的结果，我们认为我们的方法对于20个基因的基因网络是可行的，即使没有进行约束并且使用像BNRC分值这样的复杂评分方案。对于离散分值BDe和MDL(其能够计算的更加快)，甚至超过20个基因的网络也能够没有约束的被最优推导。
当亲本的数量被限定为约6(推论2.1)或，可选的，预选约20个候选亲本的集合(推论2.2)，使用BNRC分值甚至超过30个基因的基因网络也能够被最优的推导。然而，像现在这样的方法将不能允许估计超过约40个基因的网络。
尽管推论2.2中所给出的方法的理论时间复杂度低于推论2.1中所给出的方法的时间复杂度，我们指出后者可能实际上更加重要。第一，通过常数限定亲本的数量能够被容易的进行并且是生物学合理的，而为每个基因选择一组候选亲本需要基因选择的方法，其能够对计算结果有潜在的偏离。第二，必须考虑到在函数F的计算中的每个动态规划步骤都需要计算一个分值，然而函数M的一个动态规划步骤只需要寻找一些先前的结果。当亲本的数量被如推论2.1中那样限定，所需的分值计算的数量成为多项式，其使得本方法在实际应用中更快，虽然推论2.2中的方法在理论上更好。
我们已经提供了可以最优的推导20-40个基因的基因网络的方法，其取决于所使用的概率分布以及是否进行了另外的假设。这使得可以比较不同的评分方案，以评估特定评分方案的最佳参数，以及评估特定微阵列数据的有效性，因为获得了最优解决方案。同样，在研究者关注于一群确定的基因并且希望最大程度地利用这些基因的基因表达测定值的设定情况(setting)下该方法特别有用。
与启发式方法相反，如果结果不满意或与生物学知识相矛盾，能够推断该统计学模型是不正确的或数据是不充分的。甚至对于20个基因的网络，了解来自特定的大搜索范围的最佳网络仍是很大量的信息。
我们注意到方法不依赖于某种评分方案或某种类型的基因表达测定值。它能够应用于任何设定中，其中提供了如第2.2部分所定义的分值。例如，当序列信息(见参考文献19)、蛋白相互作用数据(见参考文献10)或其他知识被结合于评估函数中，该方法也能够被应用。
为了寻找具有超过40个基因的基因网络，两个方向可以被使用。第一，如果子集集合的一部分(其中算法进行实时搜索)能够被剪除，对可行性的限制将被增加。第二，基因网络的划分(见参考文献18)可以被用于分解更大的网络为更小的部分，并且最优的推导每部分网络。
参考文献1.D.M.Chickering.Learning Bayesian networks is NP-complete.In D.Fisher and H.-J.Lenz，editors，Learning from DataArtificial Intelligence andStatistics V，Springer-Verlag，1996.
2.G.F.Cooper，E.Herskovits.A Bayesian method for the induction ofprobabilistic networks from data.Machine Learning，9309--347，1992.
3.N.Friedman，M.Goldszmidt.Learning Bayesian networks with localstructure.Jordan，M.I.(ed.)，Kluwer Academic Publishers，pp.421--459，1998.
4.N.Friedman，M.Linial，I.Nachman，D.Pe′er.Using Bayesian networks toanalyze expression data.Journal of Computational Biology，7601--620，2000.
5.A.P.Gasch，et.al.Genomic expression programs in the response of yeastcells to environmental changes.Molecular Biology of the Cell，114241--4257，2000.
6.J.R.Glover，S.Lindquist.Hsp104，Hsp70，and Hsp40a novel chaperonesystem that rescues previously aggregated proteins.Cell，9473--82，1998.
7.A.J.Hartemink，D.K.Gifford，T.S.Jaakkola，R.A.Young.Using graphicalmodels and genomic expression data to statistically validate models of geneticregulatory networks.Pacific Symposium on Biocomputing，6422--433，2001.
8.A.J.Hartemink，D.K.Gifford，T.S.Jaakkola，R.A.Young.Combininglocation and expression data for principled discovery of genetic regulatorynetwork models.Pacific Symposium on Biocomputing，7437--449，2002.
9.S.Imoto，T.Goto，S.Miyano.Estimation of genetic networks andfunctional structures between genes by using Bayesian networks andnonparametric regression.Pacific Symposium on Biocomputing，7175--186，2002.
10.T.Ito，T.Chiba，R.Ozawa，M.Yoshida，M.Hattori，Y.Sakaki.Acomprehensive two-hybrid analysis to explore the yeast protein interactome.Proceedings of the National Academy of Sciences USA，974569--4574，2001.
11.S.Imoto，S.Kim，T.Goto，S.Aburatani，K.Tashiro，S.Kuhara，S.Miyano.Bayesian network and nonparametric heteroscedastic regression fornonlinear modeling of genetic network.Journal of Bioinformatics andComputational Biology，in press，2003；U.S.Patent Application No10/259,723.
12.T.I.Lee，N.J.Rinaldi，F.Robert，et.al.Transcriptional regulatorynetworks in Saccharomyces cerevisiae.Science，298799--804，2002.
13.I.M.Ong，J.D.Glasner，D.Page.Modelling regulatory pathways in E.coli from time series expression profiles.Bioinformatics，18241--248，2002.
14.D.Pe′er，A.Regev，G.Elidan，N.Friedman.Inferring subnetworks fromperturbed expression profiles.Bioinformatics，17215--224，2001.
15.R.W.Robinson.Counting labeled acyclic digraphs.New Directionsin the Theory of Graphs，pp.239--273，1973.
16.Y.Shi，D.D.Mosser，R.I.Morimoto.Molecular chaperones asHSF1-specific transcriptional repressors.Genes& Development，12654--666，1998.
17.V.A.Smith，E.D.Jarvis，A.J.Hartemink.Evaluating functional networkinference using simulations of complex biological systems.Bioinformatics，18216--224，2002.
18.E.P.van Someren，L.F.A.Wessels，E.Backer，M.J.T.Reinders.Geneticnetwork modeling.Pharmacogenomics，3(4)507--525，2002.
19.Y.Tamada，S.Kim，H.Bannai，S.Imoto，K.Tashiro，S.Kuhara，S.Miyano.Estimating gene networks from gene expression data by combiningBayesian network model with promoter element detection.Bioinformatics，inpress，2003.
实施例3比最优更好的基因网络3.1介绍从表达水平测量估计基因网络是近年来生物信息学研究的焦点之一(1，5，8，13，31)。从表达数据了解基因网络的结构可以深化我们对于细胞对环境的反应的了解和我们关于基因网络的构建原则的知识(10，29)，其在几个学科中都可能应用。
模拟基因网络的常用方法是贝叶斯网络(3，4，9，13，15，26，30)，其将基因的表达水平模拟为随机变量，将基因网络模拟为联合概率分布。这些分布通过使用我们称为网络的定向非循环图进行分解。给定数据，则使用基于网络可能性的评分函数对网络进行评分。
最近的结果显示，关于约30个或更少(25)基因的基因网络的一些数据，能够了解最优基因网络模型。该结果支持任意评分函数sG×2G，其设定基因g和g的一组亲本的分值。然而，虽然最优的基因网络模型是最可能的模型，仍存在对于一些数据具有大约同样的可能性的非常不同的模型，尤其是由于基因表达测量值通常仅提供基因网络的部分信息。另外，即使对于仅10个基因的基因网络，存在大约4.17*1018种可能的网络模型(17)。因此，即使最优的基因网络模型通常也不匹配目标基因网络。
我们解决该问题的努力是三重的。第一，我们提供了用于以可能性顺序列举最优和次优网络的理论基础和算法。第二，我们严格比较了最优基因网络模型和有关现有基因网络的已知知识。第三，我们提供了获取最优基因网络模型的可靠部分的方法并应用该方法于枯草芽孢杆菌和大肠杆菌的已知数据从而比较相关物种的基因网络模型。
我们的结果显示我们的方法识别的部分网络明显比最优网络模型本身更好的符合知识。由于我们将我们的方法以最好的n基因网络模型为基础，因而从这些估计中获得的结论比来自基于启发式方法的方法的更加可靠。
在第3.2部分给出我们的理论结果后，我们在第3.3部分中陈述和分析了估计和知识的比较结果。然后，我们在第3.4部分给出了我们的基因网络估计方法，在第3.5节评估和应用了该方法。
3.2基因网络估计和列举在该实施例中，我们假定我们给出了一组基因G和评估函数sG×2G→IR，其设定基因g和g的一组亲本的分值。对于网络N，N的分值被定义为分值(N)＝∑g0Gs(g，PN(g))，其中PN(g)表示N中g的亲本的集合。为了寻找解释给定数据的最好的联合概率分布，我们需要寻找最小化分值(N)的定向非循环图N。
由于该最优化的问题是NP-hard(2)，并且搜索范围具有超指数大小(27)，研究者经常使用启发式方法，例如贪心算法、模拟退火或遗传算法(5，13，31)。然而，最近得到了如下结果，其所得到的算法在没有约束的情况下能够估计约20个基因的基因网络，当使用现实基因网络中的基因具有限定数量的亲本的事实时，能够估计具有30个或稍多基因的网络。
定理3.1(25)最优网络能够通过使用O(n·2n)步动态规划步骤被找到。
对于更大的基因网络，经验或启发式方法能够被用于选择搜索范围的子范围。如果其如(24)中所述的进行，(25)的算法能够被用于在所选子范围中寻找最优解决方案，其产生了确定性和启发式技术的结合方法。
然而，不同网络可以产生具有相同或近似相同概率的同样的表达数据。因此，不可以假设所提供的数据含有关于基因网络的完整信息，但是不得不记住能够从给定的数据集获得的信息可能的确是部分数据。因而，即使最优网络被找到，它通常含有错误的边缘。为了克服该问题，我们已经扩展了(25)的算法并且证明了下列结果。
定理3.2假设m0IN。最佳m网络可以使用O(n·2n)动态规划步骤被找到。
3.2.1证明-列举最优基因网络3.2.1.1算法3.1我们在此显示如何扩展来自(25)的算法以使得能够列举最优和次优网络。我们首先定义了一些函数，然后显示这些函数如何能够计算相当大的基因网络。如在第3.2部分中，我们假设我们被提供了一组基因G和评分函数sG×2G！IR。对于网络N，N的分值被定义为分值(N)＝3g0Gs(g，PNg))，其中PN(g)其中PN(g)表示N中g的亲代的集合。在该部分中，我们假定具有相等分值的网络被以一些方式排序，由此产生概念“第n个最佳网络”。
定义3.1Fm假设m∈IN。我们归纳定义FmGx2GxI→→2G。首先，对于所有g∈G和AG，我们定义Fm(g，A，n)＝defarg min s(g，B).
BA (2)然后，表示所有先前解决方案{Fm(g，A，p)*p＜n}的集合为J(n)，Fm(g，A，n)＝defarg min s(g，B)BA(3)B∈J(n)对于所有1＜n≤m。
定义3.2Sm假设m0IN。我们定义SmGx2GxIN≤m→IR为
Sm(g，A，n)＝defs(Fm(g，A，n)) (4)对于所有g∈G，AG，和n∈IN≤m。
∈∈通过定义，Fm(g，A，n)是当亲本必须选自A中时基因g的亲本的第n个最佳选择，Sm(g，A，n)是该选择的分值。当m根据上下文是清楚的时候，我们分别使用F和S而不是Fm和Sm。注意到Fm和Sm是部分定义的函数，因为m可以大于A的子集的数量。
集合AfG的排序可以被描述为数列B{1，...，|A|→A。让我们使用∏A表示A的所有数列的集合。我们需要下列概念，其表示其中所有的边缘都遵循数列π所提供的指示的网络。
定义3.3π-线性假设AG和π∈∏A。假设NAxA是网络。我们认为当且仅当对于所有(g，h)∈Nπ-1(g)＜π-1(h)成立，N是π-线性的。
我们算法的基本策略是对于所有π∈∏A将集合AG的所有定向非循环图的范围分割成π-线性网络的子集合，分解寻找最优和次优网络的问题为下列两个问题1.寻找含有被搜索的(次)最优网络的搜索范围的子集。(引理2)2.在所选子范围中寻找(次)最优网络。
(引理1和引理2)我们使用π表示搜索范围的子范围。为了为给定的数列π寻找最优和次优网络，我们需要下列函数QA。对于给定的基因g，我们表示π中所有先于g的基因的集合为V(π，g)＝def{h|π-1(h)＜π-1(g)}。
定义3.4QA假设AG。对于所有π∈∏A，我们定义QA∏A×IN|A|→2A×A为QA(π，v)＝def{(h，g)|h∈F(g，V(π，g)，vπ-1，(g)，)}(5)。
在定义4中，我们已经使用向量v∈IN|A|测定特殊基因亲本选择的等级。在引理1和引理2中，显示QA(π，v)是最优或次优π-线性网络，它的等级依赖于v。然后，我们分别定义指定子范围的两个函数Mm和Dm(在子范围中可以找到(次)最优网络)，和从子范围中网络的选择。
定义3.5Mm，Dm假设m∈IN。我们在它们的第二参数之上归纳定义Mm2G×IN≤m→UAGΠA和Dm2G×IN≤m→U|G|i=0IN.]]>假设AG。首先，我们定义Dm(A，1)＝def(1，...，1)∈IN|A|(6)和Mm(A，1)＝defarg min score(QA(π，Dm(A，1))) (7)π∈∏A假设π∈IN≤m其中n＞1和假设N是关于A的网络，其具有不在{QA(Mm(A，p)，DM(A，p))|p＜n}中的网络之中的最优分值。假设π*∈∏A是数列，其使得N是π*-线性。我们定义Mm(A，n)＝defπ*(8)假设v*∈IN|A|，其使得对于每个g∈A，g的亲本的集合PN(g)等于Fm(g,V,(π*,g),v*π*-1(g))5.---(9)]]>我们定义Dm(A，n)＝defv*(10)如同Fm和Sm，Mm和Dm是部分函数。在表3.4中，我们总结了上述概念。
使用这些概念，给出m0IN，我们的算法能够如下定义，
表3.1用于定义算法的函数含义直接根据定义和引理3.1.

使用这些概念，算法3.1可以被如下定义。

算法3.1在步骤1和步骤2中为所有g∈G，AG和n≤m计算函数Fm和Sm。为了在步骤2a选择B2，仅需要更低基数(cardinality)或更低n的集合A的Fm的函数值。因此Fm和Sm可以使用动态规划(programming)进行计算。
在步骤3和步骤4中，函数Mm和Dm能够类似的使用动态规划进行计算，因为为了在步骤4a中选择三元组(triple)仅需要针对更低基数的集合A或更低基数的n的Mm和Dm的函数值以计算Mm(A，n)和Dm(A，n)。三元组(g，p，q)的含义如下。g是被搜索的数列中最后元素的候选。当g成为最后元素，余下的数列能够从多至m个的先前计算的A-{g}的数列中选择。余下的所选数列通过p指定。然后，为形成子范围中由获得的数列定义的网络，使用g的亲本的第q个最佳选择，但对于其他基因如Dm(A-{g}，p)所示地选择亲本。
由于所有四个函数Fm，Sm，Mm，和Dm都是部分定义的，对于低基数的集合A不是所有的n函数值都能够被计算。为了简单，我们没有在算法的公式中明确的提到它。
如上所规定的算法3.1计算确定数量(m)的解决方案，不论集合A的大小如何。在实际应用中，计算时间可以通过为具有大量含有l个元素的G的子集的层l计算更少量(次)最优解决方案进行优化。然后，当子集的数量由于A的高基数降低，m能够被增加。当只有m个最优解决方案在更低层已知时，不保证超过m个的(次)最优解决方案能够被产生，但是其在测试计算中起作用。
由于必须在计算过程中存储的网络的数量可能变得非常大，重要的是有效的在存储器中存储网络，以及允许快速地解码存储的网络。在我们的执行中，基于上述公式，我们存储数列π并且使用对π-线性网络中边缘的限制，其产生存储和时间高效的网络编码。
为了满足步骤4a的要求，即所选网络不同于先前所选的网络，网络必须在动态规划步骤期间被比较。然而，由于具有不同分值的网络是不同的，只有具有相等分值的网络需要被比较。因此，网络比较的数量在实际应用中可以被最小化。
当Fm，Mm，或Dm的一层被计算时，对层中的集合该计算能够以任意顺序进行，且计算仅依赖于更低层中的值。因此，上述算法是很好的可平行化的。
3.2.3正确性和复杂度我们通过N*A，n表示集合AG的第n个最佳网络。我们首先规定了两个重新用公式表示的来自(25)的引理。
引理3.1假设AG和π∈∏A。假设N*A×A是具有最小分值的π-线性网络。那么，QA(π，(1，...，1))＝分值(N*)成立。
引理3.2假设AG和m∈IN。假设g*＝arg ming∈A(Sm(g，A＝{g}，1)+N*A-{g}，n)。通过π(i)＝M(A-{g*}，1)(i)定义π∈∏A，且π|A|＝g*。那么π＝Mm(A，1).。
使用这些，我们能够证明上述定义的算法的正确性。
定理3.2假设m∈IN通过使用O(n·2n)动态规划步骤能够找到最佳m网络。
证明.通过定义，算法的输出，QG(Mm(G，i)，Dm(G，i))，i≤m，是G的最佳m网络。我们仅需要证明算法中给出的递归公式是正确的。通过Fm和Sm的定义，步骤1中给出的公式是正确的。当我们在步骤2中选择集合A？G的子集，我们基本上有两个选择整个集合A或确实的子集(ture subset)。在以前的例子中，我们能够直接计算选择的分值，在后来，我们能够使用先前计算的Fm和Sm的值，其提供了步骤2的正确性。
完成步骤2后，我们计算了Fm和Sm的所有值。使用这些值，函数Q能够直接被计算。因此，我们仅仅需要计算函数Mm和Dm以在步骤5中产生输出。通过定义，步骤3中的等式又是正确的。我们观察到通过引理3.1与定义3.5结合，下列等式成立N*a，n＝QA(Mm(A，n)，Dm(A，n))
从这个等式和引理3.2我们看到对于n＝1步骤4中的递归式是正确的。对于n＞1，我们以同样的方式计算次优数列Mm(A，n)和亲本的次优选择Dm(A，n)，其限制为先前没有被选择的网络。
一个动态规划步骤所需的时间依赖于m，但是能够被认作如在这个工作的计算中所选的m的常数。此外，使用在附录中所描述的规划技术，在实际计算中它在多数动态规划步骤中对于计算明显少于m的网络是足够的。
使用这个算法，我们能够以他们可能性的顺序列举最优和次优网络。当我们能够寻找到该网络的共同组分，共同组分的可能性是含有它的网络的可能性的总和。因此，列举的网络的共同组分，在本工作中被称为基因网络基序。因此，我们的基因网络基序的定义首先不同于例如在(21)中所使用的概念，但可能是紧密相关的。因此，网络基序分析被期望比单独的网络估计更加可靠。
3.3.评估基因网络估计的可靠性为了验证列举的最优和次优基因网络模型含有关于真实基因网络的有价值的信息，我们执行算法并将它应用于RNA微阵列数据。考虑到细菌中转录和翻译的关系，我们期望单独根据RNA数据的基因网络估计对细菌比对真核细胞更加合适，在真核细胞中转录和翻译发生于不同的地方和不同的时间。因此，我们选择枯草芽孢杆菌和大肠杆菌作为目标，它们的微阵列数据和关于基因网络的知识是已知的。
3.3.1数据和软件对于大肠杆菌，我们从Gene Expression Omnibus选择数据集GDS95--GDS100。通过色氨酸代谢干扰、UV照射和新霉素处理诱导基因表达水平的改变(17，18，19)。然后我们从RegulonDB(28)中获得有关大肠杆菌中已知转录调控的数据以与估计结果进行比较。
对于枯草芽孢杆菌，我们使用数个数据集，其包括来自不同处理下的时程实验的70个微阵列，和来自基因破坏实验(gene disruptant experiments)的99个微阵列。数据还没有公开，虽然已经证明了能够使用这些数据进行生物学有意义的估计(12)。然后我们从DBTBS(16，20)中获得已知转录调控的数据集。
我们的执行所使用的评分函数(sG×2G→IR)中(其包括MDL分值(4)、BDe分值(3，4)和BNRC分值(13))因为下列原因在本工作中我们选择BNRC分值用于计算。第一，BNRC使用不要离散化的数据，其避免了附加的参数和信息的丢失。第二，使用B-样条(B-splines)模拟基因相互关系，其允许不限于线性关系的全面模拟。第三，在根据枯草芽孢杆菌数据集的使用所有三种评分函数的预计算实验中，使用BNRC分值的估计明显更好的符合经典知识(32)。
3.3.2目标网络的选择从大肠杆菌已知调控关系的数据集中，我们选择所有提供了实验证据的关系。这产生了899个已知关系的集合。从枯草芽孢杆菌数据集中，我们选择840个具有文献证据的调控关系。为了选择这些大网络的部件，我们使用随机程序。由于我们需要以在所选基因中存在一些已知调控关系的方式选择基因，我们首先随机选择少数基因，然后迭代地选择与前面选择的基因相联系的基因，在每步中通过一个基因或至少一个边缘扩展所选的部分网络。在每次迭代中我们选择具有同样概率的部分网络的连接的组分，然后随机选择与组分中至少一个基因具有已知关系的基因，如果这样的基因存在。由于我们应该避免目标网络的无效选择，当连续选择了五个相连的基因时，我们选择不与先前所选的基因相连的基因。
选择程序产生了具有非无效结构的已知基因网络N。我们将N表示为矩阵。每对具有已知关系的基因被表示为矩阵的1项(entry)。对于大多数对来说，不具有知识的基因对被用0项表示，但是对于对(pairs)(g，h)表示为0.5项，对其下列四个条件的一个或多个成立。
1.g被h调控2.在目标网络中存在基因i，其调节g和h3.在目标网络中存在基因i，其被g(h)调节并调节h(g)4.条件2或条件3对目标网络外的基因i成立使用这些条件，几乎正确的预测从错误的预测中被区分出来。如果边缘(g，h)被预测，且(h，g)是已知的调控关系，那么至少两个基因相互作用的事实被正确的预测了。同样的，间接调控关系(很可能通过一些不包括在目标网络本身中的基因)如果被预测了，也不完全错误。
3.3.3评估结果我们随机选择了如第3.3.2部分所述的10个基因的30个集合，并且应用该方法于每个目标网络。然后我们计算在最佳的500个网络结构中每个关系被估计了多少次，并且检测估计的频率与生物学知识之间的关系。注意到我们方法获得的结果是有指向的。但是，如我们在先前的部分所提及的，具有错误定向的估计的边缘也具有信息。因此我们考虑指向的和非指向的例子。
了解每个估计的边缘正确与否是很困难的，因为可能还存在未被发现的关系。另外，已知的基因关系可以不被表达在数据中。因此，我们根据数据库信息将估计的边缘分为三群(1，0.5和0)并且检测它们的区别。我们假设好的估计能够很好的区分这些群。
图2(a)和2(b)显示枯草芽孢杆菌时程的结果及各破坏体数据(disruptantdata)。X轴显示在500个网络中出现每个边缘的百分比。Y轴显示落入间隔的相应群的边缘的分数。由于我们在x轴的两端都观察到两个峰，我们更密切的检测这两个区域。
表3.2在0到10％和90到100％周围的边缘的频率

表3.2显示了频率的百分比。从表3.2我们观察到三个群的趋势。在0％周围，群0相比其它显示超过10％更高的频率。相反地，在100％周围，频率跌至最低值。另一方面，群0.5显示距群1更近的值。比较时程和破坏体数据(disruptant data)，当使用破坏体数据时，群被很好地区分。
表3.3在0到10％和90到100％周围的边缘的频率

最后我们分析大肠杆菌微阵列数据并且结果在图3和表3.3中显示。有趣地，在这些结果中群之间的区别有些小。我们假定建议的方法不仅能够评分评估函数的好处而且还可以评估数据的准确性。
3.4基因网络基序的获取虽然在前面的部分中的评估仅仅是基于边缘的数量，我们现在考虑使用列举的网络的全部信息。
3.4.1基序获取问题定理3.2结合本实施例中描述的技术显示出列举相当大的基因网络的最优和次优网络的可能性。一个使用由列举的网络提供的信息的简单方法是计算网络中每个边缘出现的数量，仅仅选择边缘，所述边缘具有超过一些阈值的计数，并且从这些边缘组成部分网络(partial network)。这个途径在(26)中应用以从通过自助(bootstrap)法列举的网络中获取部分网络。基于边缘的置信水平相互独立的假设，显示出具有高置信水平的显著多个边缘的区域可以在基因网络中被找到。然而，这个假设不说明连接至相同基因g的所有边缘的置信水平依赖于g的表达测量这一事实。
无论如何，即使每个单独的边缘具有超过一些阈值的计数，由这些边缘组成的部分网络不必要在附近的网络中频繁出现。因此，由可能边缘组成的部分网络本身可以是不太可能的。这导致下列问题，我们表示为基因网络基序获取问题假设图表N1＝(G，E1)，...，Nm＝(G，Em)，以及k∈IR，寻找|M|＝k的集合MG×G，其最大化包括M，即MEi的图Ni的数量。
基序获取问题等同于从数据挖掘(data mining)中找寻最大频率项组集合的公知问题。找寻平衡的完整的双向子图的问题(6)能够被归纳为全部两个问题，其因此是NP-hard。然而，能够采用我们下面所描述的穷举搜索策略为实际实例解决该问题。
3.4.2基序获取运算法则(运算法则3.2)为了获取包含在给定的基因表达测定值中的有关基因网络的部分信息，我们建议下面的策略，其采用三个参数n，c，k，∈，IN.

这个运算法则使用每个在基序M中至少c次出现的边缘本身必须具有c次出现的事实。因此，具有超过c次出现的用于构成基序的候选边缘的数量在实践中可以被实质上减少。当阈值c设置的太低时，这个穷举搜索策略的运行时间成为不可实行的，但是这没有对我们的计算施加实际的限制，因为我们对具有高置信性的基因网络基序感兴趣。
3.5.基序获取结果如在第3部分中，我们选择细菌作为我们估计的预期目标。我们首先估算基序获取的可靠性，随后我们应用该方法于涉及枯草芽孢杆菌以及大肠杆菌的色氨酸代谢的基因的群。
3.5.1在细菌基因网络中寻找基序为了评估基序获取运算法则的预测能力，以及将其与单独使用最大可能基因网络模型的预测力相比较，我们将其应用于枯草芽孢杆菌的破坏体数据集(disruptant dataset)。我们以如第3.3部分所述的随机方式选择了25个的8个基因的目标基因网络Ni，列举最可能的100个网络，并且获取带有2个或更多边缘的基序，使用90作为阈值(我们算法中的参数c)。这个计算用1.9GHz的单个CPU运行大约3个小时。我们注意到如果使用超级计算机，列举可以被应用于多至约30个基因的基因网络。
我们随后随机选择最优网络模型的一个边缘，以及在带有多数边缘的基序中间带有最高分值的基序的一个边缘。我们随后使用DBTBS数据检查全部两个边缘的正确性，并且通过将具有已知调控关系的边缘的数量被Ni中可能边缘的总量除来计算随机从已知网络中正确地猜测单个边缘的概率pi。按照第3部分的结果，我们判定1-项和0.5-项是正确的。我们计算在n网络中正确地猜测至少k个单个边缘的概率的上界，P(n，k)，，其通过使用下面的不等式P(n,k)≤Σi=kn(ni)pi(1-p)n-1,---(1)]]>其中p表示maxi＝120pi。该计算的结果在表3.3中提供。我们观察到基序获取方法的结果与最优基因网络的结果相比在指向的(directed)基序的情况下具有更强的显著性，与非指向的(indirected)基序的情况具有相等的显著性。
这是由于在任何有情况下最优网络模型都不是空图，对一些目标网络而言没有超过阈值的基序。由于我们对目标网络的随机选择可能选择一些没有在数据中表达的网络，我们不能期待对任意基因集合都找到高分值基序，因此预测调控关系可能不适合一些目标网络。我们注意到我们的结果也从整体上支持预测，因为我们为评估选择了任意的边缘。
表3.3基序获取运算法则的评估结果

3.5.2比较相关物种的基因网络所描述的方法帮助了解相关物种，诸如革兰氏阳性杆状枯草芽孢杆菌和革兰氏阴性杆状大肠杆菌的基因调控网络的不同和类似处，其例如通过应用描述于微阵列数据的算法，所述数据从大肠杆菌(色氨酸代谢实验)和枯草芽孢杆菌的时程实验和已知涉及已经被详细研究的色氨酸网络的基因中获得。我们获取基于100个列举的最优网络模型的指向的基序，使用50作为阈值。由显示在数据集中找到的具有最高采样数(hits)的最大基序的输出文件产生图，基因通过节点表示并且每个边缘通过其权重标记。
获得大肠杆菌的数据集的最大基序被找到54次并且含有13个权重范围在79到100的边缘(图4)，获得枯草芽孢杆菌数据的被找到61次并且含有权重范围在77到100的15个边缘(图5)。
在两个基序中TrpA，trpB，trpC，trpD，trpE都通过高权重的边缘连接，其形成这个调控网络的核心，这符合它们在两个物种中的trp-操纵子中被置于接近彼此的位置的事实。甚至在trp-操纵子中位置的顺序能够在图表中部分被识别。但是在枯草芽孢杆菌中，七个基因编码trp生物合成途径中的酶，这不同于大肠杆菌中的五个。两个额外的基因，trpF和pabA也包含在衍生的基序中。trpF连接至trpB和trpD，其与它在trp-操纵子上大致的位置相一致。PabA，在枯草芽孢杆菌中也称为trpG，被发现与trpA和trpB紧密连接，虽然其不位于trp-操纵子中，而是在叶酸操纵子citehenner93，gollnick02中。该紧密连接可能显示出在色氨酸生物合成途径中的紧密的功能关系。在大肠杆菌中，也存在pabA基因，但是结果显示不连接至色氨酸网络中的任何一个基因，与没有被发现其被涉及的事实相一致。
另一方面，在枯草芽孢杆菌的图中mtrB与trpE和trpF紧密联系。这可以通过其基因产物，色氨酸RNA-结合衰减蛋白(TRAP)已经被发现存在于枯草芽孢杆菌的色氨酸生物合成的关键调节物中的事实(24)进行解释。有趣地，在大肠杆菌中编码色氨酸特异的转运蛋白的mtr也在图中与类似trpB和trpE的色氨酸生物合成的核心基因相关联。网络基序的结构提示大肠杆菌中mtr的功能可能与枯草芽孢杆菌中的mtrB相似。
近来，TRAP-抑制蛋白(抗-TRAP，AT)已经被发现并且已经被识别为yczA的基因产物(34)。然而，在结果中没有发现yczA和其它被检测基因之间的联系。但是这可能是由于太低的yczA-mRNA水平。或者AT可能遍布于细胞中，不需要从yczA-mRNA中新翻译。这可能主要由于AT不单独作用而是依赖于细胞中色氨酸水平并且仅结合至Trp激活的TRAP(34)。目前为止在大肠杆菌中还没有找到yczA的类似物，但是BLASTP搜索显示AT的氨基酸序列显示出与侣伴蛋白dnaJ(34)的富含半胱氨酸的区域有很高的相似性。因此dnaJ已经包含在本计算中，结果显示与大肠杆菌甚至枯草芽孢杆菌中色氨酸网络的关联。原因可能是侣伴蛋白能够与所有种类的途径相互作用。但是也可能是由于与yczA-mRNA的相互杂交所引起的。或者dnaJ实际上在trp网络中发挥作用。由于在蛋白质水平上的相似性，其甚至可能具有TRAP调节物的功能。这将解释枯草芽孢杆菌的图中dnaJ和mtrB之间的强联系。
编码色氨酰-tRNA合酶的TrpS，已经被包含在这个分析中，因为已经显示tRNATrp在两种细菌的色氨酸生物合成中具有重要，但是不同的功能。各自网络基序的结构实际不同。大肠杆菌的图显示trpS与mtr和wrbA相联系。如上面所提到的，mtr编码色氨酸特异的转运蛋白(23)并且可能涉及在tRNATrp充足时可观察到的转录终止(34)，停止于前导肽，通过trpL编码，对应图中从trpL至mtr的边缘。wrbA编码色氨酸抑制物结合蛋白，所以这个边缘对应于激活色氨酸抑制剂的Trp。然而，在编码色氨酸抑制物的基因trpS和trpR之间没有直接的联系。这可以通过阵列测量mRNA表达，以及色氨酸抑制物蛋白通过Trp结合时构象改变，而不是转录的水平来起作用的事实进行解释。而且图显示trpR与trpD的关联，其可能显示一般情况下色氨酸抑制物蛋白和色氨酸生物合成酶一起产生，尽管trpR不位于trp-操纵子中。在枯草芽孢杆菌图中，trpS表达与pabA和dnaJ相联系。如果dnaJ具有yczA的功能，这与tRNATrp已经被发现与AT钝化的TRAP的形成相关的发现(34)相一致。
虽然，这是有局限性的。微阵列数据常遭受噪音并且仅在mRNA水平进行观察。仅有少量大肠杆菌数据集(18个阵列)；对于枯草芽孢杆菌，没有来自设计用于色氨酸网络分析的试验的实验数据是发明人可以任意利用的。但是大肠杆菌试验已经被设计用于研究色氨酸代谢，并且从在各种培养基上，即在各种营养条件下生长的细菌中获得枯草芽孢杆菌数据集。所以期望在培养基中的色氨酸水平不同，其影响色氨酸基因调控网络。因此，给出这些数据，色氨酸网络为最好的目标，尽管所选的基因集合对评估所述计算方法而言可能不是最优的和合适的。
3.6.结论我们已经提供了为相当大的基因网络列举最优和次优基因网络模型的理论基础，显示了最优模型与现有知识的全面比较的结果，介绍了获取最优网络模型的可靠部分的方法，并应用了该方法。
基于技术发展水平的分值从最优估计中获取最可靠的基序本身是有价值的并且开辟了常规基因网络基序(某种意义的网络基序)分析的方法，在相关物种之间基因网络的对比，等等。
在本发明中所描述的算法方法学不依赖于某种评分方案或某种类型的基因表达测定值。因此这些技术一般在所有情况下均可应用，其中提供具有功能性sG×2G！IR的分值s(其适用于在Bayesian网络框架内所有分值的情况)以及多数其它评分函数。这对基因网络估计技术是一个重要的性质，因为与诸如序列信息(33)或蛋白相互关系数据(14，22)的已有知识结合的新分值的工作正在进行。
对于尺寸超出我们算法计算限制的基因网络，可以应用如(24)中的技术以限制于可能携带生物学上有意义的网络的一部分搜索范围，并且在所选子范围内寻找最优网络模型。我们注意到我们的列举算法能够很好的在这种情况中被应用，从而使得我们找寻基因网络基序的方法不限于小基因网络。
用已知知识严格评估基因网络估计的准确性，如我们在第3.3部分中所构建的，是建立用于比较基因网络估计方法能力的标准的有希望的方法。
本文所引用的每篇参考文献均在此全文作为参考。
参考文献1.T.Chen，H.L.He，G.M.Church.Modeling gene expression withdifferential equations.Pacific Symposium on Biocomputing，429--40，1999.
2.D.M.Chickering.Learning Bayesian networks is NP-complete.In D.Fisher and H.-J.Lenz，editors，Learning from DataArtificial Intelligence andStatistics V，Springer-Verlag，1996.
3.G.F.Cooper，E.Herskovits.A Bayesian method for the induction ofprobabilistic networks from data.Machine Learning，9309--347，1992.
4.N.Friedman，M.Goldszmidt.Learning Bayesian networks with localstructure.Jordan，M.I.(ed.)，Kluwer Academic Publishers，pp.421--459，1998.
5.N.Friedman，M.Linial，I.Nachman，D.Pe′er.Using Bayesian networksto analyze expression data.Journal of Computational Biology，7601--620，2000.
6.M.R.Garey，D.S.Johnson.Computers and intractability，W.H.Freeman and Company，1979.
7.P.Gollnick，P.Babitzke，E.Merino，C.Yanofsky.Bacillus subtilis and itsclosest relativesfrom genes to cells.A.L.Sonenshein，J.A.Hoch，R.Losick.(Eds.)，American Society for Microbiology，2002.
8.A.J.Hartemink，D.K.Gifford，T.S.Jaakkola，R.A.Young.Usinggraphical models and genomic expression data to statistically validate models ofgenetic regulato rynetworks.Pacific Symposium on Biocomputing，6422--433，2001.
9.A.J.Hartemink，D.K.Gifford，T.S.Jaakkola，R.A.Young.Combininglocation and expression data for principled discovery of genetic regulatorynetwork models.Pacific Symposium on Biocomputing，7437--449，2002.
10.L.H.Hartwell，J.J.Hopfield，S.Leibler，A.W.Murray.From molecularto modular cell biology.Nature，402C47--C52，1999.
11.D.Henner，C.Yanofsky.Bacillus subtilis and other gram-positivebacteriabiochemistry，physiology and molecular genetics.A.L.Sonenshein，J.A.Hoch，R.Losick.(Eds.)，American Society for Microbiology，pp.269--280，1993.
12.M.J.L.de Hoon，S.Ott，S.Imoto，S.Miyano.Validation of noisydynamical system models of gene regulation inferred from time-course geneexpression data at arbitrary time intervals.European Conference onComputational Biology，2003.
13.S.Imoto，T.Goto，S.Miyano.Estimation of genetic networks andfunctional structures between genes by using Bayesian networks andnonparametric regression.Pacific Symposium on Biocomputing，7175--186，2002.
14.S.Imoto，T.Higuchi，T.Goto，K.Tashiro，S.Kuhara，S.Miyano.Combining microarrays and biological knowledge for estimating gene networksvia Bayesian networks.Proc.2nd Computational Systems Bioinformatics，104--113，2003.
15.S.Imoto，S.Kim，T.Goto，S.Aburatani，K.Tashiro，S.Kuhara，S.Miyano.Bayesian network and nonparametric heteroscedastic regression fornonlinear modeling of genetic network.Journal of Bioinformatics andComputational Biology，in press，2003；U.S.PatentApplication No10/259,723.
16.T.Ishii，K.Yoshida，G.Terai，Y.Fujita，K.Nakai.DBTBSA databaseof Bacillus subtilis promoters and transcription factors.Nucleic Acids Research，29278--280，2001.
17.A.B.Khodursky，B.J.Peter，M.B.Schmid，J.DeRisi，D.Botstein，P.O.Brown.Analysis of topoisomerase function in bacterial replication fork movementUse of DNA microarrays.Proceedings of the National Academy of Sciences，979419--9424，2000.
18.A.B.Khodursky，B.J.Peter，N.R.Cozzarelli，D.Botstein，P.O.Brown，C.Yanofsky.DNA microarray analysis of gene expression in response tophysiological and genetic changes that affect tryptophan metabolism inEscherichia coli.Proceedings of the National Academy of Sciences，9712170--12175，2000.
19.J.Courcelle，A.Khodursky，B.Peter，P.O.Brown，P.C.Hanawalt.Comparative gene expression profiles following UV exposure in wild-type andSOS-deficient Escherichia coli.Genetics，15841--64，2001.
20.Y.Makita，M.Nakao，N.Ogasawara，K.Nakai.DBTBSDatabase oftranscriptional regulation in Bacillus subtilis and its contribution to comparativegenomics.submitted for publication.
21.R.Milo，S.Shen-Orr，S.Itzkovitz，N.Kashtan，D.Chklovskii，U.Alon.Network motifsSimple building blocks of complex networks.Science，298824--827，2002.
22.N.Nariai，S.Kim，S.Imoto，S.Miyano.Using protein-proteininteractions for refining gene networks estimated from microarray data byBayesian networks.Pacific Symposium on Biocomputing，in press，2004.
23.I.M.Ong，J.D.Glasner，D.Page.Modelling regulatory pathways inE.coli from time series expression profiles.Bioinformatics，18241--248，2002.
24.S.Ott，S.Miyano，Finding optimal gene networks using biologicalconstraints，submitted for publication.
25.S.Ott，S.Imoto，S.Miyano.Finding optimal models for small genenetworks.Pacific Symposium on Biocomputing，in press，2004.
26.D.Pe′er，A.Regev，G.Elidan，N.Friedman.Inferring subnetworks fromperturbed expression profiles.Bioinformatics，17215--224，2001.
27.R.W.Robinson.Counting labeled acyclic digraphs.New Directions inthe Theory of Graphs，pp.239--273，1973.
28.H.Salgado，A.Santos-Zavaleta，S.Gama-Castro，D.Millán-Zárate，E.D″EDaz-Peredo，F.Sánchez-Solano，E.Pérez-Rueda，C.Bonavides-Mart″EDnez，J.Collado-Vides.RegulonDB(version 3.2)transcriptional regulation and operonorganization in Escherichia coli K-12.Nucleic Acids Research，2972--74，2001.
29.S.S.Shen-Orr，R.Milo，S.Mangan，U.Alon.Network motifs in thetranscriptional regulation network of Escherichia coli.Nature Genetics，3164--68，2002.
30.V.A.Smith，E.D.Jarvis，A.J.Hartemink.Evaluating functional networkinference using simulations of complex biological systems.Bioinformatics，18216--224，2002.
31.E.P.van Someren，L.F.A.Wessels，E.Backer，M.J.T.Reinders.Geneticnetwork modeling.Pharmacogenomics，3(4)507--525，2002.
32.A.L.Sonenshein，J.A.Hoch，R.Losic k.Bacillus subtilis and itsclosest relativesfrom genes to cells.ASM Press，Washington，D.C.，2001.
33.Y.Tamada，S.Kim，H.Bannai，S.Imoto，K.Tashiro，S.Kuhara，S.Miyano.Estimating gene networks from gene expression data by combiningBayesian network model with promoter element detection.Bioinformatics，inpress，2003.
34.A.Valbuzzi，C.Yanofsky.Inhibition of the Bacillus subtilis regulatoryprotein TRAP by the TRAP-inhibitory protein，AT.Science，2932057--2059，2001.
35.http://www.ncbi.nlm.nih.gov/geo/
权利要求
1.推导基因网络的方法，包括(a)提供生物体的可能基因网络的推导的模型，其中包括定义搜索范围；(b)选择所述搜索范围的生物学相关的子范围；和(c)通过反复应用最优计算小基因网络的算法在所述选择的子范围中计算最优解决方案。
2.权利要求1的方法，其中所述推导模型是贝叶斯网络估计模型。
3.权利要求1的方法，其中所述生物学相关的子范围包括涉及所述生物体代谢途径的基因。
4.权利要求1的方法，其中所述算法包括下列步骤(a)对所有g∈G计算F(g，φ)＝s(g，φ)；(b)对所有AG，A≠φ和所有g∈G计算F(g，A)为min{s(g,A),mina~∈AF(g,A-{a})};]]>(c)设定M(φ)＝φ；(d)对于所有AG，A≠φ，进行下列步骤(i)计算g*＝arg mingA(F(g，A-{g})+QA-{g}(M(A-{g})))；和(ii)对所有1≤i＜|A|，设定M(A)(i)＝M(A-{g*})(i)，和M(A)(|A|)＝g*；和(e)返回QG(M(G))。
5.权利要求4的方法，其中根据下列步骤修改所述算法(a)在步骤1和步骤2中F的计算中，为所有g∈Si和所有ACg仅计算F(g，A)；和(b)替换步骤4a中的项F(g，A-{g})为F(g，(Cg＝Si)∪(Cg∩A))。
6.权利要求1的方法，其中最优网络N具有定义分值(N)＝∑g∈Gs(g，PN(g))。
7.权利要求1的方法，其中所述算法包括下列步骤(a)在G中聚类基因使得没有聚类大于c个基因；(b)根据大小递减排序聚类C1，...，Cn；(c)对于每个i∈{1，...，n}和对于每个g∈Ci，从C1∪...∪Cn中选择多至m个的候选亲本；(d)使用定理1.2计算最优基因网络模型。
8.权利要求1的方法，其中所述算法包括下列步骤(a)聚合G中的基因于群Ci中，其中|Ci|≤c，并根据生物学知识排列它们C1，...，Cn；(b)对每个i∈{1，...，n}和对每个基因g∈Ci，从C1∪...∪Ci中选出多至m个的候选亲本；和(c)使用定理2计算最优的基因网络模型。
9.权利要求1的方法，其中所述算法包括下列步骤(a)对所有g∈G计算F(g，φ)＝s(g，φ)；(b)对所有AG，A≠φ和所有g∈G计算F(g，A)为min{s(g，A)，mina~∈AF(g,A-{a})};]]>(cc)设定M(φ)＝φ；(d)对所有AG，A≠φ，进行下列两个步骤(i)计算g*＝arg mingA(F(g，A-{g})+QA-{g}(M(A-{g})))；和(ii)对于所有1≤i＜|A|，设定M(A)(i)＝M(A-{g*})(i)，和M(A)(|A|)＝g*；和(e)返回QG(M(G))。
10.权利要求1的方法，其中所述算法包括步骤(a)对所有g∈G设定Fm(g，φ，1)＝φ，Sm(g，φ，1)＝s(g，φ)；(b)对所有g∈G，所有AG，A≠N和所有n≤m进行下列两个步骤(i)从{BA|B＝AvB＝Fm(g，A-{h}，p)，h∈A，p≤m}-{Fm(g，A，p)|p＜n}中选择B*使得s(g，B*)最小化；和(ii)设定Fm(g，A，n)＝B*，Sm(g，A，n)＝s(g，B*)；(c)设定Mm(φ，1)＝φ和Dm(φ，1)＝φ；(d)对于所有AG，φ，和所有n≤m进行下列三个步骤(i)选择三元组(g，p，q)∈A×IN≤m×IN≤m使得分值(QA-{g}(Mm(A-{g}，p)，Dm(A-{g}，p)))+Sm(g，A-{g}，q)最小化并且对于r＜n，(g，p，q)诱导与QA(Mm(A，r)，Dm(A，r))不同的网络；(ii)对i＜|A|设定Mm(A，n)(i)＝Mm(A-{g}，p)(i)，且Mm(A，n)(|A|)＝g；和(iii)假设v表示Dm(A-{g}，p)。对所有I＜|A|且w|A|＝q设定w∈IN|A|为wi＝vi和设定Dm(A，n)＝w；和(e)对所有i≤m返回QG(Mm(G，i)，Dm(G，i))。
11.权利要求1-10的任意一个的方法，其中列举的基因网络的可靠性，包括步骤(a)列举最可能的基因网络模型Ni，1≤i≤n；(b)对每一个g，h∈0G，计数在网络Ni中边缘(edge)(g，h)的出现次数(occurences)；(c)选择所有至少c次出现(occurences)的边缘(g，h)；(d)对|M|＝k的所选边缘的集合的所有子集M，计算包括MI中的所有边缘的网络；和(e)返回至少c次出现的所有目的M。
1.权利要求1-11任意一个的方法，进一步含有计算选自BRNC分值、BDe分值和MDL分值的评分函数。
13.如在此所充分描述的用于检测基因网络的方法。
14.含有使用权利要求1-11任意一个的方法获得的结果的存储介质。
15.含有使用在此所充分描述的方法获得的结果的存储介质。
16.用于确定基因网络关系的系统，其含有用于提供生物体基因的定量表达数据的输入设备；适于接收所述生物体基因的定量表达数据的存储设备；适于进行所述基因间网络关系的贝叶斯网络分析，从而产生反映所述网络关系的数据集的处理器；和显示所述网络关系的所述数据集的输出设备。
17.如在此所充分描述的确定基因网络关系的系统。
全文摘要
从基因表达测定值中准确估计基因网络是生物信息学领域中是一个主要挑战。我们提供一般方法以将搜索范围减少至生物学有意义的子范围和通过使用以生物学相关信息约束的推导模型在线性时间内在子范围中寻找最优解决方案。我们列出该方案用于酵母和枯草芽孢杆菌数据的效果。并且，我们提供适合提供和存贮基因网络关系的数据和结果的系统和存贮介质。
文档编号G06F19/24GK1914510SQ200480041631
公开日2007年2月14日 申请日期2004年12月13日 优先权日2003年12月12日
发明者萨斯卡·奥特, 宫野悟 申请人:Gni株式会社, Gni美国公司