基于模型的聚类法与贝叶斯网络法的基因芯片数据分析法

摘要

本发明公开了一种结合基于模型的聚类法与贝叶斯网络法的基因芯片数据分析方法。包括如下步骤：1)利用基于模型的聚类方法，根据基因在不同芯片实验中的表达情况将差异表达的基因进行聚类；2)将聚类得到的结果与Gene Ontology注释进行比较，观察属于同一类的基因的Gene Ontology注释是否相似；3)对于在聚类中属于同一类的基因，利用贝叶斯网络对同一类的基因之间的调控关系建模，确定基因间调控网络结构。本发明克服了聚类算法不够精确的缺点以及贝叶斯网络算法速度慢的不足，并且在芯片实验数量相对基因个数显得较少的情况下采用这个方法构建基因间调控网络还是比较适合的。

说明书

技术领域

本发明涉及一种基于模型的聚类法与贝叶斯网络法的基因芯片数据分析法。

背景技术

DNA微阵列技术的出现使得在全基因组范围内观察基因的表达情况成为可能。研究人员通过微阵列杂交实验获取了大量的基因表达量的数据，这种大规模的数据为我们进一步了解生物体的基本运行机制提供了可能：从基因的调控表达，基因的功能到细胞机制；与此同时，我们需要开发新的分析方法来分析这一数据，从数据中获取信息。怎样由海量的表达轮廓挖掘出具有生物学意义的基因间相互作用的调控网络，成为目前生物信息学一个很重要的领域。如何处理大规模微阵列数据从而获取基因间调控关系网络正成为一个值得关注的课题。

在对基因芯片表达谱数据的分析方面，我们能够在获得多次实验数据的基础上，通过聚类的方法将具有相似表达特征的基因归为一组(Spellman，1998)。用这种方法可以发现被共同调控或者具有相似功能的基因。许多传统聚类方法已经被运用于基因谱表达数据的分析，例如分级聚类、自组织图、k-means以及支持向量机等。由于这些方法大多是采用随机发生的算法，在聚类时无法严格地选择正确组分数以及“好”的聚类方法。基于概率的聚类模型是一个理想的用来代替随机算法的模型。基于模型的聚类认为，数据是由有限混合概率模型产生的(如多元正态分布)。这样，对于聚类的组分数以及聚类方法的选择就衍变为对统计模型的选择，通过贝叶斯信息量(BIC)这一评价标准，我们可以在聚类时选择最优的组分数以及聚类模型。

此外，由于微阵列实验的数量和其本身模型的大小相比极其不足，聚类并不能为我们提供更深层次的信息，包括基因间关系的精细结构，基因间调控是否直接，还是通过其他基因的调控。因此，要构建合理，对生物学家有更深刻知道的代谢调控网络，我们通常采取复杂随机过程的概率模型分析基因间的依赖关系，为我们进一步推测它们之间的因果关系提供可能。在这里，我们利用贝叶斯网络构建基因间调控关系，采用这一算法有以下好处：贝叶斯网络的统计学基础以及学习算法比较成熟，并且在多个领域成功的被应用；贝叶斯网络适合在数据样本比较少，而模型比较大时的学习，这在目前的基因芯片表达谱数据分析中尤为适用(Friedman，2000)；贝叶斯网络不仅可以用来推测依赖关系，还可以用来推测因果关系(即基因间调控关系)。

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

Heckerman(1995)等人采用BDe函数评判每一个候选的贝叶斯网络建模。经过对每一个候选网络进行启发式搜索，最后在整个搜索空间中得到一个或一个等价类的最优化贝叶斯网络。对于得到的一组候选网络，采用重新抽样的方法对得到的贝叶斯网络进行特征置信度的分析，从而可以得到一个置信度较高的特征集合。

我们注意到，如果单纯利用聚类的方法分析芯片数据，运算速度较快，但只能粗略地将基因分类，无法得到基因间精确的调控关系；而单纯利用贝叶斯的方法对整个基因组的基因做调控网络分析非常耗时，并且我们只关心某些关系较密切的基因之间的调控关系，没有必要对全部基因进行调控关系分析。所以，我们决定先根据基因在不同芯片实验中的表达情况，利用基于模型的聚类方法将基因进行聚类，再利用贝叶斯网络这一基于复杂随机过程的概率模型来分析在聚类中属于同一类的基因，构建基因间调控网络。Gene Ontology(GO)数据库包含了对基因产物所参加的生物过程，细胞组成部分以及分子功能这三方面的文字描述，GO提供的描述是与物种无关的。GO已被应用于相当多的基因组数据库的注释中，如：SGD，CGD，FlyBase，MGI，TAIR，ZFIN，DictyBase，WormBase和RGD。

