应用产生自多重破坏表达文库的基因调控网络的生物发现的制作方法

文档序号:6432041阅读:1513来源:国知局
专利名称:应用产生自多重破坏表达文库的基因调控网络的生物发现的制作方法
技术领域
本发明实施方案涉及用以阐明基因之间网络关系的方法。这些方法包括布尔逻辑(boolean logic)、贝叶斯估计(Bayesian estimation)、最大似然估计(maximum likelihood analysis)、样条函数(spline function)和其它曲线拟合法,以及用于确定彼此功能相关的各组基因之间的边缘关系(edge relationship)的方法。
背景技术
从重要的基因组范围的实验数据所阐明的基因调控网络,有助于从观察已知和未知生物功能的基因之间的转录调控事件中发现新的基因功能信息和表达调控事件。
旨在阐明基因调控路径的方法先前已经被报道(1,2,3)。这些研究报道的推导网络由来自时程、细胞周期和环境变动的基因表达数据组产生(4,5)。但是,从这些数据组推导得到的控制关系不是建立在用以阐明转录相关的调控控制功能的综合试验数据上,因此受到质疑。为了严格而精确地从全新表达数据组中确定新的、复杂的基因调控网络,结合了基因组删除突变的表达试验与适当计算方法的系统、全面的策略是必需的(6-9)。
在建立推导调控网络中尤其重要的是产生推定控制关系的表达试验之间的生物相关性。产生自竞争杂交破坏试验的基因组范围的表达数据具有内部控制和定量测定一种基因存在与否对其它基因表达的直接作用的优点。选择详尽阐明控制关系的破坏试验对于产生有用的基因调控信息是有价值的。
总结本发明提供用于确定系统中基因之间或者各组基因之间的相互关系的方法,例如确定系统中的基因路径或者基因网络的方法。在一个实施方案中,通过用本发明提供的改进的布尔方法(Boolean method)、改进的贝叶斯方法(Bayesian method),或者布尔方法和贝叶斯方法结合分析基因破坏表达图谱(disruption expression profile)来构建基因网络或者基因路径。在另一个实施方案中,通过分析受到试剂如药物作用的基因表达图谱来构建基因网络或者基因路径。还有在另一个实施方案中,通过分析从基因破坏表达图谱和试剂作用表达图谱得到的基因网络或路径来构建受试剂作用的基因网络或基因路径。
总之,根据基因文库的表达图谱可以获得基因破坏表达图谱,其中各种基因各自地或者与其它基因一起被破坏,例如与相关功能的其它基因一起被破坏产生表达图谱。例如,可以选择基因文库(例如完整基因组或者根据(PCR)引物等其它标准选择的一系列基因)。无论选择怎样的基因文库,一旦选择,文库中的每个基因可以各自和/或与其它组合被破坏,产生一个文库集合,每个文库含有至少一个被破坏基因和其它未被破坏的基因。这样,如果选择100个基因,得到的文库将由至少101个不同的子文库组成,例如,一个未被破坏的基因或野生型的子文库和至少该100个基因的每个基因的一个破坏体(disruptant)子文库。因此,建立破坏基因的文库,得到各个子文库和整个文库的表达图谱。
通过向含有选择基因的系统中加入一种或多种合适的试剂和收集不同剂量和不同试剂加入时间点的基因表达图谱,获得试剂作用表达图谱。有时试剂对基因表达不起作用,而有时起抑制作用(例如降低基因表达)和另些时候试剂增加基因表达(例如,为诱导物)。
采用任何适当方法如采用微阵列,可以定量地得到基因表达图谱。在一个实施方案中,基因表达图谱被编成基因表达矩阵,如将基因分成被作用和未被作用的二元素矩阵。在另一个实施方案中,引入等价集合(equivalenceset),其中数据被归一化(normalized),从而揭示基因表达的定量关系。通过等价集合,可以建立基因间相互关联的网络信息。然后可以使用网络来确定基因之间的功能关系。然后利用推导的基因之间的功能关系来预测药物和/或生物效应。这些预测可以用如微阵列试验进行试验性检验。该过程产生所谓基因表达的“最终共同路径”,也就是说,一个基因功能的变化会产生对被改变的基因的“下游”的基因的影响。
根据本发明的实施方案,利用本发明的改进的布尔方法可以分析来自直接被作用基因的“下游”效应。根据本发明的另一个实施方案,通过应用本发明的改进的贝叶斯方法可以分析被作用基因的“上游”。本发明给贝叶斯基因网络分析提供了新方法,其中使用非线性的和/或非参数的回归模型。采用这种方法意味着不需要对因果关系进行任何先验假设,而可以推导因果关系,从而提供“上游”或“初始路径”,通过该“上游”或者“初始路径”,被药物或者处理作用的所观察的基因,将被其上游关系中的其它基因诱导从而改变其表达。本发明提供一种改进的检验数(criterion),此处术语“BNRC”用于评估基因网络的贝叶斯结构图(graph)。因此,选择使BNRC检验数最小化的基因关系。
在另一个实施方案中,提供其它的不依赖于数据中的高斯(Gaussian)或其它假定方差的新方法。而在一些实施方案中,测定数据的方差,用观察到的方差来影响BNRC检验数。利用这种具有异质误差方差(heterogeneous errorvariance)和相互作用(interaction)的非参数回归,可以优化所得到的数据的曲线拟合和预测需要多少试验来获得具有期望方差的数据。使用这些方法,可以获得具有期望精确度和可信度的并且提供上、下游效应的基因网络信息。在一些实施方案中,新的贝叶斯模型和所罚最大似然估计(penalizedmaximum likelihood estimation)(PMLE)得到相似结果。
另外还有一些实施方案中,通过一种基因表达结构图和另一种基因表达的比较分析,可以阐明这两种基因的关系。这些结构图可能是线性或者非线性的。在一些实施方案中,采用线性样条函数、B-样条函数(B-splines)、傅立叶变换(Fourier transforms)、小波变换(wavelet transforms)或者其它基本函数来表征关系。有时,适宜应用B-样条函数。
确定特定基因和系列基因对于调控其它基因是否重要是十分有用的。但是,在一些情况下,数据中的无关值(outlier)使结果解释复杂化。当无关值接近整个网络的基因群的边界时,这个问题尤其麻烦。因此,在一些实施方案中,可以阐明边界效应,在其中确定边缘强度(dege intensity)和贝叶斯因果关系方向(direction of Bayes causality)的置信度。
在另一些实施方案中,用线性样条函数确定基因表达的时程(time-course)。采用样条函数比相对不灵敏的现有技术“倍数-变换”(fold-change)能更加可靠地分析时间排序数据。
利用一种或多种上述常规方法,本发明提供有用的基因网络信息并用一些已知的酵母基因组结果进行验证。然而,这些方法可以广泛应用于任何基因网络(如“转录子组”(transcriptome))和蛋白质水平的相互作用(“蛋白质组”(proteome))。因此,利用一种或多种这样的新方法,本发明人已经确定与抗真菌治疗相关的功能性基因间的新关系。因此,应用本发明的方法,可以预测基于对基因表达的起始和最终共同路径的了解的假定治疗靶点。
附图简述用特别实施方案描述本发明,其中附

