CN1714371A - 根据时间序列基因表达数据的基因网络的非线性模拟 - Google Patents

Abstract

本发明的实施例包括新的推断方法来分析复杂的生物信息，包括基因网络。在一些实施例中，为有机体中的大量基因同时获得时间过程数据。新方法包括贝叶斯推断方法的改进以及那些方法的应用，从而确定所表达的基因之间的因果关系，并且在一些实施例中用于确定调节的基因的上游影响基因。贝叶斯模型的另外的改进包括使用自举方法和确定边缘影响，从而更准确地提供所表达的基因之间的网络信息。关于基因网络的信息可以存储在存储器设备中并可以发送到输出设备，也可以发送到远程位置。

Description

根据时间序列基因表达数据的基因网络的非线性模拟

本申请要求于2002年11月19日提交的美国临时申请序列号为60/427448的申请的优先权。该申请援引于此以供参考。

技术领域

本发明涉及根据基因表达的时间序列研究、使用具有非参数回归的贝叶斯模型来推断基因之间的网络关系。本发明尤其涉及最小化标准即BNRC_dynamic来推断最佳网络关系的方法。

背景技术

在生命科学、医学、药物发现的当前研究和开发以及在制药工业的开发中，最重要的方面之一是需要开发用于解释大量原始数据并基于这些数据得出结论的方法和设备。生物信息学有助于根本上理解生物***，并且预示了对生命***的部分之间的复杂关系提供更深的理解。特别是，随着用于快速检测所表达的基因和用于量化基因表达的新方法的出现，可以使用生物信息学来预测潜在的治疗目标而不需要明确了解特定基因在生物体(biology oforganism)中所起的确定作用。

基因***的模拟是生物***的中心主题。由于可以基于生物知识进行模拟，因此通过预测或推断预先未知的关系，网络估算方法可以支持生物模拟。

特别是，微阵列技术的开发已经允许各种有机体的大量基因表达的研究。可以从有机体的大量基因获得大量的原始数据，并且通过突变、疾病或药物的干预来研究基因表达。发现特定基因表达在特定疾病中增加，或响应于特定干预可以使人相信该基因直接参与疾病过程或药物反应。然而，在生物有机体中，基因很少独立地由任何这样的干预调节，其中特定的干预可以影响许多基因。由于可以这样影响大量不同的基因，在这样的研究中理解基因之间的因果关系是很难的。因此，正在花费大量的努力来开发用于确定基因之间的因果关系、哪些基因是生物现象的中心以及哪个(些)基因的表达对于研究中的生物过程非重要的方法。虽然这些非重要基因表达作为生物或病理生理状态的标记可能是有用的，但如果这样的基因对于生理或病理生理状态是不重要的，也可能不值得努力来开发基于这样的基因的药物。相反，对于那些识别为对过程来说是中心的基因，对于与已改变的基因表达相关的状态，开发药物或其他干预可能对于开发治疗方法是至关重要的。

用于从微阵列基因表达数据中估算基因网络的贝叶斯网络分析的开发已经引起了广泛关注，并且已经报告了许多成功的研究(Friedman et al[13]；Imoto et al[14]；Pe’er et al[18]以及我们自己的成果[作为参考完全结合在这里的美国专利申请号：10/259723])。

然而，传统的贝叶斯网络模型的缺点是它们不能构造循环网络，而确定的真实基因调节机制具有循环部分。近来，已经提出了用于构造具有循环调节部分的基因网络的动态贝叶斯网络模型(Bilmes et a.[3]；Friedman et al[12]；Someren et al[19])。动态贝叶斯网络基于时间序列数据，并且数据通常可以离散(discritized)为几类。因此，动态网络模型可以依赖于用于离散过程的阈值的设置，并且不幸的是该离散化将导致信息丢失。Imoto et al.[14，15]提出基于贝叶斯网络和非线性回归的网络估算方法以便解避免离散化，捕捉基因间的非线性关系。然而，贝叶斯网络和非参数回归模型[14，15]不能充分解决具有循环调节部分的网络。

发明内容

在特定实施例中，本发明包括在具有非参数回归的贝叶斯网络模型中使用时间序列表达数据。利用时间序列表达数据，我们可以识别循环调节部分。在其他实施例中，可以将时间延迟信息合并到贝叶斯/非线性回归模型中，然后这种模型甚至可以提取基因之间的非线性关系。在确定的这些实施例中，可以使用常微分方程作为替代。从贝叶斯统计观点看，我们也已经开发了用于选择最佳网络的新标准。这样的标准可以基于具有噪声的数据优化网络结构。

附图说明

通过结合附图对本发明的优选实施例进行详细描述，本发明的上述目的和优点将会变得更加清楚，其中：

图1描述基因表达的时间动态的示意说明；

图2a和2b描述以KEGG编译的、包含在酵母中的细胞周期调节的基因的网络关系的图表；

图2a描述在循环依赖的蛋白质激酶路径(cyclin-dependent protein kinasepathway)中的基因；

图2b描述图2a中描述的、包含在调节的循环依赖的蛋白质激酶中的基因之间的网络关系；

图3a-3c描述包含在新陈代谢路径中的酵母基因的网络关系的图表。

图3a描述包含在新陈代谢路径中的几个基因；

图3b描述从贝叶斯/非参数回归模型中导出的、在图3a中描述的基因之间的网络关系；以及

图3c描述从动态贝叶斯/非参数回归模型中导出的、在图3a中描述的基因之间的网络关系。

具体实施方式