需要注意的是GO的结构是一个DAG(Directed Acylic Graph)。我们得到的关于某个基因的注释信息在DAG中的位置可能位于DAG的分枝上端也可能位于分枝末端：我们对这个基因的了解程度越深入，那么注释信息越靠近分枝末端；并且，GO注释并不完善，特别是它的关键词库必须根据其他数据库中信息的更新而进行相应的更新，因此，对于新序列的基因注释能力要视其他数据库对该序列及其同源序列的注释情况而定。由于DNA数据库中的信息要远多于已被注释了基因的信息，所以目前，用户通过GO查询不到任何与自己感兴趣的基因有关的信息也是很有可能的；此外，GO的注释是通过生物学实验或者BLAST等方法得到，用这些方法得到的注释也有可能包含错误信息。所以在我们的本发明中，我们观察聚类结果与GO注释之间的相关性，而不是将GO作为评判聚类结果好坏的标准。

我们采用贝叶斯网络的方法来获取基因间调控关系。对DNA微阵列数据的分析和调控网络的构建主要可以分为两大类：基于传统统计学的聚类分析(Quackenbush，2001)和基于复杂随机过程的概率模型(Pearl，1988)，如贝叶斯网络。

利用贝叶斯网络分析基因表达具有以下好处：1)贝叶斯网络的统计学基础以及学习算法比较成熟，并且在多个领域成功的被应用；2)贝叶斯网络适合在数据样本比较少，而模型比较大时的学习，这在目前的DNA微阵列数据分析中尤为适用(Friedman，2000)；3)贝叶斯网络除了可以用来推测依赖关系，还可以用来推测因果关系。

通过机器学习实现从大规模芯片实验数据获取多个基因间调控关系，这个结果是传统生物学试验无法得到的。此外，我们的结果也可为今后的生物学实验提供理论指导。

我们提出了一整套达到这个目的的方法。具体说，我们利用基于模型的聚类方法，根据基因在不同芯片实验中的表达情况将它们进行聚类；将聚类得到的结果与GO注释进行比较，观察属于同一类的基因的GO注释是否相似；最后，对于在聚类中属于同一类的基因，我们利用贝叶斯网络这一基于复杂随机过程的概率模型来分析它们，并且构建基因间调控网络。我们用这套方法处理了一组酵母芯片数据，并从结果中得到一些启示。

发明内容

本发明的目的是提供一种结合基于模型的聚类法与贝叶斯网络法的基因芯片数据分析方法。

包括如下步骤：

1)利用基于模型的聚类方法，根据基因在不同芯片实验中的表达情况将差异表达的基因进行聚类；

2)将聚类得到的结果与Gene Ontology注释进行比较，观察属于同一类的基因的Gene Ontology注释是否相似；

3)对于在聚类中属于同一类的基因，利用贝叶斯网络这一基于复杂的概率模型对同一类的基因之间的调控关系建模，确定基因间调控网络结构。

所述的利用基于模型的聚类方法，根据基因在不同芯片实验中的表达情况将它们进行聚类：基于模型的聚类假设每个基因的表达谱是独立同分布的随机变量，其分布是由多个不同的分布函数混合而成，每个基因所属的类别由其属于某个类别的后验分布概率决定，每个基因聚类结果的类标被看成不能直接观察到的隐变量，这些隐变量通过Expectation Maximization算法来估计，聚类最优的类别个数的选择通过模型选择的贝叶斯信息量标准决定。

所述的将聚类得到的结果与Gene Ontology注释进行比较，观察属于同一类的基因的Gene Ontology注释是否相似：Gene Ontology数据库包含了对基因产物所参加的生物过程，细胞组成部分以及分子功能这三方面的文字描述，基因表达谱相似的基因通常属于同一生物学通路，且在生物过程和分子功能上会有相似处，利用Perl语言编程从网页上提取了酵母基因的Gene Ontology注释，然后比较在聚类时被归如同一类的基因的注释信息是否相同或者有上，下级关系，以此验证聚类结果。

所述的对于在聚类中属于同一类的基因，利用贝叶斯网络这一基于复杂的概率模型对同一类的基因之间的调控关系建模，学习基因间调控网络结构：首先，对于每个网络结构，基于所做的模型假设，得到数据和每个网络联合分布的封闭表达式，这个数值正比于这个网络的后验分布概率，利用这个表达式来打分评判所有可能网络的后验概率，其次，搜索网络拓扑空间，采用爬山算法搜索得到打分较高的网络，寻找到一个局部最优解，通过随机多次重启，得到多个局部最优解，确定最优的网络，最后，通过用马尔可夫链蒙特卡罗方法，来抽样模拟网络结构的后验分布，来比较爬山算法的搜索结果。