图1描述本发明基于模型的相互作用的生物发现的原理图。
附图2描述用本发明的布尔方法构建基因调控子网络模型。附图2A描述基因表达矩阵,其中整合了大量表达图谱。该矩阵的每个元素代表一列中的一个基因与一行中的一个基因的基因表达的比例。附图2B描述用于生成二元素矩阵的基因之间的二元素关系。如果基因‘G1’被删除,结果基因‘G2’的强度显著地改变,那么基因‘G1’作用于基因‘G2’。附图2C描述确定循环调控基因的邻接矩阵。如果基因‘G3’和基因‘G4’相相互作用用,形成了循环(紧密连接成分)调控。附图2D描述将基因划分成等价集合的步骤,其中从理论上将循环视为“虚拟基因”。等价集合是一组相相互作用用的基因或者作用一个离散基因群。附图2E描述虚拟基因之间的基干联系。应当选择最短路径关系来建立分级连接。附图2F描述从基干矩阵中形成的调控路径。
附图3描述本发明在酵母中的基于破坏的研究和分析的结果,其中552个节点表示包括的基因和2953条这些基因之间的假定的调控连接。
附图4描述按照酵母中的细胞功能角色(CFR)分类的本发明的转录因子调控网络模型。在该模型中,根据酵母蛋白质组数据库提供的信息,98个转录因子按细胞功能角色分组。该循环中的基因分到给定的细胞功能类别中。用彩线描述调控控制关系。颜色表示基因类别,控制关系从该基因发出。蓝色碳水化合物代谢,蓝紫色染色质/染色体结构,棕色能量产生,深绿色其它代谢,灰色DNA修复,绿色脂类、脂肪酸代谢,浅绿色氨基酸代谢,橙色细胞应激(cell stress),红色减数分裂/配对反应,粉红色分化,紫色细胞周期。加粗的红线表示减数分裂/配对反应组中的基因对细胞周期相关的转录因子的调控控制。加粗的蓝线表示碳水化合物代谢组发出的作用于脂肪酸代谢组基因的控制关系。描述多条从“碳水化合物代谢”发出的控制线。明显不同的,只有两个基因SKN7和HMS2单独地影响细胞周期相关的基因表达。
附图5描述了从破坏突变酵母菌株的基因表达试验的布尔分析中重建的转录因子的基因调控子网络的详细图。黑线表示调控关系,箭头表示表达影响的方向。节点的颜色和形状表示按照YPD中的描述的基因产物的细胞功能的常规分类。与细胞分裂机制相关的基因用三角形节点表示,与DNA修复和染色体结构相关的基因用正方形表示。所阐明的网络表明了与减数分裂、配对反应和DNA结构和修复机制相关的基因之间经由UME6和MET28的新拓扑控制关系。与减数分裂和配对反应相关的基因位于受INO2调控的级联(cascade)的下游,而且更低一级的与DNA修复和结构相关的基因分级出现在MET28的下游。
附图6描述表明用贝叶斯推导和BNRC检验数最小化来估计建立基因网络关系的本发明方法的流程图。
附图7a和7b描述用于分析药物反应表达数据的常规方法。附图7a包括酵母基因的列表,所述基因的表达水平显著受到暴露于100mg灰黄霉素1分钟的影响(20个试验之一)。附图7b描述来自药物反应和基因破坏试验的表达数据的分级分组。
发明详述下文描述包括来自对酵母(Sarchomyces)研究的特定实施方案。用于分析酵母基因之间关系的方法同样地适用于包括真核、原核和病毒的不同物种的基因间关系的分析。因此,下面的描述和实施方案是说明性的,而不是对本发明范围的限制。
I.产生自多重破坏全基因组表达文库的酵母基因调控网络的生物发现从重要的基因组范围试验数据所阐明的基因调控网络,有助于从对已知和未知生物功能的基因间的转录调控事件的观察中发现新的基因功能信息和表达调控事件。在建立推导调控网络中尤其重要的是从中获得假定控制关系的表达试验之间的生物相关性。产生于竞争杂交破坏试验的基因组范围的表达数据具有内部控制和定量测定一个基因存在与否对其它基因表达的直接作用的优点。选择破坏试验来详尽地阐明控制关系对于产生有用的基因调控信息是重要的。产生用于阐明对已知参与转录调控的120个基因的转录调控的可靠、广泛的数据。本发明人报道了用该表达文库所发现的已知转录因子和先前未知生物学功能的其它基因之间的新的调控关系。
如附图1所示,本发明人已经对酵母基因组实施了一种系统、相互作用的方法,该方法用基因调控推导结合全基因组生物表达试验。在一些实施方案中,以建立大量全基因组表达试验的文库开始,每个试验有一个基因表达被破坏。用计算技术从这些数据中推导基因表达调控关系的近似值。然后,用计算可视化和模拟软件检验调控关系的生物相关性,并通过其它数据库以及进一步通过试验,包括组合破坏试验,在新的生物学上感兴趣的子网络上验证本发明。
用酵母全基因组c-DNA微阵列构建基因表达数据文库。该文库包括120个酵母菌株的表达试验,每个试验有一个基因通过同源重组被破坏。该文库的每种基因被选择用于试验是因为它们在酵母蛋白质组数据库(PYD)中被报道为参与转录调控的因子。先前报道的基因调控网络表明,基因自身之间及其与其它调控基因之间会发生相相互作用用(10)。为了从表达文库重建分级调控关系,在一些实施方案中,本发明人开发了一种新布尔算法,该算法适合一般的、成环的(looped)调控关系。如附图2所示,通过该方法建模的基因调控关系,可以被描绘成各个表达试验测定的5871个基因中的两个特定基因之间的基因表达的上调或者下调关系的定向结构图(directed graph)。本发明人构建了552个基因成员的调控控制关系模型,然后进一步限定了(constrain)含有98个公知转录因子的子网络模型。在附图3中表示结果模型,含有总共552个节点代表被包含的基因和2953条这些基因之间的推定调控连接。
在一些实施方案中,根据YPD定义的细胞功能角色(CFR)对网络模型中的转录因子进行分类。附图4表明网络中的分类的转录因子间的控制关系。只有SKN7和HMS2直接影响细胞周期相关基因的表达。确定多条从“碳水化合物代谢”基因发出到所有其它功能基因群的控制线,该发现与许多细胞过程和代谢路径的能量依赖特性是一致的。
如附图4所示,一个显著特点是脂肪酸代谢转录因子的表达水平专一地受碳水化合物转录因子控制。这种关系被证明,参与磷脂合成路径的蛋白质之间的大量相互作用已经被报道,磷脂合成路径具有葡萄糖反应路径、脂质信号路径和其它脂质合成路径(11)。
本发明人的新方法被用于进一步从所有基因表达实验数据中搜索表达调控基因之间的详细关系,例如转录因子与调节和非调节基因之间。可以根据未报道生物功能的基因受到和/或对已知功能的基因的表达控制来表征其调控作用。附图5表示新阐明的转录因子间的控制关系的实例,这些因子参与细胞分化调控和DNA复制/修复调控。本发明人发现在对应于细胞分化调控和DNA复制/修复的子网络中,有两个离散的功能性分支通过UME6和MET28相连,表明了这两个转录因子在协调这些相互依赖的调控路径的协调表达调控中的重要作用。MET28,如其名字所暗示,先前被表征为与甲硫氨酸代谢相关的转录因子(12)。在配对型双杂交分析(mating-type two-hybridassays)中,Met28p作为较大连锁的染色体分离蛋白的一部分在新推定的染色体分离调控中的作用得到已报道的Met28p与已知的染色体分离成分Smclp的相互作用(13)的支持。
通过上述子网络中的基因的编码序列和上游区域的序列分析,评价转录因子与其目标基因以及DNA结合序列之间的序列水平控制机制。在UME6被认为是许多减数分裂基因及其控制系统的总转录调控(14,15)的情况下,我们对本发明模型中受UME6控制的34个基因的上游区进行了基序引发的多重期望-最大化(Multiple Expectation-Maximization for motif Elicitation)(MEME)分析(16)。发现两个共有序列,TAGCCGCCGA(SEQ ID NO1)和TGGGCGGCTA(SEQ ID NO2),分别以14.7%和32.4%存在于34个基因中,根据MEME搜索,P值显著。依照TRANSFAC数据库,TAGCCGCCGA被定义为Ume6p的结合位点(14),TCGGCGGCA被报道为CAR1的抑制子的结合位点(17),CAR1被一个含有Ume6p(Ume6p、Sin3p和Rpd3p)的三成分复合物抑制(18)。除了Ume6p相关结合基序,在11个减数分裂相关基因的上游不存在其它MEME共有序列,结果表明这11个基因专一地受UME6调节并且Ume6p直接影响它们的表达。仅两个其它基因具有推定的结合序列,但是它们在试验中不表现出由UME6产生的表达影响。
从表达控网络模型试验获得的发现可使人们获得有关基因调控路径的特定生物学认识,这些基因调控路径不能容易地从可以得到的生物学资料中重建,也不存在于已收集的路径数据库中。本申请已经表明这样一种系统在发现新基因功能信息以及新调控机制中是有用的。用此方法或者相似策略从全基因组表达文库来阐明分级调控路径,将实现迅速了解可用于理论药物发现和农业化学打靶的转录调控。
方法布尔网络推导算法从一组基因破坏试验建立基因表达矩阵E。矩阵元素(element)E(a,b)的值表示基因‘b’在删除基因‘a’引起的基因‘a’破坏体中与正常条件下的表达比。
(1)利用基因表达矩阵E,如果基因‘a’的破坏,导致基因‘b’的强度变得高于特定阈值θ或者变得低于特定阈值1/θ,那么定义基因‘a’直接或间接地作用于基因‘b’,二元素矩阵R中的元素(a,b)的值设为1;R(a,b)=1。这样,通过以阈值(θ或1/θ)分割基因表达矩阵E中的每个元素,建立二元素矩阵R。
(2)在二元素矩阵R中,如果存在基因‘a’和‘b’相互影响的关系,也就是R(a,b)=R(b,a)=1,无法确定哪个基因位于上游。这就是该方法的局限性或缺点,但是,引入等价集合,得到一个由相互影响的基因组成的组,该组被假定为一个基因。
(3)排序基因(拓扑分选(sort))的等价集合具有半排序关系,而且可以以半排序(拓扑分选)编排等价集合来推导网络。
(4)基干(skeleton matrix)矩阵;等价集合之间的半排序可达性矩阵(semi-ordered accessibility matrix)包含间接作用。为了除去它们和构建基干矩阵,设定每个等价集合的行定义如下行1的等价集合对其它等价集合不提供间接作用。行3的等价集合直接作用于行2的等价集合而间接作用于行1的等价集合。在给每个等价集合设定行之后,从半排序可达性矩阵中除去所有间接作用。
(5)绘制多水平有向图(diagraph);根据基干矩阵的各个元素的值在节点之间连线。
微阵列试验可以用任何本领域已知的常规方法获得基因表达信息。例如可以采用微阵列,其中将每个要测定表达的基因的特有互补DNAs(cDNA)被放在基质(substrate)上。得自样品的RNA通过与相应的cDNA杂交来分析,并且可用多种方法检测。在一些实施方案中,期望使用2002年5月23日申请的编号为60/382,669的美国临时专利申请描述的微阵列探测物和/或方法,在此全文引作参考。
用全基因组酵母c-DNA微阵列(13)收集基因表达数据。以BY4741(MATa,HIS3Dl,LEU2D0,METl5D0,URA3D0)作为野生型菌株。菌株BY4741的基因破坏物购自Research Genetics,Inc。细胞被接种并在YPD(1%酵母提取液,1%细菌蛋白胨,2%葡萄糖)30℃下生长,直到在对数生长阶段,OD600达到1.0,收获来分离用于基因表达分析的mRNA。亲本菌株用作每个破坏体菌株的对照。
数据归一化通过cDNA微阵列分析来测定来自155个破坏体的5871种mRNA种类的量。Cy3和Cy5之间荧光强度的不同引起表达量比例的偏差。归一化每个表达图谱的表达量比例。在各个标记区段(spotted block)中比例偏差具有固定趋势,因此计算线性回归将每个区段的平均值比例归一化为1.0。比例的对数值用于表示标准表达水平,因此,发现了比例的对数值并计算这些对数值的平均数和标准偏差(见表1)。来自UME6(YDR207C)破坏体表达阵列的所有标记基因表达水平的标准偏差(SD)为0.4931,UME6被定义为YPD中破坏体的“总调控子”(33),因此,在阵列数据中可能有不可接受量的误差,其总SD大于0.5。因此,在一些实施方案中,可以从分析中排除这些试验。用不同方式选择数据的SD是有价值的。选择SD小于0.5,可以采用相对传统的方法,其中上述类型的误差相对不会影响整个基因网络推导。然而,可以选择SD小于0.4、小于0.3、小于0.2、小于0.1小于约0.05、小于约0.01、小于约0.005或者小于约0.001。
选择基因在YPD中,有314个基因被定义为“转录因子”,其中98个已经进行过调控机制的研究。包括该98个“转录因子”所控制的基因的552个基因的表达图谱数据选自5871个图谱。从而从表达图谱数据集来构建基因调控网络,这些数据集建立在120个基因破坏试验中的552个基因的表达值上。
参考文献1.Pilpel,Y.,Sudarsanam,P.& Church,G.M.通过启动子元素件的组合分析来确定调控网络,自然遗传学(Nature Genetics,29,153-159(2001)。
2.Friedman,N.,Linial,M.,Nachman,I.& Pe`er,D.用贝叶斯网络分析表达数据,计算生物学杂志(J.Comput Biol),7,601-620(2000)。
3.Somogyi,R.& Sniegoski,C.A.给基因网络的复杂性建模了解多基因和多向调控,复杂性(Complexity),1,45-63(1996)。
4.Spellman,P.T.等,通过微阵列杂交来综合确定酵母属酿酒酵母的细胞周期调控基因,分子细胞生物学(Mol Cell Biol),9.3273-3297(1998)。
5.Gasch,A.P.等,酵母细胞对环境变化的反应中的基因组表达程序.分子细胞生物学(Mol Cell Biol),17.2099-2106(1997)。
6.Akutsu,T.,Miyano,S和Kuhara,S.基于矩阵乘法和指纹功能的用于确定布尔网络和相关生物学网络的算法,计算生物学杂志(J.ComputBiol).7,331(2000)。
7.Akutsu,T.,Miyano,S & Kuhara,S.推导基因网络和代谢路径之间的数量关系,生物信息学(Bioinformatics),16,727(2000)。
8.Brown,M.p.等,用支持向量设计来基于知识分析微阵列基因表达数据,Proc.Natl.Acad.Sci.97,262-267(2000)。
9.Ideker,T等,被系统地扰动的代谢网络的综合基因组和蛋白质的分析,科学(Science).292,929-934(2001)。
10.Shoudan,L.REVEAL,一种通用的推导基因网络构建的反转工程算法,Pac.Sym.Bio.3,18-29(1998)。
11.Hiesinger,M.,Roth,S.Meissner,E.和Schuller,H.J.Cat8和Sip4通过碳源反应元素件对酵母葡糖异生基因的转录激活作用。遗传学最新进展(Curr Genet),39,68-76(2001)。
12.Kuras,L.,Cherest,H.,Surdin-Kerjan,Y.,& Thomas,D.含有着丝粒编码因子1和两个主要亮氨酸拉链因子Met4和Met28的异质复合物协调酵母硫代谢转录活性,Embo Joumal 15,2519-2529(1996)。
13.Newman,J.R.S.,Wolf,E.& Kim,P.S.从头开始确定来自酿酒酵母的相互作用卷曲螺旋的计算导向屏幕,Proc.Natl.Acad.Sci.97,13203(2000)。
14.Bowdish,K.S.,Yuan,H.E.& Mitchell,A.P.分子细胞生物学(MolCell Biol),15,2955(1995)。
15.Rubin-Bejerano,I.,Mandel,S.,Robzyk.,K.& Kassir,Y.通过与转录激活Imel联合调控将转录抑制子Ume6反转为正向调节物来诱导酿酒酵母减数分裂,分子细胞生物学(Mol Cell Biol),16,2518(1996)。
16.Bailey,T.L.& Elkan,C.通过期望值最大化揭示生物高聚物中的基序来确定复杂的模型,分子生物学智能系统第二次国际会议会议记录(PSICISMB),pp28-36,AAAI Press,Menlo Park,加利福尼亚,1994。
17.Luche,R.M.,R.Surnrada,M.& Cooper,T.G.多基因中的反作用元素件作为酵母CAR1基因的抑制子蛋白结合位点,分子细胞生物学(Mol CellBiol),10,3884(1990)。
18.Messenguy,F.,Vierendeels,F.,Scherens,B.& Dubois,E.在酿酒酵母中,对外源氮利用响应的精氨酸代谢基因CARl和CAR2的表达受到Ume6(CargR I)-Sin3(CargR II)-Rpd3(CargRIII)复合物调节,细菌杂志(J.Bacteriol)182,3158(2000)。
II.贝叶斯网络和非参数回归估计基因之间的基因网络和功能性结构关系1.引言遗传工程的发展产生大量有价值的数据,如微阵列基因表达数据。对基因间的关系分析引起人们对分子生物学和生物信息学领域的极大关注。但是,由于数据的维度和复杂性的原因,发掘隐藏在噪声中的结构关系不是一件容易的事。为了从生物数据中提取有效信息,从统计学角度考虑,需要发展理论和方法学。本发明的目的在于建立更清楚地了解基因间关系的新方法。下述实施例1提供这些方法的数学原理。
在本发明的一些实施方案中,利用来自图像-理论方法的微阵列基因表达数据,采用贝叶斯网络构建基因网络。Friedman和Goldszmidt(1998)提出了一种用贝叶斯网络构建遗传连接的有益方法。他们离散基因的表达值,假定符合基于多项分布的模型。但是仍存在选择离散用的阈值的问题(不仅限于试验)。阈值的确说明结果作用的实质变化,而不适当的阈值产生不正确结果。另一方面,最近Friedman等(2000)指出离散会导致信息丢失。为了将表达数据用作连续值,他们采用基于线性回归的高斯模型。但是,这种模型仅仅能够检测线性相关,而不能完整说明关系。本发明人利用贝叶斯网络开发了用于构建基因网络的新方法。用具有高斯噪声的非参数回归模型既捕获基因之间的线性关系又捕获非线性结构关系。非参数回归被开发用于在预先缺乏功能性关系知识的情况下探索所期望反应的复杂的非线性形式。鉴于贝叶斯网络的新结构,需要一个估计模型的合适的检验数。因此,本发明包括来自贝叶斯统计的新检验数。采用这些方法,克服了已有方法的缺点并且获得更多信息。此外,本发明的方法包括作为特例的前述方法。通过分析酿酒酵母细胞周期数据来评价本发明的方法。
本发明人发现按照这些方法采用贝叶斯分析,可以确定基因之间的上游因果关系,从而确定潜在的治疗靶点。
2.贝叶斯网络和非参数回归利用贝叶斯网络结构,将基因视为随机变量,将联合概率分解成条件概率的乘积。例如,如果有一系列随机向量的观察值,可以根据条件概率的密度来表示获得特定观察值的概率。在一些实施方案中,用非参数回归模型来捕获变量间的关系。可采用各种绘图工具来阐明上述关系。例如,多项式(polynomials)、傅立叶级数(Fourier series)、回归样条函数(regression splinebases)、B-样条函数(B-splines)、小波函数(wavelet basese)等可以被用于定义基因关系结构图。选择恰当的结构图的困难在于恰当地评价系统中的变量和噪声。
3.选择恰当结构图的检验数在一些实施方案中,使π(θG|λ)为具有超级参数向量(hyper parametervector)λ的未知参数θG的先验分布(prior distribution),并且logπ(θG|λ)=0(n)。在参数范围内求积分得到数据Xn的边缘概率(marginal probability),选择具有最大后验概率(posterior probability)的结构图G。Fiedman和Goldszmidt(1998)将多项式分布作为贝叶斯网络模型,并且对参数θG提出Dirichlet先验概率(prior)。此时,Dirichlet先验概率是共轭先验概率(conjugate prior),而后验分布属于相同分布类型。如此得到(4)中的积分的闭式解,他们称之为选择结构图的BDe分值。Bde分值限于多项式模型,而本发明人提出了在更普遍和更多种情形下的选择结构图的检验数。
构建基于上述模型的检验数的问题之一是如何计算积分。可以采用包括马尔可夫链(Markov chain)和蒙特卡罗(Monte Carlo)方法等方法。在本发明的一些实施方案中,采用拉普拉斯逼近(Laplace approximation)来求积分。因此,选择最佳结构图,使得实施例1方程(5)中的检验数BNRC最小化。
该检验数得自log(θG|λ)=0(n)。如果log(θG|λ)=0(1),众数(mode)θG等价于最大似然估计、MLE和从贝叶斯信息标准得到的检验数,被称作为除去较高阶限制(higher order terms)0(n-j)(j≥0)的BIC。Konishi(2000)给出根据Kullback-Leibler信息和贝叶斯方法构建模型选择检验数的一般框架。从下述实施例1可以看出,由于以非循环构建结构图,所以可以在BNRC最小时选择最终结构图,而不必最小化每个局部分值BNRCj。
4.用BNRC检验数估计结构图和相关结构例如,可以以图的方式表示本发明的方法。本发明方法的关键点在于使用非参数回归和从贝叶斯统计中得到用于选择结构图的新检验数。对于本实施例1第2节的非参数回归,以B-样条函数作为基本函数。实施例1图1是具有等距点t1...,t10的3次(degree)B-样条函数的一个例子。当给定βjk值时,采用修合(backfitting)算法获得模型。修合算法在下面的实施例1中说明。
所述模型依赖于超级参数βjk,应选择最佳βik值。在本发明方法中,当BNRCj最小时选择最佳βjk值。
通过最大化实施例1方程(6)来估计B-样条函数系数向量。方程(6)的众数与所罚似然估计相同,可以将超级参数λjk或βjk视为所罚似然估计中的修匀参数(smoothing parameter)。因此,超级参数在将曲线与数据拟合中起重要作用。
5.计算试验采用蒙特卡罗模拟来检验本发明方法的特性。数据产生于变量之间的虚拟(artificial)结构图和结构(实施例1图2),结果总结如下检验数BNRC可以检测数据间的线性和非线性结构关系。但是BNRC分值可能具有扩大结构图变化的倾向。因此,采用称作AIC的Akaike信息检验数并使用这两种方法。AIC最初作为评价由最大似然方法估计的模型的检验数被引入。该方法的估计与最大所罚似然估计相同,但与MLE不同。Sjk轨迹表明拟合曲线的自由度,这十分有用。也就是说,如果trSjk接近2,可以认为相关是线性的。使用BNRC和AIC来决定是否增加亲本变量。
分析Spellan等和Friedman等论述的酿酒酵母细胞周期数据来评价这些方法。这些数据从800个基因和77个试验收集得到。由于缺少关于正确结构图规格的信息,设定先验概率πG为恒量。由20个B-样条函数组成非参数回归。B-样条函数的数量可以根据需要改变。超级参数决定拟合曲线的光滑度,期望的是通过选择超级参数和B-样条函数的个数来用最少B-样条函数获得好的数据拟合。
分析结果总结如下实施例1图3表示用一个基因预测CLN2、CDC5和SVS1时的BNRC分值。给出较小BNRC分值的基因对目标基因产生更可靠的作用。观察哪些基因与目标基因相关,找出强烈依赖于目标基因的表达的基因组群。事实上,可以用这些信息构建简明网络。可以将最佳结构图视为考虑相互作用效应的简明网络的修订版。如果基因之间线性相关,当亲本-子代关系被反转时,BNRC分值是好的。因此,特别是当相关接近线性时,结构图中的因果关系方向不精确。一些基因,如MCD1、CSI2和YOX1等,构成了Friedman等的结果。这些基因中的大多数被报道为具有重要作用。这些基因之间的大量关系几乎是线性的。但是,本发明人发现了线性模型极少能发现的一些非线性相关。
实施例1图5表示与基因相关的所估计的结构图,这些基因根据其进入细胞周期的过程及其邻近基因进行分类。在图5中省略了一些分支,但是表明了重要信息。对于本发明人和Friedman等给出的网络,确定亲本-子代关系,并且可以看出这两个网络是彼此相似的。特别地,本发明的网络包括了Friedman等报道的典型关系。至于这两个网络之间的不同处在于SVS1的亲本基因。Friedman等以CLN2和CDC5作为SVS1的亲本基因。一方面,本发明结果是为SVS1提供了CSI2和YKR090W。本发明的候选亲本基因比Friedman等的候选亲本基因更恰当。本发明的模型与实施例1图4的两种情形恰当地吻合。特别地,可以得出与Spellman数据中的其它基因相比,CDC5对SVS1的作用较弱的结论(又见实施例1图3)。事实上,作为SVS1的亲本基因,CDC5的BNRC分值的排序为第247位。考虑到上述环境因素,本发明方法提供可以理解的有用形式的有价值信息。
附图6图示说明本发明的用于确定基因网络关系的方法。
6.讨论应用贝叶斯网络和非参数回归从微阵列基因表达数据中评估基因网络的新方法产生更高的精确度。本发明人理论性地得到用于选择结构图的新检验数,并通过分析细胞周期数据证实其有效性。本发明方法的优点包括(1)将表达数据用作连续值;(2)可以检测非线性结构关系和可视化其功能结构关系使之容易理解;和(3)自动搜索实现最佳结构图的创建。
Friedman等的方法保留了已知参数,如用于离散的阈值和Dirichlet先验中的超级参数,这些参数通过反复试验来选择并且从某种意义说没有被优化。另一方面,本发明方法可以根据具有良好理论基础的检验数自动而恰当地估计任何参数。
在其它实施方案中得到更加普遍情形时的检验数BNRC。通过这些实施方案,可以根据其它统计模型建立结构图选择检验数。
III.通过贝叶斯网络和具有异质误差方差及相互作用的非参数回归建立非线性基因网络模型在其它实施方案中,用不同的统计学方法从基于贝叶斯网络的微阵列基因表达数据中构建基因网络。在这些实施方案中,估计每个随机变量的条件分布。考虑用异质误差方差及相互作用拟合非参数回归来捕获基因之间的非线性结构关系。一旦用贝叶斯网络和非参数回归来建立结构图,仍需解决选择最能代表基因系统的最佳结构图的问题。得自贝叶斯方法的用于选择结构图的新检验数包括用贝叶斯方法估计网络的前述方法。通过分析最新从破坏100个基因得到的酿酒酵母基因表达数据证实了所提出的方法的有效性。
1.引言在通过大量变量的联合分布给现象建模方面,使用贝叶斯网络是一种有效方法。Friedman和Goldszmidt离散表达值,并假定多项式分布为候选统计模型。Peter等研究了用于离散的阈值。Friedman等指出离散可能会丢失一些考虑拟合线性回归模型的信息,而所述模型将数据视为连续的来分析。但是,亲本基因线性地依赖于目标基因的假设并不总是能得到保证。Imoto等描述采用非参数附加回归模型来捕获基因间的线性关系以及非线性关系。
在一些实施方案中,使用贝叶斯网络和非参数异方差回归,其更加不受无关值的影响,并且还能捕获亲本基因间的相互作用。
在确定结构图后,评价其优点以及与完全未知的实际结构图的接近程度。建立合适的检验数变得重要。Friedman和Goldszmidt得到根据多项式模型和Dirichlet先验函数选择结构图的检验数BDe。但是,Dirichlet先验函数中的未知超级参数仍然存在,并且只能通过试验确定其值。本发明人从贝叶斯方法中得到用于选择结构图的新检验数。该检验数自动地优化所有参数模型并产生最佳结构图。此外,本发明的方法包括通过贝叶斯网络构建基因网络的前述方法。本发明通过破坏100个基因,分析酿酒酵母基因表达数据来证明这种新方法的有效性。
2.贝叶斯网络和具有相互作用的非参数异方差回归模型假定得到p维随机可变向量X=(X1...,Xp)T的n组数据{x1,....,xn},其中xi=(xi1...,xip)T和xT表示x的转置矩阵(transpose)。在微阵列基因表达数据中,n和p代表阵列和基因的数量。在贝叶斯网络结构中,我们考虑在节点间定向非循环结构图G和马尔可夫假想。联合密度函数分解成条件密度函数。假定条件密度函数通过参数向量θj参数化,从这些概率模型中提取有用的信息。
然后用非参数回归策略来捕获xij和Pij之间的非线性关系。多数情况下,这种方法可以很好地捕获目标关系。但是,当数据中含有特别接近域{pij}的边界的无关值时,非参数回归模型有时会导致不恰当的平滑估计,即估计的曲线显示一些由无关值作用产生的伪波形。为了避免这个问题,下面在实施例2中描述具有异质误差方差的非参数回归模型。
选择结构图的检验数一旦确定结构图,可以构建基于贝叶斯网络和非参数回归的统计模型(实施例2方程8),并可以通过合适方法进行估计。但是,由于似然值在较复杂模型中会变得比较大,为了选择最佳结构图,以对隐藏在数据中的系统给出最佳估计,不宜采用似然函数作为模型选择检验数。因此,基于概括的或预期的误差的统计方法,Kullback-Leibler信息、贝叶斯方法等(20)是合适的。因此,本发明根据来自贝叶斯方法的模型(实施例2方程8)构建用于评价结构图的检验数。
如下从贝叶斯理论方法构建评价结构图优点的检验数(1)通过结构图πG的先验概率和数据的边缘概率的乘积得到结构图的后验概率;(2)除去标准化恒量,结构图的后验概率与实施例2的方程9相称。
(3)用贝叶斯方法,选择结构图使得π(G|Xn)最大。
(4)确定计算是否包括高阶积分(high dimensional integral)(实施例2方程9)。
(5)如果是,回到(4),采用实施例2方程11,[见实施例2中关于积分的参考文献17和26]。
如此,选择结构图使得检验数BNRC最小。采用Laplace方法的优点之一在于不必考虑使用共轭先验分布。因此,实现了在更多种类模型分布和先验函数中建模。
3.评估基因网络非参数回归本发明提供根据上面第2节中描述方法构建基因网络的方法。非参数回归(实施例2方程5)含有两个成分用每个亲本基因的附加模型表示的主效应成分和相互作用成分。在附加模型中,通过B-样条函数(实施例2方程9)构建各个光滑函数mjk(·)(见参考文献18)。
在相互作用方面,我们采用高斯径向基本函数(Gaussian radial basisfunction)。在根据径向基本函数的回归建模中,采用两种方法来估计中心zij和宽度,这种方法称作为“完全监督学习”(fully supervised leaning)。可代替的方法通过预先仅采用亲本的观察数据来测定数值。可以采用后一种方法,采用k-法聚类算法(k-means clustering algorithm)来构建基本函数。在实施例2中的参考文献7、21和23中进一步描述径向基本函数的细节。超级参数控制基本函数中的重叠的数量。
在误差方差中,我们考虑异方差回归模型并推测出实施例2方程6所示结构。恒量的设计会影响捕获数据异方差的精确性。按照实施例2方程12设定权重。如果超级参数p设定为0,那么变量是齐次的。如果采用大的p值,那么接近于亲本变量区域边界的数据的误差方差也大。因此,如果存在接近边界的无关值,可以采用恰当的p值来降低其影响并得到合适光滑估计。
4.真实数据分析通过分析酿酒酵母基因表达数据,证明上述方法的有效性。观察通过基因破坏,微阵列上基因表达水平的变化。通过使用这种方法,揭示酿酒酵母基因之间的基因调控网络。从基因破坏试验中收集大量表达图谱来评价基因调控网络。备置400个以上的突变体,积累基因表达图谱。
用扫描仪(scanner)监控打点在微阵列上的5871个基因的转录水平。400个以上的破坏体的表达图谱被收集在本发明数据库中。评定微阵列上的所有基因的标准偏差(SD),SD值粗略地表示试验的误差。以0.5作为所有基因表达比的标准偏差的临界点。从400个图谱中选择107个破坏体,其包括68个转录因子被破坏的突变体。
使用100个微阵列,从上述数据中构建521个基因的基因网络。发现94个确定了调控基因的转录因子,从5871种表达图谱中选择了421个基因的表达图谱,这些基因受94个因子控制。用20个B-样条函数和20个径向基本函数构建非参数回归模型。由于无论使用多少基本函数,超级参数都控制着拟合曲线的光滑度,所以,对不同数量的基本函数的光滑估计间的差别不能从视觉上发现。对相互作用的成分应用双基因作用。因此,相互作用效应以拟合面得到,可以从视觉上被了解。
本发明说明了超级参数在先验分布和权重恒量中的作用。实施例2图1(a)表示用超级参数的3种不同值光滑估计的YGL237C和YEL071W的扫描图。光滑估计依赖于超级参数值。实施例2图1(b)描述图1(a)中的两个基因的BNRC检验数的变化。选择超级参数的最佳值作为BNRC最小值,最佳光滑估计(立体曲线)能够很好地捕获这些基因间的结构关系。这些虚的和带点的曲线分别与最大似然估计和参数线性拟合接近。实施例2图2(a)表明权重恒量wlj...,wnj的作用。如果采用同方差回归模型,得到虚曲线,其显示一些由左上角的数据作用产生的伪波形。通过调整实施例2方程12中的超级参数p,估计曲线变为立体曲线。通过最小化BNRC检验数来选择p的最佳值(见实施例2图2(b))。当恰当地得到光滑估计时,p的最佳值趋向0。
采用两步策略来用相互作用拟合非参数回归模型。首先,估计由附加B-样条函数回归值表示的主效应。接着,用相互作用成分拟合残差(residual)。实施例2图3描述YIL094与其亲本基因CYKL152和CYER055C的关系的拟合面的例子。当增加亲本基因时,这两个基因的相互作用导致过度表达。
在酿酒酵母中,GCN4基因编码氨基酸生物合成的“总控制”(generalcontrol)系统的转录激活子,该系统是一个至少12条不同生物合成路径的网络[实施例2的参考文献6]。试验说明了对由组氨酸类似物3-氨基三唑诱导的“氨基酸饥饿”信号应答的涉及GCN4p水平的总控制结果。GCN4激活30个以上基因的转录,这些基因参与11种对氨基酸饥饿或者tRNA合成酶受损响应的氨基酸的生物合成[见实施例2的参考文献24]。嘌呤的生物合成基因ADE1、ADE4、ADE5和ADE7响应氨基酸饥饿时,表现为GCN4-依赖的表达[24]。GCN4对氨基酸或者嘌呤饥饿相应时,GCN4激活氨基酸和嘌呤的生物合成基因的转录[24]。这些试验结果表明嘌呤代谢和氨基酸代谢之间通过GCN4密切联系。本发明的关系图谱与嘌呤和氨基酸代谢之间的已知关系十分吻合。
IV.应用贝叶斯网络和非参数异方差回归建立基因网络的Bootstrp非线性模型在其它的实施方案中修改上述方法。贝叶斯网络结构的相关特征在于各个随机变量的条件分布估计。将具有异质误差方差的非参数回归模型与微阵列基因表达数据拟合,来捕获基因之间的非线性结构。选择最能代表基因中的系统的最佳结构图的问题仍然有待于解决。本发明人从贝叶斯方法中得到新的一般情形下的结构图选择检验数。所提出的方法包括基于贝叶斯网络的前述方法。本发明人还采用一种方法来测定边缘强度和所评估的基因网络中的贝叶斯因果关系方向的置信度。通过分析破坏100个基因而新得到的酿酒酵母基因表达数据,证实所提出的方法的有效性。
一旦选定了结构图,可以评价其优点或其与未知的真实结构图的接近程度。对于本发明方法,需要建立测定其边缘强度的方法。采用“自举(Bootstrap)”方法(Efron,1979;Efron和Tibshirani,1993)来解决该问题。采用该方法不仅可以测定边缘强度,而且还可以测定贝叶斯因果关系的置信度。通过分析破坏100个基因而新得到的酿酒酵母基因表达数据,证实所提出的方法的有效性。
(1)用于基因网络非线性建模的贝叶斯网络和非参数异方差回归1.1非线性贝叶斯网络模型如本申请其它部分所描述,在贝叶斯网络结构中可以在节点之间假定存在定向非循环曲线G和马尔可夫假说。然后将联合密度函数分解成各个变量的条件密度函数。根据下述的实施例3的方程1,在通过贝叶斯网络统计建模中所关注的焦点是怎样构建条件密度函数(densities)fj。可以假设条件密度函数fj被参数向量参数化,从实施例3进一步描述的概率模型中提取信息。
在一些实施方案中,可采用非参数回归策略捕获xij和Pij之间的非线性关系,表明基因之间存在许多非线性关系而线性模型难以获得充分的结果。在许多情况下,这些方法可以很好地捕获目标关系。但是当数据中含有无关值时,特别是该无关值接近域的边界时,标准非参数回归模型有时会产生不合理的光滑估计,即,由于无关值的作用,估计曲线显示一些伪波形。为了避免该问题,采用异质误差方差来拟合非参数回归模型。如果模型中参数的数量大大地多于观察的数量,倾向于产生不稳定的参数估计。
由于线性回归不能评价贝叶斯因果方向或者在许多时候会产生错误方向,所以在一些情况下期望用非参数回归来代替线性回归。用一个简单例子说明与线性回归相比时该新模型的优点。假定已经获得实施例3图1(a)中基因1和基因2的数据。我们认为两个模型基因1>基因2和基因2>基因1,并且得到分别如实施例3图(b)和图(c)所示光滑估计。然后根据本发明的检验数可以确定模型(b基因1→基因2)好于模型(c基因2→基因1),该检验数在后面章节中得到(模型的分值为(b)120.6和(c)134.8)。由于这些数据是从真实结构图基因1→基因2中产生,所以本发明的方法得到正确结果。但是,如果用线性回归模型与这些数据拟合,会选择模型(c)(线性模型的分值为(b)156.0和(c)135.8)。如此,基于线性回归的方法得到错误结果,而非线性回归分析获得正确结果。
即使在基因之间的关系几乎为线性的情况下,本发明的方法和线性回归都可以恰当地与数据拟合。但是,采用线性模型难以确定贝叶斯因果关系的方向。
先前已经描述了检验数BNRC。本发明人得到了用于选择在一般网络结构中的结构图的检验数BNRChetero。利用方程(8),累加局部分值BNRChetero可以得到结构图的分值BNRChetero。选择最佳结构图,使得检验数BNRChetero(实施例3方程7)最小。拉普拉斯(Laplace)方法的优点在于不需要考虑使用共轭先验分布。因此实现在多种模型分布和先验分布中建模。按实施例3中进一步描述的来评估基因网络。但是考虑到误差方差,我们评价了一个异方差回归模型并且采用了实施例3方程3所示结构。
网络学习(learning Network)在贝叶斯网络内容中,确定最佳网络是一个NP-难题。为了解决该问题,可以采用如下的“贪心爬山”(greedy hill-climbing)算法(1)步骤1生成分值矩阵,其第(i,j)个元素是结构图基因i→基因j的BNRChetero分值;(2)步骤2对每个基因实施用于边缘(edge)的3种处理之一增加、去除和反转(reverse),得到最小BNRChetero;(3)步骤3重复步骤2直到BNRCherero不再减小;和(4)改变基因的计算顺序,在步骤3中得到多个候选学习顺序。
在如下情形中期望进行步骤4,即贪心爬山算法产生许多局部最小值,而该结果依赖于变量的计算顺序。网络学习中的另一个问题是,当基因数目多时,亲本基因的搜索范围很宽。在这种情形下,可以基于步骤1所给出的分值矩阵限定候选亲本基因。
可以将该学习策略应用于基因网络学习,并且可以通过蒙特卡罗模拟方法来证明该方法的有效性。本发明人也通过相同的蒙特卡罗模拟方法来检验本发明新模型的有效性,并发现由于应用了非参数异方差回归模型而产生了进步。在后面的章节中说明非参数异方差回归模型的有效性。
超级参数考虑实施例3方程4所定义的非参数回归模型。估计值θj是logπ(θj|Xn)的众数,依赖于所用超级参数。用20个B-样条函数构建非参数回归模型。我们确认,对应于多个不同基本函数的光滑估计值间的差别难以从视觉上被发现。当我们使用稍大量的基本函数时,超级函数控制拟合曲线的光滑度。(实施例3图3(a)显示YGL237C和YEL071W的离散图和超级参数的3个不同值的光滑估计值。在后面章节中说明数据的详细内容。光滑估计强烈地依赖于参数值。实施例3图3(b)描述图3(a)中的两个基因的BNRChetero检验数的变化。可以选择最佳超级参数值作为最小的BNRChetero,最佳光滑估计(图3(a)中的实心曲线)能够很好地捕获这些基因之间的结构关系。虚曲线和带点的曲线分别与最大似然估计和参数线性拟合接近。
权重恒量wij,...wnj的作用在实施例3图4(a)中表明。如果采用单一方差回归模型,得到具有由左上角的数据作用产生的伪波形的曲线。通过调整实施例3方程9中的超级参数pj,估计曲线产生实心曲线。参见实施例3图4(b),通过最小化BNRChetero检验数也可选择pj最佳值。当得到合适的光滑估计时,pj最佳值趋向于零。
最后,在实施例3的章节3中提供一种用于估计光滑曲线和优化超级参数的算法。
步骤1确定超级参数pj,步骤2初始化(initializing)γjk=0,k-1,...,qj;步骤3重复步骤3-1和步骤3-2寻找最佳βjk;步骤3-1计算
γjk=(BTjkWjkBjk+nβjkKjk)-1BTjkWjkx(x(j)-∑Bjk·γjk·)k′≠k来确定βjk,步骤3-2评价对应于βjk的候选值重复步骤3-1,并选择最佳βjk值,其使BNRChetero最小化。
步骤4收敛重复步骤3以k=1,...,qj,1,...,qj,1,...直到满足合适的收敛检验数。
步骤5对应于βjk的候选值重复步骤1到4,并选择使BNRChetero最小化的最佳βjk值。
自举边缘强度和贝叶斯因果关系的置信度通过自举方法测定所评估的基因网络中的边缘强度和贝叶斯因果方向的置信度。算法表达如下(1)通过随机取样n次,从原始基因表达数据{x1,...,xn}置换生成自举(bootstrap)基因表达矩阵X*n={x*1,...x*n}T;(2)根据所提出的方法从X*n评估建立基因网络;和(3)重复步骤1和步骤2T次。
从该算法得到T个基因网络。如下确定自举边缘强度和贝叶斯因果关系的方向。
边缘强度如果基因i→基因i和基因j→基因i边缘在T个网络中分别出现t1次和t2次,定义基因I和基因j之间的自举边缘强度为(t1+t2)/T。
贝叶斯因果关系的置信度如果t1>t2,则选择基因i→基因j的方向并且定义因果关系的置信度为ti/(t1+t2)。但是,可以使用某一阈值。例如,设定一个确定的数值δ,如果ti/(t1+t2)>δ,则确定基因1→基因2。
至少有两种方法可以用于说明生成的网络。第一,期望确定在原始基因网络中的边缘强度和贝叶斯因果方向的置信度。可以给原始网络的各个边缘增加强度。从该网络确定原始网络的可信度。第二,将自举网络与原始网络叠加。但是,叠加的网络含有强度小的边缘。因此,可设定特定的阈值,去除强度小于阈值的边缘。然而设定阈值存在问题,即可视化问题。叠加的网络可能不支持非循环假设,但是存在许多有效信息。
真实数据分析通过分析酿酒酵母基因表达数据来说明本发明方法的有效性。用扫描仪监控打点在微阵列上的5871个基因的转录水平。在本发明数据库中存储了400个以上的破坏体的表达图谱。评价微阵列上所有基因的水平的标准偏差(SD)。SD值粗略地表示试验误差。在本发明的数据中,设定数值0.5作为试验精确性的临界点。根据所有基因的表达比例的标准偏差来评价那些表达图谱的精确性。从400中表达图谱中选择107个包括68个转录因子被破坏的突变体的破坏体。
其它SD值对于精确性来说也是关键的。例如,SD值可以为0.4、0.3、0.2、0.1、约0.05、约0.01、约0.005、约0.001,或者设定用以产生具有期望可信度的信息的其它任何值。
采用100个微阵列,从上述数据中构建521个基因的基因网络。发现94个其调控基因已经被清楚确定的转录因子。由该94个因子控制的521个基因的表达图谱选自5871个图谱。
实施例3表1给出具有高自举边缘强度的基因对。“inte”和“Dire”分别表示边缘强度和贝叶斯因果关系方向的置信度。“F”为亲本基因的功能,即“+”和“-”分别为诱导和抑制。如果不能确定该功能是诱导还是抑制,在“F”中输入“0”。实施例3图5表示通过叠加100个自举网络生成的部分网络。用线条宽度代表边缘强度,线条后的数值表示贝叶斯因果关系方向的置信度。
从表1可以得出如下结论表1中60%以上的基因对符合生物学知识。边缘强度等于1的大部分基因组(set)具有YPD数据库中的“相关蛋白”关系。其它一些基因组,如ARO基因和PDR基因具有相同调控系统。这些结果表明在先前未知的基因之间存在关系,并且首次被本发明的方法揭示。还发现其它9个基因组具有“功能基因组学”关系。“功能基因组学”关系表明了一些已经用已有信息揭示的一些关系,这些信息来自旨在揭示大群体蛋白质特性的大范围、高通量试验,如微阵列分析。
对于酿酒酵母的嘌呤生物合成路径调控的研究表明,细胞外嘌呤的存在会在转录水平上抑制所有编码AMP重新生物合成所需酶的基因。称为Bas1p和Bas2p的两个转录因子,对于调节ADE基因(Daigman-Fornier和Fink,1992)以及一些组氨酸生物合成基因(Denis等,1998)的活化是必需的。那些嘌呤生物合成基因,如ADE1、ADE4、ADE5,7和ADE8,对氨基酸饥饿应答表现为GCN4依赖(Rolfes和Hinnebusch,1993)。嘌呤饥饿通过与氨基酸饥饿细胞中运行的相同机制刺激GCN4翻译,导致GCN4转录活化的靶点之一HIS4表达的大副度增加。图5表明这些ADE基因和组氨酸生物合成基因与BAS1及GCN4相关。
本发明的方法的一个合适特点包括使用非参数异方差回归来捕获基因之间的非线性关系和表达数据的异方差性。因此,本发明方法在如人类基因组,来自其它真核有机体、原核有机体和病毒的基因组的未知系统的分析中是有用的。
V.应用线性样条函数统计分析小组时间-排序基因表达数据最近,采用cDNA微阵列技术,通过在少数时间点上测定基因表达水平研究了基因对环境变化的瞬时反应。传统的时序分析方法不适用于这种短的连续时间排序数据。因此,基因表达数据的分析通常被限定为倍数-变化分析,代替系统分析方法。
采用最大似然方法以及Akaike的信息检验数来将线性样条函数与小组时间-排序基因表达数据拟合,目的在于从测定值中推导具有统计意义的信息。用Student t-检验评定测定的基因表达数据的显著性。
用线性样条函数重新分析已有的藻青菌(Synechocystis sp.)PCC6803的基因表达测定值。确定许多基因的瞬时反应,这些反应被倍数-变化分析遗漏。根据本发明的统计分析,发现在各个时间点需要约4个或者更多个基因表达的测定值。在本文下述实施例4中将给出这些结论和进一步描述。
1.引言近些年来,已经进行了许多cDNA微阵列试验,测定不同条件下的基因表达水平。被测定的基因表达数据可以广泛地从公共数据库中获得,如KEGG数据库(Nakao等,1999)。
这些试验中的一些试验,在多种环境条件下测定稳态基因表达水平。例如,在不同温度下,测定藻青菌PCC6803和突变体的基因表达水平,结果确定基因Hik33为该藻青菌的潜在冷感受器(Suzuki等,2001)。
在其它试验中,通过在大量时间点上测定基因表达水平来估计基因表达的瞬时模式(pattern)。例如,在酿酒酵母细胞周期中测定周期性变化的基因表达水平(Spellman等,1998)。同一酵母种属的基因表达水平在从发酵到呼吸的代谢变换过程中被测定(DeRisi等,1997)。在这些试验中,环境条件随时间缓慢变化。相反地,可以测定基因对突然变化的环境的反应。作为一个例子,在从低光照到强光照突然变化后的多个时间点上测定藻青菌PCC6803的基因表达水平(Hihara,2001)。
在cDNA微阵列试验中,在少数时间点上代表性地测定基因表达水平。常规的时序分析方法,如傅立叶分析或者自回归或者移动平均数建模,都不适合于这种少量数据点。替代地,基因表达数据通常通过聚类(clustering)技术或者通过只考虑基因表达水平中的相关变化来分析。这种“倍数-变化”分析可能会遗漏基因表达水平中的显著变化,而这可能无意中给受噪声支配的测定值带来显著性。此外,倍数-变化分析可能不能确定瞬时基因表达反应的重要特征。
多种分析基因表达数据的技术在过去已经被使用,例如派生的布尔或贝叶斯网络(Liang等,1998;Akutsu等,2000;Friedman等,2000)。期望按照调控网络描述基因相互作用,而在大量时间点上获得基因表达数据有利于开发网络模型。但是,目前尚未获得大量基因在大量时间点的数据。虽然在特定有机体中,基因的数目可能接近数千,但是通常仅在5个或10个时间点上测定基因表达水平。
迄今,尚缺少一种系统的方法来从少量时间-排序数据统计分析基因表达测定。本发明人开发一种基于用最大似然方法和Akaike信息检验数(Akaike,1971)将线性样条函数与时间-排序数据(time-ordered data)拟合的策略。用Student t-检验评定基因表达测定值的显著性。这使得在考虑数据的统计学显著性的同时,可以从基因表达测定值推导信息。这种分析可以被视为建立基因调控网络的附加步骤。作为一个例子,本发明人重新分析了藻青菌PCC6803的基因表达测定值(Hihara,2001)。采用线性样条函数,可以推导信息以及测定被仅用倍数-变化方法所遗漏的数据。通过用可以得到的数据的子集重复本发明的分析,确定在各个时间点上需要多少测定值才能合理地评估线性样条函数。
2.方法Student t-检验评价测定的基因表达比是否显著地区别于1(unity)。如果对于特定基因来说,得出在所有时间点上所测定的表达比与1的差别都不显著的结论,可从以后的分析中排除该基因。可以对各个时间点分别使用Student t-检验来确定显著性水平。由于对每个基因进行多重比较,应当谨慎选择显著性水平值α。
定义H0(i)为对于特定基因在特定时间点ti的表达比为1的假设,而H0为对于特定基因在所有时间点的表达比都等于1的假设。如果用α表示否定假设H0的显著性水平,用α`表示否定假设H0(i)的显著性水平,那么α`与α以1-α=(1-α`)a相关,其中a为测定基因表达比的时间点的数量。通过扩展一阶泰勒级数的右边,将该方程转换为Bonferoni方法来调节显著性水平(又见Anderson和Finn,1996)。
以α`为显著性水平,对每个基因的各个时间点进行Student t-检验,结果无论H0(i)还是H0都应被否定。如果H0未被否定,结果是该基因不会显著地受试验处理影响,从而在进一步分析中不应当包括该基因。如果对于特定基因,虚(null)12假设H0被否定,结论是该基因显著地受试验处理影响。
用线性样条函数分析时间-排序数据接着,分析被发现显著地受影响的基因的瞬时基因表达反应。测定的基因表达比构成了小组时间-排序数据,用线性样条函数拟合这些数据。线性样条函数是由分段线性函数组成的连续函数,这些分段线性函数通过节点相互连接(Fiedman和Silverman,1989;Higuchi,1999)。然而立方样条函数更加常用,对少量数据点采用线性样条函数更为合适。实施例4图1给出具有节点t*1、t*2、t*3、t*4的线性样条函数的概念上的例子。本发明人想用形式为xj=g(tj)+εj的非参数回归模型拟合这些数据,其中g为线性样条函数,εj是平均数为零、方差为σ2的正态分布的随机自变量。
用最大似然方法来评估线性样条函数。特定tj的数据点xj的概率分布在实施例4方程3中表示。实施例4方程4给出该n个数据点的log-似然函数。通过最大化log-似然函数的σ2得到方差σ2的最大似然估计。从而得到实施例4方程5。
拟合模型依赖于节点q的数量。可以使用称作为AIC的Akaike信息检验数(Akaike,1971;以及Priestly,1994)来选择节点数量。对于每个q值,在如实施例4所述拟合线性样条函数后计算AIC的值,选择最小化AIC值的q值。q=2的情况对应线性回归。对于特定情形q=1,直线有效地拟合了该数据。如果发现对于某一特定基因,得到恒定函数(q=1)的最小AIC,可以认为该基因的表达水平没有受到试验操作的影响。基因表达数据通常是以表达比的形式给出。定义在零时的表达比为1。通过修改实施例4方程7可以容易地将该固定点纳入本发明的方法学中。通过选择实施例4方程17所示线性样条函数可以获得σ2最小值。
3.结果Student t-检验通过重新分析藻青菌PCC6803在突然暴露在强光(HL)下后所测定的基因表达图谱(Hihara等,2001),在此说明应用Student t-检验和拟合线性样条函数。对暴露于HL的藻青菌以及保留于低光照(LL)条件下的藻青菌,在零时、15分钟、1小时、6小时和15小时测定3079个ORF的表达水平。实施例4表1说明每个时间点的测定值的数目。来自cDNA表达测定的数据从KEGG数据库获得(Nakao等,1999)。
用作原始分析的数据(Hihara,2001)可能不同于被提交到KEGG的原始数据(Hihara,个人信息(personal communication))。此外,KEGG数据库中缺少在t=15分钟时间点的6个测定组中的两组。重新从原始数据中计算基因表达比得到与先前公开的结果接近的数值。
在从HL和LL的原始数据中减去背景信号强度后,运用整体归一化,计算HL与LL信号强度的比值来确定相对于对照条件(LL)的基因表达水平的相对变化。在倍数-变化分析中,如果基因表达水平被两个或多个因子之一改变,该基因被认为受HL作用。通过假定这些测定值标准偏差的大小,启发式地估计这种变化的统计学显著性(Hihara,2000)。
对各个基因的表达比分别进行Student t-检验的结果在表2中给出。在显著性水平α=0.001时,发现167个基因显著地受HL条件影响。在这167个基因中可以预期大约3个1类错误。与此相比,在最初分析中发现164个ORF受HL条件影响(Hihara,2001)。
通过倍数-变化分析psbD2基因(slr0927),结果表明该基因未显著地受HL诱导(Hihara,2001)。由于该基因已经被报道在藻青菌PCC7942中受HL诱导(Bustos和Golden,1992;Anandan和Golden,1997),所以上述结论不寻常。但是对psbD2基因在t=6小时的基因表达数据进行Student t-检验,得到p=3.3×10-5,表明该基因确实受HL影响。
应用线性样条函数分析接着用线性样条函数与测定的基因表达比拟合。节点数q在1到5之间,具有在零时间点时等于1的固定节点。当q=3和q=4时,线性样条函数的线性片段之间的节点的定位存在3种可能性。在实施例4图2中说明这些内容以及q=1、q=2和q=5的情况。表明可能的节点定位的数量随着最大节点数qmax指数性地增加。
在倍数-变化分析中,瞬时基因表达模式分成6类(Hihara,2001),在实施例4表3中列出。将线性样条函数与测定的基因表达数据拟合,提供一种比分类更灵活的方法来描述数据。此外,基因表达反应模式(pattern)的数字化描述是获得基因调控网络的重要步骤。
应用线性样条函数分析表达数据接着用一个例子说明AIC的使用。在倍数-变化分析中,发现苏氨酸合成基因thrC(sll1688)大约在1小时后被抑制。表4中列出不同组节点的AIC的计算值。在0、15分钟、1小时和15小时时获得节点的最小AIC。实施例4图3给出测定的基因表达水平,以及与数据拟合的线性样条函数。
对所有基因表达测定实施该处理,以基于不同基因对HL的时间-依赖反应其进行分类。至于哪些基因被包括在该分析中可以进行多种选择。在原始分析中,如果基因的表达水平处于3079个ORF中最低表达的2000个中(Hihara,2001),其从计算中被排除。另一种情况,如果Student t-检验表明基因未显著地受HL影响,也排除该基因。表5表示基因的数量,这些基因的测定的表达水平对应于这些不同情形的各个图谱。采用显著性水平q=0.001进行Student t-检验。
比较通过Student t-检验确定为显著地受HL影响的、具有倍数-变化分析结果的167个基因(Hihara,2001),其中164个ORF被确定。首先,将数据中存在无关值的基因从分析中排除。定义无关值为偏离特定时间点的数据的平均值两个以上标准偏差的数据点。只有一个基因被发现其测定的表达数据含有无关值;与表达数据拟合的线性样条函数是一条直线。其它基因表达水平都不能用直线描述,这与Student t-检验的结果一致。
接着,排除表达水平在最低的2000个之中的基因来避免使用受噪声支配的数据。同样的处理被用于倍数-变化分析(Hihara,2001)。在排除了这些基因后,保留显著受HL影响的107个基因。
该107个基因中的42个未在倍数-变化分析中被确定(Hihara,2001)。实施例4表6中列出这些基因,以及发现的各个基因的节点的定位。对于每个线性样条函数,计算所说明方差的百分比,作为对拟合程度的检测。作为一个例子,图4显示测定的基因表达比以及基因syIR(slr0329)的拟合线性样条函数,具有零时、15分钟、1小时和15小时的4个节点。在该42个基因中,基因xyIR(slr0329)具有所说明方差的最大百分比(98.7%)。在164个由倍数-变化分析确定的ORF中,按照q=0.001的Student t-检验,其中39个没有显著受HL影响。这些ORF在实施例4表7中列出。
最后,本发明人确定在各个时间点的测定值的数量是否足以可靠地确定线性样条函数的节点的定位。为此,重复线性样条函数的估计。为此,用所测定数据的子集重复线性样条函数的估计。然后,如果用数据的子集替代整组数据,计算多少基因的所估计的节点位点发生改变。表8给出每个时间点数值的4个、3个和2个数据点的平均数和标准偏差。
即便仅有2个数据点(在1小时的时间点)被排除并且每个时间点采用4个数据点,15%的情形估计的节点位置发生改变。因此,在一些实施方案中,每个时间点需要4个或者更多个数据点来可靠地从基因表达测定值中推导信息。
讨论本发明已经描述一种基于最大似然方法的策略来分析一组时间-排序测定值。对测定的基因表达数据进行Student t-检验,首先确定被测基因中的哪些显著地受试验处理影响。然后通过拟合线性样条函数来描述那些基因的表达反应。用Kike信息检验数(AIC)确定用于线性样条函数的节点数量。
在描述测定的基因表达中,采用线性样条函数比名称上的分类更具灵活性。而且,为了建立基因调控网络,期望的是从基因表达测定值中确定的基因反应以数字化形式获得。最后,节点的位置清楚表明基因表达水平显著变化的那些时间点,从而期望确定其生物功能。基于节点位置的基因表达水平的分类,可以通过建立考虑位点的线性样条函数大小的亚类来细化。例如,可以对具有3个节点的线性样条函数建立6个亚类,其中基因表达水平用(平缓,上升)、(平缓、下降)、(上升,平缓)、(下降、平缓)、(上升、下降)或(下降、上升)来描述。
应用线性样条函数技术来测定基因表达数据,可以确定显著地受试验处理影响的基因的瞬时表达反应模式。那些基因中42个基因的反应在较早的表达数据的倍数-变化分析中未被提出。而且表明164个基因中的33基因,在倍数-变化分析中发现的表达反应水平没有被Student t-检验确定为显著。在一些实施方案中,基因表达数据可能具有噪声,受到无关值干扰。然而,此处描述的Student t-检验和最大似然方法考虑噪声数据的统计学显著性,无关值的问题需要分别地解决。作为排除无关值的简单方法,计算各个时间点的数据的平均值和标准偏差,排除那些偏离平均数超过大约2个标准偏差的数据。
最后,用于可靠地拟合线性样条函数的各个时间点上所需的表达测定值的数量,通过排除一些数据点和重新拟合线性样条函数来确定。已经发现,如果每个时间点采用4个数据点,大约15%的情形不能可靠地估计节点位置。因此,建议每个时间点进行4次以上测定。
VI.用基因网络确定和评价药物靶点本发明人描述一种用基因网络确定和评价药物靶点新方法,该基因网络从cDNA微阵列基因表达数据估计得到。建立新的基因破坏和药物反应微阵列基因表达数据文库来阐明药物靶点。用两种微阵列基因表达数据来评估基因网络,然后确定药物靶点。评估的基因网络在分析药物反应数据中是有用的,这些信息不能从聚类分析方法中获得,聚类分析是基因表达分析的标准方法。在一些实施方案的建立基因网络中,在分析的不同步骤中采用布尔和贝叶斯网络模型,并利用它们的相对强度。用一个分析酿酒酵母基因表达和药物反应数据的实例来评价本发明的将基因网络信息应用于药物发现的策略。
1.引言聚类方法已经成为分析微阵列基因表达数据的标准工具。但是,无论从理论上还是从实践上看,这些方法都不能提供足以确定药物靶点的信息。本发明人提供了确定如何将评估的基因网络用于鉴定和评价目标基因来推定和开发新的治疗方法。为了发明目的,期望得到基因调控路径信息,并且采用布尔的和贝叶斯的网络建模方法来从基因表达图谱中推导基因网络。确定药物靶点的过程分成两部分。首先,确定药物作用基因。其次,搜索“目标”基因,在基因网络中该基因通常位于药物作用基因的上游。通过采用本申请描述的“虚拟基因(virtual gene)”基因技术,布尔网络模型可用于确定药物作用基因,并且贝叶斯网络模型用于搜索与被阐明的作用基因相关的可药化(druggable)基因靶点。将本发明的方法应用于新的酿酒酵母表达数据,这些数据来自120个基因破坏和对药物的多种剂量和时间应答的表达试验。
2.确定药物靶点的基因网络A.聚类方法聚类方法如分级聚类(hierarchical clustering)和自组织图(self-organizaingmaps)是生物信息学领域中用于基因表达数据分析的标准工具。Eisen专注于分级聚类,并且提供用于基因表达数据聚类分析的软件Cluster/Tree Vew。DeHoon等改进了该软件,特别是在K-均值聚类算法方面。
聚类方法仅通过表达模式的相似性提供基因群的信息。但是,期望还有额外的分级路径信息,不仅仅是聚类信息用于检测受试剂作用的药物靶点。通过实际数据分析说明用于药物靶向目的的聚类技术的局限性。
采用两种新方法从基因表达数据中评估基因网络。在该部分中,简单介绍这两种方法。对于这些算法的详细讨论,请参见参考文献部分的文章。
B.布尔网络为了评估布尔网络模型,可以将基因表达值分成两个级别,0(未被表达)和1(被表达)。假定u1,...uk为节点v的输入节点。V的态通过布尔函数fu(ψ(u1),...,ψ(uk))确定,其中ψ(u1)是u1的态(state)或者表达模式。如果具有时序(time series)基因表达数据,该态依赖于时间t,以及时间t的节点的态依赖于时间t-1时的输入值的态(state)。另一方面,假定具有得自基因破坏的基因表达数据。Akutsu等提出评估没有延时的布尔网络模型的理论和方法学。Maki等提供名为AIGNET的系统,用于评估基于布尔网络和S-系统的基因网络。本发明采用AIGNET系统评估布尔网络模型。
采用布尔网络模型的优点包括a)该模型简单,容易被理解。当这些数据具有充分精确度和信息时,布尔网络模型能够正确地测定亲本-子代关系和b)所评估的布尔网络模型可以直接应用于基因组对象网(Genome Object Net),用于生物路径模拟软件工具。布尔方法的缺点在于数据必须被化分成两个级别,而该量化丢失信息。而且用于离散的阈值为一个参数,并且需要用合适的检验数选择。
C.贝叶斯网络贝叶斯网络是大量随机变量的复杂关系的图形表示。在贝叶斯网络内容中,用节点的马尔可夫关系考虑定向非循环结构图。通过条件概率替代随机变量的联合概率来描述复杂现象。在本申请实施例5中有进一步描述。
Friedman等提出一种从基因表达图谱评估基因网络的方法。他们将表达值离散成3个值,并且以多项式分布作为贝叶斯网络的条件分布。但是,该方法没有解决为离散选择阈值的问题。
为选择最佳结构图,本发明人开发一种非参数回归模型,该模型提供不要求量化的解决方法和称为BNRC的新检验数。BNRC定义为用拉普拉斯逼近积分的结构图的后验概率的近似值。将该提出的方法应用于酿酒酵母基因表达数据,评估基因网络。这种方法的优点在于a)可以将微阵列数据作为连续数据集分析,b)该模型不仅可以检测线性结构关系而且还可以检测基因间的非线性相关性和c)被提出的检验数能够自动地优化模型中的参数和网络结构关系。
贝叶斯网络具有一些基于推导数学的优点,我们采用数学推导方法来构建贝叶斯网络。通过在分析中结合贝叶斯和布尔网络,从产生于逻辑联合来自破坏体和药物反应试验的表达数据的数据中构建调节作用的环状调节和多水平定向模型。因此,同时使用布尔网络和贝叶斯网络能够彼此克服缺陷,获得更可靠的信息。
3.应用于微阵列数据建立两个来自酿酒酵母基因表达图谱的微阵列数据文库。一个通过破坏120基因获得,另一个由对口服抗真菌剂的反应组成。(4种浓度和5个时间点)。从酵母基因组中选择735个基因用于确定药物靶点。在YPD中,314个基因定义为“转录因子”,其中98个已经研究了它们的控制机制。选择用于分析的735个基因的表达图谱数据包括受这98个“转录因子”控制的基因,这些来自除了核受体基因之外的5871个被测的基因种类,所述核受体基因在基因表达调控中发挥作用,是常见的药物靶点。
至于药物反应微阵列基因表达数据,将酵母培养物接种到抗真菌药物剂量为10、25、50、100mg的培养基中,在试剂加入后的5个时间点(0、15、30、45和60分钟)取出培养物的等分试样。此处时间0表示刚暴露于药物后观察的起始点。然后从这些试验中提取总RNA,用cy5标记RNA,用来自未处理细胞的cy3标记的RNA与之杂交,并将它们应用于全基因组cDNA微阵列,从而产生20个微阵列数据集作为药物反应数据。在本申请中,应用基因网络以140个微阵列阐明药物靶点。
4.结果A.聚类分析在确定药物靶点中,一种常用但存在问题的现有技术策略是采用聚类分析,通常还有基础干扰控制数据(base perturbation contral data)文库作为比较(0)。本发明具有两类微阵列数据,基因破坏和药物反应,允许将药物反应图谱与由破坏导致的基因表达模式比较。在聚类分析中,如果单个破坏体或者破坏体群的表达模式与特定药物反应微阵列之间显著而强烈地相似,表明药物与被破坏的基因可能起相同作用。而且,如果这种被破坏的基因具有已知功能性作用,可以获得更多关于对药物反应的信息。
不幸的是,即使总是具有这种试验结果的情况,也不能从聚类分析数据中获得直截了当的结果。图1表示本发明微阵列数据的相关矩阵的图像表示。结合两类数据,生成矩阵Z=(X∶Y),其中X和Y分别是药物反应和基因破坏微阵列数据的矩阵。在此,每一列表示从一个2微阵列中得到的表达模式,并且每一列用平均数0和方差1进行归一化。因此,实施例5图1描述了相关矩阵R=ZTz/p的信息,其中p是基因的数量。为了说明本发明方法,集中于735个基因来评估基因网络和确定药物靶点。
在实施例5图1中,浅颜色和深颜色分别表示正和负的高相关性。药物反应微阵列之间彼此高度相关,而与任何基因破坏微阵列之间的相关性低。在这种情况下,难以从聚类分析中确定基因破坏和药物反应之间的相互作用,而所述聚类分析在药物反应分析中是有意义的。进一步进行药物反应微阵列的分级聚类(cluster),但是这产生一个聚类,而不能从该分析中获取更多任何关于真正药物靶点的信息。当采用其它的距离测定和聚类技术时,该结果没有实质性改变。因此,从聚类技术中获得确定针对性的药物靶点的信息是困难的。
B.布尔网络分析为了克服现有技术聚类方法的缺点,用微阵列数据Z来评估基因网络,微阵列数据Z是通过组合基因破坏微阵列和药物反应微阵列产生。假设药物反应数据的条件为“虚拟基因”,例如条件100mg/ml和30分钟被设定为基因YEXP100mg30分钟。通过应用布尔网络模型,确定这些虚拟基因的子代基因,药物以世代顺序作用于这些子代基因。将5个或5个以上虚拟基因为亲本基因的基因视作为假定的药物作用基因,也就是,直接受虚拟基因(药物作用)影响的基因。但是,仅只有一个虚拟基因作为亲本基因的基因可能是最初药物作用基因,这依赖于给定药物的作用方式。在初步筛选受药物诱导表达影响的基因中,虚拟基因技术使布尔网络模型的使用比贝叶斯网络模型更加突出。
此外,倍数-变化分析能够提供与所提出的虚拟基因技术相似的信息。在特定实验条件下通过倍数-变化分析可以确定被作用的基因。但是,本发明的虚拟基因技术能够改进倍数-变化分析的结果。假设通过倍数-变化分析发现基因A和基因B受药物作用。倍数-变化分析不能考虑基因A和基因B之间的基线相互作用。也就是说,基因A和基因B之间存在调控路径基因A→基因B,基因B可能不是直接受药物作用。更确切地说,药物对基因A的作用可能导致对基因B的间接作用。虚拟基因技术通过利用基因破坏数据的信息考虑了这种相互作用,从而将搜索点减少到可能性更大的目标基因,这些基因产生可以获得的相互作用数据。
不能保证受试剂作用最大的基因就是那些被该试剂“药化”(drugged)的基因,也不能保证药化的靶点即代表用于发明新药的最具生物利用度和优势的分子靶。因此,即便确定了作用的可能分子模型,还期望寻找到位于调控网络中的药物作用基因上游的最可药化的(druggable)目标基因,然后筛选对这些靶点有药物活性的低分子量化合物。在评估的布尔网络中,可以将虚拟基因置于网络的顶端。因此,在该评估的布尔网络中发掘药物作用基因的上游信息是困难的,有时是不可能的。在这种情况下,本发明人采用贝叶斯网络模型有效地搜索药物作用基因的上游区域。
C.贝叶斯网络分析本发明人发现可以通过贝叶斯网络和非参数回归方法以及BNRC优化的策略来评估基因网络。我们使用了本文所描述得自破坏120个基因的酿酒酵母微阵列基因表达数据。从布尔网络分析中有效地发现候选的药物作用基因。可药化的基因是与这些药物作用基因相关的药物靶点,这些可药化的基因是为了开发新的先导药物要确定的。通过贝叶斯网络方法在评估的基因网络中搜索位于药物作用基因上游的可药化基因。采用从敲除表达文库中获得的网络调控数据的贝叶斯模型,搜索药物作用目标基因的上游,其具有已知的对药物作用靶点表达的调控关系。例如,将核受体基因视为可药化基因,因为a)已知核受体蛋白是有用的药物靶点,总共代表当前市场上药物的靶点的20%以上。b)核受体参与转录调控作用,这些转录作用在cDNA微阵列试验中被直接测定。
实施例5图3表明部分网络,包括药物作用基因(底层)、可药化基因(顶层)和媒介基因(中间)。当然,如果发现更多媒介基因,会发现更多从可药化基因到药物作用基因的路径。由于使用贝叶斯网络模型,可以确定边缘强度和选择更可靠的路径。贝叶斯网络模型的优点之一在于寻找适当的可药化基因。在实施例5图3中,在圈中的可药化基因直接与药物作用基因相连,其它每个可药化基因具有一个媒介基因。从图3可以确定每个药物作用基因的可药化基因,例如,发现了实施例5表1所述MAL33和CDC6的可药化基因。
5.讨论本发明描述了用基因网络的计算模型确定和评价药物靶点的新策略。布尔网络和贝叶斯网络对于从微阵列基因表达数据评估基因网络是有用的。同时使用这两种方法能够比上述任何一种方法获得更多可靠信息。布尔网络适合于通过使用在此提及的虚拟基因技术来确定药物作用基因。贝叶斯网络模型能够提供关于药物作用基因上游的信息,从而获得一组候选可药化基因。基于精确地结合使用这两种网络方法,建立了本发明的新策略。在该策略中可以清楚地看出各个网络方法的优势,而且该综合方法在药物靶点的鉴定和确认中为用于基因网络推导的生物信息学技术的实际应用提供了方法学基础。
VII.药物靶点发现以及用基因调控网络评价用来自基因干扰变体细胞系(gene perturbation variant cell lines)的全基因组表达文库开发的基因调控网络可以被用于快速和有效地确定药物或者先导化合物分子作用的分子机制。本发明人开发由500个以上的酵母菌株的全基因组cDNA阵列数据组成的扩展酵母基因表达文库,每个菌株具有单个基因破坏。应用这些数据,结合口服抗真菌剂灰黄霉素的剂量和时程表达试验,用布尔和贝叶斯网络发现技术来确定其表达受该药物影响最大的基因。灰黄霉素的精确的分子靶点先前未知,它妨碍酵母中有丝分裂纺锤体的形成。本发明的系统用于直接揭示CIK1为存在灰黄霉素时的首先(primarily)被作用的目标基因。删除CIK1,通过基因网络发现所确定的首先被作用的靶点对有丝分裂纺锤体形成的形态学作用与药物所产生的相似。采用来自表达文库的基础分级数据(base bierarchical data)和非参数贝叶斯网络建模,可确定可替代的配基依赖转录因子和CIK1上游的其它蛋白,它们可以作为潜在的可替代分子靶点来诱导与灰黄霉素相同的分子反应。该基于网络的药物发现方法可以显著地减少作出合理的(rational)药物靶向决定所需的时间和资源。
1.引言合理的药物设计方法学以前集中于优化对应于预先确定的分子靶点的小分子。随机化与靶点结合的先导物进行表型筛选来选择靶点是当前药物发现中的盛行方式,但是即便出现了大规模可获得的基因组信息和高通量筛选方法,这些方式也没有产生更加有效和精确的靶点选择方法(1,2,3)。基因组序列和全基因组微阵列的可获得性以及基因网络推导计算方法的最新进展为药物靶点选择提供了新的合理方法,这些方法考虑了基因组中分子间的全网络调控相互作用。采用各种计算推导技术(6,7),可以从基于破坏体的表达数据(4,5)中获得基因调控网络数据的精确模型。本申请说明这些基因调控信息如何被用于快速确定分子网络和受特定化合物作用的基因靶点。同样的信息用于确定和选择可替代的可药化分子靶点,这些靶点在基因表达调控级联中位于药物靶向分子的上游或下游。
本发明人开发了一种基因调控网络挖掘的迭代药物靶点(gene regulatorynetwork-driven iterative drug target)的发现方法。在该方法中,首先对单基因破坏体细胞进行大量基因表达试验。这些信息用于建立分级基因表达控制的计算推导图。分级调控信息用作评价药物反应试验和产生分子作用机制假说的基础。从参考文献以及从对阐明调控子网络的生物学试验中得到的信息被用于在选择药物靶向的候选分子之前推定和确认结果。
临床应用中最有效的药物包括阿司匹林和其它大众药物,都尚未被合理设计来与特定分子靶点作用。因此,即便用这些药物获得了预期临床效果或者表现,然而隐含的分子作用机制和药物副作用机制仍然是未知的。全基因组表达试验已经被证明在确定可替代基因和受药物作用的路径中是有用的(8,9,10),但是在缺少关于潜在药物靶点的先验信息时,用标准基因表达分析方法如聚类分析来确定多种药物的初始分子靶点是不切实际的,所述药物作用于许许多多基因。本申请举例说明使用来自全基因组表达文库的分级基因表达调控网络和基因网络建模技术以及药物反应表达试验来确定先前未知的常规抗菌剂灰黄霉素的分子靶点。
灰黄霉素是被广泛开方使用的口服抗真菌剂,最初被指定用于严重的头发和指甲的真菌感染。虽然灰黄霉素的分子作用是未知的,但是该药物破坏真菌有丝分裂的纺锤体结构,导致中期休止。
2.方法微阵列试验在本发明中使用的酵母菌株是BY4741。为了监控基因表达图谱,细菌细胞在30℃的YPD(%蛋白胨、1%酵母提取物和2%葡萄糖)中预生长到指数生长中期,然后用灰黄霉素作用,将灰黄霉素加入培养基中至浓度为0、10、25、50和100mg/ml。受作用的细菌细胞在加入灰黄霉素后的0、15、30、45和60分钟时被收获并用于RNA提取。用热苯酚法提取总RNAs。
布尔网络推导算法除了别处报道(11)的用于基因破坏试验的布尔方法之外,本发明人从一组结合了破坏体的表达数据的药物处理试验中建立基因表达矩阵E。在药物产生干扰的情况下,设置一个“虚拟”基因来代表与基因破坏相似的药物作用。这样就可以使用为基于破坏的数据而设计的标准的破坏矩阵算法。
贝叶斯网络推导算法采用本申请其它部分报道的算法和方法进行非参数调控和贝叶斯网络来确定CIK1上游的调控子网络。
数据归一化通过cDNA微阵列分析,从20个药物处理菌株中测定5871种mRNA种类的量。Cy3和cy5之间的荧光强度的差别产生表达数量比的偏差。归一化(normalize)各个表达图谱的表达数量比。在各个点样区块(spotted block),比例偏差有一个固定趋势,然后计算线性回归将各个区块比例的平均值归一化为1.0。比例的对数值用于表示标准表达水平,从而得到比例的对数值,计算这些对数值的平均数和标准偏差(见表1)。来自UME6(YDR207C)破坏体表达阵列的所有点样基因的表达水平的标准偏差(SD)为0.4931,所述UME6被确定为YPD(33)破坏体中的“总调控子”,因此,本发明人意识到在总SD超过0.5的阵列数据中存在不可接受量的误差,并从分析中排除掉这些试验。
选择用于建模的基因在YPD中,314个基因定义为“转录因子”,其中98个已经研究了它们的控制机制。包括该98个“转录因子”所控制基因的552个基因的表达图谱数据选自5871个图谱。从而从表达图谱数据集构建基因调控网络,所述数据集建立在120个基因破坏试验和20个药物处理试验的552个基因的表达值上。
3.结果将酵母培养物接种到剂量为10、50和100mg的10ml DMF中,在加入灰黄霉素后的5个时间点(0、15、30、45和60分钟)取出培养物的等分试样。然后从这些试验中提取总RNA,用cy5标记RNA,用来自未处理细胞的cy3标记的RNA与之杂交,将它们应用于全基因组cDNA微阵列。552个基因中的183个基因被下调超过2δ的阈值,这些基因在药物处理酵母和正常酵母中的表达不同(附图7a)。标准分级聚类方法(附图7b)用于组合来自药物反应和来自基因破坏试验的表达文库,将基因聚类分成两个主组;第一组为受灰黄霉素作用的基因和第二组为受破坏作用的基因。在灰黄霉素聚类中,根据剂量或时程对基因进一步分组。但是,聚类没有揭示任何基因表达模式,即没有明显指示由特定基因和药物灰黄霉素产生的相关联的调控。除了抗真菌剂影响仅一个离散基因和最小限度表达的情况,才可能得到这种结果(附图7b)。
但是,采用基因调控网络模型,结合药物干扰数据,能够得到药物与转录组中的基因的相互作用的分级基因调控图。为了产生该基因网络药物干扰数据,首先建立542个单基因破坏突变体的全基因组表达文库。从该文库中选出的120个阵列数据与产生自灰黄霉素试验的阵列矩阵合理地结合。然后设计用于基因网络说明的布尔方法(11,12)被应用于各个时程的联合表达矩阵,产生各个剂量和时间点的每个试验的分级调控图。建立各个剂量和时间点试验的联合布尔调控子网络。选择布尔算法是由于其适合于处理联合矩阵和其能够处理成环的调控过程以及容易构建具有多重调控级别的分级定向结构图。
从这些数据能够确定相对于次级级联反应(secondary cascade)的一级药物作用,所述次级级联反应由那些初始干扰调控事件启动。
通过评价来自各个时间点和剂量不同试验的布尔数据确定8个基因,这些基因在各时间点和药物浓度都一贯和显著地作为第一效应被抑制。在这些基因中,CIK1在整个实验过程中具有最大的抑制效应。CIK1编码的蛋白在酵母蛋白质组数据库(YPD)中被描述为纺锤体极体的卷曲螺旋蛋白,该蛋白参与纺锤体形成和核融合的集会步骤(胞核转移)(13)。由于已知灰黄霉素的作用是影响有丝分裂纺锤体形成,其与CIK1的功能一致(14),用具有CIK1基因被破坏的酵母菌株和受灰黄霉素作用的酵母菌株进行病理学检验。虽然用常规生理可接收量灰黄霉素处理和破坏CIK1都不是致命的,但这两种培养物都表现出相似的形态学差异和生长特性。而且镜检灰黄霉素处理的酵母和CIK1被删除的酵母中的有丝分裂纺锤体结构,都表现出纺锤体和周边组织结构的极其相似的改变(15)。
本申请描述的方法清楚地证明了组合的表达阵列和应用基因网络技术的计算方法用于快速确定和评价特定化合物对细胞作用的分子机制。采用这些方法将有助于使在后基因组中的药物开发的靶点选择过程合理化,也有助于提高发现效率和降低制药行业的开发风险。该方法可以被进一步应用于其它生物发现和农业化学打靶。本发明人的实验室目前正对人类和其它的生物系统重复实施这种发现模型。更多的描述可见2002年7月12日申请的编号为60/395,756的美国临时专利申请,在此全文引作参考。
参考文献1.Smith.A.药物发现的筛选前沿问题.自然(Nature),418,453-459(2002)。
2.Aherne,G.W,McDonald,E.和Workman,P.大海捞针高通量筛选有益健康的原因,乳癌研究(Breast Cancer Res),4,148-154(2002)。
3.Willins,D.A.,Kessler,M.,Walker,S.S.,Reyes,G.R.和Cottarel,G.抗真菌药物发现的基因组策略-从基因发现到化合物筛选,Curr PharmDes.8,1137-1154(2002)。
4.Hughes,T.R.等,借助表达图谱纲要的功能发现.细胞(cell),102,109-126(2000)。
5.Glaever,G.等,酿酒酵母基因组的功能图谱,自然(Nature),418,387-391(2002)。
6.Friedman,N.,Linial,M.,Nachman,I.和Pe`er,D.用贝叶斯网络分析表达数据,计算生物学杂志(J Comput Biol),7,601-620(200)。
7.Somogyi,R.和Sniegoski,C.A.给基因网络的复杂性建模了解多基因和多向调控,复杂性(Complexity).1,45-63(1996)。
8.Marton,M.J.等,用DNA微阵列药物靶向评价和确定次级药物靶向作用,自然医学(Nature Medicine).4,1293-1301(1998)。
9.Heller,M.J.DNA微阵列技术设备、系统和应用,生物医学工程研究年报(Ann Res Biomed Eng),4,129-153(2002)。
10.Reynolds MA.微阵列技术GEM微阵列和药物发现,印度微生物生物技术杂志(J Ind Microbiol Biotechnol),28,180-185(2002)。
11.Akutsu,T.,Miyano,S & Kuhara,S.基于矩阵乘法和指纹功能的用于确定布尔网络和相关生物学网络的算法,计算生物学杂志(J.Comput Biol).7,331-343(2000)。
12.Akutsu,T.,Miyano,S & Kuhara,S.推导基因网络和代谢路径之间的数量关系,生物信息学(Bioinformatics),16,727-734(2000)。
13.Rose,M.D.酿酒酵母中的核融合,细胞发育生物学研究年报(Ann.Res.Cell.Dev.Biol),12,663-695(1996)。
14.Cottingham,F.R.,Gheber,L.,Miller,D.L.和Hoyt,A.M.酿酒酵母有丝分裂纺锤体马达的新作用,洛克菲勒大学出版社(The Rockefeller UniversityPress).147,335-349(1999)。
15.Manning,D.B.,Barrett,J.G.,Wallaces,J.A.,Granok,H.和Snyder,M.Kar3p的Kinesin相关蛋白通过两个伴生蛋白Ciklp和Viklp的不同调控,洛克菲勒大学出版社(The Rockefeller University Press).144,1219-1233(1999)。
16.Imoto,S.,Goto,T.和Miyano,S.用贝叶斯网络和非参数回归估计基因之间的基因网络和功能结构,Pac Symp Biocomput,175-186(2002)。
VII.用于发现网络关系的系统及其用途在其它实施方案中,提供用于评估基因网络的方法,包括本申请描述的实施例6。因此,期望的系统可以包括这样的系统,其中(1)当基因组结构关系被阐明时,人们可收集使网络能够被阐明的试验数据;(2)测定所有相关基因的数据;(3)使用能使大量实验实现的多种方法,例如基因芯片;(4)在使用干扰来获得大量归一化数据后,测定输出值;和(5)确定遗传关系分析实施例6图1给出这种系统的例子。实施例6图2示意性地描述从有机体中获得被表达基因的微阵列数据的方法。实施例6图3示意性地描述通过比较来自突变细胞(破坏体)的基因产物(如RNA)的量与正常细胞(野生型)的基因产物的量,估计和量化基因表达的方法。在突变细胞中,特定基因被破坏。
在一些实施方案中,用实施例6图4提供的技术路线实现大规模网络分析。实施时程研究,评价研究中的各个基因的表达。生成布尔网络模型和该网络的动态模型。绘制正的和负的相互作用,从而建立基因网络。多水平绘图方法被用于将破坏(或改变的表达)各个基因对其它被研究的基因的作用相关联。
参考文献的引入在本申请中引用的所有专利、专利申请和参考文献在此全文引作参考。
前述的本发明的实施方案被提供用于说明和描述的目的。而不是穷举或者将发明限定在所披露的确切情形中。许多更改或变化对本领域技术人员来说是显然的。选择和描述实施方案的目的在于最好地说明发明的原理和实际应用,从而使本领域技术人员能够理解发明和各种实施方案和适合于预期的特定使用的各种更改。这意味着本发明的保护范围由后面的权利要求及其等价物定义。
具体实施例方式
下列实施方案被提供用于说明本发明的具体实施方式
。不背离本发明精神时,可以使用本专利申请所教导的其它特定申请。可使用这些方法的其它更改,并且被认为包含在本发明的保护范围中。
实施例1通过使用贝叶斯网络和非参数回归估计基因之间的基因网络和功能结构关系本发明人提供一种用贝叶斯网络从基因表达数据中构建基因网络的新方法。采用非参数回归来捕获基因间的非线性关系并得到在普遍情形下选择网络的新检验数。从理论上讲,本发明人提出的理论和方法包括了已有的基于贝叶斯网络的方法。将本发明的方法应用于酿酒酵母细胞周期数据13,20,说明与前述方法相比本发明方法的有效性。
1.引言遗传工程的发展产生大量有价值的数据,如微阵列基因表达数据。对基因间的关系分析引起人们对分子生物学和生物信息学领域的极大关注。但是,由于数据的维度和复杂性的原因,发掘隐藏在噪声(noise)中的结构关系不是一件容易的事。为了从生物数据中提取有效信息,从统计学角度考虑,需要发展理论和方法学。本发明的目的在于建立更清楚地了解基因间关系的新方法。
在微阵列基因表达数据的分析中,贝叶斯网络被用于从图像-理论方法3,4,5中构建基因网络。Friedman和Goldszmidt(1998)12提出了一种用贝叶斯网络构建遗传连接的有益方法。他们离散基因的表达值,并且考虑到拟合基于多项分布的模型。但是仍存在选择离散(disctritize)用的阈值的问题,而不仅仅是根据试验。阈值的确说明结果的实质变化,而不适当的阈值产生不正确结果。另一方面,最近Friedman等(2000)13指出离散会导致信息丢失。可将表达数据用作连续值,因此他们采用基于线性回归的高斯模型。但是,这种模型仅仅能够检测线性相关,而不能获得充分结果。
在本申请中,本发明人提出了利用贝叶斯网络开发了用于构建基因网络的新方法。不但为了捕获基因之间的线性关系,而且为了捕获非线性结构关系,本发明人采用了具有高斯噪声11,13,21,22的非参数回归模型。非参数回归被开发用于在预先缺乏功能性关系知识的情况下探索期望反应的复杂的非线性形式。鉴于贝叶斯网络的新结构,需要一个估计模型的合适的检验数。本发明人从贝叶斯统计中得到了新的检验数。通过使用被提出的方法,克服了前述方法的缺点并且获得更多信息。此外,本发明的方法包括作为特例的前述方法。我们通过分析酿酒酵母细胞周期数据来评价本发明的方法。
2.贝叶斯网络和非参数回归设X={X1,X2...,Xp}T为一个p-维的随机向量,设G为定向非循环结构图(directed acyclic graph)。在贝叶斯网络结构中,将一个基因视为一个随机变量,并把联合概率分解成条件概率的乘积,也就是P(X1,X2,...,Xp)=P(X1|P1)P(X2|P2)×...×P(Xp|Pp),(1)其中,Pj=(P1(j),P2(j),…,Pqj(j))T]]>是结构图G中Xj的亲本变量的qj-维向量。
假设有随机向量X的n个观察值x1,...,xn,Pj的观察值用p1j,...,pnj表示,其中pij是一个具有第k个元素pik(j)的qj-维向量,k=1,...,qj。例如,设Xn是一个n×p矩阵,其中Xn=(x1,...,xn)T=(x(1),...,x(p))=(xij)i=1,...,n;j=1,...,p,xi=(xi1,...,xip)T,x(j)=(x1j,...,xnj)T,x(j)=(x1j,...,xnj)T和xiT是向量xi的转置矩阵。如果X1有一个亲本向量P1=(X2,X3)T,令p11=(x12,x13)T,...,pn1=(xn2,xn3)T。
当用密度函数f(xi1,xi2,...,xip)=f1(xi1|pi1)f2(xi2|pi2)×...×fp(xip|pip)置换(1)中的概率值P时,立刻发现方程被求解。然后所有需要作的是考虑如何构建条件密度函数fj(xij|pij)(j=1,...,p)。
在本申请中,用非参数回归模型xij=m1(pi1(j))+m2(pi2(j))+…+mqj(piqj(j))+ϵij,i=1,…,n;j=1,…,p,]]>来捕获xij与pij=(pi1(j),...,piqj(j))T之间的关系,其中mk(k=1,...,qj)是从R到R的光滑函数,并且εij(i=1,...,n)通常独立(independently and normally)依赖于平数0和方差σj2。对于函数mk,假设
mk(pik(j))=Σm=1Mjkγmk(j)bmk(j)(pik(j)),i=1,…,n;k=1,…,qj,---(2)]]>其中{b1k(j),...,bMjk(j)k}是已述组(set)的基本函数(如傅立叶级数、多项式函数(polynomial bases)、回归样条函数、B-方样函数和小波函数等),系数γ1k(j),...,γMjk(j)k为未知参数,而Mjk是基础函数的数量。
那么,非参数回归模型可以写成概率密度函数,形式为fj(xij|pij;γj,σj2)=12πσj2exp[{xij-Σk=1qjΣm=1Mjkγmk(j)bmk(j)(pik(j))}22σj2],---(3)]]>其中γj=(γj1T,…,γjqjT)T]]>是具有γjk=(γ1k(j),…,γMjkk(j))T]]>的参数向量。如果变量Xj没有亲本向量,我们考虑模型基于平均数为μj和方差为σj2的正态分布函数上。
最后,使贝叶斯网络模型建立在具有高斯噪声f(xi;θG)=Πj=1pfj(xij|pij;θj),i=1,…,n,]]>的非参数回归模型上,其中θG=(θ1T,...,θpT)T是包含在结构图G中的参数向量,θj是模型fj,即θj=(γjT,σj2)T或θj=(μj,σj2)T中的参数向量。
3.用于选择结构图的所提出的检验数设π(θG|λ)为具有超级参数向量λ的未知参数向量θG的先验分布,令logπ(θG|λ)=0(n)。通过在参数区域上进行积分来获得数据Xn的边缘概率,选择具有最大后验概率的结构图GπG∫Πi=1nf(xi;θG)π(θG|λ)dθG,---(4)]]>其中πG是G的先验概率。Friedman和Goldszmidt(1998)12认为多项式分布为贝叶斯网络模型f(xi;θG),还提出参数θG具有Dirichlet先验分布。如此,Dirichlet先验分布是一种联合先验分布,后验分布属于相同类型的分布。得到了(4)中积分的闭式解,他们称之为选择结构图的BDe分值6,16。由于BDe分值限于多项式模型,而本发明人提出了在更普遍和更多种情形下的选择结构图的检验数。
构建基于(4)的检验数的关键问题是如何计算积分。可以考虑采用一些方法计算积分,如马尔可夫链蒙特卡罗(Markov chain Mote Carlo),本申请采用拉普拉斯逼近来求积分7,17,23。Xn的边缘概率的拉普拉斯逼近为∫Πi=1nf(xi;θG)π(θG|λ)dθG=∫exp{nlλ(θG|Xn)}dθG]]>=(2π/n)r/2|Jλ(θ^G)|1/2-exp{nlλ(θ^G|Xn)}{1+Op(n-1)}]]>其中r是θG的维数,lλ(θG|Xn)=1nΣi=1nlogf(xi;θG)+1nlogπ(θG|λ),]]>Jλ(θG)=-∂2{lλ(θG|Xn)}/∂θG∂θGT]]>而 是1λ(θG|Xn)的众数。得到选择结构图的检验数BNRC。
BNRC(G)=-2log{πG∫Πi=1nf(xi;θG)π(θG|λ)dθG}]]>=-2logπG-rlog(2π/n)+log|Jλ(θ^G)|-2nlλ(θ^G|Xn)---(5)]]>选择最佳结构图以使检验数BNRC(5)是最小的。
该检验数来自logπ(θG|λ)=0(n)。如果logπ(θG|λ)=0(1),众数 等价于最大似然估计,MLE,而该检验数产生于贝叶斯信息的检验数,称作除去较高阶限制(higher order terms)0(n-j)(j≥0)的BIC19。Konishi(2000)18提出根据Kullback-Leibler信息和贝叶斯方法构建模型选择检验数的一般框架。
如果先验密度π(θG|λ)被分解成先验密度与θj的乘积,π(θG|λ)=π1(θ1|λ1)×...×πp(θp|λp)。因此,(5)中的 和 分别为Σj=1plλ(j)(θ^j|Xn)]]>和Σj=1plog|-∂2lλ(j)(θj|Xn)∂θj∂θjT|,]]>其中
lλ(j)(θj|Xn)=1nΣi=1nlogfj(xij|pij;θj)+1nlogπj(θj|λj).---(6)]]>如下通过结构图的局部分值(local score)来获得BNRC(5)定义第j个变量Xj的局部BNRC为BNRCj=-2log{πLj∫Πi=1nfj(xij|pij;θj)πj(θj|λj)dθj},---(7)]]>其中πLj是与Xj相关的局部结构关系的先验概率,对BNRCj采用拉普拉斯方法,通过BNRC=-2logπG+Σj=1p{BNRCj+2logπLj}]]>得到BNRC。
由于结构图是以非环状(acyclic)方式建立的,所以选择最终结构图为BNRC的最小值,且不需要最小化每个局部BNRCj。
4.用BNRC估计曲线和相关结构关系在本章节中以更详细方式说明本发明方法。本发明人提出的方法的关键点是采用非参数回归和来自贝叶斯统计的用于选择结构图的新检验数。
对于第2章节的非参数回归,采用B-样条函数作为(2)中的基本函数。图1是具有等距节点t1,...,t10的3次B-样条函数的例子。设置将区域[mini(p(ik)(j)),maxi(p(ik)(j))]分成Mjk-3个等距间隔的节点,并设3次Mjk的B-样条函数。
图1 图16个3次B-样条函数的例子。t1,...,t10称为节点。这些节点是等间距的。假设参数向量θj的先验分布为
πj(θj|λj)=Πk=1qjπjk(γjk|λjk)]]>每个先验分布πjk(γjk|λjk)是奇异Mjk变量(singular Mjk variate)的正态分布,表示为πjk(γjk|λjk)=(2πnλjk)-(Mjk-2)/2|Kjk|+3/2exp(-nλjk2γjkTKjkγjk),---(8)]]>其中λjk是超级参数,Kjk是Mjk×Mjk的矩阵,γjkTKjkγjk=Σl=3Mjk(γlk(j)-2γl-1,k(j)+γl-2,k(j))2,]]>而|Kjk|+是Mjk-2个非零的Kjk特征值(eigenvalue)。分值BNRCj(7)可以如下获得BNRCj=-2logπLj-2Σi=1nlogf(xij|pij;θ^j)-2Σk=1qjlogπk(γ^jk|λjk)]]>+log|-∂2lλ(j)(θ^j|Xn)∂θjθjT|-(Σk=1qjMjk+1)log(2πn-1),---(9)]]>其中而θ^j=(γ^jT,σ^j2)T]]>是(6)中定义的1λ(j)(θj|Xn)相对于特定λjk的众数。在计算方面,通过下式逼近(9)中的黑塞(Hessian)矩阵的确定值的对数值Σk=1qj{log|BjkTBjk+nσ^j2λjkKjk|-Mjklog(nσ^j2)}-log(2σ^j4)]]>其中Bjk是一个用bjk(pjk(j))=(bjk(j)(pjk(j)),...,bMjk(j)k(pjk(j)))T的通过Bjk=(bjk(p1k(j)),...,bjk(pnk(j))T定义的n×Mjk矩阵。因此,结合(3),(8)和(9),BNRCj的结果为BNRCj=Cj+(n-2qj-2)logσ^j2]]>+Σk=1qj{nβjkσ^j2γ^jkTKjkγ^jk+log|Λjk|-(Mjk-2)log}βjk]]>其中βjk=σj2λjk为超级参数,
Cj=-2logπLj+(n+M-j-2qj)log(2π)+n-log2]]>-2(M-j-qj)logn-Σk=1qjlog|Kjλ|+,]]>Λjk=BjkTBjk+nβjkKjk,]]>M-j=Σk=1qjMjk]]>采用修合算法(backfitting algorithm)15,当给定βjk值时,可以得出众数 修合算法表达如下步骤1 初始化γjk=0,k=1,...,qj,步骤2循环k=1,...,qj,1,...,qj,1,...
γjk=(BjkTBjk+nβjkKjk)-1BjkT(x(j)-Σk′≠kBjk′γjk′)]]>步骤3 继续步骤2,直到得到满足合适的收敛检验数。众数 由σ^j2=||x(j)-Σk=1qjBjkγ^jk||2/n]]>得到。
另外,众数 和 依赖于超级参数βjk,必须选择最佳βjk值。在本发明的方法中,自然地会选择最佳βjk值作为BNRCj的最小值。
由于B-样条函数系数向量γjk通过最大化(6)估计。(6)的众数与所罚似然估计是相同的21,24,可以将超级参数λjk或βjk视为所罚似然中的修匀(smoothing)参数。因此,超级参数在将曲线与数据拟合中起重要作用。
X1=X22+2sin(X5)-2X7+ϵ1]]>X2={1+cxp(-4X3)}-1+ε2X3=ε3,X6=ε6,X9=ε9X4=X52/3+ϵ4,X5=X3-X62+ϵ5]]> X8=cxp(-X4-1)/2+ε9X10=cos(X9)+ε10-图2蒙特卡罗模拟(左边)真实,(右边)估计。
5.计算试验在分析真实数据前,采用蒙特卡罗模拟来检验本发明方法的特性。数据来自变量之间的虚拟(artificial)结构图和结构关系(图2),结论总结如下提出的检验数BNRC能够极好地检测数据间的线性和非线性结构关系。但是BNRC分值倾向于扩大结构图。然后在分析真实数据时,考虑采用称作为AIC1,2的Akaike信息检验数并同时使用两种方法。AIC最初作为评价由最大似然方法估计的模型的检验数被引入。但是本发明方法的估计与最大所罚似然估计相同,与MLE不同。如此,修正过的AIC10为AIC=-2Σi=1nlogf(xij|pij;γ^j,σ^j2)+2(Σk=1qjtrSjk+1),]]>其中Sjk=Bjk(BjkTBjk+nβ^jkKjk)-1BjkT.]]>Sjk轨迹表明拟合曲线的自由度,这十分有用。也就是说,如果trSjk接近2,可以认为相关是线性的。使用BNRC和AIC来决定是否增加亲本变量。采用这种模型,所估计的结构图和结构关系与真实模型接近。
本发明分析Spellan等(1998)20和Friedman等(2000)13论述的啤酒焦母细胞周期数据。这些数据从800个基因和77个试验收集得到。
图3.CNL2、CDC5和SVS1的BNRC分值。
由于缺少认为大结构图是不可接受的理由和关于真实结构图规格的信息,设定先验概率πG为恒量。由20个B-样条函数组成非参数回归。事实上,B-样条函数的数量也是一个参数。然而,我们使用了稍大数量的B-样条函数,超级参数控制着拟合曲线的光滑度,对于对应不同数量B-样条函数的拟合曲线间的差别不能从视觉上发现。
分析结果总结如下图3表示通过一个基因预测CLN2、CDC5和SVS1时的BNRC分值。产生较小BNRC分值的基因给出了目标基因的较好说明。观察哪些基因与目标基因相关,找出其所强烈依赖的基因组群。事实上,可以用这些信息构建简明网络。可以将最佳结构图视为考虑相互作用效应的简明网络的修订版。如果基因之间线性相关,当亲本-子代关系被反转时,BNRC分值仍是好的。因此,特别是当依赖关系接近线性时,结构图中的因果关系方向不精确(strict)。本发明的结果基本上与Friedman等(2000)13的结果一致。但是,存在部分不同之处。一些基因如MCD1、CSI2和YOX1等构成了Friedman等的结果。这些基因中的大部分被报道其重要作用。基因间的大量关系近乎线性。但是本发明发现了线性模型难以发现的非线性相关。图5表示与基因相关的估计结构图,这些基因根据其进入细胞周期的过程及其邻近基因进行分类。在图5中省略了一些分支,但是差不多给出了重要信息。就本发明人和Friedman等13给出的网络我们确定了亲本-子代关系,并且可以看出这两个网络是相似的。特别地,本发明的网络包括了Friedman13等报道的典型关系。至于这两个网络之间的不同是我们关注SVS1的亲本基因。Friedman13等以CLN2和CDC5作为SVS1的亲本基因。一方面,本发明结果给SVS1提供CSI2和YKR090W。本发明检验了这两种结果的不同。毕竟,从BNRC和AIC上看,本发明的候选亲本基因比Friedman13等的候选亲本基因更恰当。原因可能是离散的影响,因为本发明的模型与图4的两种情形恰当地吻合。特别地,可以得出与Spellman数据中的其它基因相比,CDC5对SVS1的作用较弱的结论(又见图3)。事实上,作为SVS1的亲本基因,CDC5的BNRC分值的排序为第247位。考虑到上述影响因素,本发明方法提供可理解和有用形式的有价值信息。
图4.细胞周期数据和光滑估计(a)和(b)Friedman等(2000)13,BNRC=57.71,AIC=167.96;(c)和(d)本发明方法,BNRC=32.53,AIC=140.16。
6.讨论本发明提出通过应用贝叶斯网络和非参数回归从微阵列基因表达数据中评估基因网络的新方法。本发明人理论性地得到用于选择结构图的新检验数,并通过分析细胞周期数据证实其有效性。本发明方法的优点主要为可以将表达数据用作连续值。不仅检测线性相关,还检测非线性结构关系和可视化其功能结构关系使之容易理解。完全自动的搜索可获得最佳结构图的制备。
本发明还指出Friedman等13的方法保留了未知参数如用于离散的阈值和Dirichlet先验函数(priors)的超级参数,所述函数通过反复试验来选择并且从某种意义说没有被优化。另一方面,本发明方法可以根据所提出的具有良好理论基础的检验数自动而恰当地估计任何参数。此外,本发明方法包括了Friedman等13的方法作为特例。
将来工作中的问题(1)本发明采用基于高斯分布的统计模型,但是本发明得到适合更加普遍情形的检验数BNRC。事实上,可以构建基于其它统计模型的结构图选择检验数。(2)无关值导致奇异的结果是可能的。因此,开发强大的检测无关值的方法和技术很重要。(3)可能可以用自举方法9测定联合的强度(intensities of the unions)。本发明人将在以后的文章中探讨这些问题。
在其它实施方案中得到更加普遍情形时的检验数BNRC。通过这些实施方案,可以根据其它统计模型建立结构图选择检验数。
图5细胞周期数据结果。
参考文献1.H.Akaike,in Petrov,B.N.和Csaki,F.eds.,2nd Inter,信息理论讨论(Symp.on Information Theory),Akademiai Kiado,Budapest,267(1973)。
2.H.Akaike,IEEE iwans.Autom.Contr.,AG19,716(1974)。
3.T.Akutsu,S.Miyano and S.Kuhara,太平洋生物计算研讨会(PacificSymposium on Biocomputing),17,(1999)。
4.T.Akutsu,S.Miyano和S.Kuhara,生物信息学(Bioinformotics),16,727(2000)。
5.T.Akutsu,S.Miyano和S.Kuhara,计算生物学(J.Comp.Biol.),7,331(2000)。
6.G.F.Cooper和E.Herskovits,机器学习(Machine Learning),9,309(1992)。
7.A.C.Davison,Biometrika,73,323(1986)。
8.C.de Boor,.样条函数实用手册(A Practical Guide to Splines).Springer,柏林.(1978)。
9.B.Efron,统计年报(Ann.Stat.),7,1(1979)。
10.P.H.C.Eilers and B.Marx,Statistical Science,11,89(1996)。
11.R.L.Eubank,样条修匀函数和非参数回归,Marcel Dekker,纽约,(1988)。
12.N.Friedman和M.Goldszmidt,in M.I.Jordan ed.,图形模型学习,KluwerAcademic Publisher,421(1998)。
13.N.Friedman,M.Linial,I.Nachman和D.Pe′er,计算生物学杂志,7,601(2000)。
14.P.J.Green和B.W.Silverman,非参数回归和广义线性模型,Chapman &Hall,伦敦.(1994)。
15.T.Hastie和R.Tibshirani,广义加性模型,Chapman & Hall,伦敦.(1990)。
16.D.Heckerman,D.Geiger和D.M.Chickering,机器学习,20,274(1995)。
17.D.Heckerman,in M.I.Jordan ed.,图形模型学习,301,Kluwer AcademicPublisher.(1998)。
18.S.Konishi,(日文),Sugaku,52,128(2000)。
19.G.Schwarz,统计年报(Ann.Stat.),6,461(1978)。
20.P.Spellman,G.Sherlock,M.Zhang,V.Iyer,K.Anders,M.Eisen,P.Brown,D.Botstein和B.Futcher,细胞分子生物学(Mol.BioL Cell),9,3273(1998)。
21.B.W.Silverman,J.R.Stat.Soc.Series B,47,1(1985)。
22.J.S.Simonoff,统计学重的修匀方法,Springer,纽约.(1996)。
23.L.Tinerey和J.B.Kadane,美国统计协会杂志(J.Amer.Statist.Assoc.),81,82(1986)。
24.G.Wahba,J.R.Stat.Soc.Series B,40,364(1978)。
实施例2通过贝叶斯网络和具有异质误差方差及相互作用的非参数回归非线性地建立基因网络模型本发明提出通过基于贝叶斯网络从微阵列基因表达数据中构建基因网络的新的统计学方法。在贝叶斯网络中,网络构建的关键点在于估计各个随机变量的条件分布。本发明考虑用异质误差方差及相互作用拟合非参数回归模型来捕获基因间的非线性结构关系。一旦用贝叶斯网络和非参数回归选定结构图,仍需解决选择最能代表存在于基因之间的系统的最佳结构图的问题。本发明人理论地从贝叶斯方法中得到在普遍情形下选择结构图的新检验数。所提出的方法包括用贝叶斯方法评估基因网络的前述方法。通过分析从破坏100个基因最新得到的酿酒酵母基因表达数据证实了所提出的方法的有效性。
关键词微阵列基因表达数据;贝叶斯网络;非参数回归;异方差性;相互作用;后验概率1.引言随着微阵列技术的发展,现在可以同时观察到大量的基因表达数据。在基因表达数据分析中,构建基因网络在分子生物学和生物信息学领域引起极大关注(见[3],[4],[5],[12],[13],[18]和[22])。但是,由于数据的维度和复杂性的原因,影响了微阵列基因表达数据分析的发展。也就是说,我们所需的信息隐藏在具有噪声的大量数据中。在本申请中,本发明人提出了新的构建基因网络的统计学方法,其甚至能够更清晰地捕获基因间的非线性关系。
贝叶斯网络([19])在通过大量随机变量的联合分布给现象建模中是一种有效方法。近些年来,在通过贝叶斯网络从微阵列基因表达数据中构建基因网络方面已经进行了一些有意义的工作。Friedman和Goldszmidt[12]离散表达值并采用多项式分布(multinominal distribution)为候选的统计模型。Pe`er等[22]研究了用于离散的阈值。另一方面,Friedman等[13]指出离散可能会丢失数据的信息,并且考虑拟合线性回归模型来分析连续的数据。但是,亲本基因线性地依赖于目标基因的假设并不总是能得到保证。Imoto等[18]提出采用非参数附加回归模型来捕获基因间的线性关系以及非线性关系。在本申请中,本发明人提出用贝叶斯网络和非参数异方差回归模型构建基因网络的方法,其更加不受无关值的影响,并且还能捕获亲本基因间的相互作用效果。
一旦确立结构图,必须评价其优点以及与完全未知的实际结构图的接近程度。因此,构建合适的检验数成为统计学基因网络建模的关注焦点。Friedman和Goldszmidt[12]得到根据多项式模型和Dirichlet先验函数选择结构图的检验数BDe。但是,他们保留了Dirichlet先验函数中的未知超级参数并且只能通过试验确定其值。本发明人将结构图选择问题以统计学模型的选择和评价问题来研究,并理论性地从贝叶斯方法([8])中得到用于选择结构图的新检验数。该检验数自动地优化模型中的所有参数并产生最佳结构图。此外,本发明的方法包括通过贝叶斯网络构建基因网络的前述方法。本发明通过破坏100个基因,分析酿酒酵母基因表达数据来说明这种新方法的有效性。
2.贝叶斯网络和具有相互作用的非参数异方差回归模型2.1贝叶斯网络中的非线性模型假定有p维随机可变向量X=(X1...,Xp)T的n组数据{x1,....,xn},其中xi=(xi1...,xip)T和xT表示x的转置矩阵。在微阵列基因表达数据中,n和p代表阵列和基因的数量。在贝叶斯网络结构中,在节点间假定定向非循环结构图G和马尔可夫假想。联合密度函数分解成各个变量的条件密度函数(见[10]),即,f(xi1,...,xip)=f1(xi1|pi1)×...×fp(xip|pip) (1)其中,pij=(pi1(j),...,piqj(j))T是结构图G中的xij的qj-维亲本观察向量。当P1=(X2,X3)T是X1的亲本变量向量,则pi1=(xi2,xi3)T,(i=1,...,n)。通过公式(1),贝叶斯网络统计建模的着重点在于如何构建条件密度函数fj。我们采用条件密度函数fj被参数向量θj参数化,从这些概率模型中提取有用的信息。
Imoto等[18]提出采用非参数回归策略来捕获xij和Pij之间的非线性关系。多数情况下,这种方法可以很好地捕获目标关系。但是,当数据中含有特别接近域{pij}的边界的无关值时,非参数回归模型有时会导致不恰当的光滑估计,即估计的曲线显示一些由无关值作用产生的伪波形。由于所估计的是一个生命系统,太多复杂关系是不合适的。事实上,这种不恰当的情况有时也存在于真实数据分析中。为了避免这个问题,考虑用异质误差方差拟合非参数回归模型。
xij=mj(pij)+ϵij,ϵij~N(0,σij2),---(2)]]>其中mj(·)是一个从Rqj到R的光滑函数,而N(μ,σ2)表示具有平均数μ和方差σ2的高斯分布。此处R代表一组真实数据。该模型包括了Imoto等[18]的模型,而且显然线性回归模型是其中的特例。一种构建系统成分mj(pij)的可能方法是非参数附加回归[18]mj(pij)=mj1(pi1(j))+…+mjqj(piqj(j)),---(3)]]>其是多重线性回归的直接(straightforward)扩展。一般每个光滑函数mjk(·)的特征为n个值mjk(p1k(j)),...,mjk(pnk(j)),系统(3)含有n×qj个参数。于是该模型中的参数个数大大多于观察值的个数,倾向于不稳定的参数估计。在本申请中,通过基础函数方法(basis function approach)mjk(pik(j))=Σm=1Mjkγmk(j)bmk(j)(pik(j)),k=1,…,qj,]]>构建光滑函数mjk(·),其中γ1k(j),...,γMjk(j)k是未知的系数参数,而b1k(j)(·),...,bMjk(j)k(·)为基本函数。从这些表示中,n个参数mjk(p1k(j)),...,mjk(pnk(j))被Mjk的系数参数γ1k(j),...,γMjk(j)k参数化。另一方面,对于相互作用的影响,可以构建相互作用模型mj(pij)=Σl=1sjψjl(pil(j)),---(4)]]>其中pi1(j)是pij的子向量(subvector),也就是说,如果考虑pi1(j)和pi2(j)之间的相互作用,得到pi1(j)=(pi1(j),pi2(j))T。还可以用
ψjl(pil(j))=Σk=1Lj1ξlk(j)c1k(j)(pil(j));l=1,…sj,]]>构建函数ψj1(·)。其中c1k(j)(pi1(j))为基础函数,ξ1k(j)为参数。结合(3)和(4),具有相互作用的非参数回归通常表示为mj(pij)=Σk=1qjΣm=1Mjkγmk(j)bmk(j)(pik(j))]]>+Σl=1sjΣt=1Ljiξil(j)cil(j)(pil(j)).---(5)]]>对于误差方差σij2,采用如下结构σij2=ωij-1σj2,i=1,…,n;j=1,…,p,---(6)]]>其中ω1j,...,ωnj为恒量,σj2是未知参数。通过确立反映误差方差特征的恒量ω1j,...,ωnj,可以表示数据的异方差性。结合(2)、(5)和(6),得到具有异质误差方差和相互作用的非参数回归模型。
fj(xij|pij;γj,ξj,σj2)]]>=(ωij2πσj2)1/2exp[-ωij2σj2{xij-Σk=1υγjkTbjk(pik(j))]]>-Σl=1sjξjlTcjl(pil(j))}2],---(7)]]>其中γjk和bjk(pik(j))是Mjk-维向量,分别表示为γjk=(γ1k(j),…,γMjkk(j))T,]]>bjk(pik(j))=(b1k(j)(pik(j)),…,bMjkk(j)(pik(j)))T,,]]>而ξjk和cjk(pik(j))是Ljk-维向量,分别表示为ξjl=(ξ1l(j),…,ξLj1l(j))T,]]>cjl(pil(j))=(c1l(j)(pil(j)),…,cLj1l(j)(pil(j)))T.]]>因此,模型(7)中的参数为γj=(γj1T,…,γjqjT)T,]]>ξj=(ξj1T,…,ξjsjT)T]]>和σj2。如果变量Xj在结构图中没有亲本变量,指定模型为基于平均数μj和方差σj2的正态分布。因此,贝叶斯网络模型表示为f(xi;θG)=Πj=1pfj(xij|pij;θj),---(8)]]>其中θG=(θ1T,...,θpT)T是包含在结构图G中的参数向量,而θj是条件密度函数fj的参数向量,也就是说,θj=(γjT,ξjT,σj2)T或θj=(μj,σj2)T。
2.2选择结构图的检验数一旦确定结构图,可以构建基于贝叶斯网络和非参数回归的统计模型(8),并可以通过合适方法估计。但是,仍然有待于解决如何选择最近似于隐藏在数据中的系统的最佳结构图的问题。由于似然值在较复杂模型中变大,不能采用似然函数作为模型选择检验数。因此,有必要考虑基于概括的或预期的误差的统计方法、Kullback-Leibler信息、贝叶斯方法等([20])。在本章节中,本发明人从贝叶斯方法中构建用于评估基于本发明模型(8)的结构图的检验数。
当设定了一个结构图,可如下从贝叶斯理论方法构建用于评价该结构图的优点的检验数通过结构图πG的先验概率和数据的边缘概率的乘积得到结构图的后验概率。除去标准化恒量,结构图的后验概率与π(G|Xn)=πG∫Πi=1nf(xi;θG)π(θG|λ)dθG,---(9)]]>相称,其中Xn=(x1,...,xn)T是n×p个基因表达图谱矩阵。π(θG|λ)是满足logπ(θG|λ)=0(n)的参数θG的先验分布,λ为超级参数向量。根据贝叶斯方法,可以选择最佳结构图,例如π(G|Xn)最大化。构建基于结构图的后验概率的检验数的关键问题是高阶积分的计算(9)。Friedman和Goldszmidt[12]采用联合先验分布(conjugate priors)来积分并得到闭式解。为了计算该高阶积分,本发明采用拉普拉斯逼近([11]、[17]和[26])进行积分。
∫Πi=1nf(xi;θG)π(θG|λ)dθG]]>=(2π/n)r/2|Jλ(θ^G)|1/2exp{nlλ(θ^G|Xn)}{1+Op(n-1)}]]>其中r是θG的维数,
lλ(θG|Xn)=1nΣi=1nlogf(xi;θG)+1nlogπ(θG|λ),]]>Jλ(θG)=-∂2{lλ(θG|Xn)}/∂θG∂θGT]]> 是1λ(θG|πn)的众数。从而得到检验数BNRC用于选择结构图BNRC(G)]]>=-2log{πGΠi=1nf(xi;θG)π(θG|λ)dθG}]]>≈-2logπG·rlog(2π/n)+log|Jλ(θ^G)|]]>-2nlλ(θ^G|Xn).---(10).]]>选择最佳结构图使得检验数BNRC(10)最小。采用拉普拉斯方法的优点在于不必考虑使用共轭先验分布。因此,实现了在更多种类模型分布和先验函数中建模。
假定参数向量θj相互独立,那么先验分布可以被分解成π(θG|λ)=Πj=1pπj(θj|λj)]]>因此,(10)中的log|Jλ(θG|Xn)|和nlλ(θG|Xn)分别为log|Jλ(θG|Xn)|=Σj=1plog|-∂2lλj(θj|Xn)∂θJ∂θjT|,,]]>nlλ(θG|Xn)=Σj=1pnlλj(θj|Xn),]]>lλj(θj|Xn)=1nlogfj(xij|pij;θj)+1nlogπj(θj|λj)]]>其中该λj是超级参数向量。因此定义BNRCj=-2log{∫πLjΠi=1nfj(xij|pij;θj)πj(θj|λj)dθj},]]>其中πLj是满足Σj=1plogπLj=logπG]]>的先验概率,BNRC分值从各个局部BNRC分值的总和中得到BNRC=Σj=1pBNRCj.---(11)]]>
通过分别用 和 替代γj和ξj得到基于非参数回归的光滑估计。假定logπ(θG|λ)=0(n),得到检验数BNRC。如果采用满足logπ(θG|λ)=0(1)的先验密度函数,BNRC检验数结果为称作为BIC或SIC的Schwarz[25]的检验数。此时,众数 等价于最大似然估计。
3.评估基因网络3.1非参数回归本章节根据上面第2节中描述方法说明构建基因网络的方法。首先提及非参数回归模型。非参数回归(5)含有两个成分用每个亲本基因的附加模型表示的主效应成分和相互作用成分。在附加模型中,通过B-样条函数([9]和[18])构建各个光滑函数mjk(·)。
在相互作用成分中,采用高斯径向基本函数。cil(j)(pil(j))=exp(-||pil(j)-zjl||22ξjlsjl2),]]>其中zjl是中心向量,sjl2是宽度参数以及ζjl是超级参数。在基于径向基本函数的回归建模中,有两种方法来估计中心zjl和宽度sjl2。首先,通过最小化或最大化适当的目标函数如平方丢失(squared loss)和似然来估计zjl和sjl2。
这种方法称为完全监督学习(fully supervised learning)。另一方面,一种可替代方法为预先仅采用亲本的观察数据pi1(j)来测定zjl和sjl2。在本发明中采用后一种方法,并且采用k-法聚类算法来构建基本函数。径向基本函数的细节在[7]、[21]和[23]中进一步说明。超级参数ζjl控制基本函数中的重叠的数量。
在误差方差方面,考虑异方差回归模型和结构(6)。恒量ω1j,...,ωnj的设计是捕获数据异方差的重要问题。设定权重为ωij=g(pij;ρ)=exp{-ρ||pij-p-j||2/2sj2},---(12)]]>其中ρ是超级参数,p-j=Σi=1npij/n]]>和sj2=Σi=1n||pij-p-j||2/nqj.]]>如果超级参数ρ=0,那么权重ω1j=...=ωnj=1,并且模型具有异质误差方差。如果采用大的ρ值,那么接近于亲本变量区域边界的数据的误差方差也大。因此,如果存在接近边界的无关值,可以通过使用恰当的ρ值来降低其影响并得到合适光滑估计。
3.2先验分布假设先验分布πj(θj|λj)被分解成因数(factorize) 其中λjk和vjl是超级参数。用奇异变量Mjk正态分布作为γjk的先验分布,πjk(γjk|λjk)=(2πnλjk)-(Mjk-2)/2|Kjk|+1/2]]>×exp(-nλjk2γjkTKjkγjk),---(13)]]>其中Kjk是一个满足γjkTKjkγjk=Σα=3Mjh(γαk(j)-2γα-1,k(j)+γα-2,k(j))2]]>的Mjk×Mjk的对称半正定矩阵(symmetric positive semidefinit matrix)。ξjl的先验分布为πjl(ξjl|vjl)=(2πnvjl)-vjl/2exp(-nvjl2||ξjl||2),---(14)]]>其中vjl是ξjl的维数。
接着假设结构图πG的先验概率。Friedman和Goldszmit[12]采用基于结构图编译的MDL的先验分布。在本发明中,数据的边缘概率与由超级参数调整的II型似然估计相似。因此,设结构图πG的先验概率, =Πj-1pexp{-(qj+2sj+1)}=Πj=1pπLj]]>该先验分布的论证(justification)基于ABIC的Akaike[2]的贝叶斯信息检验数和Akaike[1]的信息检验数AIC。
3.3检验数在章节2.2中,得到在普遍结构中选择结构图的检验数BNRC。采用方程11,通过累加局部分值BNRCj可以得到结构图的BNRC分值。结果总结成如下定理。
定理1.令f(xi;θG)=Πj=1pfj(xij|pij;θj)]]>为(8)给出的贝叶斯网络和非参数回归模型,令π(γjk|λjk)和π(ξjl|υjl)分别为(13)和(14)定义的参数γjk和ξjk的先验密度。那么,用于评价结构图的检验数表示为BNRC=Σj=1pBNRCj,]]>其中
BNRCj=2(qj+2qj+1)]]>-(Σk=1qjMjk+Σl=1sjLj1+1)log(2π/n)]]>-Σi=1nlogωij+nlog(2πσ^j2)+n]]>+Σk=1qj{log|Λjk|-Mjklog(nσ^j2)}-log(2σ^j2)]]>+Σl=1nj{log|Pj1|-Lj1log(nσ^j2)}]]>+Σk=1qj{(Mjk-2)log(2πσ^j2nβjk)-log|Kjk|+]]>+nβjkσ^j2γ^jkTKjkγ^jk}]]>+Σl=1sj{vjllog(2πσ^j2nTjl)+nTjlσ^j2||ξ^jl||2},]]>以及Λjk=BjkTWjBjk+nβjkKjk;(Mjk×Mjk),]]>Γjl=CjlTWjCjl+nTjlIjl;(Ljl×Ljl),]]>Iij=diag[1,...,1};(Ljl×Ljl),Bjk=(b1k(j)(p1k(j)),…,bMjkk(j)(pnk(j)))T;(n×Mjk),]]>Cjl=(c1l(j)(p1l(j)),…,cLjll(j)(pnl(j)))T;(n×Ljl),]]>Wj=diag(ω1j,...,ωnj);(n×n)σ^j2=Σi=1nωij{xij-Σk=1qjγjkTbjk(pik(j))]]>-Σl=1sjξjlTcjl(pil(j))}2/n.]]>通过如下逼近黑塞矩阵(Hessian matrix),log|-∂2lλj(θj|Xn)∂θj∂θjT|≈Σk=1qjlog|-∂2lλj(θj|Xn)∂γjk∂γjkT|]]>+Σl=1sjlog|-∂2lλj(θj|Xn)∂ξj1∂ξj1T|+log|-∂2lλj(θj|Xn)∂(σj2)2|.]]>3.4所罚似然估计假设(7)中定义非参数回归模型。当估计方法为最大化log似然函数的最大似然方法时,参数估计结果不稳定,并且由于模型的适应性会导致过拟合(overfitting)。在本发明方法中,估计值θ^j=(γ^jT,ξ^jT,σ^j2)T]]>是作为lλj(θj|Xn)]]>的众数得到,并且该估计方法与最大所罚似然方法([14]和[15])等价,其包括了最大似然方法作为特例。
众数 依赖于超级参数。事实上,超级参数在估计光滑曲线中起重要作用。本发明分别采用分布(13)和(14)作为γjk和ξjl的先验分布。那么lλj(θj|Xn)的结果为lλj(θj|Xn)∝Σi=1nlogfj(xij|pij;θj)]]>-Σk=1qjnλjk2γjkTKjkγjk-Σl=1sjnvjl2||ξjl||2]]>右边的第一张图是log-似然函数,第二和第三张图称为粗罚(roughnesspenalties)。超级参数λjk和υjl称为控制拟合曲线光滑度的修匀参数。也就是说,如果修匀参数小,估计曲线贴近数据,称为过拟合。另一方面,采用大的修匀参数值,估计接近线性拟合。因此,超级参数的选择对于捕获基因间的关系起重要作用。超级参数的功能将在下一章节中说明。
图1.超级参数不同值的光滑估计。
4.真实数据分析在本章节中,将通过分析酵母基因表达数据来说明所提出的方法的有效性。本研究小组已经采取了系统的实验方法,观察由于基因破坏引起的微阵列上基因表达水平的变化。通过使用这种方法,本发明人已进行一项旨在揭示酿酒酵母的5871个基因之间的基因调控网络的项目,而一些实验室也已经报到了相似的研究。本发明人已经从基因破坏试验中收集大量表达图谱来评价遗传调控网络。备置400个以上的突变体,积累基因表达图谱。
用扫描仪监控打点在微阵列上的5871个基因的转录水平。400个以上的破坏体的表达图谱被收集在本发明数据库中。评定了微阵列上所有基因的标准偏差(SD),该SD值粗略地表示试验的误差。在本发明的数据中,以0.5作为试验精确性的临界点。本发明人已经根据所有基因的表达比的标准偏差,评价了那些表达图谱的精确性。从400个图谱中选择107个破坏体,其中包括68个转录因子被破坏的突变体。
使用100个微阵列,从上述数据中构建521个基因的基因网络。发现了94个清楚地确定所调控的基因的转录因子,从5871种表达图谱中选择了受该94个因子控制的521个基因的表达图谱。在本发明模型中,用20个B-样条函数和20个径向基本函数构建非参数回归模型。确定由于无论使用多少基本函数,超级参数都控制着拟合曲线的光滑度,所以不能从视觉上发现不同数量基本函数的光滑估计间的差别。对相互作用的成分应用双基因作用。因此,相互作用效应以拟合面得到,可以从视觉上看出。
图2.权重恒量的作用。
图3.具有相互作用的非参数回归的估计表面
图4.分析521个酿酒酵母基因得到的部分网络在论述分析结果之前,先说明了超级参数在先验分布和权重恒量中的作用。图1(a)表示用超级参数的3种不同值光滑估计的YGL237C和YEL071W的扫描图。显然,光滑估计强烈地依赖于超级参数值。图1(b)描述图1(a)中的两个基因的BNRC检验数的变化。可以选择超级参数的最佳值作为BNRC最小值,而且最佳光滑估计(实心曲线)能够很好地捕获这些基因间的结构关系。这些虚的和带点的曲线分别与最大似然估计和参数线性拟合接近。图2(a)表明权重恒量ω1j...,ωnj的作用。如果采用同方差回归模型,得到虚曲线,其显示一些由左上角的数据作用产生的伪波形。通过调节(12)中的超级参数ρ,使估计曲线变为实心曲线。还可以通过最小化BNRC检验数来选择最佳ρ值(见图2(b))。当然,当恰当地得到光滑估计时,ρ的最佳值趋于0。
采用两步策略来拟合具有相互作用的非参数回归模型。首先,估计由附加B-样条函数回归表示的主效应。接着,用相互作用成分拟合残差。图3描述YIL094与其亲本基因CYKL152和CYER055C之间关系的拟合面的例子。清楚的表明当增加亲本基因时,这两个亲本基因的相互作用导致过度表达。
分析及其评价的结果如下在酿酒酵母中,GCN4基因编码氨基酸生物合成的“总控制”(general control)系统的转录激活子,该系统是一个有至少12条不同生物合成路径的网络[6]。试验说明了对由组氨酸类似物3-氨基三唑诱导的“氨基酸饥饿”信号应答的涉及GCN4p水平的总控制结果。GCN4激活30个以上基因的转录,这些基因参与11种对氨基酸饥饿或者tRNA合成酶活性受损响应的氨基酸的生物合成[24]。嘌呤的生物合成基因ADE1、ADE4、ADE5,7和ADE8响应氨基酸饥饿时,表现为GCN4-依赖性表达[24]。GCN4对氨基酸或者嘌呤饥饿响应时,GCN4激活氨基酸和嘌呤的生物合成基因的转录[24]。这些试验结果表明嘌呤代谢和氨基酸代谢之间通过GCN4密切联系。本发明的关系图谱与嘌呤和氨基酸代谢之间的已知关系十分吻合。
5.结论在本申请中,本发明人提出了用贝叶斯网络和非参数回归从微阵列基因表达数据中构建基因网络的新的统计学方法。本发明方法的关键点在于用非线性异方差回归模型来捕获基因间的非线性关系和表达数据的异方差性。网络构建的关键问题在于结构图的评价。本发明将该问题作为统计模型选择或估计的问题来研究,并从贝叶斯方法中得到了选择结构图的新的检验数。本发明方法包括采用叶斯网络构建基因网络的前述方法,并从理论上和方法学上改进了这些方法。本发明人已经通过分析酿酒酵母基因表达数据来说明本发明方法的有效性,并通过与生物学知识比较来评价得到的网络。本发明人不用生物信息就构建了基因网络,但得到的网络却包括了许多重要关系,这些关系与生物学知识吻合。因此,本发明人认为本发明方法可以用于分析完全未知的系统如人类基因组。
参考文献[1]H.Akaike,信息理论和最大似然法则的扩展,B.N.Petrov,&F.Csald,eds.,Akademiai Kiado,Budapest,267(1973)。
H.Akaike,似然与贝叶斯方法,J.M.Bernard,M.H.DeGroot,D.V.Lindley& A.F.M.Smith,cds.,Univ.Press,Valencia(1980)。
T.Akutsu,S.Miyano和S.Kuhara,根据布尔网络模型Proe从少量基因表达图谱确定基因网络,1999年太平洋生物计算研讨会(Pacific Symposiumon Biocomputing 1999),17(1999)。
T.Akutsu,S.Miyano和S.Kuhara,推导基因网络和代谢途径中的定性关系,生物信息学(Bioinformatics),16,727(2000)。
T.Akuten,S.Miyano和S.Kuhara,基于矩阵乘法和指纹功能的用于确定布尔网络和相关生物学网络的算法,计算生物学杂志(J.ComputBiol).7,331(2000)。
G.Albcrlt,H.-U.Mosch,B.Hoffmann,U.Re和G.H.Braus,追踪酿酒酵母中的Gcn4蛋白介导反应.生物化学杂志(J.Biol.Chem.),273,12696(1998)。
T.Andou,S.Imoto和S.Konishi,用径向基础函数建立非线性回归模型.省略。
J.O.Berger,统计决定理论与贝叶斯分析,Springer-Verlag纽约(1985)。
C.deBoor,样条函数使用手册,Springer-Verlag.柏林(1978)。
R.Cowell,A.Dawid,S.Lauritzen和D.Spiegelhalter,或然网络与专家系统,SpringerVerlag纽约(1999)。
A.C.Davison,近似预言似然,Biometrika,73,323(1986)。
N.Friedman和M.Goldszmidt,具有局部结构关系的学习贝叶斯网络,M.I.Jordan ed.,Kluwer Academic Publisher.421(1998)。
N.Friedman,M.Linial,L.Nachmank和D.Pe′er,用贝叶斯网络分析表达数据,计算生物学杂志,7,601(2000)。
P.J.Green和B.W.Silverman,非参数回归和归纳线性模型,Chapman &Hall(1994)。
I.J.Good和R.A.Gaskins,概率密度的非参数粗罚,Biometrika,58,255(1971)。
T.Hastie和R.Tibshirani,归纳附加模型,Chapman & Hall(1990)。
D.Heckerman,贝叶斯网络学习指南,M.I.Jordan cd.,KtnwcrAcademic Publisher.(1998)。
S.Imoto,T.Goto和S.Miyano,用贝叶斯网络和非参数回归估计基因间的基因网络和功能结构关系,2002年太平洋生物计算研讨会记录(ProcPacific Symposium on Biocomputing 2002),(2002)。
F.V.Jensen,贝叶斯网络导论.伦敦大学出版社(1996)。
S.Konishi,统计模型评价和信息检验数,S.Ghosh ed MarcelDamer(1999)。
J.Moody和C.J.Darken,局部协调的处理单元素的快速网络学习,Neur.Comp.1,281(1989)。
D.Pe′er,A.Regev,G.Elidan和N.Friedman,从被干扰的表达图谱中推导子网络,生物性信息学,17,副刊,1,215(ISMB 2001)。
B.D.Riplcy,图形识别和神经网络,剑桥大学出版社(1996)。
R.J.Rolfes和A.G.Hinnebusch,嘌呤清除刺激酵母转录激活子GCN4的翻译隐含蛋白激酶GCN2的激活,细胞分子生物学,13,5099(1993)。
Schwarz,估计模型的维数,Ann.Sts.,6,461(1978)。
L.Tinerey和J.B.Kadane,后验和边缘密度的精确逼近,美国统计学协会杂志(J.Amer.Statist.Assoc.),81,82(1986)。
实施例3用于非线性建立基因网络的贝叶斯网络和非参数异方差回归摘要本申请提出用贝叶斯网络来从微阵列基因表达数据中构建基因网络的新统计学方法。贝叶斯网络构建的关键点在于估计各个随机变量的条件分布。本发明考虑用具有异质误差方差的非参数回归模型拟合微阵列基因表达数据来捕获基因间的非线性结构关系。仍需解决选择最能代表存在于基因之间的系统的最佳结构图的问题。本发明人理论地得到在普遍情形下从贝叶斯方法中选择结构图的新检验数。所提出的方法包括基于贝叶斯网络的前述方法。通过分析从破坏100个基因最新得到的酿酒酵母基因表达数据证实了所提出的方法的有效性。
1.引言由于微阵列技术的发展,构建基因网络在分子生物学和生物信息学领域引起极大关注[3,4,5,14,15,17,22和28]。但是,数据的维度和复杂性影响了微阵列基因表达数据分析的发展。也就是说,我们所需的信息隐藏在大量具有噪声的数据中。在本申请中,本发明人提出了构建基因网络的新的统计学方法,其甚至能够更清晰地捕获基因间的非线性关系。
在通过大量随机变量的联合分布而对现象建模中,贝叶斯网络[7,23]是一种有效方法。近些年来,在通过贝叶斯网络从微阵列基因表达数据中构建基因网络方面已经进行了一些有意义的工作。Friedman和Goldszmidt[12,13,14]离散表达值并假定多项式分布为候选的统计模型。Pe`er等[28]研究了用于离散的阈值。另一方面,Friedman等[15]指出离散可能会丢失数据的信息。实际上,离散值和阈值的数量都是未知参数,必须从数据中估计。得到的网络强烈地依赖于它们的值。Friedman等[15]考虑拟合线性回归模型,以连续的方式分析数据(又见[20])。但是,亲本基因线性地依赖于目标基因的推定并不总是能得到保证。Imoto等[22]提出采用非参数附加回归模型(见[16,18])来捕获基因间的线性关系以及非线性关系。在本申请中,本发明人提出用贝叶斯网络和非参数异方差回归构建基因网络的方法,其更加不受无关值的影响。
一旦确立结构图,必须评价其优点以及与完全未知的实际结构图的接近程度。因此,建立合适的检验数成为统计学基因网络建模的关注焦点。Friedman和Goldszmidt[14]使用检验数BDe,其最初从[21]得出,用于选择结构图的。检验数BDe仅评价基于多项式分布和Dirichlet先验函数的贝叶斯网络模型。然而,Friedman和Goldszmidt[14]在Dirichlet先验函数中保留了未知超级参数并且只能通过试验确定其值。本发明人将结构图选择问题作为统计学模型的选择和评价问题来研究,并理论性地从贝叶斯方法(见[6])得到用于选择结构图的新检验数。该检验数自动地优化模型中的所有参数并产生最佳结构图。此外,本发明的方法包括前述的根据贝叶斯网络构建基因网络的方法。本发明人采用蒙特卡罗模拟方法来说明这种新方法的有效性。本发明还分析从破坏100个基因最新得到的酿酒酵母基因表达数据。
2.贝叶斯网络和非参数异方差回归模型2.1非线性贝叶斯网络模型假定有p个基因的n组数据{x1,....,xn},其中xi=(xi1...,xip)T和xT表示x的转置矩阵。在贝叶斯网络结构中,在节点间假定定向非循环结构图G和马尔可夫假设。联合密度函数分解成各个变量的条件密度函数,即,f(xi1,…,xip)=Πj=1pfj(x1j|pij),---(1)]]>其中,pij=(pi1(j),...,piqj(j))T是结构图G中的xij的qj-维亲本观察向量。当基因2和基因3是基因1的亲本基因时,则pi1=(xi2,xi3)T,(i=1,...,n)。通过公式(1),贝叶斯网络建模的着重点在于如何构建条件密度函数fi。假定条件密度函数fj被参数向量θj参数化,从这些概率模型中提取有用的信息。
Imoto等[22]提出采用非参数回归策略来捕获xij和pij之间的非线性关系,并认为基因间存在许多非线性关系而线性模型对此难以获得充足的结果。多数情况下,这种方法可以很好地捕获目标关系。但是,当数据中含有特别接近域{pij}的边界的无关值时,非参数回归模型有时会导致不恰当的光滑估计,即估计的曲线显示一些由无关值作用产生的伪波形。由于所估计是一个生命系统,太多复杂关系是不合适的。事实上,这种不恰当的情况有时也存在于真实数据分析中。为了避免这个问题,本发明人考虑用异质误差方差拟合非参数回归模型。
xij=mj1(pi1(j))+…+mjqj(piqj(j))+ϵij,---(2)]]>其中ξij单独和一般地依赖于平均数0和方差σij2,而mjK(·)是一个从R到R的光滑函数。此处R代表一组真实数据。该模型包括了Imoto等[22]的模型,而且线性回归模型显然是其中的特例。一般每个光滑函数mjk(·)的特征为n个值mjk(p1k(j)),...,mik(pnk(j)),系统(2)含有(n×qj+n)个参数。该模型中的参数个数大大多于观察值的个数,倾向于不稳定的参数估计。在本申请中,通过基础函数方法来构建光滑函数mjk(·)mjk(pik(j))=Σm=1Mjkγmk(j)bmk(j)(pik(j)),k=1,…,qj,]]>其中γ1k(j),...,γMjk(j)k是未知的系数参数,而b1k(j)(·),...,bMjk(j)k(·)为基本函数。从这些表示中,n个参数mjk(p1k(j)),...,mjk(pnk(j))被Mjk的系数参数γ1k(j),...,γMjk(j)k参数化。
由于线性回归不能确定贝叶斯因果关系的方向而且在多数情况下产生错误方向,因此,强烈地推荐用非参数回归替代现行回归。本发明通过一个简单例子说明所提出的模型与线性回归相比的优点。假设有图1(a)中基因1和基因2的数据。设定两种模型基因1→基因2和基因2→基因1,得到分别如图1(b)和图1(c)所示光滑估计。然后根据所提出的检验数可以确定模型(b基因1→基因2)好于模型(c基因2→基因1),该检验数在后面章节中得到(模型的分值为(b)120.6和(c)134.8)。由于这些数据是从真实结构图基因1→基因2中得到,所以本发明方法得到了正确结果。但是,如果用线性回归模型拟合这些数据,会选择模型(c)(线性模型的分值为(b)156.0和(c)135.8)。如此,基于线性回归的方法得出错误结果。
假设基因之间的关系几乎为线性的情况,本发明的方法和线性回归都可以恰当地与数据拟合。但是,采用线性模型显然难以确定贝叶斯因果关系的方向。此时,方向是不精确的。
图1.模拟数据真实的因果关系是基因1→基因2。(a)模拟数据的分散点。(b)结构图基因1→基因2的光滑曲线。(c)结构图基因2→基因1的光滑曲线。这些曲线通过本发明的方法得到。
在误差方差方面,我们采用结构,σij2=ωij-1σj2,i=1,…,n;j=1,…,p,---(3)]]>其中ω1j,...,ωnj为恒量,σj2是未知参数。通过确立反映误差方差特征的恒量ω1j,...,ωnj,可以表示数据的异方差性。结合(2)和(3),得到具有异质误差方差的非参数回归模型,
fj(xij|pij;γj,σj2)=(ωij2πσj2)1/2exp[-ωij2σj2{xij]]>-Σk=1qjγjkTbjk(pjk(j))}2],---(4)]]>其中γjk和bjk(pik(j))是Mjk-维向量,分别表示为γjk=(γ1k(j),…,γMjkk(j))T]]>和bjk(pik(j))=(b1k(j)(pik(j)),…,bMjkk(j)(pik(j)))T.]]>如果第j个基因在结构图中没有亲本基因,指定模型为基于平均数μj和方差σj2的正态分布。因此,定义贝叶斯网络模型为f(xi;θG)=Πj=1pfj(xij|pij;θj),---(5)]]>其中θG=(θ1T,...,θpT)T是包含在结构图G中的参数向量,而θj是条件密度函数fj的参数向量,也就是说,θj=(γjT,σj2)T或θj=(μj,σj2)T。
2.2选择结构图的检验数一旦确定结构图,可以构建基于贝叶斯网络和非参数回归的统计模型(5),并可以通过合适方法估计。但是,仍然有待于解决如何选择最近似隐藏在数据中的系统的最佳结构图的问题。由于似然值在较复杂模型中变大,不适宜采用似然函数作为模型选择检验数。因此,有必要考虑基于概括或预期的误差的统计方法、Kullback-Leibler信息、贝叶斯方法等(见[1,24,25])。在本章节中,从贝叶斯方法中构建了用于估计基于本发明模型(5)的结构图的检验数。
通过结构图πG的先验概率和数据的边缘概率的乘积得到结构图的后验概率。除去标准化恒量,结构图的后验概率与π(G|Xn)=πG∫Πi=1nf(xi;θG)π(θG|λ)dθG,---(6)]]>相称,其中Xn=(x1,...,xn)T是n×p个基因表达图谱的矩阵。π(θG|λ)是满足logπ(θG|λ)=0(n)的参数θG的先验分布,λ为超级参数。根据贝叶斯方法,可以选择最佳结构图使得π(G|Xn)最大化。构建基于结构图的后验概率的检验数的关键问题是高阶积分的计算(6)。Hekerman和Geiger[20]采用联合先验分布来求积分并得到闭式解。为了计算该高阶积分,本发明采用拉普拉斯逼近([9,19,316])进行积分,
∫Πi=1nf(xi;θG)π(θG|λ)dθG]]>=(2π/n)r/2|Jλ(θ^G)|1/2exp{nlλ(θ^G|Xn)}{1+Op(n-1)}]]>其中r是θG的维数,lλ(θG|Xn)=Σi=1nlogf(xi;θG)/n+logπ(θG|λ)/n,Jλ(θG)=-∂2{lλ(θG|Xn)}/∂θG∂θGT]]>以及 是lλ(θG|Xn)的众数。从而定义用于选择结构图的贝叶斯网络和非参数异方差回归检验数BNRChetero,BNRChetero(G)]]>=-2log{πG∫Πi=1nf(xi;θG)π(θG|λ)dθG}]]>≈-2logπG-rlog(2π/n)+log|Jλ(θ^G)|]]>-2nlλ(θ^G|Xn).---(7)]]>选择最佳结构图使得检验数BNRChetero(7)最小。采用拉普拉斯方法的优点在于不必考虑使用共轭先验分布。因此,实现了在更多种类模型分布和先验函数中建模。
假定参数向量θj相互独立,那么先验分布可以被分解成π(θG|λ)=Πj=1pπj(θj|λj).]]>因此,(7)中的log|Jλ(θG|Xn)|和nlλ(θG|Xn)分别为log|Jλ(θG|Xn)|=Σj=1plog|-∂2lλj(θj|Xn)∂θj∂θjT|,]]>lλ(θG|Xn)=Σj=1plλj(θj|Xn),]]>其中lλj(θj|Xn)=logfj(xij|pij;θj)/n+logπj(θj|λj)/n.]]>该λj是超级参数向量。因此定义BNRChetero(j)=-2log{∫πLjΠi=1nfj(xij|pij;θj)πj(θj|λj)dθj},]]>其中πLj是满足Σj=1plogπLj=logπG,]]>的先验概率,BNRChetero分值从各个局部BNRC分值的总和中得到BNRChetero=Σj=1pBNRChetero(j).---(8)]]>
图2.模拟数据的拟合曲线。细曲线使用系数加权的B-方样函数,粗曲线是通过加权B-方样函数线性组合得到的光滑估计。
通过用 替代γj得到基于非参数异方差回归的光滑估计。假定logπ(θG|λ)=0(n),得到检验数BNRChetero。如果采用满足logπ(θG|λ)=0(1)的先验密度函数,BNRChetero检验数得到称作为BIC或SIC的Schwarz[30]的检验数。此时,众数 等价于最大似然估计。
3.评估基因网络3.1非参数回归在本章节中,给出根据上面提出的方法构建基因网络的方法。首先提及非参数回归模型。非参数回归(5)含有两个成分用每个亲本基因的附加模型表示的主效应成分和相互作用成分。在附加模型中,通过B-样条函数([10,22])构建各个光滑函数mjk(-)。图2是B-样条函数修匀的曲线的例子。细曲线为系数加权的B-样条函数,粗线为通过加权的B-样条函数线性组合得到的光滑曲线。
在误差方差方面,考虑异方差回归模型和结构(3)。选择恒量ω1j,...,ωnj对于捕获数据异方差很重要。设定权重为ωij=g(pij;pj)=exp{-pj||pij-p-j||2/2sj2},---(9)]]>其中ρj是一个超级参数,p-j=Σi=1npij/n]]>和sj2=Σi=1n||pij-p-j||2/nqj.]]>如果超级参数ρj=0,那么权重ω1j=...=ωnj=1,并且模型具有异质误差方差。如果采用大的ρj值,那么接近于亲本变量区域边界的数据的误差方差也大。因此,如果存在接近边界的无关值,可以通过使用恰当的ρj值来降低其影响并得到合适光滑估计。
3.2先验分布假设先验分布πj(θj|λj)被分解为因数rj(θj|λj)=Πk=1qjπjk(γjk|λjk),]]>其中λjk是超级参数。用奇异变量Mjk正态分布作为γjk的先验分布,πjk(γjk|λjk)=(2πnλjk)-(Mjk-2)/2|Kjk|+1/2]]>×exp(-nλjk2γjkTKjkγjk),---(10)]]>其中Kjk是一个满足γjkTKjkγjk=Σα=3Mjk(γαk(j)-2γα-1,k(j)+γα-2,k(j))2]]>的Mjk×Mjk的对称半正定矩阵。
接着考虑结构图πG的先验概率。Friedman和Goldszmit[14]采用基于结构图编译的MDL的先验分布。在本发明中,数据的边缘概率等价于由超级参数调整的II型似然估计。因此,设定结构图πG的先验概率, =Πj=1pexp{-(qj+1)}=Πj=1pπLj]]>该先验分布的正确性得到称为ABIC的Akaike[2]的贝叶斯信息检验数和Akaike[1]的信息检验数AIC的证实。
3.3检验数本发明人得到在普遍结构中选择结构图的检验数BNRChetero。采用方程(11),通过累加局部分值BNRC(j)hetero可以得到结构图的BNRChetero分值。结果总结成如下定理。
定理1.令f(xi;θG)为(5)给出的贝叶斯网络和非参数异方差回归模型,令π(γjk|λjk)为(10)定义的参数γjk的先验密度。那么,用于评价结构图的检验数表示为BNRChetero=Σj=1pBNRC(j)hetero,]]>其中,BNRChetero(j)=2(qj+1)-(Σk=1qjMjk+1)log(2π/n)]]>-Σi=1nlogωij+nlog(2πσ^j2)+n]]>+Σk=1qj{log|Λjk|-Mjklog(nσ^j2)}-log(2σ^j2)]]>+Σk=1qj{(Mjk-2)log(2πσ^j2/nβjk)-log|Kjk|+]]>+nβjkγ^jkTKjkγ^jk/σ^j2},]]>
以及Λjk=BjkTWjBjk+nβjkKjk;(Mjk×Mjk),]]>Bjk=(b1k(j)(p1k(j)),…,bMjkk(j)(pnk(j)))T;(n×Mjk),]]>Wj=diag(ω1j,...,ωnj);(n×n)σ^j2=Σi=1nωij{xij-Σk=1qjγ^jkTbjk(pik(j))}2/n.]]>BNRChetero(j)=2(qj+1)-(Σk=1qjMjk+1)log(2π/n)]]>-Σi=1nlogωij+nlog(2πσ^j2)+n]]>+Σk=1qj{log|Λjk|-Mjklog(nσ^j2)}-log(2σ^j2)]]>+Σk=1qj{(Mjk-2)log(2πσ^j2/nβjk)-log|Kjk|+]]>+nβjkγ^jkTKjkγ^jk/σ^j2},]]>通过如下逼近黑塞矩阵,log|-∂2lλj(θj|Xn)∂θj∂θjT|≈Σk=1qjlog|-∂2lλj(θj|Xn)∂γjk∂γjkT|]]>+log|-∂2lλj(θj|Xn)∂(σj2)2|.]]>3.4网络学习在贝叶斯网络内容中,确定最佳网络是一个NP-难题。在本文中,对网络学习采用如下的“贪心爬山”算法步骤1生成分值矩阵,其第(i,j)个元素是结构图基因i→基因j的BNRC(j)hetero分值;步骤2对每个基因实施3种边缘处理之一“增加”、“删除”和“反转”,得到最小BNRChetero;步骤3重复步骤2直到BNRChetero不再减小。
一般地,贪心爬山算法具有许多局部最小值,而计算结果依赖于变量的计算顺序。为了解决该问题,改变基因的计算顺序,在步骤3中得到多个候选学习顺序。网络学习的另一个问题是,当基因数目多时亲本基因的搜索范围很宽。在这种情形下,可以严格基于步骤1给出的分值矩阵的候选亲本基因的设定。
图3.超级参数的不同值的光滑估计。(a1)超级参数βjk在B-样条函数的系数的先验分布中的作用。该参数控制拟合曲线的光滑度。(b1)和(c1)超级参数ρj在误差方差的参数中的作用。该参数能够捕获数据的异方差性和能够降低无关值的作用。
3.5超级参数考虑(4)所定义的非参数回归模型。估计值 是1λj(θj|Xn)的众数并且依赖于超级参数。事实上,超级参数对于估计光滑曲线起重要作用。
在本发明模型中,用20个样条函数构建非参数回归模型。因为无论使用多少基本函数,超级参数控制拟合曲线的光滑度,所以基于不同多个基本函数的光滑估计间的差别难以从视觉上发现。用超级参数的3个不同值的光滑估计时,图3(a1)显示了YGL237C和YEL071W的离散图。数据的详细内容在后面章节中说明。显然,光滑估计强烈地依赖于超级参数值。图3(a2)描述图3(a1)中的两个基因的BNRChetero检验数的变化。可以选择最佳超级参数值作为最小的BNRChetero,而最佳光滑估计(图3(a1)中的实心曲线)能够很好地捕获这些基因之间的结构关系。虚曲线和点曲线分别与最大似然估计和参数线性拟合接近。
权重恒量ωij,...ωnj的作用在图3(b1)和(c1)中表明。如果采用非参数同方差(homoscedastic)回归模型[22],得到具有由左上角的数据(b1)作用产生的伪波形的虚曲线。通过调整(9)中的超级参数ρj,估计曲线为实心曲线。最小化BNRChetero检验数来获得ρj的最佳值(参图3(b2))。当然,当得到合适的光滑估计时,ρj的最佳值趋向于零。
最后,给出用于估计光滑曲线和优化超级参数的算法。
步骤1确定超级参数ρj,步骤2初始化γjk=0,k-1,...,qj;步骤3重复步骤3-1和步骤3-2寻找最佳βjk;3-1计算γjk=(BjkTWjkBjk+nβjkKjk)-1BjkTWjk]]>×(x(j)-Σk′≠kBjk′γjk′),]]>来确定βjk,步骤3-1计算步骤3-2评价用候选的βjk值重复步骤3-1,选择最佳βjk值,其使BNRC(j)hetero最小化。
步骤4收敛以k=1,...,qj,1,...,qj,1,...重复步骤3直到满足合适的收敛检验数。
步骤5用候选的ρj值重复步骤1到4,选择使BNRC(j)hetero最小化的最佳ρj值。
4.计算试验4.1蒙特卡罗模拟用蒙特卡罗模拟方法来证明本发明方法的有效性。数据来源于图4(a)的虚拟网络(artificial network),各节点之间具有如下功能结构关系X1=X22+2sin(X5)-2X7+ϵ1,ϵ1~N(0,(4s)2),]]>X2={1+exp(-4X3)}-1+ε2,ε2~N(0,s2),X3=ε3,,ε3~N(0,1),X4=X52/3+ϵ4,ϵ4~N(0,(4s)2),]]>X5=X3-X62+ϵ5,ϵ5~N(0,(4s)2),]]>X6=ε6,ε6~N(0,1), X8=exp(-X4-1)/2+ε8,ε8~N(0,(2s)2),X9=ε9,ε9~N(0,1),X10=cos(X9)+ε10,ε10~N(0,(4s)2),其中s为恒量。在将亲本变量的观察值变换成平均数0和方差1后,得到子代变量的观察值。
从该虚拟网络中得到100个观察值,目的在于从这些模拟数据中重建图4(a)的网络。对噪声变量进行两种不同的设定,其一设定s=0.2,另一个为s=0.1。来自噪声设定s=0.2的观察值通过实验发现与真实的微阵列数据相似。重复蒙特卡罗模拟1000次,重点关注的是正确的估计的数量。图4(b)和(c)分别是s=0.2和s=0.1时的蒙特卡罗模拟结果。
s=0.2x1-X415 87%X6- X332581%X673 99%X919053%X972 99%X7- X212 58%X10 10 70%X332 97%x2-X420 65%X424 71%X671 87%X527 85%X865 78%X612 98%X912092%X968 96%X10 27 59%X10 3 67%x3-X131 55%X8 -X11 100%X420 90%X613987%X866 52%X914665%X914970%X10-X371 87%x4-X635 51%X48 75%X913795%X613995%x5-X29 100% X883 94%X828 97%X9108100%X10 6 100%s=0.1X1-X416 100% X6- X330874%X653 83%X44 50%X968 100% X7- X21 100%X10 1 100% X328 100%X2-X44 100% X420 95%
X6 5791% X518 89%X8 3197% X693 99%X9 8996% X943 98%X101080% X8-X11 100%X3-X1 2357% X615481%X4 5 80% X912377%X8 5852% X9-X616951%X9 157 69% X10 -X381 88%X4-X9 172 99% X49 67%X5-X2 8 75% X6156100%X8 27100% X71 100%X9 106 100% X877 96%X106 100%表1.蒙特卡罗模拟中的假阳性。节点名称后面所附的数字是缺少方向信息的估计关系的数量,百分比是方向信息。例如在s=0.2中,在1000个蒙特卡罗模拟中,本发明方法估计关系“x1→x4”或“x4→x1”15次,15次中的87%反映从左到右的方向(“x1→x4”)。
蒙特卡罗模拟的结果总结如下在设定噪声变量s=0.2中,本发明的模型能够很好地重建目标网络。表1说明了蒙特卡罗模拟中的假阳性,可以看出假阳性的百分比几乎小于10%。由于设定s=0.2中的估计数据与真实的微阵列数据相似,可以认为本发明的网络估计方法在真实数据分析中也是有效的。在图4(b)和(c)以及表1中,阴性数量远远少于假阳性数量。这种趋势在研究数据分析中是优选的。与s=0.2的结果相比,在设定s=0.1中,本发明方法可以更精确地重建目标网络并且减少假阳性的数量。
图4.蒙特卡罗模拟的结果。(a)真实网络。(b)s=0.2的结果。(c)s=0.1的结果。边缘旁边的数字表示从1000个蒙特卡罗试验所估计的相关的数量。百分比包含了边缘方向的信息。例如,在s=0.2中,1000个蒙特卡罗试验出现x5和x1之间相关958次,其中97%是正确的方向(从x5到x1)。
4.2真实数据分析在本章节中,通过分析酿酒酵母基因表达数据来说明本发明方法的有效性,这些数据通过破坏100个基因新得到。本研究小组已经采取了系统的实验方法,观察由于基因破坏引起微阵列上基因表达水平的变化。通过使用这种方法,本发明人已进行一项旨在揭示酿酒酵母的5871个基因之间的基因调控网络的项目,而一些实验室也已经报到了相似的研究。本发明人已经从基因破坏试验中收集大量表达图谱来评价遗传调控网络。备置400个以上的突变体,并积累基因表达图谱。
用扫描仪监控打点在微阵列上的5871个基因的转录水平。400个以上的破坏体的表达图谱被收集在本发明数据库中。评定微阵列上所有基因表达水平的标准偏差(SD),SD值粗略地表示试验的误差。在本发明的数据中,以0.5作为试验精确性的临界点。本发明人已经根据所有基因的表达比的标准偏差,评价了那些表达图谱的精确性。从400个图谱中选择107个破坏体,其包括68个转录因子被破坏的突变体。
采用100个微阵列,从上述数据中构建521个基因的基因网络。发现94个其调控基因已经被清楚确定的转录因子。由该94个因子控制的521个基因的表达图谱选自5871个表达图谱。
Bas1p和Bas2p还激活组氨酸生物合成路径中的3个基因的表达。在gcn4的背景中,消除BAS1或BAS2的功能的突变体会产生组氨酸营养缺陷。先前的研究表明Bas1p和Bas2p是HIS4以及如GCN4等ADE基因[8,11,29]转录所需的DNA结合蛋白。在本申请中说明了这两个基因关系。图4表明这些ADE基因和组氨酸生物合成基因与BASI之间比GCN4更直接相关。嘌呤核苷酸的核糖成分来自核糖5-P,其是戊糖磷酸化周期的中间体。碱基部分的分子由一些化合物提供。他们被逐步加入来形成核糖。其与组氨酸合成路径存在显著的相互影响。
对酿酒酵母中的嘌呤生物合成路径的调控的研究表明所有编码AMP重新生物合成所需酶的基因在转录水平上被细胞外嘌呤的存在所抑制。ADE基因和一些组氨酸生物合成基因在转录水平上被激活。特别地,HIS4的表达与ADE基因相关的事实是已知的。在本发明的调控网络中,HIS4与一些ADE基因密切相关,而一些HIS基因如HIS4基因与ADE基因相关。在酿酒酵母中表现的必需氨基酸组氨酸的生物合成与嘌呤代谢密切相关,本发明的结果符合该事实。
图5.分析521个酿酒酵母基因得到的部分网络5.结果本申请提出了一种用贝叶斯网络和非参数回归从微阵列基因表达数据中评估基因网络的新统计学方法。本发明方法的关键点在于采用非参数异方差回归模型来捕获基因间的非线性关系和表达数据的异方差性。如果具有代表基因间因果关系的网络,可以在计算机上模拟该基因网络,例如基因组对象网(Genomic Object Net)[26,27]。其中要求基因间的关系被恰当地估计。从这一点上看,由于已有的模型有时会产生对系统的不恰当估计,因此所提出的异方差模型能够带来关键的改进。本发明人将生物系统模拟作为未来的工作。
网络构建的关键问题在于结构图的评价。本发明将该问题作为统计模型选择或估计的问题来研究,并从贝叶斯方法中得到了新的选择结构图的检验数。本发明方法包括了采用贝叶斯网络构建基因网络的前述方法,并从理论上和方法学上改进了这些方法。所提出的方法成功地提取了有用的信息,并且可以在得到的基因网络中观察到这些信息。但是,这种算法需要许多时间来确定最佳结构图。因此,开发更好的算法是非常重要的问题之一,将在以后的文章中讨论。
本发明人已经通过蒙特卡罗模拟和分析酿酒酵母基因表达数据来说明本发明方法的有效性,并通过与生物学知识比较来评价得到的网络。本发明人不用生物信息就构建了基因网络,而且,得到的网络包括了许多重要相关关系,这些相关关系与生物学知识吻合。因此,本发明人认为本发明方法可以用于分析完全未知的系统如人类基因组。
参考文献[1]H.Akaike,信息理论和最大似然法则的扩展,B.N.Petrov,& F.Csald,eds.,Akademiai Kiado,Budapest,267,1973。
H.Akaike,似然与贝叶斯方法,J.M.Bernard,M.H.DeGroot,D.V.Lindley & A.F.M.Smith,cds.,Univ.Press,Valencia,41-166,,1980。
T.Akutsu,S.Miyano和S.Kuhara,根据布尔网络模型Proe从少量基因表达图谱确定基因网络,1999年太平洋生物计算研讨会(Pacific Symposiumon Biocomputing 1999),4,17-28,(1999)。
T.Akutsu,S.Miyano和S.Kuhara,推导基因网络和代谢途径中的定性关系,生物信息学,16,727-734,2000。
T.Akuten,S.Miyano和S.Kuhara,基于矩阵乘法和指纹功能的用于确定布尔网络和相关生物学网络的算法,计算生物学杂志(J.ComputBiol).7,331-344,2000。
J.O.Berger,统计决定理论与贝叶斯分析,Springer-Verlag纽约,1985。
R.Cowell,A.Dawid,S.Lauritzen和D.Spiegelhalter,或然网络与专家系统,Springer Verlag纽约,1999。
B.Daignan-Fornier和GR.Fink,转录激活子BAS1和BAS2共调控嘌呤和组氨酸生物合成,美国国家科学院院报(Proc.Natl.Acad.Sci.USA),89,6746-6750,1992。
A.C.Davison,近似预言似然,Biometrika,73,323-332,1986。
C.deBoor,样条函数使用手册,Springer-Verlag.柏林,1978。
V.Denis,H.Boucherie,C.Monribot和B.Daignan-Fornier,Myb-样蛋白BasIp在酿酒酵母中的作用蛋白质组分析.分子微生物学,30,556-566,1998。
N.Friedman和M.Goldszmidt,在贝叶斯网络学习中离散连续特征值.第13届机器学习国际会议记录,157-165,1996。
N.Friedman和M.Goldszmidt,具有局部结构关系的学习贝叶斯网络.第12届人工智能中的不确定性的讨论会会议记录,252-262,1996。
N.Friedman和M.Goldszmidt,具有局部结构关系的学习贝叶斯网络,M.l.Jordan ed.,Kluwer Academic Publisher,1998。
N.Friedman,M.Linial,I.Nachman和D.Pe′er.N.Friedman,M.Linial,L.Nachmank和D.Pe′er,用贝叶斯网络分析表达数据.计算生物学杂志,7,601-620,2000。
P.J.Green和B.W.Silverman,非参数回归和归纳线性模型,Chapman &Hall,1994。
A.J.Hartemink,D.K.Gifford,T.S.Jaakkola和R.A.Young(2002).结合测定和表达数据来原则地发现遗传调控网络模型.太平洋生物计算研讨会记录,7,437-449,2002。
T.Hastie和R.Tibshirani,归纳附加模型,Chapman & Hall(1990).1990。
D.Heckerman,贝叶斯网络学习指南,M.I.Jordan cd.,KtnwcrAcademic Publisher.1998。
D.Heckerman和D.Geiger.贝叶斯网络学习离散与高斯结构域的结合,第11届人工智能中的不确定性的讨论会会议记录,274-284,1995。
D.Heckerman,D.Geiger和D.M.Chickering,贝叶斯网络学习知识与统计数据的结合,机器学习,20,197-243,1995。
S.Imoto,T.Goto和S.Miyano,用贝叶斯网络和非参数回归估计基因间的基因网络和功能结构关系,太平洋生物计算研讨会记录,7,175-186,2002。
F.V.Jensen,贝叶斯网络导论.伦敦大学出版社,1996。
S.Konishi,统计模型评价和信息检验数,S.Ghosh ed MarcelDamer,1999。
S.Konishi和GKitagawa.模型选择中的归纳信息检验数,Biomerrika,83,875-890,1996,[26]H.Matsuno,A-Doi,M.Nagasaki和S.Miyano,基因调控网络的杂合Petri网表示,太平洋生物计算研讨会记录,5,338-349,2000。
R.Matsuno,A.Doi,Y Hirata和S.Miyano,生物路径XML资料及其在基因组目标网中的模拟,基因组信息学,12,54-62,2001。
D.Pe′er,A.Regev,G.Elidan和N.Friedman,从被干扰的表达图谱中推导子网络,生物性信息学,17,副刊,1,(ISMB 2001),215-224,2001。
R.J.Rolfes和A.G.Hinnebusch,嘌呤清除刺激酵母转录激活子GCN4的翻译隐含蛋白激酶GCN2的激活,细胞分子生物学,13,5099-5111,1993。
Schwarz,估计模型的维数,Ann.Sts.,6,461-464,1978。
L.Tinerey和J.B.Kadane,后验和边缘密度的精确逼近,美国统计学协会杂志,81,82-86,1986。
实施例4应用线性样条函数统计分析小组时间-排序基因表达数据摘要目的近来采用cDNA微阵列技术通过在少数时间点上测定基因表达水平研究了基因对环境变化的瞬时反应。传统的时序分析(time series analysis)方法不适用于这种短串时间-排序数据。因此,基因表达数据的分析通常被限定于倍数-变化分析,代替系统统计方法。
方法本发明人采用最大似然方法以及Akaike的信息检验数来将线性样条函数与小组时间-排序基因表达数据拟合,目的在于从测定值中推导具有统计意义的信息。用Student t-检验评定测定的基因表达数据的显著性。
结果用线性样条函数重新分析已有的藻青菌PCC6803的基因表达测定值。确定许多被倍数-变化分析遗漏的基因瞬时反应。根据本发明的统计分析,发现在各个时间点需要4个或者更多个基因表达的测定值。
联系mdehoon@ims.u-tokyo.ac.jp1.引言近些年来,已经进行了许多cDNA微阵列试验,测定不同条件下基因的表达水平。被测定的基因表达数据可以广泛地从公共数据库中获得,如KEGG数据库(Nakao等,1999)。
这些试验中的一些试验,在多种环境条件下测定稳态基因表达水平。例如,在不同温度下,测定藻青菌PCC6803和突变体的基因表达水平,结果确定基因Hik33为该藻青菌的潜在冷感受器(Suzuki等,2001)。
在其它试验中,通过在有限时间点上测定基因表达水平来估计基因表达的瞬时模式。例如,在酿酒酵母细胞周期中测定周期性变化的基因表达水平(Spellman等,1998)。同一酵母种属的基因表达水平在从发酵到呼吸的代谢变换过程中被测定(DeRisi等,1997)。在这些试验中,环境条件随时间缓慢变化。相反地,可以测定基因对突然变化的环境的反应。作为一个例子,在从低光照到强光照突然变化后,在多个时间点上测定藻青菌PCC6803的基因表达水平(Hihara,2001)。
在cDNA微阵列试验中,在少数时间点上代表性地测定基因表达水平。常规的时序分析方法,如傅立叶分析或者自回归或者移动平均数建模(moving average modelling),适合于这种少量数据点。此外,基因表达数据通常通过聚类技术或者通过只考虑基因表达水平中的相关变化来分析。这种“倍数-变化”分析可能会遗漏基因表达水平中的显著变化,而这可能无意中给受噪声支配的测定值带来显著性。此外,倍数-变化分析可能不能确定瞬时基因表达反应的重要特征。
多种分析基因表达数据的技术在过去已经被报道,例如获得的布尔或贝叶斯网络(Liang等,1998;Akutsu等,2000;Friedman等,2000)。然而按照调控网络描述基因相互作用是非常重要的,获得网络模型需要在大量时间点上的基因表达数据,而目前通常不能获得这些数据。应当指出基因的数目可达数千个次序,而通常仅在5个或10时间点上测定基因表达水平。
迄今,尚缺少统计分析来自少量时间-排序数据的基因表达测定值的系统方法。在本申请中,本发明人将略述基于用最大似然方法和Akaike信息检验数(Akaike,1971)将线性样条函数与时间-排序数据拟合的策略。用Studentt-检验评定基因表达测定值的显著性。这使得在考虑数据的统计学显著性的同时,可以从基因表达测定值推导信息。这种分析可以被看成是建立基因调控网络的第一步。作为一个例子,本发明人重新分析了藻青菌PCC6803的基因表达测定值(Hihara,2001)。表明可以从仅考虑倍数-变化时遗漏的测定数据中推导信息。通过重复分析可以得到的数据的子集,可以确定在各个时间点上需要多少个测定值才能合理地评估线性样条函数。
2.方法Student t-检验首先,应评价测定的基因表达比是否显著区别于1。如果对于特定基因来说,得出在所有时间点上所测定的表达比与1的差别都不显著的结论,可以在以后的分析中排除该基因。可以对各个时间点分别使用Student t-检验来确定显著性水平。由于对每个基因进行多重比较,应当谨慎选择显著性水平值α。
定义H0(j)为对于特定基因在特定时间点tj的表达比为1的假设,而H0为特定基因在所有时间点的表达比都等于1的假设。如果用α表示显著性水平用于否定假设H0,用α`表示否定假设H0(j)的显著性水平,那么α`与α相关,1-α=(1-α`)a(1)其中a为测定基因表达比的时间点的数量。通过扩展一阶泰勒级数(firstorder Tayler series)的右边,将该方程降为Bonferoni方法来调节显著性水平(又见Anderson和Finn,1996)。
以α`为显著性水平,对每个基因在各个时间点上进行Studentt-检验,结果无论H0(j)还是H0都将不成立。如果H0成立,则可得结论是该基因不会显著地受试验处理影响,该基因不应当被包括在进一步分析中。如果对于特定基因,虚假设H0不成立,结论是该基因显著地受试验处理影响。
2.2用线性样条函数分析时间-排序数据接着分析被发现显著地受影响的基因的瞬时基因表达反应。测定的基因表达比构成了小组时间-排序数据,可以用线性样条函数拟合这些数据。线性样条函数是由分段线性函数组成的连续函数,这些分段线性函数通过节点相互连接(Fiedman和Silverman,1989;Higuchi,1999)。然而立方样条函数更加常用,对少量数据点采用线性样条函数处理更为合适。图1给出具有节点t*1、t*2、t*3、t*4的线性样条函数的概念性例子。
考虑一组数据点(tj,Xj),j∈{1,...,n}。将用非参数回归模型xj=g(tj)+ξj(2)拟合这些数据,其中g为线性样条函数,εj是平均数为零、方差为σ2的正态分布的随机自变量(independent random variable)。
用最大似然方法来估计线性样条函数g。特定tj的数据点xj的概率分布为
f(xj|tj;g,σ2)=12πσ2exp{-(xj-s(tj))22σ2}.---(3)]]>该n个数据点的log-似然函数为L(g,σ2)=-n2ln[2πσ2]-12σ2Σj=1n(xj-g(tj))2.---(4)]]>通过最大化log-似然函数的σ2得到方差σ2的最大似然估计。从而得到σ^2=1nΣj=1n(xj-g(tj))2.---(5).]]> 图1.拟合测量数据的线性样条函数草图。
然后log-似然函数可以写成如下形式L(g,σ2=σ^2)=-n2ln[2πσ^2]-n2.---(6)]]>通过最小化ε2可以确定线性样条函数g的最大似然估计 如果选择线性样条函数,将会得到ε2的最小值,以使A‾g^‾=b‾,---(7)]]>其中 是含有在节点ti*的线性样条函数的q值的向量。 是如下的三对角线对称矩阵(tridiagonal symmetric matrix) for1<i<q,---(9)]]> for1≤i<q,---(11)]]>而 是如下得出的向量
for1<i<q;---(13)]]> 拟合模型依赖于节点q的数量。可以使用称为AIC的Akaike信息检验数(Akaike,1971;以及Priestly,1994)来选择节点数量。
AIC=-2L(g^,σ^2)+2(q+1),---(15)]]>其中q+1是被估计参数(ε2和估计向量 )中的q条目(entries)的数量)。从方程(6)置换估计log似然函数,发现AIC=nln[2πσ^2]+n+2q+2,---(16)]]>其中ε2在线性样条函数g被最大似然估计 置换后从方程(5)中得到。
对于每个q值,在如上所述拟合线性样条函数后计算AIC的值,选择得到AIC最小值的q值。q=2的情况对应线性回归。对于特定情形q=1,有效地用直线拟合了这些数据。如果发现对于某一特定基因,从恒定函数(q=1)得到最小AIC,那么可以认为该基因的表达水平没有受到试验操作的影响。
基因表达数据通常是以表达比的形式给出。定义零时的表达比为1。通过修改(7)可以容易地将该固定点包纳入本发明的方法中。通过选择如下所示线性样条函数可以获得ε2最小值A22Λ23000A32A33Λ34000A43A4400000Aq-1,q-1Aq-1,q000Aq,q-1Aq,qg^2g^3g^4g^q-1g^q=b2′b3b4bq-1bq,---(17).]]>其中b2′=b2-A21,以及g^1=1.]]>3.结果3.1 Student t-检验通过重新分析藻青菌PCC6803在突然暴露在强光(HL)后所测定基因的表达图谱(Hihara等,2001),在此说明应用Student t-检验和拟合线性样条函数。对暴露于HL的藻青菌以及保留于低光照(LL)条件下的藻青菌,在零时、15分钟、1小时、6小时和15小时测定3079个ORF的表达水平。表1说明每个时间点的测定值的数目。来自cDNA表达测定的数据从KEGG数据库获得(Nakao等,1999)。
表1.作为时间-排序数据的基因表达测定。
时间点 测定次数t=15分钟4(*)t=1小时 6t=6小时 4t=15小时4(*)在t=15分钟的时间点上,KECG数据库中缺少的两组测定值表2.基因表达测定值的Student t-检验(共3079个ORF)。
显著性水平 ORF数量p<0.0003980.0003<p<0.001 690.001<p<0.01 3200.01<p<0.05517p>0.05 2075应当指出,用作原始分析的数据(Hihara,2001)可能不同于被提交到KEGG的原始数据(Hihara,个人信息)。此外,KEGG数据库中缺少在t=15分钟时间点的6个测定组中的两组。重新从原始数据中计算基因表达比得到与先前公开的结果接近的数值。
在从HL和LL的原始数据中减去背景信号强度后,实行整体归一化,计算HL数值与LL信号强度的比值来确定相对于对照条件(LL)的基因表达水平的相对变化。在倍数-变化分析中,如果基因表达水平被两个或多个因子之一改变,该基因被认为受HL作用。通过假定这些测定值标准偏差的大小,启发式地估计这种变化的统计学显著性(Hihara,2000)。
对各个基因的表达比分别进行Student t-检验的结果在表2中给出。在显著性水平α=0.001时,发现167个基因显著地受HL条件影响。在这167个基因中可以预期大约3个1类错误。与此相比,在最初分析中发现164个ORF受HL条件影响(Hihara,2001)。
通过倍数-变化分析psbD2基因(slr0927),结果表明该基因未显著地被HL诱导(Hihara,2001)。由于该基因已经被报道在藻青菌PCC7942中受HL诱导(Bustos和Golden,1992;Anandan和Golden,1997),所以上述结果值得关注。但是对psbD2基因在t=6小时的基因表达数据进行Student t-检验,结果p=3.3×10-5,表明该基因确实受HL影响。
3.2应用线性样条函数分析接着用线性样条函数拟合测定的基因表达比。节点数q在1到5之间,具有在零时等于1的固定节点。当q=3和q=4时,线性样条函数的线性片段之间的节点的定位存在3种可能。这在图2中说明,以及q=1、q=2和q=5的情况。表明可能节点的定位的数量随着最大节点数qmax以1+2qmax-2指数性地增加。
在倍数-变化分析中,瞬时基因表达模式分成6类(Hihara,2001),在表3中列出。将线性样条函数与测定的基因表达数据拟合,提供一种比分类更灵活的方法来描述数据。此外,基因表达反应图谱的数字化描述是获得基因调控网络的重要起始步骤。
图2时间排序组的测量值在5个时间点上的节点的可能定位。
表3在倍数-变化分析中用于分类基因瞬时表达模式的六种分类(Hihata,2001)。
Type 1 15分钟内被诱导,然后下降Type 2 持续高水平被诱导
Type 3约1小时时被诱导Type 415分钟内被抑制,然后上升Type 5持续低水平被抑制Type 6约1小时时被抑制举例说明AIC的使用。在倍数-变化分析中,发现苏氨酸合成基因thrC(sll1688)大约在1小时时被抑制。表4中列出不同组节点的AIC的计算值。在0、15分钟、1小时和15小时时获得节点的最小AIC。图3表明测定的基因表达水平以及与数据拟合的线性样条函数。
对所有基因表达测定实施该处理,根据对HL的时间-依赖反应来对不同基因进行分类。至于哪些基因被包括在该分析中可以有多种选择。在最初分析中,如果基因的表达水平在3079个ORF中的2000个最低中(Hihara,2001),则从计算中被排除。或者,如果Student t-检验表明基因未显著地受HL影响,也排除该基因。表5表示基因的数量,这些基因的测定的表达水平对应于这些不同情形的各个应答模式。采用显著性水平α=0.001进行Student t-检验。
表4.苏氨酸合成基因thrC(sll1688)的不同组节点的计算的AIC节点的定位 AIC直线 62.55(0,15小时) 64.54(0,15分钟,15小时) 64.54(0,1小时,15小时) 66.24(0,6小时,15小时) 66.38(0,15分钟,1小时,I5h) 25.24(0,15分钟,6小时,15小时) 64.69(0,1小时,6小时,15小时)68.23(0,15分钟,1小时,6小时,15小时)26.45
图3.苏氨酸合成酶基因thrC(sll1688)的测量的基因表达水平,以及拟合的线性样条函数。
表5.各个反应模式的基因数量基因数量节点定位 保留所有基排除最低的用Student t-检验因2000个基因排除的基因直线(0,15小时)(0,15分钟,15小时) 927 249 1(*)(0,1小时,15小时) 280 693(0,6小时,15小时) 663 217 61(0,15分钟,1小时,15386 178 7小时)160 643(0,15分钟,6小时,15296 119 27小时)178 8339(0,1小时,6小时,15 103 4912小时)865114(0,15分钟,1小时,6小时,15小时)(*)由于存在无关值被排除比较Student t-检验确定为显著地受HL影响的、具有倍数-变化分析结果(Hihara,2001)的167个基因,其中164个ORF被确定。首先,将数据中存在无关值的基因从分析中排除。无关值定义为偏离特定时间点的数据的平均值两个以上标准偏差的数据点。只有一个基因被发现其测定的表达数据中含有无关值;与该表达数据拟合的线性样条函数是一条直线。其它基因表达水平都不能用直线描述,这与Student t-检验的结果一致。
接着排除表达水平在最低的2000个之中的基因,以避免使用受噪声支配的数据。同样的处理被用于倍数-变化分析(Hihara,2001)。在排除了这些基因后,保留显著受HL影响的107个基因。
在该107个基因中,有42个未在倍数-变化分析中被确定(Hihara,2001)。表6中列出这些基因以及为各个基因确定的节点定位。对于每个线性样条函数,计算所说明方差的百分比,作为对拟合程度的检测。作为一个例子,图4显示测定的基因表达比以及基因sylR(slr0329)的拟合线性样条函数,具有零时、15分钟、1小时和15小时的4个节点。在这42个基因中,基因xylR(slr0329)具有所说明方差的最大百分比(98.7%)。
在倍数-变化分析确定的164个ORF中,根据根据显著性水平α=0.001的Student t-检验,其中39个受HL影响不显著。这些ORF在表7中列出。
表6.显著受HL影响的ORF,但是未在倍数-变化分析中被鉴定(Hihara,2001)。
ORF 基因 节点定位 %所说明方差sll1330(0,15分钟,15小时)91.8sll1797 yct21 (0,15分钟,I5h) 61.4slr0121(0,15分钟,15小时)78.7slr0343 petD (0,15分钟,15小时)93.5slr0442(0,15分钟,15小时)81.1slr0769(0,15分钟,15小时)64.7slr0906 psbB (0,15分钟,15小时)94.6slr1128(0,15分钟,15小时)63.7slr1152(0,15分钟,15小时)60.7slr1277 gspD (0,15分钟,15小时)40.5slr1610(0,15分钟,15小时)70.2slr18I3(0,15分钟,15小时)84.3slr2052(0,15分钟,15小时)87.1sll0144 pyrH或smbA(0,1小时,15小时) 88.5
sll1327 atpC (0,1小时,15小时) 70.3sll1496(0,1小时,15小时) 77.4slr1367 glgP或glgY(0,1小时,15小时) 88.2slr1793 talB (0,1小时,15小时) 87.0sll1878(0,6小时,15小时) 65.8slr1600(0,6小时,15小时) 46.6sll0851 psbC (0,15分钟,1小时,15小时) 91.9s1l1471 cpcG (0,15分钟,1小时,15小时) 96.6sll1812 rFs5 (0,15分钟,1小时,15小时) 60.3slr0076(0,15分钟,1小时,15小时) 85.8slr0329 xylR (0,15分钟,1小时,15小时) 98.7slr0927 psbD2 (0,15分钟,1小时,15小时) 48.9slr1513(0,15分钟,1小时,15小时) 94.8sll0547(0,15分钟,1小时,15小时) 79.3sll1528(0,15分钟,1小时,15小时) 90.9slr0056 G4(0,15分钟,1小时,15小时) 76.6slr0161 pilT (0,15分钟,6小时,15小时) 96.2slr1239 pntA (0,15分钟,6小时,15小时) 96.7slr1545(0,15分钟,6小时,15小时) 96.7slr1799(0,15分钟,6小时,15小时) 95.5sll0617 hr30 (0,1,6小时,15小时) 84.9Sll0921(0,1,6小时,15小时) 46.9sll0985(0,1,6小时,15小时) 89.2slr1237 codA (0,1,6小时,15小时) 80.7sll0108 amt1 (0,15分钟,1小时,6小时,15小时) 84.4sll0185(0,15分钟,1小时,6小时,15小时) 94.9sll1873(0,15分钟,1小时,6小时,15小时) 92.0slr0116 hemD (0,15分钟,1小时,6小时,15小时) 89.2
图4.基因xylR(slr0329)的测量的基因表达水平,以及拟合的线性样条函数。在倍数-变化分析中,该基因被认为未受HL的作用(Hihara,2001)。
表7.在倍数-变化分析中确定的ORF(Hihara,2001),但是Student t-检验无显著性。
Slr1311 sll1867 sll1732sll1733 sll1734 slr1280sll0170 sll0430 sll1514slr1641 slr1350 slr1516slr0228 slr7604 slr1675ssl3044 stl1483 sll1541sll2012 slr0394 slr0426sIr1351 slrI594 sll0217sll0218 sll0219 sll0528sll0668 sll0814 sll0846sll1884 slr0007 slr0476slr0740 slr0959 slr1544slr1674 slr1686 slr1963表8.线性样条函数拟合处理作为所用数据数量的函数的可靠性。仅考虑那些根据Student t-检验确定为显著受HL影响的基因,它们的数据中不含有任何无关值。
在15分钟、6小时和15小时用4估计166个基因的线性样条函个数据点,数在1小时用6个数据点各个时间点用4个数据点 区别地估计25±8各个时间点用3个数据点 区别地估计53±10各个时间点用2个数据点 区别地估计79±17最后确定在各个时间点的测定值的数量是否足以可靠地确定线性样条函数的节点定位。为此,重复线性样条函数的估计。为此,用测定数据的子集重复线性样条函数的估计。然后计算如果用数据的子集替代整组数据,所估计的节点位置发生改变的基因数量。表8给出每个时间点数值的4个、3个和2个数据点的平均数和标准偏差。
即便仅有2个数据点(在1小时时间点)被排除,每个时间点采用4个数据点,15%的情形估计的节点位置发生改变。这表明每个时间点需要4个或者更多数据点来可靠地从基因表达测定值中推导信息。
4.讨论本发明已经描述一种基于最大似然方法的策略来分析一组时间-排序测定值。对测定的基因表达数据进行Student t-检验,首先确定被测基因中的哪些基因显著地受试验处理影响。然后通过拟合线性样条函数来描述那些基因的表达反应。用Akaike信息检验数(AIC)确定用于线性样条函数的节点数量。
在描述测定的基因表达中,采用线性样条函数比标准分类更具灵活性。而且,为了建立基因调控网络,从基因表达测定值中确定的基因反应以数字化形式获得是十分重要的。最后,节点的位置清楚表明基因表达水平显著变化的那些时间点,这对于确定其生物学功能是重要的。作为后续步骤,根据节点位置的基因表达水平的分类,可以通过建立亚类来细化,所述亚类考虑了节点处线性样条函数数的大小(magnitude)。例如,可以对具有3个节点的线性样条函数建立6个亚类,其中基因表达水平用(平缓,上升)、(平缓、下降)、(上升,平缓)、(下降、平缓)、(上升、下降)或(下降、上升)来描述。
应用线性样条函数方法来测定基因表达数据,可以确定显著地受试验处理影响的基因的瞬时表达反应模式。那些基因中42个基因的反应在较早的表达数据的倍数-变化分析中未被提出。而且,在倍数-变化分析中确定的164个基因中的33基因的表达反应水平,Student t-检验确定为不显著。
基因表达数据趋向于为噪声并且经常受到无关值干扰。然而本申请描述的Student t-检验和最大似然方法考虑噪声数据的统计学显著性,无关值的问题需要分别地处理。作为排除无关值的简单方法,计算各个时间点的数据的平均值和标准偏差,排除那些偏离平均数超过大约2个标准偏差的数据。
最后,在各个时间点上用于可靠地拟合线性样条函数的表达测定值的数量,通过排除一些数据点和重新拟合线性样条函数来确定。已经发现,如果每个时间点采用4个数据点,大约15%的情形不能可靠地估计节点定位。因此,建议每个时间点进行4次以上测定。
参考文献Akaike,H.(1971)信息理论和最大似然法则的扩展,第46号研究备忘录,统计数学协会,东京,Petrov,B.N.,Csaki,F.(eds.),第2届信息理论国际讨论会会议记录,Akademiai Kiado,布达佩斯(1973),第267-281页。
Akutsu,T.,Miyano,S.和Kuhara,S.(2000)推导基因网络和代谢途径中的定性关系,生物信息学,16,727-734。
Anandan,S.和Golden,S.S.(1997)酿酒酵母菌株PCC 7942中psbDII基因的光反应所需的反作用序列,细菌学杂志179,6865-6870。
Anderson,T.W.和Finn,J.D.(1996)新的数据统计分析,Springer-Verlag,纽约。
Bustos,S.A.和Golden,S.S.(1992)酿酒酵母菌株PCC 7942中的psbD基因家族的光调控表达复制的psbD基因在藻青菌中的作用的证据.分子遗传学,232,221-230。
DeRisi,J.L.,Iyer,V.R.和Brown,P.O.(1997)基因组水平上探讨基因表达的代谢和遗传控制,科学,278,680-686。
Friedman,J.H.和Silverman.B.W.(1989)灵活收敛的光滑和附加建模(讨论中),技术计量学,31,3-39。
Friedman,N.,Linial,M.,Nachman,I.和Peter,D.(200D)用贝叶斯网络分析表达数据,计算生物学,7,601-620。
Higuchi,T.(1999)自动确定大规模试验排列现有系统作为从大数据库中知识发现的例子(日文),统计数学协会的学报,47,291-306。
Hihara,Y,Kamei,A.,Kanehisa,M.,Kaplan.A.和Ikeuchi,M.(2001)在适应强光照过程中,藻青菌基因表达的DNA微阵列分析,植物细胞,13,793-806。
Liang,S.,Fuhrman,S.和Somogyi,R.(1998)展示用于推导基因网络体系的通用反向工程算法,太平洋计算生物学研讨会会议记录,3,18-29。
Imoto,S Goto,T.和Miyano,S.(2002)用贝叶斯网络和非参数回归估计基因间的基因网络和功能结构关系,太平洋生物计算研讨会记录。
Nakao,M.,Bono,H.,Kawashima,S.,Kamiya,T,SatoasK,Goto,S.和Kanehisa,M.(1999)KEGG中的基因组水平的基因表达分析和路径重组,Dunker,A.K.,Konagaya,A.,Miyano,S.,Tiakagi,T.(eds),基因组信息学第10次讨论会记录,大学学术出版社,东京,日本,第94-103页。
Priesdey,M.B.(1994)光谱分析和时序.学术出版社,伦敦。
Spellman,P.,Sherlock,G.,Zhang,M.Q.,Iyer,V.R,Anders,K.,Eisen,M.B.,Brown,P.O.,Botstein,D.和Futcher.B.(1998)微阵列杂交综合确定酿酒酵母中的细胞周期调控基因,细胞分子生物学,9,3273-3297。
Suzuki,I.,Kanesaki,Y.,Mikami,K.,Kanehisa,M.和Murata,N.(2001)Synechocystis中的冷感受器Hik33控制的冷调控基因,分子微生物学,40,235-244。
关键词基因表达数据,线性样条函数,AIC,时间-排序数据,最大似然方法。
实施方案5应用基因网络确定和评价药物靶点本发明人提出了一种用基因网络确定和评价药物靶点新方法,该基因网络从cDNA微阵列基因表达数据估计得到。本发明建立新的基因破坏和药物反应微阵列基因表达数据文库来阐明药物靶点。用两种微阵列基因表达数据来评估基因网络,然后确定药物靶点。估计的基因网络在理解药物反应数据中起重要作用,这些信息不能从聚类分析方法中获得,聚类分析是基因表达分析的标准方法。在基因网络的构建中,在分析的不同步骤中采用布尔和贝叶斯网络模型,捕获它们的相对强度。用分析酿酒酵母基因表达和药物反应数据的实例来评价本发明的将基因网络信息应用于药物发现的具体策略。
1.引言近些年来,微阵列技术的发展产生了大量各种试验条件下的基因表达数据。从基因表达数据中构建基因网络已经引起极大关注,并且在计算机科学和统计学领域多种用于基因网络推导的方法已被报道,例如见参考文献部分的文章1,2,3,4,9,11,13,14,16,20,22。构建基因网络的新方法的发展十分重要并且应当被进一步发展。另一方面,估计的基因网络实际应用于解决医学和生物学难题近来已成为关注的焦点。在本申请中,将介绍确定和评价基于基因表达数据的药物靶点的方法学。
聚类方法6,7,23已经成为分析微阵列基因表达数据的标准方法。但是,无论从理论上还是从实践上看,这些方法都不能提供足以确定药物靶点的信息。本申请将说明如何将估计的基因网络用于确定和评价目标基因的方法,推定和开发新的治疗药物。基因调控路径信息对于实现本发明目的是关键的,采用布尔和贝叶斯网络作为从基因表达图谱中推导基因网络的方法。确定药物靶点的过程分成两部分。首先,确定药物作用基因。其次,搜索可药化基因,在基因网络中该基因通常位于药物作用基因的上游。通过采用本申请提出的“虚拟基因”技术,布尔网络模型对于确定药物作用基因是最有用的,而贝叶斯网络模型对于搜索与阐明的作用基因相关的可药化基因靶点十分有效。将本发明的方法应用于新的酿酒酵母表达数据,这些数据包括来自120个基因破坏和多种剂量和对药物响应的时间的表达试验。通过本申请研究来说明用于确定药物靶向的基因的具体策略。
2.确定药物靶点的基因网络2.1聚类方法聚类方法如分级聚类和自组织图23(self organizing map)是生物信息学领域中常用的用于基因表达数据分析的标准工具。Eisen8专注于分级聚类,并且提供用于基因表达数据聚类分析的软件Cluster/TreeVew。De Hoon等5,6特别是在K-均值聚类算法方面改进了该软件。聚类算法在分析的早期阶段是有效的。在基因表达数据的聚类分析中,假设具有相似表达图谱的基因在细胞中起相同的作用。因此,聚类方法常被用于预测基因功能。
本发明人认可聚类方法的有用性,但是聚类方法不能提供用于本发明目的的充足信息。聚类方法仅通过表达图谱的相似性提供基因群(gene group)的信息。事实上,为了检测受试剂作用的药物靶点,需要额外的分级路径信息,而不仅仅是聚类信息。在第3章节中,将通过实际数据分析来说明聚类技术用于药物靶向目的的局限性。
2.2网络方法本发明采用两种方法从基因表达数据中评估基因网络。在该分部中,简单介绍这两种方法。对于这些算法的详细讨论,请参见参考文献部分的文章1,2,3,4,9,11,13,14,16,20,22布尔网络许多研究小组1,2,3,4,15,22已经提出以布尔网络模型作为基因网络构建方法。为了估计布尔网络模型,必须将基因表达值离散成两个级别,0(未被表达)和1(被表达)。假定u1,...uk为节点v的输入节点。V的态通过布尔函数fv(ψ(u1),...,ψ(uk))确定,其中ψ(ui)是ui的态或者表达模式。如果有时序基因表达数据,该态依赖于时间t,以及t时的节点的态依赖于时间t-1时的输入值的态。另一方面,假定具有得自基因破坏的基因表达数据。Akutsu等1提出估计没有延时的布尔网络模型的理论和方法学。Maki等16提供名为AIGNET的系统,用于估计基于布尔网络和S-系统21的基因网络。本发明采用AIGNET系统估计布尔网络模型。
采用布尔网络模型的优点包括a)该模型简单,容易被理解。b)当这些数据足够精确和信息时,布尔网络模型能够正确地测定亲本-子代关系和c)估计的布尔网络模型可以直接应用于基因组对象网10,16,19(Genomic ObjectNet),一种用于生物路径模拟软件工具。布尔方法的缺点在于必须将数据离散成两个级别,而该量化可能导致信息丢失。而且用于离散的阈值是一个参数,需要用合适的检验数来选择。
贝叶斯网络贝叶斯网络15是大量随机变量的复杂关系的图形表示。在贝叶斯网络中,假定定向非循环结构图具有节点间的马尔可夫关系。通过条件概率替代随机变量的联合概率来描述复杂现象。也就是,假设有第i个阵列的第j个基因的基因表达数据xij,i=1,...,n和j=1,...,p。令f(xi1,xi2,...,xip)=f1(xi1|pi1)f2(xi2|pi2)×...×fp(xip|pip),其中pij=(pil(j),...,piqj(j))T是xij的qj-维亲本观察值向量。
Friedman等9提出一种从基因表达图谱评估基因网络的有益方法。他们将表达值离散成3个值,并且以多项式分布作为贝叶斯网络的条件分布。但是,该方法没有解决为离散选择阈值的问题。最近,Imoto等13,14采用非参数回归模型提供一种解决办法,不需要量子化xij=mj1(pi1(j))+…+mjqj(piqj(j))+ϵij,]]>其中εij单独地和一般地依赖于平均数0和方差σij2,而mjk(·)是由B-样条函数构建的光滑函数,表示为mjk(pik(j))=Σm=1Mjkγmk(j)bmk(j)(pik(j)),]]>k=1,...,qj。此处γ1k(j),...,γMjk(j)k是未知系数,而b1k(j)(·),...,bMik(j)k(·)为B-样条函数。Imoto等13,14提出非线性贝叶斯网络模型f(xi;θ)=Πj=1pexp[-{xij-Σk=1qjΣm=1Mjkγmk(j)bmk(j)(pik(j))}2/2σj2]/2πσij2,]]>以及用于选择最佳结构图的新检验数BNRC。BNRC定义为用拉普拉斯逼近求积分的结构图的后验概率近似值。Imoto等14将所提出的方法应用于酿酒酵母基因表达数据评估基因网络。这种方法的优点在于a)可以将微阵列数据作为连续数据集分析,b)该模型不仅可以检测线性结构关系而且还可以检测基因间的非线性相关性和c)被提出的检验数能够自动地优化模型中的参数和网络结构关系。
在本申请中,贝叶斯网络具有一些基于推导数学的理论优点,所以我们采用了Imoto等13,14提出的方法来构建贝叶斯网络。不幸的是,贝叶斯网络不能构建环状调控,对于从数据中构建调控作用的多水平方向模型(multilevel directional model)是无用的,所述数据通过逻辑联合来自破坏体和药物反应试验的表达数据来建立。但是,可以在分析中结合布尔网络和贝叶斯网络来解决该问题。因此,同时使用布尔网络和贝叶斯网络能够逐个地克服这些缺陷,而且使用这两种网络方法能够获得更可靠的信息。
3.应用3.1微阵列数据从酿酒酵母基因表达图谱建立两个微阵列数据文库。一个通过破坏120基因获得,另一个含有对口服抗真菌剂的反应(4种浓度和5个时间点)。从酵母基因组中选择735个基因用于确定药物靶点。在YPD中,314个基因被定义为“转录因子”,其中98个基因的控制机制已经被研究。选择用于分析的735个基因的表达图谱数据包括受这98个“转录因子”控制的基因,这些基因来自除了核受体基因之外的5871个被测的基因种类,所述的核受体基因在基因表达调控中发挥作用,是常见的药物靶点。从在超过120个基因破坏的条件下的735个基因的数据集来构建基因网络模型。破坏的数据的细节还在Imoto等14中描述。
至于药物反应微阵列基因表达数据,将酵母培养物接种到抗真菌药物的剂量为10、25、50、100mg的培养基中,在试剂加入后的5个时间点(0、15、30、45和60分钟)取出培养物的等分试样。此处时间0表示刚被暴露于药物后观察的起始点。然后从这些试验中提取总RNA,用cy5标记RNA,用cy3标记的未处理细胞的RNA与之杂交,将它们应用于全基因组cDNA微阵列,从而产生20个微阵列数据集作为药物反应数据。在本申请中,应用基因网络以该140个微阵列阐明药物靶点。
3.2结果聚类分析在确定药物靶点中,一种常用但存在问题的策略是采用聚类分析12,17,通常还有基础干扰控制数据文库(library of base pertubation control data)以基于相关关系与药物反应数据比较。本发明具有两类微阵列数据,基因破坏和药物反应,可以将药物反应图谱与由破坏导致的基因表达图谱比较。在聚类分析中,如果单个破坏体或者破坏体群的表达图谱与特定药物反应微阵列之间显著而强烈地相似,则可以得出药物可能与被破坏的基因起相同作用的结论。而且,如果这种被破坏的基因具有已知功能作用,则可能获得更多关于对药物反应的信息。
图1.组合基因表达数据的相关矩阵R的图像。
不幸的是,情况往往是即使有这种试验结果,也不能从聚类分析数据中获得直截了当的结论。图1表示本发明微阵列数据的相关矩阵的图像表示。结合两类数据,生成矩阵Z=(X∶Y),其中X和Y分别是药物反应和基因破坏微阵列数据的矩阵。在此,每一列表示由一个微阵列得到的表达图谱,并且每一列用平均数0和方差1进行归一化。因此,图1是相关矩阵R=ZTZ/p的信息,其中p是基因的数量,为735。
在图1中,浅颜色和深颜色分别表示正和负的高相关性。显然,药物反应微阵列之间高度相关,而与任何基因破坏微阵列之间的相关性低。在这种情况下,从聚类分析中不能发现基因破坏和药物反应之间的任何相互作用,聚类分析对于确定强烈被作用的基因是有用的,而这些基因在药物反应中有意义。进一步对基因破坏和药物反应微阵列进行分级聚类,但是药物反应微阵列产生一个聚类,而从该分析中不能获取任何关于真正药物靶点的更多信息。当采用其它的距离测定和聚类技术时,该结果没有实质性改变。因此,不能够采用聚类方法从本发明数据中获得用于确定有意义药物靶点的信息。
布尔网络分析用微阵列数据Z来评估基因网络,微阵列数据Z是通过组合基因破坏微阵列和药物反应微阵列建立的。本发明人引入称为虚拟基因技术的方法。假设药物反应数据的条件为“虚拟基因”,例如条件100mg/ml和30min被给定一个赋值为基因YEXP100mg30min。通过应用布尔网络模型,可以确定这些虚拟基因的子代基因,药物以世代顺序作用于这些子代基因。将以5个或5个以上虚拟基因为亲本基因的基因视为假定的药物作用基因,也就是,直接受虚拟基因(药物作用)影响的基因。但是,仅只有一个虚拟基因作为亲本基因的基因可能是最初的药物作用基因,这依赖于特定药物的作用方式,这需要逐个地分析。在初步筛选受药物诱导表达影响的基因中,虚拟基因技术使布尔网络模型的使用比贝叶斯网络模型更加优越。
此外,倍数-变化分析能够提供与所提出的虚拟基因技术相似的信息。事实上,在特定实验条件下,可以通过倍数-变化分析确定被作用的基因。但是,本发明的虚拟基因技术能够改进倍数-变化分析的结果。假设通过倍数-变化分析发现基因A和基因B受药物作用。倍数-变化分析不能考虑基因A和基因B之间的基线相互作用。也就是说,基因A和基因B之间存在调控路径基因A→基因B,基因B可能不是直接受药物作用。虚拟基因技术通过利用基因破坏数据的信息能够考虑这种相互作用,从而将搜索点减少到可能性更大的目标基因,这些基因产生可以获得的相互作用数据。
关于最受试剂作用的基因是那些被该试剂“药化(drugged)”的基因的观点,以及关于被药化的靶点表示新药干扰的最具生物利用度和优势的分子靶点的观点都是毫无保障的。因此,即便确定了作用的可能分子模型,我们应当寻找位于调控网络中的药物作用基因上游的最可药化的(druggable)基因,然后筛选低分子量化合物对这些靶位点的药物活性。在估计的布尔网络中,应将虚拟基因置于网络的顶层。因此,在该估计的布尔网络中发掘药物作用基因的上游信息是困难的,有时是不可能的。此时,可以采用贝叶斯网络模型有效地搜索药物作用基因的上游区域。
图2.虚拟基因YEXP100mg30min的下游路径。
贝叶斯网络结果通过贝叶斯网络和非参数回归方法以及BNRC优化策略来评估基因网络13,14。使用得自破坏120个基因的酿酒酵母微阵列基因表达数据。从布尔网络分析中,可以有效地发现药物作用基因的候选组。可药化的基因是与这些药物作用靶点相关的药物靶点,其是为了新的先导药物的开发确定的。可以通过贝叶斯网络方法在估计的基因网络中搜索位于药物作用基因上游的可药化基因。采用从敲除表达文库中获得的网络调控数据的贝叶斯模型,搜索药物作用基因的上游来获得基因,这些基因具有已知的对药物作用靶点表达的调控控制关系。例如,将核受体基因视为可药化基因,因为a)已知核受体蛋白是有用的药物作用靶点,并且总共代表当前市场上药物的20%以上靶点。b)核受体参与转录调控作用,这些转录作用在cDNA微阵列试验中被直接测定。
图3.在可药化基因(顶层)、药物作用基因(底层)、和媒介基因(中间)中的部分网络。
图3表示部分网络,包括药物作用基因(底层)、可药化基因(顶层)和媒介基因(中间)。当然,如果容许更多媒介基因,会发现更多从可药化基因到药物作用基因的路径。由于使用贝叶斯网络模型,可以确定边缘强度和选择更可靠的路径。这是贝叶斯网络模型在搜索理想的可药化靶点的一个优点。在图3中,圈中的可药化基因直接与药物作用基因相连,其它每个可药化基因具有一个媒介基因。从图3可以确定每个药物作用基因的可药化基因,例如,发现了表1所述MAL33和CDC6的可药化基因。
表1.MAL33和CDC6的可药化基因。“亲本基因”是指这些基因与药物作用基因直接相连。“祖代基因(Grandparent)”是指在这些基因与药物作用基因之间有一个媒介基因。
药物作用的 MAL33(YBR297W)麦芽糖发酵调控蛋白可药化的亲本基因HSP82(YPL240C)热休克蛋白SRB4(YER022W)DNA指导的RNA聚合酶II全酶和Kornberg介质子(mediator)(SRB)亚复合体亚基祖代基因BAR1(YILO15W)抑胃蛋白酶(Barrierpepsin)前体GPA1(YHR005C)信息素路径的GTP-结合蛋白α-亚基KAR2(YJL034W)核融合蛋白药物作用的 CDC5(YJL194W)细胞分裂控制蛋白可药化的亲本基因ARP7(YPR034W)SWl-SNF总转录活化剂复合物成分和RSC染色质重构(remodeling)复合物BAR1(YIL015W)抑胃蛋白酶前体祖代基因GAL11(YOL051W)DNA指导的RNA聚合酶II全酶和Kornberg介质子(SRB)亚复合体亚基FAR1(YJL157C)细胞周期依赖的蛋白激酶抑制剂(CKI)SLA2(YNL243W)细胞骨架装配控制蛋白4.讨论本申请提出了用基因网络计算模型确定和评价药物靶点的新策略。说明了聚类分析方法不能提供足够信息,尚需通过网络方法产生分级相互作用数据。本发明关注于布尔网络和贝叶斯网络用于从微阵列基因表达数据评估基因网络。这两种方法最初是根据不同观点被提出的,但是,他们都具有相对优势和缺点,通过在本发明分析中同时使用这两种方法,我们能够获得更加可靠的信息。布尔网络理想地通过使用上述虚拟基因技术来确定药物作用基因。本申请还描述了虚拟基因技术相对于从单一倍数-变化分析中获得的结果的理论优势。将较简单的布尔网络用于特定目的是一种积极的做法。但是,本发明人说明当搜索位于估计网络中的药物作用基因的上游的可药化基因时,布尔网络不能产生可信的信息。通过应用贝叶斯网络模型我们解决了该问题。贝叶斯网络模型能够有效地提供关于药物作用基因上游区域的信息,从而获得各个药物作用基因的一组候选可药化基因。基于精确地结合使用这两种网络方法,建立了本发明的策略。在该策略中明显地体现各个网络的优势,而且所提出的综合方法在药物靶点的确定和评价中为用于基因网络推导的生物信息学技术的实际应用提供了方法学基础。
参考文献1.T.Akutsu,S.Kuhara,O.Maruyama和S.Miyano,从由基因破坏和过度表达产生的基因表达图谱中确定基因网络的系统,基因组信息学,9,151,(1998)。
2.T.Akutsu,S.Miyano和S.Kuhara,根据布尔网络模型从少量基因表达图谱确定基因网络,太平洋计算生物学研讨会会议记录,4,17,(1999)。
3.T.Akutsu,S.Miyano和S.Kuhara,推导基因网络和代谢途径中的定性关系,生物信息学,16,727,(2000)。
4.T.Akutsu,S.Miyano和S.Kuhara,基于矩阵乘法和指纹功能的用于确定布尔网络和相关生物学网络的算法,计算生物学杂志,7,331,(2000)。
5.M.J.Lde Boon,C聚类文库手册,(2002)。
6.M.J.L.de Hoon,S.Imoto和S.Miyano,生物信息学。
7.M.B.Eisen,P.T.Spellman,P.O.Brown和D.Botstein,基因组范围的表达图谱的聚类分析和展示,美国国家科学院院报,95,14863,(1998)。
8.M.B.Eisen,http//rana.1bl.gov/EisenSoftwareSource.htm(2000)。
9.N.Friedman,M.Linial,I.Nachman和D.Pe′er,用贝叶斯网络分析表达数据,计算生物学杂志,7,601,(2000)。
10.基因组对象网,http//www.genomicobject.net/public/index.html(2002)。
11,A.J.Hartemink,D.K.Gifford,T.S.Jaakkola和R.A.Young,结合测定和表达数据来原则地发现遗传调控网络模型,太平洋生物计算研讨会记录,7,437,(2002)。
12.T.R.Hughes等,通过表达图谱概略的功能发现,细胞,102,109,(2000)。
13.S.Imoto,T.Goto和S.Miyano,用贝叶斯网络和非参数回归估计基因间的基因网络和功能结构关系,太平洋生物计算研讨会记录,7,175-186,(2002)。
14.S.Imoto,S.Kim,T.Goto,S.Aburatani,K.Tashiro,S-Kuhara和S.Miyano,基因网络非线性建模的贝叶斯网络和非参数异方差回归,IEEE计算机协会生物信息学研讨会会议记录,(2002)。
15.F.V.Jensen,贝叶斯网络导论,伦敦大学出版社,(1996)。
16.Y.Maki,D.Tominaga,M.Okamoto,S.Watanabe和Y.Eguchi,推导大规模基因网络的系统的开发,太平洋生物计算研讨会记录,6,446,(2001)。
17.M.J.Marton等,用表达图谱的DNA微阵列评价药物靶点和确定次级药物靶向作用,Nat.Med,4,1293,(1998)。
18.H.Matsuno,A.Doi,M.Nagasaki和S.Miyano,基因调控网络的杂合Petri网表示,太平洋生物计算研讨会记录,5,338-349,(2000)。
19.H.Matsuno,A.Doi,Y Hirata和S.Miyano.生物路径XML资料及其在基因组目标网中的模拟.基因组信息学,12,54-62,(2001)。
20.D.Pe′er,A.Regev,G.Elidan和N.Friedman,从被干扰的表达图谱中推导子网络,生物性信息学,17,副刊.1,215,(ISMB 2001)。
21.M.A.Savageau,生物化学系统分析分子生物学中功能和设计的研究,Addison-Wesley,Reading,(1976)。
22.I.Shmulevich,E.R.Dougherty,S.Him和W.Zhang,或然布尔网络基因调控网络的规则基础的不确定模型,生物信息学,18,261,(2002)。
23.P.Tamayo,D.Slonim,J.Mesirov,Q.Zhu,S.Kitareewan,E.Dmitrovsky,E.S.Lander和T.R.Gloub,用自组图解释基因表达图谱方法和在造血细胞分化中的应用,美国国家科学院院报,96,2907,(1999)。
实施例6用于发现网络关系的系统及其应用1)目标在众多领域中极其期望确定网络关系,这些领域包括药物发现、基因组学、蛋白质组学、计算网络、人类网络和其中元素件以复杂方式相相互作用用的其它系统。
药物开发中的要点之一是探索先导化合物的作用位点。药物筛选包括过程1)合成先导化合物并筛选,2)寻找作用位点和指导疾病处理,3)改进为更有效的药物和4)改进为副作用更少的药物。在药物开发中至今尚无重要技术是利用基因组的分析辅助系统。虽然人类基因组已经被阐明,但是目前还难以利用人类基因组来构建一个辅助系统。潜在的困难基于以下事实,即利用基因敲入/敲除来构建模型十分困难而且由于每个器官有其特定表达方式,应当针对各个器官建立模型。
在这些情形下,把酵母用作真核微生物来确定药物靶点,并构建提炼有关靶向位点的网络的特征的系统被认为是用于最有效地开发具有较低副作用的药物的有效系统,所述的酵母容易被建模并且对其已经积累了大量数据。实际上,通过各个基因的敲入/敲除作用积累了基因间的表达差异,并且从这些数据建立了表达控制网络。然后,先导化合物被加到运用该网络的酵母的培养液中,并且通过比较基因表达与已积累的数据之间的不同以估计先导化合物的任何作用位点。GNI已经开发了一种基因组分析辅助系统,提供分析环境和分析服务,帮助制药公司更加精确和有效地开发药物。
2)构建网络的技术基因组项目被认为是20世纪最伟大的分子生物学项目之一,而许多动物模型的基因组结构已经被阐明是基因组项目成就之一,并且对甚至人类基因组结构的分析也在进行中。随着基因组结构分析的迅速发展,阐明基因组功能将成为可能,而这不可能通过传统的生物化学和遗传学方法来实现。最近的结果包括用于研究所有基因的表达的转录组和用于研究所有蛋白质表达的蛋白质组。除了这些结果,通过分析这些数据研究构成该机体的所有基因的表达网络,能够从实践上阐明整个细胞功能。基于这些事实,确定基因组结构的下一步工作是阐明各个机体包含的基因组功能。根据各个基因组本身特性进行分类和系统化已经开始了对控制机体所有功能的研究。因而,一个新时代一开始,在这个时代观察到的数据最终将被组合并与基因组功能阐明相联系。在这些情形下,几乎不可能通过分析单个基因的传统方法来收集和分析多样的信息。
需要新的分析方法来实现至今所描述的目标。与通过试验确定各个关系而构建的传统模型不同,应当考虑从这些大规模数据完整地构建模型。例如,可以通过分析在各种培养条件下机体所有基因的表达图谱和时程,或者通过分析基因被破坏的突变体的表达图谱来构建基因表达控制网络。如图1所示的研究可能是从系统的角度考虑,并且从输入输出数据和干扰来假定系统内在结构的研究已经成为可能。当将干扰视为对基因的破坏(基因不表达)时,从系统的角度考虑,问题是基因表达系统能否被推测,假定输入指所有基因的正常表达,输出指当部分基因被破坏时所有基因的表达。
从系统的角度来确定系统当基因组结构被阐明时,有可能收集足以使网络被假定的试验数据。
—可以测定所有相关基因的数据。
—很多工具使许多试验被实施—例如已经开发了芯片当通过应用干扰来获得许多归一化数据来测定输出值,那么
图1.分子网络使表达图谱能够被分析的起点是微阵列1,2)或者基因芯片3)方法的出现,所述的微阵列能够同时测定多种DNA或RNA的变化(表达图谱)。“自然遗传学,1月,1999”出版了关于基因芯片的专刊4)。现在这种方法开始被广泛用于生物学和医学领域,并且世界上已经建立许多企业,开始销售整个系统,包括芯片和盘(见,如www.inmcite.com和www.affimetrix.com)。在日本,DNA芯片实验室公司、Takara酿酒公司和Hokkaido系统科学公司也已经涉足芯片生产。由于期望从这些方法获得大量试验结果,需要新的假说和分析方法。这些方法将使新模型能够替代通过试验确定各个关系而构建的传统模型,而新模型可一举从大规模数据构建。例如,通过分析在各种环境下的全部基因的表达图谱可能构建基因表达的控制网络。综合方法搜索能够探测到导致疾病的基因并建立疾病模型。
分析基因表达网络(1)基本阵列技术在阵列试验中基本上采用与传统Southern印迹方法相同的原理。促进阵列技术迅速发展的技术进步包括测定仪器的开发,采用这些仪器通过在具有良好光学特性的玻璃片上结合高浓度DNA来精确测定样品和对照之间的数量比。最近,采用美国斯坦福大学的布朗博士等人开发仪器(阵列仪(arrayer)或打点器)能够成功地将10,000或更多cDNA点在玻璃片上。通过标记了两种以上荧光物质的样品和参考DNA与玻璃片上的cDNA杂交,可以测定基因表达。在不同方法中,用工业上更可获得的照相平版印刷技术已经生产了基因芯片,这种方法在玻璃的表面上直接合成高浓度寡核苷酸。这种方法现在在实际中使用3)。此外,正在开发一种利用喷墨机制的方法,这种方法有望能够较高速较高密度地形成阵列。
3)微阵列的应用下文将描述应用试验的例子。在本发明人实验室中,通过破坏或者过度表达芽殖酵母中的一些基因来综合收集基因的直接和间接作用,以研究所有基因的表达图谱。微阵列试验过程将以图2和3在本文描述。
(I)首先制备将被研究的微生物的所有基因。一般地,制备一组所有基因的PCR(聚合酶链反应)引物(长20个碱基),对基因组或者克隆进行PCR反应来获得所有基因的扩增片段。
(II)用称为打点器或者阵列仪的机器人将各个待扩增的基因的DNA点阵在玻璃片上,将基因固定。每个点的直径为150μ,每平方厘米上点阵2,500个基因。
(III)之后,制备样品。在芽殖酵母的情况下,在相同条件下培养基因组分析后选择的菌株和特定基因被破坏的突变菌株。从这两种菌株中提取mRNAs,用两类荧光物质(两者之间具有不同激发波长和荧光波长)同时标记。混合等量的从各个菌株制备的cDNA。将混合物加到微阵列上,cDNA通过互补键与固定在玻璃表面的探针结合。清洗玻璃片以除去未特异性地结合的cDNAs,从而在一个试验中获得总共约6,000个基因的表达图谱。
(IV)用双-波长荧光计来计数微阵列结合的样品,计算样品与参考基因的比率。由于点阵的尺寸很小,要求荧光计具有高灵敏度。
图2制备微矩阵的简略图
图3.用微阵列测定表达图谱图中的突变细胞是指用化学药品处理的细胞。基因表达图谱是为了图谱间的比较而准备,接着确定这些图谱在网络中的位置。
3.用信息技术分析基因表达网络(1)数据类型在选定时间获得数据,或者在培养过程中的特定时间间隔对细胞进行取样得到时程此时的表达模型数据。此外,在不同试验条件下,在时间轴的一个点上进行试验,尽可能多地收集表达图谱变化。表1给出部分基因表达矩阵,其中汇集了基因破坏试验的数据。在表格的列中以与野生型菌株的比率显示6,000个基因的基因破坏试验结果,其中收集了各个被破坏的菌株的数据。