通常，使用任意适合的、用于确定基因表达的方法可以获得动态贝叶斯网络模型。在特定实施例中，由于可以根据应用到阵列的单个样本研究大量的基因，因此期望微阵列实验，使得很容易确定基因表达中的相对差异。可以期望通过从反映真实基因表达的信号减去背景信号，和/或纠正用于度量基因表达(如cy3/cy5)的标签中的内在差异，来改善微阵列方法的精确度。

我们利用贝叶斯网络框架，将基因考虑为随机变量，并将联合概率分解为条件概率的乘积。例如，如果我们有一系列随机向量的观测，则可以表示获得给定观测的概率可以依赖于条件概率密度。在确定实施例中，可以使用非参数回归模型来捕捉变量间的关系。可以使用各种图形工具来说明这些关系。例如，可以使用多项式、傅立叶级数、回归样条基、B样条基、小波基等来定义基因关系的图形。在美国专利申请序列号10/259723中公开了说明网络关系的确定方法，将其援引于此以供参考。选择合适的图形的困难之处在于合适地评估***中的方差和噪声。

在本发明的一些实施例中，可以使用具有非参数回归的贝叶斯估算构造网络，其中贝叶斯估算使用来自时间序列研究的数据。在许多基因网络中，干预导致某些基因的表达在观察到其它基因改变之前发生改变。可以推断在时间上后发生的干预后的某些基因表达可能与其表达是早期的基因因果地相关。时间序列信息有助于定义“早期”或基因和“后期”或基因。后期基因的表达的改变不太可能是其基因比后期基因在时间上改变得更早的早期基因的表达的改变的起因。虽然该假定不能应用到所有情况中，但是很有可能的是，在网络中早期基因比后期基因更可能是“上游”的，而后期基因更可能是“下游”基因。因此，基因表达的时间关系可以有助于改进贝叶斯估算和非参数回归来提供可靠的网络解。

在本发明的方面中，我们将贝叶斯网络和非参数回归模型扩展为动态贝叶斯网络模型，当具有时间序列基因表达数据时，该模型可以用来构造循环关系。关于基因表达的变化之间的时间延迟的信息可以轻易包括在模型中，并且该模型甚至可以轻易提取非线性关系。

在某些实施例中，为了构造具有循环调节部分的基因网络，可以使用常微分方程模型(Chen et al.[5]；de Hoon et al.[8])。然而，该模型基于线性***，并且可能不适合捕捉复杂的现象。从贝叶斯统计的观点看来，我们已经获得了用于选择最佳网络的新标准。该标准可以优化网络结构，其给出了由具有噪声的数据描述的基因交互的最佳表达。新标准在这里用术语BNRC_dynamic表示。

如图1所示，可以使用一阶马尔可夫关系评估BNRC_dynamic。在这样的关系中，上游基因X_i被描述为对一个或多个下游基因X₂具有影响(右箭头)，X₂对X₃有影响(未示出)等，直到在X_n上观察到影响。在其中X₁不具有其自己的“上游”基因的情况中，在网络中X₁称为“父”基因。受父基因影响的基因称为“目标”基因。应当注意的是，在本上下文中使用的“目标基因”不应与作为关于诸如由潜在的药物之类的干预的对象的基因相混淆。实际上，父基因可以是用于治疗干预的对象。在这种方案中，直到已经引出对X₁、X₂等的影响为止，不能观察到对X_n的影响。应当注意的是，图1图解了“串行(series)”因/果关系，而不存在并行或反馈***，但是在许多基因***中存在串行影响和“并行”影响，其中两个或多个基因可以由上游基因影响，和/或它们可以影响下游基因。此外，可以存在循环影响(“反馈”)，其中X_a可以影响另一个基因X_b，X_b可以影响X_c，而X_c自己可以影响X_a(或X_b)。此外，这样的反馈可以是其中X_c刺激X_a的正反馈，或者是其中X_c抑制X_a的负反馈。在串行、并行、正反馈和负反馈都存在的情况中出现更进一步的复杂度。

通常，时间点之间的关系可能是任意的，但是在一些情况下，基于研究中的基因的生物影响的知识和它们的表达动态，来使用预选的时间点是有利的。在一阶状态下，联合概率可以分解为下面的例1中的方程(1)那样。然后，利用例1中的方程(2)，可以将条件概率分解为条件概率的乘积。方程(1)和(2)可以成立，并且密度函数可以使用概率测度替代。因此，例如根据例1的方程(3)，使用在例1中描述的密度可以表示动态贝叶斯网络，从而得出基因及其父基因的本地网络结构。

例如，如例2所述，可以应用具有非参数回归的动态贝叶斯模型。一旦收集到实验数据，网络的求解可以当作统计模型选择问题。在确定实施例中，我们可以使用贝叶斯方法解决该问题，并导出用于评估动态贝叶斯网络和非参数回归方法的好处的标准。假设可以根据例2中方程(4)确定先验分布(priordistribution)、边际似然性(marginal likelihood)和后验概率。如方程(4)所述，基因网络的后续构造包含高次积分计算。在一些实施例中，例如可以使用积分的拉普拉斯方法来近似该积分。因此，可以求解在例2中的方程(5)所示的标准BRNC_dynamic。

为了将BRNC_dynamic应用到实验***，例如，可以在影响基因***之后在许多时间点上实验地获得cDNA微阵列数据。为了平滑曲线，我们可以使用样条函数，如在例3中描述的B样条。可以根据例3的方程(6)分解BRNC_dynamic。当最小化BRNC_dynamic时获得最佳网络关系。

