CN101105841A - 由大规模基因芯片表达谱数据构建基因调控亚网络的方法 - Google Patents

Abstract

本发明公开了一种由大规模基因芯片表达谱数据构建基因调控亚网络的方法。包括如下步骤：1)基因芯片表达谱缺失数据值估计；2)利用贝叶斯网络构建基因间的调控网络；3)对原始的表达数据重新抽样，重复步骤2)，得到一个可能的贝叶斯网络的集合；4)统计分析网络特征，重构显著性的亚网络模块。本发明实现从大规模基因芯片实验数据中获取多个基因间调控关系，这个结果是传统生物学实验无法得到的；它在一定程度上弥补了芯片数据不足而导致贝页斯网络学习噪音比较大的缺点。通过构建小的亚网络，在一定程度上来说，紧密调控基因之间的关系可以通过一致图的方法来对那些不是直接调控的基因间关系去噪。得到的亚网络为下步的生物学实验提供更好指导。

Description

由大规模基因芯片表达谱数据构建基因调控亚网络的方法

技术领域

本发明涉及一种由大规模基因芯片表达谱数据构建基因调控亚网络的方法。

背景技术

微阵列分析技术的出现，使得我们可以同时监控整个基因组范围内的基因的表达谱。研究人员通过微阵列杂交实验获取了大量的基因表达量的数据；与此同时，我们需要开发新的分析方法来分析这一数据，从数据中获取信息。怎样由海量的表达轮廓挖掘出具有生物学意义的基因间相互作用的调控网络，成为目前生物信息学一个很重要的领域。

在对基因芯片表达谱数据的分析和调控网络的构建方面，传统的聚类方法能够在多次实验的基础上，将具有相似表达特征的基因归为一组(Spellman，1998)。用这种方法可以发现被共同调控或者具有相似功能的基因。然而，由于微阵列实验的数量和其本身模型的大小相比极其不足，聚类并不能为我们提供更深层次的信息，包括基因间关系的精细结构，基因间调控是否直接，还是通过其他基因的调控。因此，要构建合理，对生物学家有更深刻知道的代谢调控网络，我们通常采取复杂随机过程的概率模型分析基因间的依赖关系，为我们进一步推测它们之间的因果关系提供可能。

我们工作主要的任务是利用贝叶斯网络这一基于复杂随机过程的概率模型来对基因表达数据建模，构建调控网络以及重建具有生物学意义的功能模块。利用贝叶斯网络分析基因表达具有以下好处：1)贝叶斯网络的统计学基础以及学习算法比较成熟，并且在多个领域成功的被应用；2)贝叶斯网络适合在数据样本比较少，而模型比较大时的学习，这在目前的基因芯片表达谱数据分析中尤为适用(Friedman，2000)；3)贝叶斯网络不仅可以用来推测依赖关系，还可以用来推测因果关系(即基因间调控关系)。我们以离散型随机变量为例。贝叶斯网络概率图模型是一个有向无环图，图中每一个节点对应一个随机变量，Heckerman(1995)等人采用BDe函数评判每一个候选的贝叶斯网络建模。经过对每一个候选网络进行启发式搜索，最后在整个搜索空间中得到一个或一个等价类的最优化贝叶斯网络。对于得到的一组候选网络，采用重新抽样的方法对得到的贝叶斯网络进行特征置信度的分析，从而可以得到一个置信度较高的特征集合。

另外，利用基因网络的2阶图对由聚类得到的多重网络进行分析，可以有效地鉴定出具有生物学意义的功能模块(Hu，2005)。然而，基于Pearson相关系数的聚类分析，这个无向网络图只能粗略的构建出功能模块的概况，尽管它可以为生物学家的研究提供进一步的方向。而且，这个功能模块图不能反应基因间的直接相互作用关系，只能反应它们的表达具有显著的相关性。因此，在这里我们提出基于贝叶斯网络构建的模型的有向的功能模块图的重构。

发明内容

本发明的目的是提供一种由大规模基因芯片表达谱数据构建基因调控亚网络的方法。