本发明实现从大规模基因芯片实验数据中获取多个基因间调控关系，这个结果是传统生物学实验无法得到的；本方法结合基于模型的聚类法与贝叶斯网络法的基因芯片数据分析方法来探讨生物节点(基因/蛋白质)之间的调控关系，克服了聚类算法不够精确的缺点以及贝叶斯网络算法速度慢的不足，并且在芯片实验数量相对基因个数显得较少的情况下采用这个方法构建基因间调控网络是比较适合的。同时，本方法将两种算法相结合，避免了一些基因间的调控关系很可能由于采用了聚类算法而丢失。本方法为开发更具生物学意义芯片数据分析算法显得尤为重要。

附图说明

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

图2是DAG G的网络结构示意图；

图3是DAG G*的网络结构示意图；

图4是DAG G*的网络结构示意图；

图5是标准化之前对YMR095C和YMR096W基因在300次实验中的表达谱线图；

图6是标准化之后对YMR095C和YMR096W基因在300次实验中的表达谱线图；

图7是采用各种聚类个数和聚类方法时的BIC值示意图；

图8是在20个类中基因归属的不确定性示意图；

图9是类的大小类中基因归属不确定性的关系示意图；

图10是对含有9个基因的类中各基因在300次实验中的表达情况示意图；

图11是对含有20个基因的类中各基因在300次实验中的表达情况示意图；

图12是YDL248W的MF的树形图注释图；

图13是YML132W的MF的树形图注释图；

图14是MCMC抽样接收率示意图；

图15是用贝叶斯网络法分析含有20个基因的类，得到基因之间调控关系网络图(每个节点代表一个基因，带有方向的箭头表示基因之间的调控关系)；

图16是A基因和B基因结合后共同调控C基因示意图；

图17是将YOL106W和YPR065W的表达情况标准化之后示意图；

图18是将YOL106W和YPR065W的表达情况标准化并取绝对值后示意图。

具体实施方式

结合基于模型的聚类法与贝叶斯网络法的基因芯片数据分析方法，包括如下步骤：

1)利用基于模型的聚类方法，根据基因在不同芯片实验中的表达情况将差异表达的基因进行聚类；

2)将聚类得到的结果与Gene Ontology注释进行比较，观察属于同一类的基因的Gene Ontology注释是否相似；

3)对于在聚类中属于同一类的基因，利用贝叶斯网络这一基于复杂的概率模型对同一类的基因之间的调控关系建模，确定基因间调控网络结构。

发明技术方案

1.基于模型的聚类

基于模型的聚类假设每个基因的表达谱是独立同分布的随机变量，其分布是由多个不同的分布函数混合而成，每个基因所属的类别由其属于某个类别的后验分布概率决定，每个基因聚类结果的类标被看成不能直接观察到的隐变量，这些隐变量通过Expectation Maximization算法来估计，聚类最优的类别个数的选择通过模型选择的贝叶斯信息量标准决定。

1.1基于模型聚类方法的优势

许多传统聚类方法已经被运用于基因谱表达数据的分析，例如分级聚类、自组织图、k-means以及支持向量机等。在进行基于距离的聚类时，我们首先定义一个相似度的度量。通常，我们选择两种距离：欧拉距离和Pearson相关系数。在定义了一个相似度之后，我们可以使用许多聚类的算法，如：完全连接聚类，最小连接聚类，平均连接聚类；以及权重成对分组聚类等等(Sokal，1963)。大多数聚类分析方法是阶梯的，聚类的结果象一个系统发育进化树。非阶梯聚类算法同样存在，比如k-means法，这种方法把所有的对象归结到几个预先定义好的分区中。许多传统聚类方法已经被运用于基因谱表达数据的分析，但是问题在于，由于这些方法大多是采用随机发生的算法，他们在聚类时无法严格地选择正确组分数以及“好”的聚类方法。这些方法的运用取得了一定的成功，但是至今没有一个方法被认为是处理基因谱数据的最好方法。

基于概率的聚类模型是一个理想的用来代替随机算法的模型。特别的，基于模型的聚类认为，数据是由有限混合概率模型产生的(如多元正态分布)。这样，对于聚类的组分数以及聚类方法的选择就衍变为对统计模型的选择(Dasguptaand Raftery，1998；Fraly and Raftery，1998)。

1.2基于模型聚类的框架

在基于模型的聚类中，我们认为，每个观察值x都是p维的向量，它是从一个有限的混合分布中得到的。

$L_{\mathrm{mix}}({\theta }_1,...,{\theta }_G;{\pi }_1,...,{\pi }_G|x)=\underset{j=1}{\overset{n}{\Pi }}\underset{i=1}{\overset{G}{\Sigma }}{\pi }_i{f}_i\left(x_j\right|{\theta }_i)$

其中π_i是先验概率，f_i是每个组分所特有的参数为θ_i的概率分布。混合分布中的每个组分代表了一个类。类的个数g需在对数据进行分析时确定。