利用具有非参数回归的动态贝叶斯网络模型和标准BRNC_dynamic，我们可以阐述网络学习过程。然而，当考虑所有可能的基因组合和关系时，确定哪些基因是父基因和哪些是目标基因是非常耗时的。为了降低所需的分析量，我们可以选择候选父基因。然后可以使用贪婪登山(greedy hill-clibming)算法。计算BRNC_dynamic然后添加或删除额外的父基因，并且根据例3的步骤2重新评估BRNC_dyamic。重复这些处理直到找到合适的收敛性。然后，置换计算顺序并重新评估BRNC_dynamic。最佳网络提供最小的BRNC_dynamic。

在图2a和2b中的例4中显示了以上方法的特别说明。通过来自Saccharomyces cerevisiae的基因表达数据的分析，显示这些方法的效率。图2a描述了涉及细胞周期调节的一组S.cerevisiae基因。这些基因根据所包含的整体新陈代谢路径编组进行描述，并且集中在循环依赖蛋白质激酶基因(YBR160w)。应当注意的是，基于图2a，父/目标基因网络关系是未知的。相反，使用本发明的方法，可以在图2b中评估并描述那些基因的网络关系。

在图3a-3c中描述了另一个例子，图3a描述包含在新陈代谢路径中的基因。图3a没有示出基因网络关系。图3b描述使用具有非参数回归的贝叶斯网络分析获得的网络解，但是不考虑BRNC_dynamic。图3c描述通过最小化BRNC_dynamic获得的网络解。应当注意的是，在图3b中，网络关系更简单，并且与图3b描述的那些相比，存在更少假阳性关系(fewer false positiverelationships)(“x”)。

使用本领域中已知的方法例如自举(bootstrap)算法可以确定网络中基因组之间的边界。这样的方法包括使用下列步骤确定边缘的强度：

(1)通过从原始基因库表达数据随机地以移位方式采样多次，以提供自举基因表达矩阵；

(2)对基因gene_i和基因gene_j估算基因网络；

(3)重复步骤(1)和(2)T次，由此产生T个基因网络；及

(4)计算基因gene_i和基因gene_j之间的自举边缘强度为(t₁+t₂)/T。

与诸如贝叶斯或布尔网络之类的其它估算方法相比，新方法的优点包括：(1)可以容易地结合时间信息；(2)不需要诸如离散化的额外数据预处理，可以将微阵列数据作为连续数据分析；和(3)发现更少的假阳性关系。甚至可以由本发明实施例检测并模拟非线性关系。本发明的方法有助于分析基因网络，并有助于开发把控制重要基因的基因表达的特定基因作为目标的新药物。

在下面的例子中描述本发明的其它方面。

例子

下面列出的例子代表本发明的特定实施例。不需要不适当的实验，通过本领域普通技术人员可以开发本发明的其它方面。所有这些实施例被认为是本发明的一部分。

例1：贝叶斯网络和非线性回归

假设我们有n×p微阵列基因表达数据矩阵X，其中n和p分别是微阵列和基因的数量。通常，基因的数量p远远大于微阵列的数量n。在基于贝叶斯网络的基因网络估算中，将基因当作随机变量。当我们通过使用由密度或概率函数描述的统计模型模拟基因网络时，统计模型应当包括p个随机变量。然而，我们仅有n个采样，并且n通常远远小于p。在这种情况下，由于模型具有许多参数而且采样的数量不够估算这些参数，该模式的推断是非常困难或者是不可能的。已经在这样的模拟中建议贝叶斯模型。

在动态贝叶斯网络的上下文中，我们考虑时间序列数据，并且X的第i列向量x_i对应于p个基因在时间i上的状态。至于时间依赖性，我们考虑图1中描述的一阶马尔可夫关系。在这种条件下，联合概率可以被分解如下：

P(X₁₁，...，X_np)＝P(X₁)P(X₂|X₁)×…×P(X_n|X_n-1)，      (1)

其中X_i＝(X_i1，...，X_ip)是在时间i的p个基因的随机变量向量。随机概率P(X_i|X_i-1)也可以分解为条件概率的乘积，其格式为：

P(X_i|X_i-1)＝P(X_i1|P_i-1，1)×…×P(X_ip|P_i-1，p)，       (2)

其中，P_i-1，j是在时间i-1的第j基因的父基因的状态向量。当我们使用密度函数而不是概率测度时，方程(1)和(2)仍成立。因此，动态贝叶斯网络可以通过使用密度表示如下：

$f\left(x_{i1,...,}{x}_{\mathrm{np}}\right)=f_1\left(x_1\right)f_2\left(x_2\right|x_1)\×\·\·\·\×f_n\left(x_n\right|x_{n-1})$

$=f_1\left(x_1\right)\underset{i=2}{\overset{n}{\mathrm{\Π}}}{g}_1\left(x_{i1}\right|p_{i-\mathrm{1,1}})\×\·\·\·\×g_p\left(x_{\mathrm{ip}}\right|p_{i-1,p})$

$=f_1\left(x_1\right)\underset{j=1}{\overset{p}{\mathrm{\Π}}}\left\{\underset{i=2}{\overset{n}{\mathrm{\Π}}}{g}_j\left(x_{\mathrm{ij}}\right|p_{i-1,f})\right\}.$

这里我们具有来自(2)的分解

f_i(X_i|X_i-1)＝g₁(X_i1|P_i-1，1)×…×g_p(X_ip|P_i-1，p)，

其中 $P_{i-1,j}=(p_{i-\mathrm{1,1}}^{\left(j\right)},....p_{i-1,q_j}^{\left(j\right)})$ 是父基因的q_j维观察向量。

为了模拟在xij和P_i-1，j之间的关系，我们使用如下非参数加法回归：