由大规模基因芯片表达谱数据构建基因调控亚网络的方法包括如下步骤：

1)基因芯片表达谱缺失数据值估计；

2)利用贝叶斯网络构建基因间的调控网络；

3)对原始的表达数据重新抽样，重复步骤2)，得到一个可能的贝叶斯网络的集合；

4)统计分析网络特征，重构显著性的亚网络模块。

所述的基因芯片表达谱缺失数据值估计：利用Least Local Squares算法，通过计算Euclidean距离获知基因表达量是否相似，缺失的基因表达量数据由表达情况与其最相似的k个基因的表达量数据通过线性组合得到。

利用贝叶斯网络构建基因间的调控网络：选取表达量有差异的基因，用它们来构建贝叶斯网络；利用Gaussian-inverse Gamma分布来处理连续变量，即基因表达量数据；在参数独立的假设下，从训练数据和用户定义的网络中得到一个初始的贝叶斯网络。

对原始的表达数据重新抽样，重复步骤2)，得到一个可能的贝叶斯网络的集合：利用随机重启的贪婪算法，得到贝叶斯分值最高的一个优化的网络；在优化的网络中，得到基因之间的条件独立关系；如果两个基因之间有直接的相互作用，代表它们的两个结点会由一条带有方向的边连接；如果他们存在v-structure中，那么他们之间有调控关系；如果不在一个v-structure中，那么这两个基因作为结合物行使生物学功能；用Efron非参数重新抽样法对芯片数据进行有放回的重新抽样，然后利用每次抽样的数据构建多个优化的贝叶斯网络；通过选取那些置信度高于阈值的边，得到较为准确的基因间的相互关系；将重新抽样法得到的一阶网络构建成二阶网络。

构建成二阶网络：采取构建二阶图的方法将二阶图中的一个点代表一阶图中的一条边，二阶图中的一条边代表了一阶图中两条边的同时出现。

统计分析网络特征，重构显著性的亚网络模块：运用扩充CODENSE算法计算有向网络xCODENSE内部调用HCS，HCS是一个利用聚类来挖掘具有高连接度的亚网络的算法。扩充CODENSE算法：xCODENSE构建有向的概括子图，作为下一步挖掘子图的候选图，二阶图不变，因为它是代表一阶图中的连接关系，然后通过概括图不对称的领接矩阵，将二阶图转换为一阶有向子图。

本发明实现从大规模基因芯片实验数据中获取多个基因间调控关系，这个结果是传统生物学实验无法得到的；本方法综合了贝叶斯学习和图论两方面的证据来探讨生物节点(基因/蛋白质)之间的调控关系。它在一定程度上弥补了芯片数据不足而导致贝页斯网络学习噪音比较大的缺点。通过构建小的亚网络，我们可以在小范围内对其行驶的生物学功能进行更精确的探讨。在一定程度上来说，紧密调控基因之间的关系可以通过一致图的方法来对那些不是直接调控的基因间关系去噪。这样得到的亚网络为下步的生物学实验提供更好指导的可能性。

附图说明

图1是一个简单的贝叶斯网络；

图2是重构一致图的流程；

图3是处理基因芯片数据的总体流程图；

图4是构建显著亚网络的流程图；

图5是由原始数据构成的一个最优的贝叶斯网络；

图6是重新抽样得到的置信度大于t(t＝0.75)的概括图；

图7是置信度对应边的数目对相应得置信度作图；

图8是对产生的贝叶斯网络进行一致亚网络的挖掘。

具体实施方式

由大规模基因芯片表达谱数据构建基因调控亚网络的方法包括如下步骤：

1.基因芯片表达谱缺失数据值估计；

2.利用贝叶斯网络构建基因间的调控网络；

3.对原始的表达数据重新抽样，重复步骤2，得到一个可能的贝叶斯网络的集合；

4.统计分本析网络特征，重构显著性的亚网络模块。

发明技术方案

1贝叶斯网络

1.1学习贝叶斯网络