特别的，有限正态混合分布被广泛应用，即每个组分f_i都是多元正态分布，我们在本文中就是应用这个模型。我们有：$f_i(x;{\theta }_i)=\frac{1}{{\left(2\pi \right)}^{p/2}{\left|V\right|}^{1/2}}\exp {(-\frac{1}{2}(x-{\mu }_i)}^{\prime }{V}^{-1}(x-{\mu }_i)),$其中$V=\mathrm{diag}({\sigma }_1^2,{\sigma }_2^2,...{\sigma }_p^2),$$\left|V\right|={\Pi }_{q=1}^p{\sigma }_p^2.$

当得到一个数据集x_j，j＝1，...n，在基于模型的聚类中，EM算法(Dempster etal.，1977)被反复利用从而估计参数^Θ(McLachlan and Peel，2002；Fraley andRaftery，2002)。在EM算法中，我们认为每个变量x都是从一个组分中得到的，并认为指示某个变量是属于哪个组分的指示变量是缺失的。在E步骤中，我们根据当前的参数估计指示变量，在M步骤中，我们根据更新组分的均值向量、方差协方差矩阵以及组分的权重。我们用Θ^(m)表示在第m次重复后估计所得的参数，EM的算法如下：

M-Step：

${\hat{\mu }}_i^{(m+1)}=\underset{j}{\Sigma }{\tau }_{\mathrm{ij}}^{(m+1)}{x}_j/\underset{j}{\Sigma }{\tau }_{\mathrm{ij}}^{(m+1)},$

${\hat{\sigma }}_k^{2,(m+1)}=\underset{i=1}{\overset{g}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}{\tau }_{\mathrm{ij}}^{\left(m\right)}{(x_{\mathrm{jk}}-{\mu }_{\mathrm{ik}}^{\left(m\right)})}^2/n,$

${\hat{\pi }}_i^{(m+1)}=\underset{j=1}{\overset{n}{\Sigma }}{\tau }_{\mathrm{ij}}^{\left(m\right)}/n$

E-Step：

${\tau }_{\mathrm{ij}}^{\left(m\right)}=\frac{{\pi }_i^{\left(m\right)}{f}_i(x_j;{\theta }_i^{\left(m\right)})}{{\Sigma }_{i=1}^g{\pi }_i^{\left(m\right)}{f}_i(x_j;{\theta }_i^{\left(m\right)})}$

以上的算法被重复，直至观察数据的似然度函数收敛，这样就得到了最大似然估计(MLE)。我们通过τ_ij获知x_j来自组分i的后验概率，并且在分类时把x_j分配到概率最大的那一类中。

1.3模型的选择

在应用中，我们需要决定组分的个数g。这是通过最优化一组包含不同数目组分的模型，然后选择最优模型来实现的。在基于模型的聚类中，通常使用贝叶斯信息标准(BIC)(Schwartz，1978)，定义如下：$\mathrm{BIC}=-2\log L\left(\hat{\Theta }\right)+\log \left(n\right)d,$其中d是有效参数的总个数(Fraley and Raftery，1998)。在标准模型中，由于我们有三组参数，π_i，σ_q，μ_iq，我们有限制$\underset{i=1}{\overset{g}{\Sigma }}{\pi }_i=1,$所以d＝p+gp+g-1。

2.GO注释聚类比较

Gene Ontology数据库包含了对基因产物所参加的生物过程，细胞组成部分以及分子功能这三方面的文字描述。基因表达谱相似的基因在聚类时会被归入同一类中，并且，基因表达谱相似的基因通常属于同一生物学通路，即它们在生物过程和分子功能上会有相似处。因此我们决定将聚类结果与GO注释进行比较。对于每一方面的注释信息，在GO中的结构都是以DAG(Directed AcylicGraph)形式存在，在DAG中，越靠近根部的注释信息越概括化，越靠近分支末端的注释信息越细节化。因此我们得到的关于某个基因的一个方面的GO注释信息，例如生物过程，其在DAG中的位置可能位于DAG的分枝上端也可能位于分枝末端：我们对这个基因的了解程度越深入，则注释信息越靠近分枝末端；此外，拥有相同父结点的子结点共享父结点的属性。我们利用Perl语言编程从网页HTML源文件上提取酵母基因在GO中的三方面注释，然后比较在聚类时被归入同一类的基因的注释信息是否相同或者有上，下级关系(既基因的注释信息在DAG中是否属于同一结点或者邻近结点)，以此验聚类结果。

3.基于贝叶斯网络的基因调控网络的构建

首先，对于每个网络结构，基于所做的模型假设，得到数据和每个网络联合分布的封闭表达式，这个数值正比于这个网络的后验分布概率，利用这个表达式来打分评判所有可能网络的后验概率，其次，搜索网络拓扑空间，采用爬山算法搜索得到打分较高的网络，寻找到一个局部最优解，通过随机多次重启，得到多个局部最优解，确定最优的网络，最后，通过用马尔可夫链蒙特卡罗方法，来抽样模拟网络结构的后验分布，来比较爬山算法的搜索结果。