表1基因表达矩阵各列表示被破坏的基因的名称,各行指基因的名称。各图以正常表达水平为1描绘相对表达量。
(2)分析模型从阵列数据确定网络按照根据分子生物学的中心规则预测的微生物的功能表达是由基因表达信息确定,用于阐明表达网络已经被研究的微生物的功能的信息有望通过分析大规模表达数据得到。一般而言,所用的方法包括从采用了聚类方法6, 7,8)的方案产生的基因表达图谱来阐明具有相同基因表达图谱的基因组群。除了这些方法之外,全世界都已经开始推测基因表达控制网络的研究。其基本目标在于开发计算方法来从基因表达矩阵或时程数据产生如图4所示基因表达控制网络的构建草图。
分析大规模网络的导向图 图4.基因网络鉴定的概述图采用这种模型的网络分析大体上分为两组。一组采用布尔网络,基因间关系用双元素状态表示,其中一个基因被或不被另一个基因控制9,10)。不同模型涉及实际数量的表达,其中大部分模型建立在于连续值相关的动力学的不同方程上11,12。但是,由于在针对整个基因组进行的表达控制网络分析中应当优化模型中包含的反应恒量13),建立在各种方程上的模型需要一些装置。除了这些方法之外,还引入了采用贝叶斯网络的分析方法14)。
2)本发明试验(1)表达图谱构建程序由于当前已经可能在环境中进行基因组水平研究,这使得多维数据容易被测定,为了阐明基因表达控制网络,有必要收集大量基因被破坏和基因过度表达的菌株的表达图谱。本发明人还构建了许多基因破坏的基因突变体,这些突变体主要含有表2所示的与转录相关的基因、与信号传导相关的基因和与细胞周期相关的基因,并试图收集其表达图谱。至今已经构建了主要涉及转移因子的173个突变体,并且今年将通过揭示表达图谱数据的项目构建它们的表达图谱。
表2.基因破坏菌株的构建方案