贝叶斯网络是定义在一个在集台 $\stackrel{\→}{X}=\{X_1,...,X_n\}$ 域上节点间的条件独立关系，又称为马尔可夫独立性。贝叶斯网络为一个有向无环图(DAG)G，G的定点对应随机变量X₁，...X_n，以及刻画给定变量双亲的每个变量的条件分布参数。用贝叶斯网络来描述基因表达，优势在于贝叶斯网络不仅能反映基因间的依赖关系，还能反映基因间的调控关系。一个简单的贝叶斯网络如图1。

通常，我们考虑的节点有离散型和连续型两种分布。对于不同类型的节点(变量)，我们分别假设它们具有多项分布和正态分布，相应的，它们的参数分别取Dirichlet分布和Inverse Gamma分布。再又参数独立性和参数模块性的假设，我们可以由一个样本数据库来学习贝叶斯网络，通过对网络的打分搜索得到与数据吻合最好的一个或一个等价类的贝叶斯网络G。

1.2贝叶斯网络的搜索

在搜索得分高的贝叶斯网络时，我们理论上能计算所有可能DAG的得分，然后选择最高得分的DAG。但是所有可能的DAG数目随着节点的增长呈超指数分布，一般情况下，搜索最高得分的网络是NP-完全的。因此，遍历所有可能的网络是不可行的，我们需要使用其它方法，如随机重启的贪婪搜索算法。

2基于贝叶斯网络的基因调控网络的重建

2.1特征置信度

我们假设网络G正确地表示了某个域的相互作用。这个假设是否合理呢？如果我们有足够大量的取样，我们能很确定我们学习的这个网络是吻合数据的好的模型。然而，通常对于DNA芯片数据来说，我们只有少量的训练实例，可能有许多模型可以很好解释这些数据。但是这些模型又可以拥有性质不同的结构。我们不能确定某一个网络是这个生物学领域的准确描述。

因此，我们需要检验由这些给定数据产生的网络特征的后验概率。形式上，我们考虑网络特征的分布。一个网络的特征是其的一个性质，比如：X→Y在网络中，或者在网络中

d-分割X和Y等。我们用一个指示函数f(G)定义这种特征。

f(G)＝1当G满足这个特征，否则f(G)＝0。一个特征的后验概率为：

$P\left(f\left(G\right)\right|D)=\underset{G}{\mathrm{\Σ}}f\left(G\right)P\left(G\right|D)$

这个概率反应了特征f的置信度。

一个计算等式简单的方法是枚举所有得分较高的网络。不幸的是，这些网络的数目是网络节点数目的指数方，因此精确地计算这个后验概率是不实际的。但是我们可以通过抽样代表性的网络来估计这个后验概率，然后并且再去估计含有此特征的这部分代表性网络所占的比例。我们能通过非参数重新抽样的方法产生这样的网络，或者利用更精确但是耗时的MCMC模拟方法。Friedman等利用和基因表达数据分布吻合的模拟数据评估了重新抽样的方法。他们发现假阴性的比率很高(Friedman，1999)。因此，如果我们不能检测出一个特征的高置信度，我们不能认为这个特征不存在，而只能认为这些数据并不是有力的支持这个特征。

2.2显著亚网络的重构

由上述模型训练出来的最优贝叶斯网络在适当的模型规模，数据足够大时，可以很好的解释实际现象。但是在DNA芯片实验中，这种情况往往得不到满足。到目前为止，可以做的实验数小于要建模的模型大小，或者实验数据不足以训练出可信度很高的模型。这样，就要求我们利用有限的数据挖掘出比较可靠的一些网络特征。通常，生物网络是以一个个模块为调控单位，这就要求我们对已经建模的贝叶斯候选网络重构显著性的亚网络。Zhou(2005)引入了通过从每一个1阶图取出它们的元信息来构建2阶图。归结重构一致图的流程如图2。处理缺失数据：