3.1学习贝叶斯网络

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

给定从未知分布中独立抽样的训练集$D=\left\{\overrightarrow{x}\right[1],...,\overrightarrow{x}[M\left]\right\},$我们想由一个网络G估计这个分布。解决这个问题通常的方法是引入一个有统计学依据的打分函数来评估每个网络，相关联训练数据，以及通过得分来搜寻优化网络。一个常用基于贝叶斯推理的评估，对候选可能网络图G打分，通过它们的在给定数据后的后验概率。我们定义打分函数S(G:D)与P(G|D)成正比。这个打分函数的重要特征是当数据完全时(没有丢失值)，打分函数是可以分解的：

$S(G:D)=\underset{i}{\Sigma }{S}_{\mathrm{local}}(X_i,{\mathrm{Pa}}_i^G:D)$

每个变量X_i对整个分数的贡献仅仅与X_i和Pa_i^G有关。

$S_{\mathrm{local}}(X_i,U:D)=\log P({\mathrm{Pa}}_i=U)+\log \int \underset{m}{\Pi }P\left(X_i\right[m\left]\right|U_i\left[m\right],\theta )\mathrm{dP}\left(\theta \right)$

第一项是选择集合U作为X_i双亲的先验概率。第二项衡量数据的可能性，我们对所有可能的条件分布参数θ积分。这些每个变量的局部贡献可以用一个封闭式等式计算。

对于离散型变量，我们做如下几个假设：

1)所有数据都是从多项分布中抽样产生(Multinomial Sample)

2)参数符合Dirichlet分布

3)统计学习的数据是完全的

由上面3个假设，我们可以得到一个完全图的离散型网络的BD(BayesianDirichlet)的分数，表达式如下：

$S(G:D)=\log P({\mathrm{Pa}}_i=U)+\underset{i=1}{\overset{n}{\Sigma }}\underset{j=1}{\overset{q_i}{\Sigma }}\log (\frac{\Gamma \left(N_{\mathrm{ij}}^{\prime }\right)}{\Gamma (N_{\mathrm{ij}}^{\prime }+N_{\mathrm{ij}})}+\underset{k=1}{\overset{r_i}{\Sigma }}\log \left(\frac{\Gamma N_{\mathrm{ijk}}^{\prime }+N_{\mathrm{ijk}}}{\Gamma \left(N_{\mathrm{ijk}}^{\prime }\right)}\right))$

再由假设：

4)参数的独立性

5)参数的模块性

我们可以由所有完全图的等价类得到一个完全的连续型网络等价类的BDe(Bayesian Dirichlet equivalent)的分数。

对于连续型变量，我们做如下几个假设：

1)所有数据都是从以未知均值和未知精确度矩阵W的多元正态分布中产生

2)参数的先验分布为正态-威沙特(normal-Wishart)分布

3)所有数据是完全的

由上面3个假设，我们可以得到一个完全的连续型网络的BG(BayesianGaussian)的分数，表达式如下：

$\rho \left(D\right|B_{s_c}^e,\xi )={\left(2\pi \right)}^{-\frac{\mathrm{nm}}{2}}{\left(\frac{v}{v+m}\right)}^{\frac{n}{2}}\frac{c(n,\alpha )}{c(n,\alpha +m)}{\left|T_0\right|}^{\frac{\alpha }{2}}{\left|T_m\right|}^{-\frac{\alpha +m}{2}}$

其中：

$c(n,\alpha )={\left[2^{\frac{\mathrm{\alpha n}}{2}}{\pi }^{\frac{n(n-1)}{4}}\underset{i=1}{\overset{n}{\Pi }}\Gamma \left(\frac{\alpha +1-i}{2}\right)\right]}^{-1}$

再由假设：

4)参数的独立性

5)参数的模块性

我们可以得到一个完全的连续型网络等价类的BGe(Bayesian Guassianequivalent)的分数，它们与完全网络的关系如下：

$\rho \left(D\right|B_s^e,\xi )=\underset{i=1}{\overset{n}{\Pi }}\frac{\rho \left(D^{x_i{\Pi }_i}\right|B_{s_c}^e,\xi )}{\rho \left(D^{{\Pi }_i}\right|B_{s_c}^e,\xi )}$

对于混合型网络，我们同样做如上几个假设，并且约定：

$\overrightarrow{x}=\Delta \cup \Gamma$

其中Δ中的每一个变量都是离散的，Γ中每一个变量都是连续的，那么我们得到：