$x_{\mathrm{ij}}=m_{j1}\left(p_{i-\mathrm{1,1}}^{\left(j\right)}\right)+\·\·\·+m_{j{q}_j}\left(p_{i-1,q_j}^{\left(j\right)}\right)+{\mathrm{\ϵ}}_{\mathrm{ij}},$

其中ε_ij通常独立正常地依赖于均值0和方差σ_j ²。这里m_jk(·)是从R到R的光滑函数，并且可以通过使用基本函数的线性组合表示为

$m_{\mathrm{jk}}\left(p_{i-1,k}^{\left(j\right)}\right)=\underset{m=1}{\overset{M_{\mathrm{jk}}}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{mk}}^{\left(j\right)}{b}_{\mathrm{mk}}^{\left(j\right)}\left(p_{i-1,k}^{\left(j\right)}\right),k=1,...,q_j,$

其中

是未知系数参数，并且

是规定的基本函数组。然后，我们定义动态贝叶斯网络和非参数回归的形式为

$f(x_{11},...,x_{\mathrm{np}};{\mathrm{\θ}}_G)$

$f_1\left(x_1\right)\underset{j=1}{\overset{p}{\mathrm{\Π}}}\left[\underset{i=2}{\overset{n}{\mathrm{\Π}}}\frac{1}{\sqrt{{2\mathrm{\π\σ}}_j^2}}\exp \right\{-\frac{{(x_{\mathrm{ij}}-\mathrm{\μ}\left(p_{i-1,j}\right))}^2}{{2\mathrm{\σ}}_j^2}\left\}\right],$

其中 $\mathrm{\μ}\left(p_{i-1,j}\right)=m_{j1}\left(p_{i-\mathrm{1,1}}^{\left(j\right)}\right)+...+m_{{\mathrm{jq}}_j}\left(p_{i-1,q_j}^{\left(j\right)}\right).$ 当第j基因没有父基因时，μ(p_i-1，j)产生常数μ_j。

我们假设f₁(x₁)＝g₁(x₁₁)×...×g₁(x_1p)，那么联合密度f(X₁₁，...X_np；θ_G)可以重写为

$f(x_{11},...,x_{\mathrm{np}};{\mathrm{\θ}}_G)=\underset{j=1}{\overset{P}{\mathrm{\Π}}}\left[g_1\left(x_{\mathrm{ij}}\right)\underset{i=2}{\overset{n}{\mathrm{\Π}}}\frac{1}{\sqrt{{2\mathrm{\π\σ}}_j^2}}\exp \right\{-\frac{{(x_{\mathrm{ij}}-\mathrm{\μ}\left(p_{i-1,j}\right))}^2}{{2\mathrm{\σ}}_j^2}\left\}\right]$

$=\underset{j=1}{\overset{p}{\mathrm{\Π}}}\underset{i=1}{\overset{n}{\mathrm{\Π}}}{g}_j\left(x_{\mathrm{ij}}\right|p_{i-1,j};{\mathrm{\θ}}_j),---\left(3\right)$

其中 $p_{0_j}=\mathrm{\φ}.$ 因此，g_j(X_ij|P_i-1，j；θ_j)表示第j基因和其父基因的本地结构。

例2：获得用于选择网络的标准

当我们固定网络结构时，可以构造在前一节引入的动态贝叶斯网络和非参数回归模型，并且由合适的过程进行估算。然而，基因网络通常是未知的，并且我们应该根据数据估算最佳网络。这个问题可以看作是统计模型选择问题(见Akaike[1]；Konishi and Kitagawa[17]；Burnham and Anderson[4]；Konishi[16])。我们从贝叶斯统计方法解决该问题，并获得用于评估动态贝叶斯网络和非参数回归模型的好处的标准。

假设在动态贝叶斯网络和非参数回归模型中π(θ_G|λ)为关于θ_G参数的先验分布，并且假设π(θ_G|λ)＝O(n)。该边际似然性可以表示为

∫f(X₁₁，...，X_np;θ_G)π(θ_G|λ)dθ_G.

因此，当给定数据时，网络G的后验概率是

${\mathrm{\π}}_{\mathrm{post}}\left(G\right|X)\frac{{\mathrm{\π}}_{\mathrm{prior}}\left(G\right)\∫f(x_{11},...,x_{\mathrm{np}};{\mathrm{\θ}}_G)\mathrm{\π}\left({\mathrm{\θ}}_G\right|\mathrm{\λ})d{\mathrm{\θ}}_G}{{\mathrm{\Σ}}_G\{{\mathrm{\π}}_{\mathrm{prior}}\left(G\right)\∫f(x_{11},...,x_{\mathrm{np}};{\mathrm{\θ}}_G)\mathrm{\π}\left({\mathrm{\θ}}_G\right|\mathrm{\λ})d{\mathrm{\θ}}_G\}},---\left(4\right)$

其中π_prior(G)是网络G的先验概率。方程(4)的分母与模型评估无关。因此网络的评估依赖于分子的大小。因此，我们可以作为下式的最大值选择最佳网络π_prior(G)∫f(x₁₁，...，x_np；θ_G)π(θ_G|λ)dθ_G.

非常清楚的是，构造网络选择标准的基本点是怎样计算高次积分。Imotoet al.[14 15]使用用于积分的拉普拉斯近似(同样见Tinerey and Kadane[21]；Davison[6])，并且我们可以将该技术直接应用到动态贝叶斯网络模型和非参数回归中。因此，我们拥有一个标准即所谓BRNC_dynamic，其形式为

${\mathrm{BNRC}}_{\mathrm{dynamic}}\left(G\right)$