利用Least Local Squares算法估计缺失数据。通过计算Euclidean距离可以获知基因表达量是否相似，缺失的基因表达量数据可以由表达情况与其最相似的k个基因的表达量数据通过线性组合得到。例如，在对g1基因进行的5次芯片实验中，表达量数据有一个缺失值，我们选取k个表达量数据完整的并且与g1基因表达情况最接近的基因，然后按下图所示构建矩阵A，向量b和ω，以及缺失信α。

$\left(\begin{array}{cc}a& {\mathrm{\ω}}^T\\ b& A\end{array}\right)=\left(\begin{array}{ccccc}\mathrm{\α}& {\mathrm{\ω}}_1& {\mathrm{\ω}}_2& {\mathrm{\ω}}_3& {\mathrm{\ω}}_4\\ b_1& A_{\mathrm{1,1}}& A_{\mathrm{1,2}}& A_{\mathrm{1,3}}& A_{\mathrm{1,4}}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ b_k& A_{k,1}& A_{k,2}& A_{k,3}& A_{k,4}\end{array}\right)$

α是g1基因表达量数据在第一次实验中的缺失值，ω代表可获取的g1基因4个表达量数据。b中的元素是k个基因中第一次实验的数据，α由b₁，b₂…b_k线性求和得到。

构建贝叶斯网络：

我们选取表达量有差异的基因，用它们来构建贝叶斯网络。我们利用Gaussian-inverse Gamma分布来处理连续变量，即基因表达量数据。(在BNArray中，我们认为基因表达量符合多元正态分布)。在参数独立的假设下，我们从训练数据和用户定义的网络中得到一个初始的贝叶斯网络。

为了避免局部最优，我们再利用随即重启的贪婪算法，得到贝叶斯分值最高的一个优化的网络。最后得到的网络中，我们可以得到基因之间的条件独立关系。如果两个基因之间有直接的相互作用，代表它们的两个结点会由一条带有方向的边连接。如果他们存在v-structure中，那么他们之间有调控关系；如果不在一个V-structure中，那么我们认为这两个基因作为结合物行使生物学功能。对数据进行重新抽样构建多个贝叶斯网络：

为了重复利用我们有限的实验数据，我们用Efron非参数重新抽样法对芯片数据进行有放回的重新抽样，然后利用每次抽样的数据构建多个优化的贝叶斯网络。这就提供了一个有效的检验网络中边的可靠性的方法：我们可以知道两个基因间的相互作用关系或者调控关系是否能用计算学方法加以保证。通过选取那些置信度高于阈值的边，我们得到较为准确的基因间的相互关系。此外，我们将重新抽样法得到的一阶网络构建成二阶网络，二阶图可以为我们提供更多关于网络的信息。

具有显著性、一致性亚网络的挖掘：

我们扩充了CODENSE算法，将它运用于有向网络。xCODENSE内部调用HCS。HCS是一个利用聚类来挖掘具有高连接度的亚网络的算法。

从做bootstrap训练出来的一个最优网络集合中，首先取出它们共有的网络特征，即概括子图。在这一步中，我们先计算每个可能边的支持度，然后可以通过基于连接度的算法来挖掘稠密的边的集合。这些稠密边的集合作为下一步提取一致子网络的特征的候选一阶图。在提取原网络集合的概括图后(可能大于一个)，我们将起转化为二阶的元信息图。在二阶的元信息图中，一个节点代表一阶图中一条支持度大于阈值的边，二阶图中一条边表示一阶图中两条高支持度并且相连接的两条边。通过对二阶图的挖掘(同样运用挖掘一阶稠密子图的方法)，我们可以得到二阶的稠密子图(可能大于一个)。这些子图中的边，代表一阶图中稠密一致的一个一阶子图，最后，通过将二阶图转化为一阶子图，我们可以得到一个(一组)一致的亚网络。xCODENSE构建有向的概括子图，作为下一步挖掘子图的候选图。二阶图不变，因为它是代表一阶图中的连接关系。然后通过概括图不对称的领接矩阵，我们可以将二阶图转换为一阶有向子图。