$\rho \left(D\right|B_S^e,\xi )=p\left(D^{\Delta }\right|B_S^e,\xi )\underset{l=1}{\overset{m}{\Pi }}\rho \left(C_l^{\Gamma }\right|C_l^{\Delta },C_1^{\Gamma },...,C_{l-1}^{\Gamma },B_S^e,\xi )$

其中

$\underset{l=1}{\overset{m}{\Pi }}\rho \left(C_l^{\Gamma }\right|C_l^{\Delta },C_1^{\Gamma },...,C_{l-1}^{\Gamma },B_S^e,\xi )=\underset{i=1}{\overset{\gamma }{\Pi }}\underset{j=1}{\overset{s_i}{\Pi }}\rho \left(D^{{\Gamma }_i,j}\right|{\Delta }_i=j,B_S^e,\xi )$

3.2贝叶斯网络搜索

在搜索得分高的贝叶斯网络时，我们理论上能计算所有可能DAG的得分，然后选择最高得分的DAG。Robinson给出了计算含有n个结点的所有可能DAG数目的递归公式(Robinson，1977)，

$f\left(n\right)=\underset{i=1}{\overset{n}{\Sigma }}{(-1)}^{i+1}\left(\begin{array}{c}n\\ i\end{array}\right)2^{i(n-i)}f(n-i)$

那么在我们的一般混合型网络中(不允许离散型结点有连续型的父结点)，所有可能的混合DAG数目为：

f(|Δ|，|Γ|)＝f(|Δ|)×f(|Γ|)×2^|Δ|*|Γ|

其中f(|Δ|)和f(|Γ|)分别为离散型，连续型变量的所有可能的DAG的数目，2^|Δ|*|Γ|代表从离散到连续变量所有可能的箭头组合。

从上面我们可以看出，所有可能的DAG数目随着节点的增长呈超指数分布(表1)，一般情况下，搜索最高得分的网络是NP-完全的(表1)。

表1.含有n个结点的贝叶斯网络的可能个数。

  Number of  nodes  Number of DAGs  1  1  2  39116  3  39441  4  144-543  5  4800-29281  6  320000-3781503  7  ≈56*10E6-10E9  8  10E10-10E11  9  10E13-10E15  10  10E16-10E18

由于网络结构的可能数目随节点数增加呈超指数增长，我们不能够完全的搜索整个DAG空间，因此我们可以利用局部优化搜索方法(如：K2，多次重启的爬山算法)或者全局搜索方法(如：MCMC)。在这个项目里，我们采用全贝叶斯框架下的MCMC方法来逼近模拟DAG结构的后验分布。这种方法，如其它许多方法一样，使用贝叶斯因子来比较两个不同DAG的网络得分。

3.3贝叶斯因子

为了比较两个网络(G和G^*)的得分，我们计算验后加权值，定义如下：

$\frac{p\left(G\right|D)}{p\left(G^{*}\right|D)}=\frac{p(G,D)}{p(G^{*},D)}=\frac{p\left(G\right)}{p\left(G^{*}\right)}\times \frac{p\left(D\right|G)}{p\left(D\right|G^{*})}$

其中：

为验前加权，为贝叶斯因子。

对于只有一个边不同的两个网络，由于可分解性，贝叶斯因子的计算很简单。在只有一个边不同的条件下，又分为两种情况：

两个DAG只有一个边的方向不同

令在G中，存在边v←w，在G^*中存在边v→w，G与G^*仅有v，w之间边的方向不同。并且设w在G中的父结点为U_w，v在G^*中的父结点为Uv，那么G和G^*的网络结构如图2和图3所示。

由似然函数的可分解性，贝叶斯因子可分解为：

$\frac{p\left(D\right|G)}{p\left(D\right|G^{*})}=\frac{p\left(v\right|U_{v,}w,G)p\left(w\right|U_w,G)}{p\left(w\right|U_w,v,G^{*})p\left(v\right|U_v,G^{*})}$

$=\frac{\int p\left(x_v\right|x_{w\cup \mathrm{Uv},}{\Theta }_{v|w\cup U_v},G)p\left({\Theta }_{v|w\cup U_v}\right|G)d{\Theta }_{v|w\cup U_v}}{\int p\left(x_w\right|x_{v\cup U_{w,}}{\Theta }_{w|v\cup U_w},G^{*})p\left({\Theta }_{w|v\cup U_w}\right|G^{*})d{\Theta }_{w|v\cup U_w}}\times \frac{\int p\left(x_w\right|x_{U_w},{\Theta }_{w|\mathrm{Uw}},G)p\left({\Theta }_{w|\mathrm{Uw}}\right|G)d{\Theta }_{w|\mathrm{Uw}}}{\int p\left(x_v\right|x_{U_v},{\Theta }_{v|U_v},G^{*})p\left({\Theta }_{v|U_v}\right|G^{*})d{\Theta }_{v|U_v}}---\left(11\right)$