(2)鉴定网络的试验下面将介绍本发明尝试的多水平绘图方法。通过图5所示的方法确定网络,包括如下步骤。(1)从所收集的基因表达图谱数据构建基因表达矩阵GE;(2)一个基因对其它任何基因的作用如果存在的话,通过对表达比已经显著改变的基因赋值为1和赋予其它基因为0的操作进行探明,从而形成二元素矩阵。(3)探寻矩阵形成环的部分,收集该部分并定义为新的基因组群B’。由于形成环形控制结构的处理很复杂,把这些部分收集在一起再定义为新基因组群。(4)替换矩阵的行和列来排列矩阵右上角的相关基因。这使得相互控制的基因的网络拓扑关系的大部分变得清楚。(5)最后,可以通过除去间接联系部分进一步定义基因表达网络。通过这种方法从本发明人构建的70个基因被破坏的菌株的表达图谱确定的网络示于图6。
多水平绘图方法从表达模式假定表达控制网络 图5.多水平绘图方法的略图 图6.由70个被破坏菌株确定的网络的一部分4.将来的问题最大的问题是包含在数据中的误差。鉴于当前可以获得的试验技术可能引起约30%-40%的误差,有必要开发基于这些误差的信息学(imformatex)技术。但是,目前还不足以获得用于处理这么大误差的方法。因此,期望引入试验性地减少误差的方法或者在一定程度上估计误差的方法。
结果的显示方法的发展对于处理大量基因来说是十分重要的。应当为选择用于分析的转录因子展示6,000个基因之间的网络。本发明人构建初始的552个基因成员的调控控制关系的模型,然后构建含有98个已知转录因子的子网络模型。在构建这些模型中,采用Shoudan等X报道的布尔计算算法。图7显示本发明方法的略图。
多水平绘图方法用基因表达矩阵重构按基因和一些等价设定分类的网络 图7.揭示等价集合和基干网络之间的部分排序,并定义基因之间的布尔函数图7.用布尔方法构建基因调控子网络模型。首先,从一组基因破坏试验建立基因表达矩阵。矩阵的每个元素表示基因表达比,双元素关系矩阵表明如果基因‘a’被删除基因‘b’的强度显著地变化,那么基因‘a’作用于基因‘b’。如果基因‘a’和基因‘b’相相互作用用,形成环(强烈联系的成分)。引入等价集合,逻辑地将所述的环视为一个基因。一个等价集合是相相互作用用的或者仅作用于一个基因的一组基因。接着将基因划分成等价集合。等价集合之间的这些半排序可达性矩阵将作用描述成一个等价集合组,即使具有该设定的各个基因并不直接地相相互作用用。然后除去较低的行间还包含的关系以建立网络中的分级关系,生成控制关系中的基因的网络拓扑关系。
通过表达试验的布尔模型构建的网络模型产生包含在有丝分裂和减数分裂路径中的非常明确的调控转录因子子网络。图8详细地说明了转录因子子网络。从本图中可以看出,该子网络描述了包含在配对、减数分裂和DNA结构和修复机制中的转录因子的特定调控关系。例如HAP2,其功能尚未被报道,其通过特定功能同样未知的中间产物THI2对已知的配对反应基因Alpha2起间接调控作用。而且,被认为与已知的减数分裂路径无关的各种与染色体结构和DNA修复相关的元素件,受MET28和CIN5作用。这两个结果具有生物学意义;成功的有性生殖和减数分裂与配对间的关系对无误差DNA复制的依赖是显然的。
为了理解转录因子的各种功能组中的广义调控关系,将网络模型中的转录因子根据其“细胞角色”进行分类。本文的“细胞角色”定义为涉及YPD中的蛋白质的主要生物学过程。在本发明的转录因子网络中定义了10个组用于分类。图8a图示本发明网络中被分类的转录因子之间的关系,图8b表示这些组间的大量关系。
从图8a可以看出,细胞周期转录因子仅受10组转录因子中的两组、总共仅3个基因的作用。单独地作用于细胞周期转录因子的这3个基因与减数分裂、配对反应和细胞应激相关,由于功能性原因,这些过程有必要将相互作用限定在细胞周期过程。众所周知,真核生物采用谨慎和高度保守的细胞周期关卡(checkpoint)来确保在DNA复制和修复时抑制核分裂(35-37)。
同样地,显然脂肪酸代谢转录因子受碳水化合物代谢转录因子影响。本发明数据证明,例如,INO4,UDFHGDHG的表达受RTG1,JJHGDGGH的影响。还有,PDR1,BIG HORSE受CAT8,SMALL DOG表达的影响。磷脂合成路径中包含的蛋白质与葡萄糖反应路径、脂质信号路径和其它脂质合成路径的基因之间的其它这种相互作用已经在文献中充分描述(26)。众所周知,编码许多磷脂生物合成酶的基因在其启动子中都含有UASINO的变体,它们在对生长阶段和营养缺陷反应时被调控。(G.M.Carman(1991))。特别地,Snf1p/Snf4p和Irelp被报道为与UASINO序列结合的重要分子,与葡萄糖反应相关,是抑制含有UASINO的基因必需的。
为了证实本发明模型在基因组范围内的有效性,在本发明的552个基因成员的调控模型中搜索是否存在先前在文献中描述的调控关系。在98个已知转录因子中,其中26个因子的调控网络上游和下游可以从文献中已经报道的相互作用中阐明。检验有多少已知基因的这种路径在本发明模型中被重建。结果发现,这些关系中大约27%直观地从本发明表达数据中得出。
采用基于基因破坏体的表达图谱数据以及布尔计算模型的本发明的基因调控网络构建已经阐明了大量基因调控关系,这使得特定预测与路径生物学相关。从新的生物学视角对这些模型的进一步调整和数据与其它如蛋白质表达试验的表达数据的相关性,将发现更加确定的特异性路径相互作用。此外,本发明人目前正在研究其它计算方法如贝叶斯建模,来检验表达模型数据中的随机关系。采用这种以及相似的策略从统计的表达数据中阐明动力学调控路径将开始新的阶段,有成本效益地、高通量地研究生物学路径和机制。
图8为细胞功能的基因调控网络。图8a为按照YPD数据库提供的信息通过细胞功能分组的98个转录因子的调控子网络。圈中的基因是被分入特定细胞功能分类的基因。调控控制中的分级关系的线路用彩线描述,线条具有相同颜色为一类基因,,从这些基因发出控制关系到它们施加调控作用的基因。粗的红色线条表示来源于减数分裂/配对反应组基因的对与细胞周期相关转录因子的调控控制。深蓝色线条表示发自碳水化合物代谢组的施加于脂肪酸代谢组基因的控制关系。这些结果得到文献支持,也是常规知识。这种及其保守基因调控关系保证如果食用过多甜食就会发胖。b.转录因子的功能分类之间的控制关系频率。每行的标题表示起始于该类的基因控制反应,列的标题代表所控制的基因的分类。数字表示从一个分类中一个基因发出的控制关系的数量,并将该控制施加于其它分类的一个基因。粗的红色线条突出了相对少的对细胞周期转录因子施加作用的控制关系。粗的蓝色线条的数量显示相对少的对脂肪酸代谢施加作用的控制关系,这两种控制关系都发自碳水化合物代谢分类。
图9.从对破坏体突变的酵母菌株的基因表达试验的布尔分析构建的转录因子基因调控子网络的详图。黑色线条表示调控关系,箭头表示表达作用的分级方向。节点的颜色和形状表示按照在YPD数据库中描述的基因产物的细胞功能的一般分类。与细胞分裂机制相关的基因用三角形节点表示,与DNA修复和染色体结构相关的基因用正方形描述。所阐明的网络表明与减数分裂、配对反应和DNA结构及修复机制相关的基因之间的新的拓扑控制关系。与减数分裂和配对反应基因相关的基因位于受INO2调控的级联的下游,而与DNA修复和结构相关的更低一级的子组基因分级出现在UME6的下游。
图10描述用于发现本发明网络的系统的示意图。细胞(第一行“A”)生长并进行或者不进行干扰,来产生反映缺少其时处理的mRNA。信使RNA被反转录成cDNA并点阵成阵列(从上而下第二行)。收集关于表达水平的信息并储存在数据库中。按照本发明的网络建模分析来自数据库的信息,生成不同基因间的关系网络(盒子中)。然后来自系统的结果显示在监视器或者其它可视化工具上。
上述确定网络中的方法可以被应用于各种线性和非线性分析方法中,包括但不限于布尔分析、贝叶斯分析、图形建模(graphical modeling)、聚类分析、最大似然估计、最大所罚似然估计和其它类型分析方法。
参考文献1)Schena,M.等,用互补DNA微阵列定量追踪基因表达图谱,科学270,467-470(1995)。
2)DeRisi,J.L.,在基因组水平上搜索对基因表达的代谢和遗传控制,科学278680-686(1997)。
3)Lockhart,D.J.等,通过与高浓度挂核苷酸阵列杂交的表达追踪,自然生物技术14.1675-1680(1996)。
4)芯片预测.自然遗传学增刊.21(1)1999。
5)田代康介,牟田滋,原哲,用酵母基因发现分析DNAマイクロアレイを,细胞工程学,18.1050-1056(1999)。
6)Tamayo P.等,用自组图方法说明基因表达(geno oxpression)图谱方法和应用于造血细胞分化,美国国家科学院院报.96,2907-2912(1999)。
7)Toronen,P.等,用自组图方法分析基因表达数据,FEBSleft.451,142-146,1999。
8)Arkin,A.等,从测量值相关度量构建反应路径的实例,科学277.1257-1279(1997)。
9)Liang S.,等,展示用于推导基因网络体系的通用反向工程算法,太平洋计算生物学研讨会会议记录,18-29(1998)。
10)Akutsu,T.等,根据布尔网络模型从少量基因表达图谱确定基因网络.太平洋计算生物学研讨会会议记录,17-28(1999)。
11)R.Somogyi,等,基因表达矩阵为了提取基因网络结构,第二次非线性分析世界大会(WCNA)(1996)。
12)D.Weaver等,用加权矩阵建立调控网络模型,生物计算太平洋研讨会,112-123(1999)。
13)冈本正宏,S-系统用于推定基因相相互作用用,生物信息学,165-188(2000)。
14)Friedman,N.,等,用贝叶斯网络分析表达数据,RECOMB2000,127-135(2000)。
图8a表3.各个“细胞角色”的关系分析 图8b
图9<基因网络的商业模型>
主要商业目的 目标基因网络模型的例子1)新的药物靶点信息信号传导基因网络2)新的毒副作用信息有丝分裂基因网络3)新的路径/相互作用信息 细胞周期基因网络4)新的组合分子信息人类同种组织基因网络人类组织基因网络 图10
工业实用性用于确定基因之间网络关系的方法用于在制药中开发先导化合物、卫生保健和其它期望得到基因间关系的工业中。本发明的推导方法在任何期望得到各组数据间的复杂关系的统计分析领域中应用。这种领域包括工程学、经济学和生物学。
序列表<110>基因网络公司(Gene Networks,Inc.)<120>应用产生自多重破坏表达文库的基因调控网络的生物发现<130>GENN1000CN0-DBB<150>60/325,016<151>2001-09-26<150>60/334,230<151>2001-11-29<150>60/370,824<151>2002-04-08<150>60/334,372<151>2001-11-29<150>60/334,255<151>2001-11-29<150>60/397,458<151>2002-07-19<160>2<170>PatentIn vetsion 3.2<210>1<211>9<212>DNA<213>酿酒酵母(Saccharomyces cerevisiae)<400>1tagccgcca9<210>2<211>10<212>DNA<213>酿酒酵母(Saccharomyces cerevisiae)<400>2tgggcggcta 10
权利要求
1.一种构建基因网络的方法,包括步骤(a)给机体的一组基因提供定量破坏体数据文库,所述文库包括基于破坏所述组基因中的每个基因的表达结果,量化平均效应和每种破坏对每个其它所述基因的可变性测定结果;(b)从所述文库构建基因表达矩阵;(c)生成所述基因之间的网络关系;和(d)确定是否一组或多组基因表达区别于其它所述组基因的表达。
2.权利要求1的方法,还包括步骤(e)提供一种贝叶斯计算模型,其中所述贝叶斯模型包括最小化BNRC检验数。
3.权利要求2的方法,其中所述最小化BNRC检验数的步骤包括采用非线性曲线拟合方法,该方法选自多项式函数、傅立叶级数、小波函数、回归样条函数或B-样条函数。
4.权利要求1的方法,其中所述数据文库是用药物改变基因表达而建立的。
5.权利要求2的方法,其中所述最小化所述BNRC检验数的步骤还包括用修合算法选择贝叶斯模型。
6.权利要求2的方法,其中所述最小化所述BNRC检验数的步骤还包括使用Akaike信息检验数。
7.权利要求2的方法,其中所述最小化所述BNRC检验数的步骤还包括使用最大似然估计。
8.权利要求1的方法,其中所述基因与细胞周期有关。
9.权利要求2的方法,其中所述可变性测定结果是变量。
10.权利要求3的方法,其中所述非线性曲线拟合方法是一种非参数方法。
11.权利要求10的方法,其中所述用于最小化BNRC检验数的非参数方法包括使用异质误差方差。
12.权利要求11的方法,其中所述最小化BNRC检验数的步骤还包括以下步骤(1)生成分值矩阵,其第(i,j)个元素是结构图基因i→基因j的BNRCjhetero的分值;(2)实施一个或多个增加、去除和反转而产生最小的BNRChetero;和(3)重复步骤2直到BNRChetero不再进一步减小。
13.权利要求11的方法,其中所述最小化BNRC检验数的步骤还包括使用爬山算法来最小化BNRC(j)hetero。
14.权利要求11的方法,其中用自举方法来确定边缘强度。
15.权利要求14的方法,其中所述自举方法包括步骤(1)通过从原始基因文库表达数据中多次随机取样并置换,产生自举基因表达矩阵;(2)评估基因i和基因j的基因网络;(3)重复步骤(1)和(2)T次,从而生成T个基因网络;和(4)用(t1+t2)/T计算基因i和基因j之间的自举边缘强度。
16.一种阐明基因网络的方法,包括步骤(a)给机体众多基因提供时程基因表达数据的原始数据库;(b)从所述原始数据库中减去背景信号强度;(c)对所述众多基因中的每一个计算基因表达的相对变化;(d)用Student t-检验分析所述基因表达的相对变化的统计学显著性;(e)将基因表达中的所述变化拟合成线性样条函数。
17.权利要求16的方法,还包括从考虑中去除其表达水平极低以至于主要由噪声确定的那些基因的步骤。
18.权利要求1的方法,其中所述分组步骤包括将所述基因分成一个或多个等价集合。
19.一种评估基因网络关系和优化所述关系的超级参数的方法,包括步骤(1)确定超级参数pj;(2)初始化γjk=0,k-1,...,qj;(3)重复步骤3-1和步骤3-2寻找最佳βjk;(3-1)计算γjk=(BTjkWjkBjk+nβjkKjk)-1BTjkWjkx(x(j)-∑Bjk′γjk′)k′≠k来确定βjk,(3-2)针对候选的βjk值重复步骤3-1,选择使BNRChetero最小化的最佳βjk值;(4)以k=1,...,qj,1,...,qj,1,...重复步骤3直到满足合适的收敛检验数;和(5)针对候选的pj值重复步骤1到4,并选择使BNRChetero最小化的最佳pj值。
20.一种构建含基因网络的系统的基因网络模型的方法,包括使用贝叶斯计算模型,其中贝叶斯计算模型包括最小化BNRC检验数。
21.权利要求20的方法,其中最小化BNRC检验数包括采用非线性曲线拟合方法,该方法选自多项式函数、傅立叶级数、小波函数、回归样条函数或B-样条函数。
22.权利要求20的方法,其中最小化BNRC检验数包括用修合算法选择贝叶斯模型。
23.权利要求20的方法,其中最小化BNRC检验数包括使用Akaike信息检验数。
24.权利要求20的方法,其中最小化BNRC检验数包括采用最大似然估计。
25.权利要求20的方法,其中最小化BNRC检验数包括采用非线性曲线拟合方法,其中该非线性曲线拟合方法是非参数方法。
26.权利要求25的方法,其中非参数方法包括使用异质误差方差。
27.权利要求26的方法,其中最小化BNRC检验数还包括步骤(1)生成分值矩阵,其第(i,j)个元素是结构图基因i→基因j的BNRCjhetero的分值;(2)实施一个或多个增加、去除和反转而产生最小的BNRChetero;和(3)重复步骤2直到BNRChetero不再进一步减小。
28.权利要求26的方法,其中最小化BNRC检验数还包括使用爬山算法来最小化BNRCjhetero的步骤。
29.权利要求26的方法,其中用自举方法确定边缘强度。
30.权利要求29的方法,其中所述自举方法包括步骤(1)通过从原始基因文库表达数据中多次随机取样和置换,提供自举基因表达矩阵;(2)评估基因i和基因j间的基因网络;(3)重复步骤(1)和(2)T次,从而生成T个基因网络;和(4)用(t1+t2)/T计算基因i和基因j之间的自举边缘强度。
31.权利要求20的方法,其中贝叶斯计算模型被用于分析该系统的基因表达图谱。
32.权利要求31的方法,其中基因表达图谱包括该系统中的各个基因的基因表达水平。
33.权利要求32的方法,其中该系统中至少一个基因被破坏。
34.权利要求32的方法,其中基因表达图谱包括亚基因表达图谱,并且其中亚基因表达图谱包括当该系统中至少一个基因被破坏时该系统中的各个基因的基因表达水平。
35.权利要求34的方法,其中基因表达图谱包括至少两个不同的亚基因表达图谱。
36.权利要求32的方法,其中该系统用一种试剂处理。
37.一种构建含基因网络的系统的基因网络模型的方法,包括使用贝叶斯计算模型和布尔方法。
38.权利要求37的方法,其中贝叶斯计算模型包括最小化BNRC检验数。
39.权利要求37的方法,其中贝叶斯计算模型和布尔方法被用于分析该系统的基因表达图谱。
40.权利要求39的方法,其中基因表达图谱包括该系统中的各个基因的基因表达水平。
41.权利要求40的方法,其中该系统中至少一个基因被破坏。
42.权利要求40的方法,其中基因表达图谱包括亚基因表达图谱,并且其中亚基因表达图谱包括当该系统中至少一个基因被破坏时该系统中的各个基因的基因表达水平。
43.权利要求42的方法,其中基因表达图谱包括至少两个不同的亚基因表达图谱。
44.权利要求40的方法,其中该系统用一种试剂处理。
45.一个数据文件,包含通过权利要求20的方法构建的基因网络模型。
46.计算机可读形式的权利要求45的数据文件。
47.可以远程获取的权利要求45的数据文件。
48可以从因特网网址获得的权利要求45的数据文件。
49.一种在含基因网络的系统中确定一种试剂的靶向基因的方法,包括(a)用贝叶斯计算模型构建第一和第二基因网络模型,其中所述贝叶斯计算模型包括最小化BNRC检验数,其中该第一基因网络模型是通过分析第一基因表达图谱获得和该第二基因网络模型是通过分析第二基因表达图谱获得,其中该第一基因表达图谱从没有用试剂处理的系统中得到和该第二基因表达图谱从用试剂处理的系统中获得,和(b)用所述贝叶斯计算模型分析该第一和第二基因网络模型,其中该试剂被视为该系统中的一个基因,和其中该试剂的靶向基因被确定。
50.权利要求49的方法,该靶向基因是直接受所述试剂影响的基因。
51权利要求49的方法,该靶向基因是间接受所述试剂影响的基因。
52.含有按照权利要求49的方法确定的该试剂的一个或多个靶向基因资料的数据文件。
53.计算机可读形式的权利要求52的数据文件。
54.可以远程获取的权利要求52的数据文件。
55.可以从因特网网址获得的权利要求52的数据文件。
56.一种提供服务的方法,包括获得来自用户的试剂,和按照权利要求49的方法为该用户鉴定该试剂的靶向基因。
57.权利要求56的方法,其中所述获得试剂包括获得该试剂的资料。
58.一种提供服务的方法,包括获得来自用户的试剂,和采用按照权利要求20的方法构建的基因网络模型为该用户鉴定该试剂的靶向基因。
全文摘要
本发明的实施方案包括应用新的推导方法来分析复杂的生物信息,包括基因网络。在一些实施方案中,同时获得生物体中的大量基因的破坏体的数据和/或药物诱导/抑制数据。新方法包括布尔推导方法的改进以及采用这些方法来确定生物体中被表达的基因之间的关系。另外的新方法包括贝叶斯推导方法的改进以及采用那些方法来确定被表达的基因之间的因果关系,以及在一些实施方案中确定被调控基因的上游效应子。贝叶斯方法的其它改进包括使用异质方差和不同的曲线拟合方法,包括样条函数,来提高对被表达基因的网络的结构图的估计。其它实施方案包括使用自举方法和确定边缘效应来更精确地提供被表达基因之间的网络信息。用从先前研究以及新近进行的基因表达研究中获得的信息评价了本发明的方法。
文档编号G06F19/24GK1592852SQ02823466
公开日2005年3月9日 申请日期2002年9月26日 优先权日2001年9月26日
发明者井元清哉, 后藤多嘉绪, 官野悟, 田代康介, 米歇尔·德胡恩, 克里斯托弗·J·萨瓦, 久原哲 申请人:Gni株式会社
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1