从那些高置信度的边中我们可以得到经常出现的亚网络，这为我们提供了可能的基因间的调控关系。但是，生物网络往往以“模块”出现，就是说模块中的边是同时出现或同时不出现的。因此，我们期望从那些一阶贝叶斯网络中挖掘出具有显著性、一致性的亚网络。我们采取构建二阶图的方法：二阶图中的一个点代表一阶图中的一条边，二阶图中的一条边代表了一阶图中两条边的同时出现。由于定义了阈值，我们相信在一致性亚网络中出现的边是高度可信的。

应用实例：

对已有模式生物(酵母)芯片数据进行分析。采用Spellman等人做的关于酿酒酵母细胞周期的表达谱实验。具体数据可参见http://genome.www-standford.edu的Saccharomyces Genome database。我们实验选取了包含4个时间序列的78张酵母细胞周期的mRNA表达微阵列。其中每个微阵列含有6179个基因。通过聚类分析，我们选取其中799个差异表达的基因做贝叶斯网络的建模分析。

根据已有的芯片注释，我们以参与DNA损伤修复有关的所有基因构建贝叶斯网络，并对构建出的网络特征进行统计分析，得到显著性的亚调控模块。S.cerevisiae中参与DNA损伤有关的基因共有17个，由原始的78张芯片表达数据，我们构建了这17个基因间的调控网络，如图5。

从图5我们可以看出该贝叶斯网络趋向于完全图。这是因为我们的训练数据集相对于我们的模型比较小，不能很好地除去置信度不高的边。故通过定义一个置信度阈值，做重新抽样，我们可以得到密度超过预定阈值的概括图，如图6。从图6我们可以看出，许多置信度较低的边，通过网络的特征分析被除去。此图中，我们可以看出某些基因可以作为其他许多基因的“控制基因”，如YKL113C。我们把这些马尔科夫特征对依照它们的置信度排序，取出最显著的16个网络特征，即基因之间边指向关系，结果如表1。

表1.边特征最显著的前16个边及它们的相关结点

开始基因	结束基因	置信度
			YDR097CYKL113CYDR097CYKL113CYILL066CYDR097CYKL113CYDR097CYML060WYML061CYLR383WYDR097C	YDL101CYDL101CYER095WYER095WYGL163CYIL066CYIL066CYLR032WYLR032WYLR032WYML021CYNL082W	1.01.01.01.01.00.90.90.90.90.90.90.9

表2.由重构的调控亚模块来确定的基因间调控关系

开始基因	结束基因
		YML061CYML061CYML061CYML061CYML061CYML061CYML061CYML060WYML060WYLR032WYLR032WYLR032W	YDL101CYER095WYGL021WYGL163CYIL066CYKL113CYLR032WYDL101CYDR097CYDL101CYDR097CYER095W

YKL113CYDR097CYNL312WYKL113CYDR097C

YNL082WYOL090WYOL090WYOR033CYDL101C

0.90.90.90.90.9

YLR032WYLR032WYKL113CYKL113CYKL113CYKL113CYIL066CYIL066CYIL066CYGL163CYGL021WYER095W

YGL163CYIL066CYDR097CYGL021WYGL163CYIL066CYGL163CYER095WYDL101CYGL021WYER095WYDL101C

同时，我们对网络特征的置信度和它们所对应的边的数目作图，如图7。从图7我们可以看出随着置信度的上升，相应边的个数减少。其中，当置信度为0-0.3时，下降最为剧烈。可见，我们做出来的贝叶斯网络大部分的边可信度是不高的，我们应该取出其中置信度很高边的网络做下一步分析。

注意到边的置信度很高不一定代表参与这些边的所有结点是一致出现或不出现的。这样，我们需要对上面做出来的稠密概括图做显著性的一致性亚图的挖掘。利用BNArray，我们对上面数据进一步重构显著性的亚网络，结果如图8，网络中边的特征如表2