因此，计算两个网络G和G^*的贝叶斯因子，我们只需要考虑v，w的条件分布。

注意，当U_v＝U_w时，G和G^*为等价的网络，并且贝叶斯因子等于1。

两个DAG只有一个边的是否存在不同

令在G中，存在边v←w，在G^*中不存在边v→w。并且设w在G中的父结点为U_w，v在G^*中的父结点为Uv，那么G^*的网络结构如图4所示。

由似然函数的可分解性，贝叶斯因子可分解为：

$\frac{p\left(D\right|G)}{p\left(D\right|G^{*})}=\frac{p\left(v\right|U_{v,}w,G)}{p\left(v\right|U_v,G^{*})}$

$=\frac{\int p\left(x_v\right|x_{w\cup \mathrm{Uv},}{\Theta }_{v|w\cup U_v},G)p\left({\Theta }_{v|w\cup U_v}\right|G)d{\Theta }_{v|w\cup U_v}}{\int p\left(x_v\right|x_{U_{v,}}{\Theta }_{v|U_v},G^{*})p\left({\Theta }_{v|U_v}\right|G^{*})d{\Theta }_{v|U_v}}$

因此，对于只有一条边不同的2个贝叶斯网络，我们只需要比较它们不同边连接的2个节点的局部贝叶斯因子即可以得到2个网络的验后加权值。

4应用实例：

4.1实验数据

我们采用的数据是Rosetta Inpharmatics Compendium(Timothy R.Hughes，2000)，它是一个由300个全基因表达轮廓组成的参考数据库。这些表达轮廓从276个删除突变，11个可由四环素调节的关键等位基因，以及13个化学处理过的S.cerevisiae的培养体，每一组都与野生型或模拟培养组的基线做对比。为了方便聚类以及构建贝叶斯网络，我们先对表达数据进行缺失值的估计处理后再进一步运算。

4.2实验步骤

·从http://www.rii.com/tech/pubs/cell_hughes.htm下载得到芯片数据表格

·用KNN法补全表格中的缺失数据

·用基于模型的聚类方法将酵母基因根据表达谱分类

·用Perl提取酵母基因的GO注释信息，并与聚类结果进行比较

·通过贝叶斯网络得到属于同一类的基因的调控关系

4.3结果与讨论

4.3.1标准化每个基因的表达谱

从文献中我们得知，将基因表达量进行标准化之后更有利于对基因进行聚类分析。所谓标准化就是将一个基因在每个芯片实验中的表达量减去这个基因表达量的均值，除以标准差。这样，每个基因表达量的均值被调为0，方差被调为1。通过标准化后，那些表达量不同但是表达趋势相同或相似的基因的表达谱会更加靠近，而那些表达情况本来就十分相似的基因的表达谱也还是吻合得很好。例如，图5和图6分别是标准化之前和标准化之后对YMR095C和YMR096W两个基因在300次实验中的表达谱线作的matplot，我们发现标准化前后，这两个基因的表达谱线都吻合得很好。

我们对数据进行处理之前，也先把它进行了标准化。同时，我们发现，标准化之后再进行聚类，各模型的BIC值都大于对未标准化时得到的BIC值。

4.3.2通过BIC值选择模型

我们观察了聚类个数分别是5，10，......，80时，各种聚类方法的BIC值，如图7所示。之后的步骤中，我们采用其中BIC值最大的方法，即类的个数为20的VVI聚类法。

当模型的参数调整完毕后，对每个基因我们都能得到它由某个类产生的概率，我们将基因归入最有可能产生它的那一类中。表2显示了采用类的个数为20的VVI聚类法进行聚类时，每个类所含的基因数。

表2.每个类中所含的基因数。

  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  884  317  575  582  228  628  266  816  225  339  20  23  9  67  110  186  75  952  9  14

我们又观察了各个类中的基因属于这个类的不确定性，结果如图8所示。

进一步观察，我们发现，如果一个类越大，即类中所含基因数越多，那么这一类中基因归属的不确定性的均值越大；反之，如果一个类越小，这个类中基因归属的不确定性的均值越小。从图9我们可以发现这一趋势。

另外，我们还发现，如果一个类中的基因数较少，那么对这个类中所有基因的表达谱线作图，它们将吻合得非常好；如果一个类中的基因数较多，那么对这个类中所有基因的表达谱线作图，它们互相吻合程度就差一些。图10和图11分别是取含有9个、20个基因的类来作matplot的结果。

4.3.3将聚类结果与GO注释进行比较

我们利用Perl从网页上提取了酵母基因的GO注释，并将它与基因归类的表格关联起来。这样，我们就可以观察属于同一类的基因的GO注释是否有相关性。我们选取第11和第13个类进行观察，这个类中基因的GO注释如表3所示。