$=-2\log \{{\mathrm{\π}}_{\mathrm{prior}}\left(G\right)\∫f(x_{11},...,x_{\mathrm{np}};{\mathrm{\θ}}_G)\mathrm{\π}\left({\mathrm{\θ}}_G\right|\mathrm{\λ}){\mathrm{d\θ}}_G\}$

$\≈-2{\log \mathrm{\π}}_{\mathrm{prior}}\left(G\right)-r\log (2\mathrm{\π}/n)+\log \left|J_{\mathrm{\λ}}\left({\widehat{\mathrm{\θ}}}_G\right)\right|-2{\mathrm{nl}}_{\mathrm{\λ}}\left({\widehat{\mathrm{\θ}}}_G\right|X_n).---\left(5\right)$

其中r是θ_G的维数

$l_{\mathrm{\λ}}\left({\mathrm{\θ}}_G\right|X_n)=\log f(x_{11},...,x_{\mathrm{np}};{\mathrm{\θ}}_G)/n+\log \mathrm{\π}\left({\mathrm{\θ}}_G\right|\mathrm{\λ})/n,$

$J_{\mathrm{\λ}}\left({\mathrm{\θ}}_G\right)=-{\∂}^2\left\{l_{\mathrm{\λ}}\left({\mathrm{\θ}}_G\right|X_n)\right\}/{\∂\mathrm{\θ}}_G{\∂\mathrm{\θ}}_G^T$

并且 是l_λ(θ_G|X_n)。选择最佳图形使得标准BRNC_dynamic(5)最小。

例3：基因网络的估算

在这一节中，我们示出用于从cDNA微阵列时间序列基因表达数据中估算基因网络的具体策略。

3.1非线性回归

我们使用基本函数方法来构造第2节中描述的光滑函数m_jk(·)。在本文中我们使用B-样条(de Boor[7])作为基本函数。De Boor的算法(de Boor[7]，第10章，130页(3))是用于计算任意阶B-样条的有用方法。我们使用20个具有等距节点(knot)的B-样条(同样见Dierckx[10]；Eilter and Marx[11]来得到B-样条的详细内容)。

3.2关于模型中参数的先验分布

对于参数θG的先验分布，假设参数向量θ_j是相互独立的，那么先验分布可以分解为 ${\mathrm{\π}}_j\left({\mathrm{\θ}}_G\right|\mathrm{\λ})={\mathrm{\Π}}_{j=1}^p{\mathrm{\π}}_j\left({\mathrm{\θ}}_j\right|{\mathrm{\λ}}_j).$ 假设先验分布π_j(θ_j|λ_j)被因式分解为 ${\mathrm{\π}}_j\left({\mathrm{\θ}}_j\right|{\mathrm{\λ}}_j)={\mathrm{\Π}}_{j=1}^{q_j}{\mathrm{\π}}_{\mathrm{jk}}\left({\mathrm{\γ}}_{\mathrm{jk}}\right|{\mathrm{\λ}}_{\mathrm{jk}}),$ 其中λj_k是多下标(hyper)参数。我们使用奇异M_jk变量正态分布作为γ_jk的先验分布，

${\mathrm{\π}}_{\mathrm{jk}}\left({\mathrm{\γ}}_{\mathrm{jk}}\right|{\mathrm{\λ}}_{\mathrm{jk}})={\left(\frac{2\mathrm{\π}}{{\mathrm{n\λ}}_{\mathrm{jk}}}\right)}^{-(M_{\mathrm{jk}}-2)/2}{\left|K_{\mathrm{jk}}\right|}_{+}^{1/2}\exp (-\frac{{\mathrm{n\λ}}_{\mathrm{jk}}}{2}{\mathrm{\γ}}_{\mathrm{jk}}^T{K}_{\mathrm{jk}}{\mathrm{\γ}}_{\mathrm{jk}}),$

其中K_jk是满足 ${\mathrm{\γ}}_{\mathrm{jk}}^T{K}_{\mathrm{jk}}{\mathrm{\γ}}_{\mathrm{jk}}={\mathrm{\Σ}}_{\mathrm{\α}=3}^{M_{\mathrm{jk}}}{({\mathrm{\γ}}_{\mathrm{\αk}}^{\left(j\right)}-{2\mathrm{\γ}}_{\mathrm{\α}-1,k}^{\left(j\right)}+{\mathrm{\γ}}_{\mathrm{\α}-2,k}^{\left(j\right)})}^2$ 的M_jk×M_jk对称半正定矩阵。关于θ_G的该先验分布的这种设置与Imoto et al.[14，15]相同，并且详细内容存在于那些文章中。

3.3建议的标准

通过使用在笫4.2节中描述的先验分布，BRNC_dynamic可以分解为如下：

${\mathrm{BNRC}}_{\mathrm{dynamic}}=\underset{j=1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{BNRC}}_{\mathrm{dynamic}}^{\left(j\right)},---\left(6\right)$

其中BRNC_dynamic ^(j)是第j基因的标准分数，并且由下式定义

${\mathrm{BNRC}}_{\mathrm{dynamic}}^{\left(j\right)}$

$=-2\log \left\{{\∫\mathrm{\π}}_{\mathrm{prior}}\left(L_j\right)\underset{i=1}{\overset{n}{\mathrm{\Π}}}{g}_j\left(x_{\mathrm{ij}}\right|p_{i-1,j};{\mathrm{\θ}}_j){\mathrm{\π}}_j\left({\mathrm{\θ}}_j\right|{\mathrm{\λ}}_j){\mathrm{d\θ}}_j\right\}$