从图8中我们可以看出，参与DNA repair的17个基因中有9个基因形成一个调控的亚模块。注意到，在亚模块中出现的边不一定要在概括图或者每一张图中出现，它是对可能贝叶斯网络信息的重新挖掘得到的。因此，我们可以结合概括图的置信度分析(表1)和挖掘的亚网络(表2)来分析预测基因间的调控关系，结果如表3。可见两种网络特征结合起来分析，可以出去很多小于阈值的网络特征边。

表3.结合置信度分析和一致性亚网络分析，得到2种方法一致结果的网络特征的基因

开始基因	结束基因	置信度
			YIL066CYKL113CYML061CYLR032W	YGL163CYIL066CYLR032WYDR097C	1.00.90.90.9(方向逆转)

由上述分析，我们可以看出结合两种分析方法得到的结果正确率很高，但同时得到的结果是很苛刻的，即很多不确定的边的关系被出去，这就意味着，没有在结果里出现的边不一定是不存在的。也就是说，假阴性的比例比较高。

Claims (7)

1.一种由大规模基因芯片表达谱数据构建基因调控亚网络的方法，其特征在于，包括如下步骤：

1)基因芯片表达谱缺失数据值估计；

2)利用贝叶斯网络构建基因间的调控网络；

3)对原始的表达数据重新抽样，重复步骤2)，得到一个可能的贝叶斯网络的集合；

4)统计分析网络特征，重构显著性的亚网络模块。

2.如权利要求1所述的一种由大规模基因芯片表达谱数据构建基因调控亚网络的方法，其特征在于，所述的基因芯片表达谱缺失数据值估计：利用LeastLocal Squares算法，通过计算Euclidean距离获知基因表达量是否相似，缺失的基因表达量数据由表达情况与其最相似的k个基因的表达量数据通过线性组合得到。

3.如权利要求1所述的一种由大规模基因芯片表达谱数据构建基因调控亚网络的方法，其特征在于，所述的利用贝叶斯网络构建基因间的调控网络：选取表达量有差异的基因，用它们来构建贝叶斯网络；利用Gaussian-inverseGamma分布来处理连续变量，即基因表达量数据；在参数独立的假设下，从训练数据和用户定义的网络中得到一个初始的贝叶斯网络。

4.如权利要求1所述的一种由大规模基因芯片表达谱数据构建基因调控亚网络的方法，其特征在于，所述的对原始的表达数据重新抽样，重复步骤2)，得到一个可能的贝叶斯网络的集合：利用随机重启的贪婪算法，得到贝叶斯分值最高的一个优化的网络；在优化的网络中，得到基因之间的条件独立关系；如果两个基因之间有直接的相互作用，代表它们的两个结点会由一条带有方向的边连接；如果他们存在v-structure中，那么他们之间有调控关系；如果不在一个v-structure中，那么这两个基因作为结合物行使生物学功能；用Efron非参数重新抽样法对芯片数据进行有放回的重新抽样，然后利用每次抽样的数据构建多个优化的贝叶斯网络；通过选取那些置信度高于阈值的边，得到较为准确的基因间的相互关系；将重新抽样法得到的一阶网络构建成二阶网络。

5.如权利要求1所述的一种由大规模基因芯片表达谱数据构建基因调控亚网络的方法，其特征在于，所述的构建成二阶网络：采取构建二阶图的方法将二阶图中的一个点代表一阶图中的一条边，二阶图中的一条边代表了一阶图中两条边的同时出现。

6.如权利要求1所述的一种由大规模基因芯片表达谱数据构建基因调控亚网络的方法，其特征在于，所述的统计分析网络特征，重构显著性的亚网络模块：运用扩充CODENSE算法计算有向网络xCODENSE内部调用HCS，HCS是一个利用聚类来挖掘具有高连接度的亚网络的算法。

7.如权利要求6所述的一种由大规模基因芯片表达谱数据构建基因调控亚网络的方法，其特征在于，所述的扩充CODENSE算法：xCODENSE构建有向的概括子图，作为下一步挖掘子图的候选图，二阶图不变，因为它是代表一阶图中的连接关系，然后通过概括图不对称的领接矩阵，将二阶图转换为一阶有向子图。

Images