表3.属于同一类的9个基因在GO中的注释。

从表格中，我们发现，这九个基因的注释有很多空缺，所以我们对那些已有的注释进行观察和分析。

在MF中，YDL248W和YML132W的注释并不相同。我们再深入观察两个基因在注释结构上是否有上下级的关系(如图12，图13所示)。

通过对树形图的观察发现，两个基因的molecular function在GO注释中非上下分枝关系。此外，我们发现，GO的注释还标记了YML132W有receptor-associated protein activity，这说明虽然属性图上YML132W和YDL248W并非上下集关系，但它们的MF很可能有一定的相关性。

接着，我们观察YHL048W和YML132W的BP树形图，我们发现这两个基因在树形图上的注释也没有上下级关系。

最后，这九个基因中，八个有CP注释，但不尽相同。我们通过树形图观察它们的CP是否存在相关性，发现，这些基因的产物有的在细胞质内，有的在线粒体上，而有的在细胞核核膜上，这些位置并没有上下级的关系。我们认为，具有相同或相似生物学功能的基因不一定处于同一位置。

在芯片实验中表达谱相似的基因，在聚类分析时会被聚为一类。我们发现，在聚类中被分入同一类的基因，它们的GO注释可能一样，也有可能不一样。我们分析，主要有三种导致注释与分类结果间差异的原因：1)通过Perl我们只能得到GO提供的对基因最细化一个字面注释。由于人们对基因了解程度不同，一个基因的注释可能位于DAG图的上端也可能位于DAG图的末端。两个在参加的生物过程，细胞组成部分以及分子功能都非常相似的基因，可能由于人们对它们了解程度的深浅而拥有位于DAG的分枝末端或分枝上端的注释。若需要了解两个基因是否真正相关，必需细观察具体的DAG图，而不是GO提供的对基因最细化的那一个字面注释；2)GO注释不完善，即人们通过传统生物学实验对基因进行分析，可能忽略了基因的许多特性，导致注释信息与基因在细胞内的表达情况不相符。例如，两个基因虽然行使不同的生物学功能，但他们之间存在间接但是精确的调控关系，这个关系可能是通过GO注释无法得到的；或者两个基因有各自不同的生物学功能，但他们通常行使一个相同的功能，由于生物学知识的空缺，GO注释并没有体现这一相同的功能；3)GO注释错误，GO注释主要来源于过程复杂的传统生物实验和通过BLAST进行的相似性简单推断，这两种方法有可能出现错误，所以GO注释出错也是可能的。

4.3.4用贝叶斯网络法构建基因间调控关系网络

我们选取了一个含有20个基因的类，这20个基因的基因名见表4，值得一提的是这20个基因中的17个具有GO注释，并且注释内容相同。

表4.用来构建调控网络的20个基因。

  1  YAR009C  2  YBL101W-B  3  YBR012W-B  4  YCL019W  5  YER138C  6  YER160C  7  YFL002W-A  8  YHR214C-B  9  YIL060W  10  YJR027W  11  YJR029W  12  YLR035C-A  13  YLR334C  14  YML039W  15  YML045W  16  YMR045C  17  YMR050C  18  YMR158C-B  19  YOL106W  20  YPR002C-A

用贝叶斯网络法构建这20个基因间的调控关系。我们对MCMC抽样接收率的作图14。

我们可以大致看出，从第500次开始，接受率显著下降，经过7000循环过后，接收率基本稳定在40％左右，这个是MCMC理想的抽样接收率。

我们针对这20个基因构建的调控网络如图15所示。

结论

我们提出了将基于模型的聚类方法与贝叶斯网络法相结合，用以处理生物芯片数据、发现基因之间的调控关系网络。我们把这个方法用在处理Rosetta公司发布的酵母双色芯片数据上，并将实验结果与现有的基因注释进行了比较。我们的方法克服了聚类算法不够精确的缺点以及贝叶斯网络算法速度慢的不足，并且我们认为在芯片实验数量相对基因个数显得较少的情况下采用这个方法构建基因间调控网络还是比较适合的。同时，我们也注意到，如果将两种算法相结合，一些基因间的调控关系很可能由于采用了聚类算法而丢失。

总之，细胞内部的生物学过程是非常复杂的，生物芯片这项技术能很好得辅助我们研究细胞内各基因的表达情况以及它们之间的相互作用。在今后，开发更具生物学意义芯片数据分析算法显得尤为重要。

附.中英文对照表

  缩略词 英文名称 中文名称  mclust model-based cluster 基于模型的聚类  BN Bayesian network 贝叶斯网络  DAG Directed Acyclic Graph 有向无环图  BDe Bayesian Dirichlet equivalent 贝叶斯Dirichlet等价  BG Bayesian Gaussian 贝叶斯高斯  GO Gene Ontology 基因注释