$\≈-2{\log \mathrm{\π}}_{\mathrm{prior}}\left(L_j\right)-r_j\log (2\mathrm{\π}/n)+\log \left|J_{{\mathrm{\λ}}_j}^{\left(j\right)}\left({\widehat{\mathrm{\θ}}}_j\right)\right|-2{\mathrm{nl}}_{{\mathrm{\λ}}_j}^{\left(j\right)}\left({\widehat{\mathrm{\θ}}}_j\right|X),$

其中r_j是θ_j的维数，

$l_{{\mathrm{\λ}}_j}^{\left(j\right)}\left({\widehat{\mathrm{\θ}}}_j\right|X)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\log g_h\left(x_{\mathrm{ij}}\right|p_{i-1,j};{\mathrm{\θ}}_j)/n+\log \mathrm{\π}\left({\mathrm{\θ}}_j\right|{\mathrm{\λ}}_j)/n,$

$J_{{\mathrm{\λ}}_j}^{\left(j\right)}\left({\widehat{\mathrm{\θ}}}_j\right)=-{\∂}^2\left\{l_{{\mathrm{\λ}}_j}^{\left(j\right)}\left({\widehat{\mathrm{\θ}}}_j\right|X)\right\}/{\∂\mathrm{\θ}}_j{\∂\mathrm{\θ}}_j^T$

并且

是l_λj ^(j)(θ_j|X)的模。这里π_prior(L_i)是满足 ${\mathrm{\Σ}}_{j=1}^p\log {\mathrm{\π}}_{\mathrm{prior}}\left(L_j\right)=\log {\mathrm{\π}}_{\mathrm{prior}}\left(G\right)$ 的先验概率。我们设置本地结构的先验概率π_prior(L_j) 为

通过使用动态贝叶斯网络和非线性回归模型以及建议的标准BNRC_dynamic，我们可以阐明网络学习过程如下：从(3)和(6)很清楚地得出，网络结构的优化等效于调节目标基因的父基因的选择。然而当考虑多有基因组合作为父基因时，这是个耗时的任务。因此，我们通过选择候选父基因来消减学习空间。在这个步骤之后，采用贪婪登山算法来找出更好的网络。我们的算法可以表示为如下步骤：

步骤1：预处理阶段

我们生成其第(i，j)元素是图形“gene_i→gene_j”的BNRC分数的p×p矩阵，并且定义给出小BNRC分数的gene_j的父基因的候选组。我们设置父基因10的候选组的元素数。

步骤2：学习阶段

对于贪婪登山算法，我们从空网络开始并且重复下面的步骤：

步骤2-1：对于gene_i，实现添加父基因，删除父基因两个过程之一，所述父基因给出较小BNRC_dynamic分数。

步骤2-2：对基因的规定的计算顺序，重复步骤2-1直到满足合适的收敛标准。

步骤2-9：置换计算顺序来找出更好的解并重复步骤2-1和2-2。

步骤2-4：我们选择给出最小BNRC_dynamic分数的最佳网络。

例4：计算实验

我们通过对由Spellman et al.[20]收集的Saccharomyces cerevisiae细胞周期基因表达数据的分析，演示了本发明的一个实施例。该数据包括两个短时间序列(两个时间点；cln3、clb2)和四个中型时间序列(18、24、17和14个时间点；alpha、cdc15、cdc28和elu)。在基因网络的估算中，我们使用四个中型时间序列。为了合并四个中型时间序列，当我们拟合非参数回归模型时，忽略关于每时间序列中的目标基因的第一观察和父基因的最后一观察。

首先，我们关注汇编在KEGG数据库[22]中的细胞周期路径。目标网络在环绕CDC28(YBR16Ow；循环依赖蛋白质激酶)。该网络包含45个基因，并且登记在KEGG中的路径如图2(a)所示，以及估算的网络在图2(b)中。虚圆中的边缘可以当作是正确的边缘。因此，我们模拟了一些正确的关系。我们通过与边缘相邻的圆表示正确的估算。三角形表示正确方向的反向或跳过。“x”符号表示不正确的关系。

用于演示我们的方法的第二个例子是由DeRisi et al.[9]报告的新陈代谢路径。该网络包含57个基因，并且在图3(a)中示出了目标路径。

我们将贝叶斯网络和非参数回归模型[14，15]应用到该数据，并在图3(b)中描述的产生的网络。通过动态贝叶斯网络和非参数回归模型获得了图3(c)的网络。很难从cDNA微阵列数据中估算新陈代谢路径。然而，我们的模型检测到了基因间的正确关系。与贝叶斯网络和非参数回归相比，在图3(c)中描述的该方法的假阳性(false positive)的数量远远小于在图3(b)中由“x”符号表示的假阳性的数量。

所有在这里引用的参考文献作为其整体合并在此。

参考文献：

1.Akaike，J.：Information theory and an extension of the maximumlikelihood principle.In：Petrov，B.N.，Csaki，F.(eds.)：2nd InternationalSymposium on Information Theory.Akademiai Kiado，Budapest pp：267-281(1973).

2.Berger，J.O.：Statistical Decision theory and Bayesian analysis.Springer-Verlag New York(1985).

3.Bihnes，J.A.：Dynamic Bayesian multinets.Proc.16th Conference onUncertainty in Artificialhitelligence.pp：38-45(2000).

4.Burnham，K.P.，Anderson，D.R.：Model selection and inference，apractical information-theoretical approach.Springer-Verlag New York(1988).

5.Chen，Tl.，He，H.L.，Church，G.M.：Modeling gene expression withdifferential equations.Proc.Pacific Symposium on Biocomputing 4：29-40(1999).

