本申请要求2015年8月27日提交的美国临时申请第62/210502号的优先权，通过引用将其整体专门并入本文。
本发明涉及用于通过构建基因-基因调控影响网络并且与多组学患者特异性测量结果进行比较来提供患者特异性基因表达预测的数据驱动的整合系统和方法，所述基因-基因调控影响网络包括团体策划的生物通路网络信息(community-curatedbiologicalpathwaynetworkinformation)和组学数据，例如，基于rnaseq的表达数据、拷贝数变异(cnv)数据和dna甲基化数据，所述多组学患者特异性测量结果包括基于rnaseq的基因表达、基于阵列的dna甲基化(表观遗传学)和基于snp阵列的体细胞拷贝数变化(scna)。更具体地，使用患者特异性基因表达预测来识别个体患者样本中基因表达水平与预期水平的显著偏差和不一致性，以提供与癌症和癌症处置相关的预测信息。
背景技术：
癌症的病理学与控制正常细胞生长和分化的自然复杂生物过程中的显著畸变相关联。然而，甚至源自于同一组织类型中的癌症内也存在显著的异质性，可能反映出正常信号传递网络可能被病理改变的多种方式。这种异质性是诊断学和治疗诊断学生物标志物开发以及肿瘤学中潜在的治疗介入所带来的重大挑战的基础，并且指出需要对癌症病因和进展的系统水平的理解。
例如，编码受体酪氨酸激酶的表皮生长因子(egf)受体家族成员并在细胞增殖中起重要作用的erbb2基因在多种癌症中，尤其是在乳腺癌、胃肠癌和卵巢癌中高度过度表达。该基因在大约20％的乳腺癌中失调，并且在大多数情况下，该基因的过度表达与拷贝数扩增相关联，并且得到以该基因开头而命名的乳腺癌的特定亚型(her2阳性乳腺癌)的定义。尽管可以获得针对乳腺癌的该特定亚型的靶向治疗介入(即，赫塞汀)，但是乳腺癌患者对这种治疗的反应率仍保持在50-55％的范围内。反应中的这种异质性表明存在肿瘤进展的其他基因调谐子。事实上，已经显示出akt/pi3k通路中的畸变，例如，肿瘤抑制基因pten的缺失和pik3ca基因中的突变导致对赫赛汀的抗性。然而，当前还没有系统水平的通路模型能够将所有这些因素整合到用于治疗抗性的单个整合生物标志物中。
尽管已知的癌症驱动基因中的特定复发突变的致瘤效应已被很好地表征，但是对于在癌症中观察到的绝大多数复发突变的功能相关性知之甚少。评估突变的功能相关性的计算方法与背景突变过程相比在很大程度上取决于估计它们对蛋白质结构的影响或基于它们发生的相对频率。为了揭示突变对下游细胞过程的潜在影响，最近的方法试图通过将在癌症样本中的多组学测量结果与团体策划的生物通路网络相整合来确定基因组畸变的功能效应。然而，这些方法中的绝大多数都倾向于忽略关键的生物学考虑因素，包括多种调控因素对下游基因转录的不等的和可能的非线性影响以及通路相互作用的组织特异性。
为了评价癌症样本中突变或基因组畸变的功能重要性，已经开发出若干计算框架。尽管基于对蛋白质结构的突变效应的推断的方法已被广泛用于团体，但是最近的工作已经集中于通过评价与背景突变过程相比的基因突变的相对频率来确定基因中的驱动突变。认识到沉默突变对于任何候选基因通常是罕见的，可能导致背景突变率估计结果不准确，因此mutsigcv试图利用与候选基因具有相似基因组属性的基因以改善背景突变率估计结果。其他方法旨在识别在给定的癌症亚型内经常被体细胞突变击中的子网。然而，这些方法并不能提供对体细胞畸变的下游失调或信号传递效应的机制见解。这些缺点已经导致基于网络的方法，其中细胞实体(例如，基因、rna、蛋白质、蛋白质复合物和mirna)之间良好策划的生物相互作用按照通路网络被并入模型。其他研究仅聚焦于癌症临床结果与分子实体的激活水平(例如，基因和蛋白质表达水平)之间的关联，但是没有明确地对癌症生物学中突变的功能效应进行建模。最近，提出了典范转移(paradigm-shift)来将通路调控网络与多组学数据进行整合，以对体细胞突变对通路中各个节点的活性的功能影响进行建模。基于从其上游调控网络曾经获得的对应节点的活性与从其下游目标节点再次获得的对应节点的活性的差异来推断任何给定蛋白质中的体细胞畸变的功能效应。
尽管在开发中有所不同，但是在这些方法中存在共同的缺陷，这就是它们绝对依赖于生物通路网络，因此这些方法的使用应当受限于良好策划的通路网络中，并且并不推荐用于部分验证的网络或从不同的组织背景导出的分子网络。更重要的是，这些技术通常假定所有的亲本基因对对应的相互作用具有同等贡献，因此忽略了网络节点之间的相互作用之间的影响强度的变化的可能性。例如，如果多个基因作为特定目标基因的转录调控因子出现，则认为它们对目标基因的表达水平具有同等贡献，这在生物学上是有问题的。实际上，相邻节点之间的配对影响可能非常不同。hotnet算法考虑了链路之间的异质性，该算法旨在通过基于网络拓扑结构定义基因对之间的配对影响度量来发现这种异质性。然而，并不能从推定的通路网络拓扑结构中完全提取由复杂的潜在调控相互作用引起的实际配对影响异质性。
由于通路水平畸变可能源于多种来源，例如，体细胞突变、拷贝数变化、表观遗传学变异以及调控基因表达变化、因此对这些变异来源联合建模对于开发肿瘤学中使用的综合的基于通路的整合预测模型至关重要。此外，利用分子生物学中低成本全基因组数据采集技术的最新进展，不同来源的变异的测量结果变得越来越可用。然而，研究团体和诊断团体两者都缺乏能够充分利用这些多组学谱中存在的信息的建模框架。因此，开发用于整合各种数据来源(包括rna表达水平、拷贝数变异、dna甲基化模式以及体细胞突变)的计算框架来实现找到在临床上有用的生物标志物的目标是肿瘤学团体中的基本需求。
最近，提出了若干整合方法将各种信息来源并入统一的框架以促进癌症早期诊断、临床结果预测以及更为相关的治疗介入。这些方法中的大多数都采用以下两种极端的观点中的一种：i)完全忽视概念生物信息并纯粹依赖数据驱动技术，或者ii)经由并入相互作用的分子实体的网络来充分信任概念生物信息。由于可能会对数据过度拟合，因此在第一种方法中忽视细胞分子实体(例如，基因和蛋白质)之间的生物相互作用使得在寻找具有显著集体预测能力的生物学相关实体子集方面效率非常低。事实上，这个问题在癌症研究中特别突出，这是因为任何给定的研究中癌症样本的数量趋向于比测量的分子特征的数量低一个数量级。另一方面，充分依赖描述性生物网络忽视了它们的局限性：通路网络通常是基于特定细胞背景中的实验证据来构建的，其可能并不总是能转化到其他组织和病理学背景。
本发明采用混合式方法，并将基于测量的多组学数据与部分可信的通路信息并入到统一的框架中以构建基因-基因影响网络，该网络能够鉴于调控网络状态来预测特定的基因表达水平。该框架不仅细化和扩展了我们对组织特异性蛋白质-蛋白质相互作用的知识，而且还提供了网络实体(例如，基因)的患者特异性预测和条件分布。然后利用这些患者特异性基因表达预测从个体患者样本中的预期水平中找到基因表达水平的显著偏差和不一致性，因此允许发现与表现型(例如，治疗反应和预后)的潜在关联。
本发明克服了将生物信息与各种分子测量数据来源整合到统一的基于网络的计算框架中的若干显著限制。这得到了揭示更为相关的患者特异性功能障碍基因和扰乱的生物过程。
例如，本发明的方法并入了生物信息并且仅报告与潜在的基于网络的预测和患者特异性测量结果显示出显著不一致的基因。因此，这种方法在识别与考虑中的表现型相关联的最为功能相关的基因中得到更高的特异性和灵敏度。
而且，当前的基于集合的方法通过首先注释基于先前的生物学知识与特定表现型或细胞/生物过程共同关联的基因集合来考虑生物信息。然而，基于集合的方法不能自适应整合，并且用户需要经由形成潜在地更为相关的基因集合来手动地包括生物信息。相反，在本发明内不需要关于癌症生物学的任何先前信息。该方法为来自通路网络注释的每个基因开发基因调控网络。得到的与表现型相关联的通路子网络提供了功能性见解以及鲁棒的生物标志物，并且因此能广泛应用于各种癌症。
当前可用的基于网络的方法(例如，典范、病理学家和spia)旨在将通路信息与测量数据进行整合，以便识别显示出与从网络获得的预测有显著偏差的干扰的通路和基因。这些方法有两个重要的缺点。首先，这些方法完全信任生物通路网络关系，而不考虑通路网络连接中潜在的组织特异性变化。第二个也是更重要的问题是，这些技术忽略了网络中相互作用链路之间功能异质性的可能性。它们假定所有直接亲本节点的影响是等同的，但是实际上一些调控亲本基因的影响可能明显高于其他亲本基因。
内部方法和系统并不完全依赖通路网络，而是通过将不同系数分配给从多组学数据学习到的网络边缘来细化影响网络。参见例如表2和表3；表示上游调控因子的网络边缘是使用针对祖先的系数而捕捉的；顺式调控影响被捕捉为cnv和甲基化系数。另外，松散连接的链路被移除。因此，我们的方法突出并发现了网络节点(例如，基因、rna、蛋白质)之间的异质关系。
相比之下，我们的方法使用生物通路和多组学测量数据两者来不仅捕捉拓扑结构，而且还捕捉以上所提及的网络中的节点之间的影响的强度。因此，它提供了网络节点之间的更为准确且现实的影响。其次，内部方法不仅限于找到经常受体细胞突变影响的通路，而且还可以找到功能障碍的节点。
为了解决这些问题，我们称之为受突变影响的信息流(“inflo-mut”)的本发明的过程包含多组学测量结果以构建基因-基因调控影响网络，所述多组学测量结果包括基于rnaseq的基因表达、基于阵列的dna甲基化(表观遗传学)和基于snp阵列的体细胞拷贝数变化(scna)以及生物通路网络信息。inflo-mut从正常组织和癌症组织的分子谱中学习调控节点对目标基因的配对影响。为了推断新样本中节点的活性，inflo-mut使用已经从训练数据集中学习到的网络系数。这是通过学习非线性贝叶斯模型来实现的，以使用其自身的scna和甲基化谱以及从生物通路网络推断的上游调控影响来预测任何给定基因的表达水平。这种方法不仅通过捕捉异质配对影响系数来解决不相等的亲本节点贡献的问题，而且能够学习节点之间的非线性关系。inflo-mut还允许评估体细胞突变与下游目标基因之间的关联，以便揭示突变对目标基因失调有较高影响的基因子集。我们通过将inflo-mut应用于乳腺癌和结肠癌上的两个大型多组学数据集来证明inflo-mut的鲁棒性和生物有效性，并且揭示这些疾病中基因上的突变在关键致癌通路中的潜在调节效应。
技术实现要素：
具体地，本发明的目的是提供一种系统和方法，其通过将策划的通路网络与多组学生物信息和各种分子测量数据来源整合到统一的基于网络的计算框架中以识别体细胞突变的影响来解决现有技术的上述问题。本发明的另一目的是提供一种系统和方法，其用于提供患者特异性基因表达预测并识别患者基因表达水平与预测水平的显著偏差和不一致性，从而识别更为相关的功能障碍基因和受干扰的生物过程。本发明的另外的目的是识别与诸如治疗反应和预后的表现型的潜在关联。本发明的又一目的是提供对现有技术的替代方案。
因此，通过提供用于识别和报告驱动失调基因的潜在体细胞畸变的系统和方法，本发明的第一方面旨在获得上述目的和若干其他目的，这样的方法包括以下步骤：
通过从良好策划的公众可获得的通路网络获得生物网络通路信息并将所述通路信息输入到被配置为接收所述通路信息的处理器上来确定针对每个特定的感兴趣目标基因的上游调控亲本基因信息的主数据集；
通过应用来确定针对每个特定目标基因的调控树，所述调控树捕捉所述基因的表达水平与所述基因自身的基因组和表观遗传学状态之间的关系以及所述基因的上游转录调控因子；感兴趣基因存在于根节点中，并且树的叶子表示潜在地通过中间信号传递伙伴直接或间接调控基因的转录的所有基因；
确定基于测量的组学数据的第二数据集，例如，rnaseq表达数据、拷贝数变化数据以及dna甲基化数据，并且将所述基于测量的组学数据输入到被配置为接收这样的数据的处理器上，
通过计算机应用计算技术，基于基因的表观遗传学信息和调控网络状态来学习针对每个感兴趣基因的非线性函数，以便将该特定基因表达水平与和调控树叶子相关联的测量结果相关；所述非线性函数的参数是使用包括新颖的深度惩罚机制的贝叶斯推断方法来估计的，所述新颖的深度惩罚机制用于捕捉更靠近所述树中的根节点的节点的潜在更强的调控影响。
通过计算机应用分析技术来预测针对每个感兴趣基因的表达水平；
确定与针对期望的目标基因的观察到的表达水平有关的患者特异性信息，并且输入所述患者特异性信息作为第三数据集，所述患者特异性信息包括新的癌症样本数据，例如，rna表达数据、cnv数据、甲基化数据以及体细胞突变数据；
使用所述患者特异性信息和预测表达水平信息来计算针对给定样本中的所述期望的目标基因的所述预测表达水平与和所述观察到的表达水平之间的相对患者特异性不一致性评分；
评价针对所有测试样本获得的激活评分和不一致性评分，以发现所述目标基因表达水平的不一致性与该特定基因的上游调控网络中的体细胞突变之间的在统计学上显著的关联。
根据本发明的第二方面，提供了一种系统，所述系统用于利用个体患者样本中的目标基因表达水平的不一致性与上游调控网络中的体细胞突变之间的在统计学上显著的关联来识别患者特异性生物标志物，这样的系统包括用于识别基因表达水平的显著偏差和不一致性的整合的、统一的网络，包括：
从良好策划的生物网络通路信息获得的针对每个特定的感兴趣目标基因的上游调控亲本基因信息的主数据集，所述主数据集被包含在被配置为接收所述通路信息的处理器上；
针对每个特定目标基因的调控树，所述调控树捕捉所述目标基因的表达水平与所述目标基因自身的基因组和表观遗传学状态之间的关系以及所述目标基因的上游转录调控因子，所述感兴趣基因存在于根节点中，并且所述树的叶子表示潜在地通过中间信号传递伙伴直接或间接调控基因的转录的所有基因，所述树是根据所述主数据集来确定的；
基于测量的组学数据的第二数据集，例如，rnaseq表达数据、拷贝数变化数据以及dna甲基化数据，所述第二数据集也位于被配置为接收这样的数据的处理器上，
根据所述目标基因的表观遗传学信息和调控网络状态确定的针对每个目标基因学习的非线性函数，所述非线性函数将特定目标基因的表达水平与和所述调控树相关联的测量结果相关；其中，所述非线性函数的参数是使用包括新颖的深度惩罚机制的贝叶斯推断方法来估计的，所述新颖的深度惩罚机制用于捕捉更靠近所述树中的根节点的节点的潜在更强的调控影响；
与针对所述目标基因的观察到的表达水平有关的患者特异性信息的第三数据集，所述患者特异性信息包括新的癌症样本数据，例如，rna表达数据、cnv数据、甲基化数据以及体细胞突变数据；
其中，利用所述非线性函数来确定所述目标基因的表达水平，并且确定在给定样本中针对所述目标基因的所述预测表达水平与所述观察到的表达水平之间的相对患者特异性不一致性评分；并且
其中，激活评分和不一致性评分被确定为与针对所述目标基因的观察到的表达水平有关的患者特异性信息的第三数据集，所述患者特异性信息包括新的癌症样本数据，例如，rna表达数据、cnv数据、甲基化数据和体细胞突变数据；
其中，利用所述非线性函数来确定所述目标基因的表达水平，并且确定在给定样本中针对所述目标基因的所述预测表达水平与所述观察到的表达水平之间的相对患者特异性不一致性评分；并且
其中，确定针对所有测试样本的激活评分和不一致性评分，由此识别所述目标基因的表达水平的不一致性与该特定基因的上游调控网络中的体细胞突变之间的在统计学上显著的关联。
附图说明
现在将参考附图来更加详细地描述根据本发明的方法。附图示出了实施本发明的方式，并且不应被解释为限制落入权利要求的范围内的其它可能的实施例。
图1是图示将基因调控和/或信号传递通路网络与基于测量的组学数据整合以提供患者特异性基因表达预测的步骤通路的内部方法的概述。本发明的这个方面的步骤是：i)提取针对每个未隔离的目标基因的调控树，ii)使用训练数据集来学习针对每个目标基因的非线性函数，iii)预测针对感兴趣目标基因的基因表达值并计算激活和一致性评分以及iv)功能突变影响分析；
图2图示了使用从针对样本基因ppp3ca的通路数据库导出的调控相互作用生成的调控树；
图3是针对基因的祖先计数的直方图，示出了针对通路网络中的所有基因的高至水平2的祖先的数量的分布，并且图示了大多数基因在某处具有10至50个上游调控因子；
图4是包括中心s形状和软阈值化以捕捉两个潜在非线性效应的非线性函数的图形：i)接近平均灵敏度和ii)接近平均忽视；x轴指代测得的拷贝数或dna甲基化水平；y轴指代测量对基因表达的影响程度。在接近平均灵敏度的情况下，测得的接近平均结果的dna甲基化的小变化引起基因表达的大偏差。然而，在接近平均忽视中，接近平均值的拷贝数的小变化不会引起基因表达的重大变化；
图5图示了针对crc正常样本和肿瘤样本的jun基因表达水平预测与观察结果的关系。与正常样本(*)相比，癌症样本(*)示出广泛的不一致性。该方法预测是依据后验平均值(o)和置信区间多达由误差棒┬呈现的3个标准偏差来提供；
图6图示了针对brc和crc肿瘤样本的所有基因的不一致性评分；
图7是总结本发明的用于根据基于网络的预测与患者特异性测量结果之间的显著不一致来识别患者特异性功能障碍的基因的方法的流程图；
图8是本发明的图示体细胞突变对结肠癌样本中的目标基因表达的影响的方法的结果的图形表示；
图9是针对基因pten的rna表达的直方图；
图10图示了针对样本基因myb、gata3、pten和erbb2的预测与观察结果的关系；
图11图示了针对基因erbb2的rna表达水平与拷贝数变化cnv的关系；并且
图12图示了pten的上游调控子网络中的体细胞突变对其基因表达不一致性的影响。
具体实施方式
本发明提供了用于将多组学生物信息和各种分子测量数据来源整合到统一的基于网络的计算方法中以用于提供患者特异性基因表达预测并根据预期水平识别基因表达水平的显著偏差和不一致性的系统和方法。下面参考图1-12进一步详细描述本发明。
根据本发明的实施例，由图1中勾画的步骤或模块来阐述呈现该方法的总体框图的流程图，该方法用于提供患者特异性基因表达预测、根据预期水平识别基因表达水平的显著偏差和不一致性并报告患者特异性生物标志物。如图1所示，该方法包括四个主要连续步骤或模块以识别和报告驱动失调基因的潜在体细胞畸变。在第一步骤、模块1中，从通路网络提取针对每个感兴趣基因的调控树，所述调控树捕捉基因的表达水平与基因自身的基因组和表观遗传学状态之间的关系以及基因的上游转录调控因子。感兴趣基因存在于树根节点中，并且树表示基因的转录的上游调控因子的网络。树的叶子表示潜在地通过中间信号传递伙伴直接或间接调控基因的转录的所有基因。我们使用术语“祖先基因”或者简称“祖先”来指代这些基因。
在第二步骤、模块2中，我们确定针对每个基因的非线性函数，以便将该特定基因表达水平与和调控树叶子相关联的测量结果相关。因此，使用每个树子网络来学习非线性函数，以根据其自身的表观遗传学信息(例如，dna甲基化和拷贝数)及其调控祖先基因表达水平来预测对应的基因表达水平。非线性函数的参数是使用包括新颖的深度惩罚机制的贝叶斯推断方法来估计的，所述新颖的深度惩罚机制用于捕捉更靠近树的根节点的节点的潜在更强的调控影响。这提供了函数库，每种功能对应于特定组织类型的背景下的特定基因。该功能数据库被学习一次，并且能够用于由模块3和模块4执行的两个后续步骤中的患者特定分析。
在第三步骤、模块3中计算给定样本中的针对所期望的目标基因的预测表达水平与观察表达水平之间的相对患者特异性不一致性评分。也就是说，模块3接收针对给定患者的信息，并且使用函数库对针对调控网络内的所有基因的基因表达水平执行预测。该模块还通过比较基因表达的实际测量结果或观察值与预测值来计算针对每个基因的一致性评分。在第四步骤、模块4中，评估针对所有测试样本获得的激活和不一致性评分，以发现目标基因表达水平与该特定基因的上游调控网络中的体细胞突变的不一致性之间的在统计学上显著的关联。因此，模块4识别表达水平与从调控网络获得的预测值显著不一致的基因。这些基因可能由于基因中的拷贝数畸变或其祖先中的体细胞突变而变得功能障碍。模块4还提供统计数据以评价可能与子基因表达水平中的不一致性相关联的祖先基因突变的显著性。模块1：并入通路网络——调控树构建
基因转录是复杂的生物过程，其通过多种相互作用的蛋白质和复合物以及dna甲基化的程度和含有(harboring)的dna节段拷贝数在不同水平上进行调控，如生物通路数据库中所注释的。通路网络广泛用于以网络格式呈现细胞内相互作用和基因调控网络。该网络建立节点和边缘的有向图。这些节点可以包括多种多样的实体，例如，基因、蛋白质、rna、mirna、蛋白质复合物、信号受体，以及甚至诸如细胞凋亡、减数分裂、有丝分裂和细胞增殖等抽象过程。网络边缘确定相互作用的节点对并且指定每个相互作用的类型。开发了若干公开可用的通路网络来对各种物种与组织类型之间的细胞内活动进行建模。
在本发明中，我们使用综合网络，其汇集来自各种良好组织的通路来源的通路，所述各种良好策划的通路来源包括nci-pid、biocarta和reactome。这个“超级通路网络”包括六个节点类型，包括：蛋白质或对应的基因、rnas、蛋白质复合物、基因家族、mirna以及抽象事物。这些节点以六种不同的方式相互作用：i)阳性转录、ii)阴性转录、iii)阳性激活、iv)阴性激活、v)基因家族成员，以及vi)成为蛋白质复合物的组分。通常，转录仅终止于由对应的蛋白质表示的基因，而激活适用于所有节点类型。
为了学习将基因的mrna表达水平与基因的表观遗传学参数(dna甲基化和拷贝数变化)相关的函数及基因的调控网络，我们从超级通路网络数据库中提取针对每个基因的调控网络，并且将所述调控网络表示为“树”(图2)。随后，我们提取一列“调控祖先基因”，其被称为调控因子或调控基因，这些基因共同捕捉形成调控树的所有节点的影响。调控因子中的一些是目标基因的直接亲本，因此直接调控其转录，而其他调控因子通过蛋白质复合物和直接调控因子的翻译后修饰而间接影响目标基因表达。
在开发针对每个基因的调控树中，我们从特定的目标基因开始，使用具有一些修改的深度优先遍历算法(例如，公知的深度优先搜索(参见下面的伪代码))在链路的相反方向上遍历上游网络，以收集所有上游节点并捕捉调控基因及其深度(其被定义为到根节点的链路数量，如图2所描绘的)，所述一些修改基于基因转录调控的生物学以及我们有兴趣使用参与调控网络的其他基因的表达来预测目标基因表达的事实。
我们一旦到达预定义的最大深度水平，我们首先终止遍历分支，其中，深度被定义为从访问节点到根节点的链路数量。然后我们消除所有不终止于基因节点的分支；因此，树的叶子总是基因。除了表示概念抽象过程的抽象节点以外，我们还通过了所有节点，以避免不必要的网络复杂化并避免包括不相关的相互作用。当到达基因节点时，我们只通过非“转录”类型的链路，这是因为已经通过考虑这个特定基因的表达水平而考虑了经由“转录”链路终止于基因节点的上游调控网络的部分。针对这个规则的唯一例外是根节点，在该处我们做如下的确切反转：
只有当连接边缘是“转录”类型时，才允许从根节点传递到根邻域的第一环中的直接邻居，从而将亲本限制为影响存在于树根中的基因的表达水平的那些基因。我们还记录从叶子到根节点的距离，这些距离还在函数学习过程中被使用；最后，如果我们经由两条不相交的路径来满足一个节点，则考虑最短路径。用于模块1选择过程的伪代码被总结如下，并且在图2中描绘了从网络提取的针对基因ppp3ca的样本上游树。
图2是使用从针对样本基因ppp3ca的通路数据库导出的调控相互作用生成的调控树的范例。子网络包括具有高至第三水平的深度1的祖先基因。形状定义节点类型：基因(椭圆形)、蛋白质复合物(矩形)、基因家族(五角形)、抽象概念(菱形)。根据其调控功能对边缘进行着色：阳性激活(黄色)、阴性激活(红色)、阳性转录(绿色)、阴性转录(蓝色)、蛋白质复合物组分(黑色)以及基因家族成员(灰色)。根节点的表观遗传学测量结果和scna测量结果(圆角矩形)(被视为额外的调控亲本)由绿色箭头连接。收集高至第三水平(d最大＝3)的调控因子。根节点ppp3ca的第一水平祖先(直接亲本)被示为经由调控基因表达水平的“转录”边缘进行连接。例如，复合物cam/ca++经由激活链路被连接到根节点，并且因此不调控基因表达水平。因此，经由图2的左侧中的复合物cam/ca++连接的所有基因被排除出最终的祖先列表。在通过其他基因时，只允许非转录的链路。例如，myb的上游子网络被限制到非转录的节点，例如，pias3基因和map3k7基因，其影响尚未经由myb表达水平来进行捕捉。通过基因myb的表达水平来隐含地考虑基因gata3和e2f1的影响。
作为范例，在图3中，以对数标度呈现在遍历根节点上游的高至7个链路时的祖先数量的经验分布。大量的基因是上游隔离的孤儿基因(orhpangene)。仅有839个基因具有祖先，这些祖先从针对23个基因的仅一个祖先到针对基因cdkn1a的1152个祖先。在通路网络中并未呈现具有零祖先的基因。
模块2：学习针对每个基因的非线性函数
本发明的方法的第二步骤是学习将存在于根节点处的基因的表达水平与基因的调控网络及其自身的表观遗传学信息(例如，dna甲基化和cnv)相关的函数。“学习”函数意味着量化调控基因的表达水平对目标基因的表达的影响。而且，内部方法训练针对目标基因的模型，该模型基于如在训练数据中观察到的它们的配对影响为亲本基因分配不同的系数(如在下面的贝叶斯模型估计中所描述的，具体为估计βg的方法)。由于多个dna甲基化探针可以与基因的编码区域或调控区域重叠，因此本发明通过包括若干代表性统计数据(例如，最小值、最大值和加权平均值)来利用甲基化测量结果，其中，为了更加准确，在计算加权平均值时，我们排除了具有小于10个探针的区域。因此，如果基因g区域与
其中，i(.)是恒等函数。
为了包括拷贝数变化，本发明使用节段平均值，所述节段平均值被提供用于含有特定基因的区域。大多数基因落入单个cnv节段中。否则，如果基因落在两个节段的边界上，我们就简单地取两个阶段的测量值的平均值。
为了学习针对每个基因的函数，模块2使用其祖先的mrna表达、体细胞拷贝数变化和针对ng个样本的dna甲基化测量结果，以形成下面的经典回归模型：
其中，yg是针对所有ng样本中的基因g的表达水平的n×1向量。
项目
这里的目标是找到经由使均方误差(mse)最小化而提供最佳预测能力的最优模型参数βi，i＝1，2，……，p。人们可以在学习阶段使用正常样本，以避免由于严重扰乱的相互作用而导致高度污染/紊乱的癌细胞模型崩溃。然而，当预测因子的数量很大或者与样本数量相当时(n<o(p))，这可能导致较差的预测能力。在大多数研究中，剖析的癌症样本的数量趋向于显著高于正常样本的数量。例如，在针对乳腺癌的tcga数据的情况下，癌症样本的数量超过正常样本的10倍。因此，排除所有癌症样本是效率地下的。另一方面，由于上述基因组事件，在训练集中包括癌症样本可能会使针对显著偏离某些样本中真正的潜在生物学功能的特定基因的模型性能恶化。因此，我们将所有正常样本和部分未受体细胞突变影响的癌症样本包括在该特定基因及其祖先中，以便学习预测功能。这种方法使得针对每个基因的训练集大小不同，但在模型预测能力方面提供了相当大的改进。
当没有关于模型参数βi的先验信息可用时，最小均方误差(lse)解决方案使针对训练集的均方误差最小化。
当有关于模型参数的先验信息时，lse解决方案不是最优的。这里，有关于能够用于增强模型准确度的模型的先验知识。首先，可能并非所有的祖先基因都可以对给定基因的表达水平产生实质性影响。因此，大量的模型参数βi可以缩小到零。因此，通过避免噪声过度拟合，施加稀疏度增强了模型泛化属性。尽管通过使用通路网络并且仅包括祖先基因而不是使用所有基因作为输入数据已经考虑了部分稀疏度，但是当祖先基因的数量增加(数十倍和数百倍)时，预计稀疏度水平会更高。
施加稀疏度的常用的基于优化的解决方案之一是规范模型参数的范数。惩罚能够应用于系数向量β＝[β1,β2,…,βp]t的lp(p≥0)范数，其被称为桥回归。这种方法的重要特例是分别针对l、l2、l0范数惩罚的套索(lasso)、岭(ridge)和子集选择。在弹性网(elasticnet)中，惩罚项是l1和l2惩罚的线性组合；
其中，λ1和λ2是用于施加稀疏度和泛化的收缩参数。基于凸优化、基追踪、lars、坐标下降、dantzig选择器、正交匹配追踪和近似消息传递的高效算法可以用于解决该问题。然而，这些方法的最大限制缺陷在于只能提供针对回归系数的点估计结果。
相反，本发明采用贝叶斯框架，其通过后验分布提供关于模型参数的更加详细的信息以用于随后的一致性核查分析。除了稀疏度以外，它还允许结合其他的先验知识，如以下所解释的。
历史上，在分析基因表达研究中，生物测量结果之间的潜在的非线性关系被忽略。为了捕捉这种非线性关系，本发明的模块2使用中心s形函数
在通过向上遍历通路网络而开发针对每个基因的祖先集合中，另一个重要的生物学考虑因素是叶子节点到根节点的距离的变化。人们可以预期，更为接近的祖先比经由中间节点的长链连接的更远的节点对后代下游基因表达水平做出更多的贡献。因此，更为接近的节点倾向于在回归模型中产生更高的系数。模块2通过贝叶斯框架中的深度惩罚机制将这一事实用于该方法中，其在下面描述的贝叶斯模型中由
这里，本发明使用贝叶斯框架经由基因自身的表观遗传学数据的非线性变换/投影以及基因调控祖先基因的表达水平来预测基因表达水平。贝叶斯框架经由模型参数的完全后验分布提供期望的统计数据(例如，中值、平均值、时刻和……)。此外，我们使用分层贝叶斯模型来并入关于模型参数的先验知识。得到的后验分布提供了畸变在通路中的功能效应的重要见解。
本发明基于祖先基因与表达正被预测的基因的距离(即，调控网络中从叶子到根的链路数量)来使用具有惩罚的全局收缩和局部收缩的想法。构建了以下模型，其中，为便于标记，省略了下标g：
以上公式扩展了正常的伽马先验构建，以便将链路深度信息并入伽玛先验构建。经由包括在模型参数的方差中的系数k来利用该信息。因此，经由设定
上述分层模型产生下面的完全联合分布：
其立即使用这样的事实提供了下面的后验分布：即，针对每个参数的完全条件后验分布仅仅是包括该变量的项与其它项的乘积，作为归一化常数，以确保所得到的概率积分到一。这种方法被称为项完成：
当n＜p时，woodbury矩阵求逆公式用于计算a-1，以获得更为稳定的结果，并且通过将p×p矩形矩阵求逆转换为n×n矩形矩阵求逆而节省了计算。我们应用了吉布斯采样器，其中进行了老化迭代1000次和计算迭代5000次，以获得模型参数βi,σ的近似后验分布。使用所有样本s∈s对所有基因g∈重复该过程，其中，g和s分别是基因id和样本id的集合。
模块3：预测针对新样本的基因水平表达并且报告针对所有基因的激活和一致性水平
为了评价目标基因g对任何给定样本的破坏，我们获得激活评分ag(新)和不一致性评分cg(新)，其中，第一项示出基因表达水平，其可能与其调控网络一致，并且第二项示出指向基因失调的与预期值的偏差(其可能与体细胞突变相关联)。
使用来自正常同生群(cohort)和癌症同生群的训练样本来执行模块2以函数库的形式提供了结果，其中，每种功能对应于特定基因。然后在模块3中使用该函数库来分析测试样本以识别潜在的不一致性。因此，该模块对所有基因执行基因表达水平预测。针对每个基因，我们提取祖先基因的表达水平以及针对所有样本的自我表观遗传学信息。然后，我们使用为该基因学习的对应函数来预测针对所有样本的该特定基因的表达水平。预测过程为该基因的表达水平提供了条件后验分布。我们使用最大后验(map)方法来获得预期的基因表达水平。
为了计算针对学习其函数的未隔离目标基因的一致性评分，我们注意到针对每个新测试样本y新的任何基因的rna表达的预测分布是通过将模型参数从针对给定输入x新(自我表观遗传学信息和祖先表达水平)的条件后验分布边缘化而获得的：
f(y新|x新)＝∫f(y新|x新,β,σ2)f(β,σ2|y,x)dβdσ2
尽管为条件分布的第一项可以采用封闭形式，但是为模型参数的后验分布的第二项不可以采用封闭形式。这个分布能够用下面的表达式来近似，其中，模型参数(β(i),σ2(i))的实现是使用吉布斯采样方法获得的。
上述分布是高斯混合模型(gmm)，它具有均值(ψ(x新)tβ(i))和方差(σ2(i))的大量均等概率分量。如果吉布斯采样器收敛，则利用协方差矩阵
其中，||.||2是矩阵诱导范数。基于这种分布，我们如下计算针对观察值的z-评分或等价的可能性：
此外，由于针对每个基因的潜在生物过程的复杂性以及不同水平的继承随机性、自然规律性和未知因素的影响，对于每个基因，学习的函数的预测能力可能显著不同。因此，我们将针对正常样本的每个基因的平均经验可预测性作为针对一致性核查的基础水平。因此，只有具有远低于正常样本的平均不一致性的一致性水平的癌症样本才被报告为不一致的样本。使用以下归一化可能性：
其中，n0和n1是正常样本和癌症样本的数量，并且α是0到1之间的调谐参数，以便推动对正常同生群和癌症同生群的不同强调。选取针对α的较低值，以便更多地强调正常癌症并补偿较低数量的正常样本。在本发明中，我们任意设定
除了一致性评分以外，使用以正态分布建模的基因表达水平分布来获得每个基因的激活评分；
其中，μ和σ是在迭代地排除异常值之后针对每个基因表达水平学习的正态分布的均值和标准偏差。后记g被省略以方便标记。类似的归一化用于激活评分。
如以上所讨论的，应用该模块是要基于目标基因表观遗传学以及在使用的调控树中起到转录调控作用的基因的表达水平来使用调控网络顶部的训练模型以预测针对给定样本的期望的目标基因表达水平。在图5中，图示性范例被示为预测包括从tcga结肠癌数据集导出的42个正常样本和42个肿瘤样本的测试样本中的基因jun表达水平。使用模块1和模块2，使用具有5折叠交叉验证的338个正常样本和368个癌症样本来训练该模型。如使用模块1所导出的，基因jun在所采用的通路网络中具有51个高至水平2的上游调控因子。在图5中，针对正常样本和肿瘤样本两者均示出了预测值以及后验均值附近的标准偏差，其是通过在模块3中采用在模块2中学习的模型所获得的。在该图中示出的置信区间的呈现是本发明的方法与点估计方法相比在预测基因表达水平方面的优点，点估计方法仅获得预测值并且没有提供关于预测置信度的统计数据。第二观察结果是基因jun在正常样本中受到严密调控，这是因为基于它的调控因子的表达水平的预测值对于正常样本比对于癌症样本更为准确。事实上，与具有相似偏差水平的14个肿瘤样本相比，仅5个正常样本经历jun表达水平偏离预测值超过3个标准偏差。
为了进一步图示在该模块中建立的基因表达水平与体细胞突变的不一致性之间的关联，图6提供了针对调控网络可用的所有基因上的brca和crc两者的全局统计分析。在这点上，对于每个基因，肿瘤样本被分成两个子集：i)其中感兴趣基因或其第一水平调控因子和第二水平调控因子中的一些产生突变；以及ii)所有监控因子都是野生型。然后，我们取针对突变子集和非突变子集两者的绝对不一致性水平的平均值(图6a、图6c)。针对两个子集的不一致性评分的直方图(图6b和图6d)揭示了在两种癌症中针对突变子集的不一致性评分显著高于非突变子集的不一致性评分。
在图6a和图6c中，每个茎对应于特定基因，其中，红色茎是针对在该目标基因或其(高至水平2的)调控网络中具有突变的样本的平均绝对不一致性，而绿色茎是所有样本上的平均绝对一致性评分的否定结果，其中感兴趣基因及其亲本是野生型的。用于具有野生型调控基因的样本的绿色茎为了便于呈现而垂直翻转。这些基因是基于它们在野生型样本中的平均不一致性水平进行分类的。图6b和图6d是针对平均不一致性评分获得的直方图。顶部行和底部行分别针对乳腺癌和结肠直肠癌。结果显示，在目标基因或其在调控网络中的亲近亲本含有体细胞突变的样本上平均不一致性具有更高水平。
模块4：体细胞突变与不一致之间的关联
基因表达水平可能由于调控网络中存在体细胞突变而偏离预测值，引起调控功能的丧失/获得。也就是说，调控基因中的任一个中的突变都可能影响其在调控基因表达中的适当作用，并对目标基因表达产生偏差。内部方法的模块4提供了评估调控基因中的体细胞突变对针对下游目标基因的不一致性评分的影响的方法。因此，本模块采用由模块3提供的激活和一致性评分，并且针对每个新的测试样本，识别显著不一致的基因并检查它们是否潜在地由当前基因或其调控子网中的cnv畸变或体细胞突变所驱动。
首先，识别由cnv畸变事件驱动的不一致性。如果不一致性是由于基因的过度表达以及基因经历拷贝数扩增(cnv>0.5)引起的，则报告cnv扩增是引起不一致性的主要原因。同样，如果拷贝数缺失(cnv<-0.5)与基因的表达降低(downexpression)相关联，则cnv缺失被认为是不一致性的驱动因素。
对于没有经历相关拷贝数畸变的基因，这种不一致性可能是由影响下游基因的转录的基因的上游调控网络中的突变引起的。调控基因越靠近下游目标基因，预期对下游基因表达水平不一致性的影响越大。因此，模块4分配全局深度惩罚参数0＜α≤1，使得具有到根节点g的di,g跳跃的突变基因i的影响按值
为了量化调控树中的突变的影响，我们针对按其绝对不一致性水平和深度惩罚因子进行缩放的癌症样本中的每个，对影响目标基因或其调控因子的所有非沉默突变进行计数。一般而言，基因h突变对基因g的表达的功能影响(由fg(h)指代)被计算如下：
其中，pg是基因g的调控祖先基因的集合(即，对应的调控树的叶子)，m(j)是在样本j中突变的基因的集合，
图7中的流程图总结了该方法中对每个样本不一致性的解读。对所有样本重复该流程并基于它们被分配的体细胞突变影响谱
图8是以图形形式图示的在模块4中生成的结果的范例。具体地，图8a显示了apc中的体细胞突变对被识别有结肠癌的基因的wnt通路目标基因表达的相对影响。标绘的是目标基因激活和不一致性与影响结肠癌样本中的apc的突变的关联的显著性的-log10(p值)。以绿色突出显示的基因受到显著影响(fdr≤15％)。在图8中，显示了pten的上游调控子网络中的体细胞突变对其基因表达不一致性的影响。深度惩罚参数被设定为α＝1/2。示出了pten的亲本中的体细胞突变的组合对其调控的调控效应，其中基因集合{pten，dyrk2，e4f1和atf2}中的突变示出与pten的表达降低的显著关联。因此，这些基因调控pten中的体细胞突变的影响。因此，dyrk2、e4f1和atf2中的突变联合影响pten的表达，因此这些突变的组合提供了pten在肿瘤中更为准确的功能状态。鉴于pten的破坏导致akt通路的致癌激活，这些基因中的突变是用于选择治疗的预后和/或生物标志物。
范例
为了说明本发明的方法的预测能力，将其性能与包括套索(lasso)、岭(ridge)和弹性网(elastic-net)回归在内的若干接近最优的现有技术的点估计器进行比较。
为了证明本发明的方法的准确度，在迭代地排除显著的异常值之后，我们首先经由最大似然方法来学习针对每个基因表达水平的高斯分布。我们首先通过在每次迭代中学习针对样本的高斯分布，然后移除不在均值的第二标准偏差附近的样本。在随后的迭代中，我们针对剩余的样本重复该过程，直到算法收敛并不再存在异常值。在图9中呈现了针对样本基因pten的经验分布和学习的正态分布。为了比较的目的，我们还学习了student-t分布。student-t分布具有对异常值鲁棒的优点，并且在排除异常值之后非常接近正态分布，如图9所示。
接下来，我们基于预定义阈值将基因表达水平分为三个状态(中性、过度表达和表达不足)。阈值是任意设定的，使得表达降低、中性和过表达状态的概率分别变为10％、80％和10％。模块3提供了针对所有839个未隔离基因的患者特异性基因表达预测。经由对所有基因和患者的状态变化事件求平均来计算状态变化率。针对每个同生群单独计算结果。如果针对样本i和基因g的观察表达状态和预测表达状态分别为
在表1中，计算了针对一些重要基因的预测误差，所述重要基因与癌症高度关联并且在全局通路网络中有一组有效的上游调控基因。可以看出，内部方法优于现有技术的稀疏度施加回归模型，并具有提供针对基因表达水平的完全后验分布的额外优势。
表1：针对内部方法与基于基准优化的稀疏度回归模型的基因状态预测误差率的比较结果。针对所有方法的模型训练和测试都是相同的。针对正常样本和癌症样本的预测准确度被单独呈现。
另一个重要的观察结果是：尽管癌症样本对模型训练的贡献较高，但是由于相对于正常样本，癌症样本的数量较大，因此正常同生群呈现更好的可预测性。这一观察结果适用于所有模型，并且揭示了正常组织中基因表达的功能状态与上游调控网络更为一致。
在图10中也观察到这样的事实：与癌症样本相比，正常样本中目标基因表达水平的预测值与观察值之间的一致性更高，其中呈现了针对样本基因myb、gata3、pten和erbb2的观察值和预测值。这里，正常样本中的基因表达水平与根据基因自我表观遗传学数据以及基因的上游转录调控网络获得的预测更为一致。该图示出了对可以源自于不同来源的癌症样本的不一致性分析的重要性，并且揭示了在仅分析基因表达水平的方法方面关于通路扰乱和基因失调的额外信息。不一致性可能因各种来源引起，例如，目标基因中的拷贝数扩增和缺失以及调控网络中的突变，其破坏调控网络作用的正常行为并因此影响存在于调控网络的根中的目标基因的表达水平。
为了更深入地了解模型系数，针对两个基因erbb2和gata3获得的模型参数被呈现在表2和表3中。每行呈现了通过不同学习方法获得的并且用于内部非线性贝叶斯方法的对应系数值。针对后验分布的标准偏差也被呈现在最后一列的括号中。结果表明，erbb2的表达水平高度依赖于影响其基因座的拷贝数畸变事件，如所提出的非线性软阈值函数的模型参数中所见。这种非线性反映出模型忽视可能是测量噪声的零附近的小扰动。因此，能够在模型中直接使用从snp阵列导出的与拷贝数相关联的对数比率值，而无需将对数比率值离散为放大/中性/删除状态。所有学习方法都有兴趣采用非线性函数的相关性。图11验证了这种相关性，其中，针对基因erbb2描绘了观察到的rna以及预测的rna与cnv之间的关系。在图11中，蓝色点和红色点对应于从模型中获得的观察值和预测值。黑色曲线是由表2中的模型参数获得的非线性rnacnv关系。
该图表明，具有从学习过程中获得的系数的非线性cnv项很好地定义了针对erbb2的rna表达水平，其中由于其他项(例如，dna甲基化和祖先基因表达水平)而具有一些小的变异性。实际上，通过套索方法和弹性网方法将dna甲基化和大部分祖先的系数从预测因子列表中明确地移除，并且值得注意的是，内部发明为dna甲基化分配了可忽略的系数。
表2：针对两个基因的模型系数：erbb2
另一方面，针对gata3的rna表达水平更加受到dna甲基化以及上游调控网络的影响。针对dna甲基化系数的预期负号能够提示针对两个基因的基因表达水平与dna甲基化之间呈负相关关系。最后，对于gata3，上游调控网络在调控该基因的表达中起关键作用，表明该基因在乳腺癌中表达的大部分变异主要是由于转录因子的活性所致。由用于由表2和表3中提供的两个基因erbb2和gata3的方法估计的回归系数揭示出，由于基因调控功能的高度异质性，回归系数对于基因可能显著不同。
表3：针对基因gata3的回归系数
不一致的一个重要来源是由于目标基因的上游调控网络的突变。需要注意，在目标基因表达水平的预测值与观察值不一致的情况下，调控基因的表达水平的影响已经被该方法捕捉到，那么我们推断调控网络不能适当地发挥其调控作用。调控网络的这种功能障碍很可能是由调控网络中的体细胞突变引起的，所述体细胞突变阻止了体细胞的基因或产物蛋白质适当执行它们的功能(复合物形成、基因转录、蛋白质活化和……)，这继而影响到下游目标基因表达水平。
作为说明性范例，体细胞突变对基因pten失调的功能影响被描绘在图12中，揭示了pten表达的不一致性与tp53、pten、pik3ca、map3k1和map2k4中的突变高度相关联。鉴于pik3ca比tp53更频繁地发生突变(分别为387个样本对333个样本)，tp53突变比pik3ca突变产生更高影响是特别有意义的。我们观察到map3k1突变和map2k4突变(其先前被示为与luminal型乳腺癌相关联)影响pten失活，因此在驱动乳腺癌的关键亚型中提供了对这些基因之间的有趣联系。我们还计算蛋白质截断和其他非同义突变对针对pten的不一致性评分的相对影响。该模型确定这两种突变在它们影响pten的任何调控基因时具有相似的影响，而pten中的蛋白质截断突变对其失调的影响更高，与ptenmrna的无意义介导降解一致。深度惩罚参数被设定为α＝1/2。