6.Davison，A.C.：Approximate predictive likelihood.Biometrika 73：323-332(1986).

7.DeBoor，C.：A pracitial guide to splines.Springer-Verlag Berlin(1978).

8.De Hoon，M.J.L.，Imoto，S.，Kobayashi，K.，Ogasawara，NH.，Miyano，S.：Inferring gene regulatory networks from time-ordered gene expression data ofBacillus subtilis using differential equations.Proc.Pacific Symposium onBiocomputing 8：2003，in press.

9.DeRisi，J.，Lyer，V.R.，Brown，P.O.：Exploring the metabolic and genecontrol of gene expression on a genonmic scale.Science 278：680-686(1997).

10.Dierckx，P.：Curve and surface fitting with splines.Oxford(1993).

11.Eiler，P.H.C.，Marx，B.：Flexible smoothing with B-splines and penalites(with discussion).Statistical Science 11：89-121(1996).

12.Friedman，N.，Murphy，K..，Russell，S.：Learning the structure ofdynamic probabilistic networks.Proc.Conf.On Uncertainty inArtificialkgtelligence pp：139-147(1998).

13.Firedman，N.，Linial，Ml，Nachman，I.，Pe′er，D.：Using Bayesian networkto analyze expression data.J.Comp.Biol.7：601-620(2000).

14.Imoto，S.，Goto，T.，Miyano，Sl.：Estimationof gnetic networks andfunctional structures between genes by using Bayesian network and noparametricregression.Proc.Pacific Symposium on Biocomputing 7：175-186(2002).

15.Imoto，S.，Kim，S.，Goto，T.，Aburatani，S.，Tashiro，K.，Kuhara，S.，Mjiyano，S.：Bayesian network and nonparametric heteroscedastic regression fornonlinear modeling of genetic network.Proc.IEEE Computer SocietyBioinformatics Conference；PP：219-227(2002).

16.Konishi，S.：Statistical model evaluation and information criteria.In：Ghosh，S.(ed.).Multivariate Analysis，Design of Experiments and SurveySampling.Marcel Dekker，New York，pp：369-399(1999).

17.Konishi，S.，Kitagawa，G.：Generalized information criteria in modelselection.Biometrika 83：875-890(1996).

18.Pe′er，D.，Regev，A.，Elidan，G.，Friedman，N.：Inferring subnetworksfrom perturbed expression profiles.Bioinformatics 17：215-224(ISBM 2001).

19.Someren，E.V.，Wessels，L.，Reinders，M.：Linear modeling of geneticnetworks from experimental data.Bioinformatics 18：355-366(ISBM 2002).

20.Spellman，P.T.，Sherlock，G.，Zhang，M.Q.，Iyer，V.R.，Anders，K.，Eisen，M.B.，Brown，P.O.，Botstein，D.，Futcher，B.：Comprehensive identification ofcell cycle-regulated genes of the yeastSaccharofnyces ceyvisiae by microarrayhybridization.Molecular Biology of the Cell 9：3273-3297(1998).

21.Tinerey，L.，Kadane，J.B.：Accurate approximations for posterior momentsand marginal densities.J.Amer.Statist.Assoc.81：82-86(1986).

Claims (44)

1.一种用于构造基因网络的方法，包括步骤：

(a)为有机体的基因组提供定量时间过程数据库，所述库包括基于所述基因组中的每个基因的表达的时间过程的表达结果，量化平均影响和关于所述基因相互之间的每个时间点的可变性的度量；

(b)根据所述库创建基因表达矩阵；

(c)产生所述基因之间的网络关系；和

(d)确定一个或多个基因组是否与其它所述基因组不同地被表达。

2.如权利要求1所述的方法，还包括步骤：

(e)提供贝叶斯计算模型，其中所述贝叶斯模型包括最小化BNRC_dynamic标准。

3.如权利要求2所述的方法，其中最小化BNRC_dynamic标准的所述步骤包括使用从由多项式基、傅立叶级数、小波基、回归样条基和B-样条组成的组中所选择的非线性曲线拟合方法。

4.如权利要求1所述的方法，其中使用时间过程研究创建所述数据库来改变基因表达。

5.如权利要求2所述的方法，其中最小化所述BNRC_dynamic标准的所述步骤还包括使用向后拟合算法选择贝叶斯模型。

6.如权利要求2所述的方法，其中最小化BNRC_dynamic标准的所述步骤还包括使用Akaike的信息标准。

7.如权利要求2所述的方法，其中最小化BNRC_dynamic标准的所述步骤还包括使用最大似然估算。

8.如权利要求1所述的方法，其中所述基因与细胞周期相关。

9.如权利要求2所述的方法，其中所述可变性度量是方差。

10.如权利要求3所述的方法，其中所述非线性曲线拟合方法是非参数方法。

11.如权利要求10所述的方法，其中用于最小化BNRC_dynamic标准的所述非参数方法包括使用非均匀误差方差(heterogeneous error variance)。

12.如权利要求11所述的方法，其中最小化BNRC_dynamic标准的所述步骤还包括步骤：

(1)产生其第(i，j)元素是图形gene_i→gene_j的BNRC_dynamic ^j分数的分数矩阵；

(2)实现提供最小BNRC_dynamic的添加、删除和反向之一或多个；和

(3)重复步骤2直到BNRC_dynamic不再减少。

13.如权利要求11所述的方法，其中最小化BNRC_dynamic标准的所述步骤还包括应用登山算法来最小化 的步骤。

14.如权利要求11所述的方法，其中使用自举方法确定边缘强度。

15.如权利要求14所述的方法，其中所述自举方法包括步骤：

(1)从原始基因库表达数据中，通过移位地、多次随机采样来提供自举基因表达矩阵；

(2)为gene_i和gene_j估算基因网络；

(3)重复步骤(1)和(2)T次，由此产生T个基因网络；和

(4)计算gene_i和gene_j之间的自举边缘强度，作为(t₁+t₂)/T。

16.一种用于说明基因网络的方法，包括步骤：

(a)为有机体的多个基因提供时间过程基因表达数据的原始数据库；

(b)从所述原始数据库中减去背景信号强度；

(c)为所述多个基因的每个计算基因表达的相对变化；

(d)使用学生t-检验分析基因表达中所述相对变化的统计显著性；和

(e)将基因表达中的所述变化拟合为线性样条函数。

17.如权利要求16所述的方法，还包括从考虑中移除那些其表达级别足够低以至于主要由噪声确定的基因。

18.如权利要求1所述的方法，包括将所述基因编组为一个或更多等效组。

19.一种用于估算基因网络关系的方法，包括步骤：

(1)产生其第(i，j)元素是图形gene_i→gene_j的BNRC分数的p×p矩阵；

(2)选择给出小BNRC分数的gene_j的父基因的候选组；

(3)选择所述父基因的计算顺序；

(4)重复下面的步骤；

(4.1)对于gene_j，添加父基因或删除父基因；

(4.2)重新计算BNRC_dynamic分数；

(4.3)重复步骤3.1和3.2直到满足合适的收敛；

(5)置换步骤(3)中的所述父基因的计算顺序；

(6)重复步骤(4)；和

(7)重复步骤(5)和(6)直到最小化BNRC_dynamic。

20.一种用于根据时间过程基因表达数据构造包含基因网络的***的基因网络模型的方法，所述方法包括使用贝叶斯计算模型，其中所述贝叶斯计算模型包括最小化BNRC_dynamic标准。

21.如权利要求20所述的方法，其中最小化BNRC_dynamic标准的所述步骤包括使用从由多项式基、傅立叶级数、小波基、回归样条基和B-样条组成的组中选择的非线性曲线拟合方法。

22.如权利要求20所述的方法，其中最小化BNRC_dynamic标准包括选择使用向后拟合算法的贝叶斯模型。

23.如权利要求20所述的方法，其中最小化BNRC_dynamic标准包括使用Akaike的信息标准。

24.如权利要求20所述的方法，其中最小化BNRC_dynamic标准包括使用最大似然估算。

25.如权利要求20所述的方法，其中最小化BNRC_dynamic标准包括使用非线性曲线拟合方法，其中非线性曲线拟合方法是非参数方法。

26.如权利要求25所述的方法，其中非参数方法包括使用非均匀误差方差。

27.如权利要求26所述的方法，其中最小化BNRC_dynamic标准还包括步骤：

(1)产生其第(i，j)元素是图形gene_i→gene_j的 分数的分数矩阵；

(2)实现提供最小BNRC_dynamic的添加、删除和反向之一或多个；和

(3)重复步骤2直到BNRC_dynamic不再减少。

28.如权利要求26所述的方法，其中最小化BNRC_dynamic标准还包括应用登山算法来最小化 的步骤。

29.如权利要求26所述的方法，其中使用自举方法确定边缘强度。

30.如权利要求29所述的方法，其中所述自举方法包括步骤：

(1)从原始基因库表达数据中通过移位地多次随机采样来提供自举基因表达矩阵；

(2)为gene_i和gene_j估算基因网络；

(3)重复步骤(1)和(2)T次，由此产生T个基因网络；和

(4)计算gene_i和gene_j之间的自举边缘强度，作为(t₁+t₂)/T。

31.一种包括由权利要求20的方法构造的基因网络模型的数据文件。

32.如权利要求31所述的数据文件，其为计算机可读格式。

33.如权利要求31所述的数据文件，其可从远程位置访问。

34.如权利要求31所述的数据文件，其可从因特网万维网位置访问。

35.一种在包含基因网络的***中识别目标基因的方法，包括：

(a)使用贝叶斯计算模型来构造第一和第二基因网络模型，

其中所述贝叶斯计算模型包括最小化BNRC_dynamic标准，其中通过分析第一基因表达文档获得第一基因网络模型，并通过分析第二基因表达文档获得第二基因网络模型，和

其中在第一时间点从***获得第一基因表达文档，并在所述第一时间点之后的第二时间点从***获得第二基因表达文档，和

(b)使用所述贝叶斯计算模型分析第一和第二基因网络模型，其中基因表达的时间过程被量化，并且父基因被识别为目标基因。

36.如权利要求35所述的方法，其中目标基因是父基因。

37.如权利要求35所述的方法，其中目标基因是父基因的下游基因。

38.一种包含根据权利要求35所述的方法获得的一个或多个目标基因的标识的数据文件。

39.如权利要求38所述的数据文件，其为计算机可读格式。

40.如权利要求38所述的数据文件，其可从远程位置访问。

41.如权利要求38所述的数据文件，其可从因特网万维网位置访问。

42.一种提供服务的方法，包括：

(1)从合作方接收数据组，所述数据组包括一组基因的时间过程表达数据，和

(2)通过最小化BNRC_dynamic标准确定所述组中的基因之间的网络关系。

43.如权利要求42所述的方法，其中接收所述数据组包括接收所述基因的至少一个标识。

44.一种提供服务的方法，包括从合作方接收代理，并且使用根据权利要求35的方法构造的基因网络模型为合作方识别目标基因。

Images