用于计算系统、生物学和药物靶标发现的预测网络建模的系统和方法

政府许可权利

本发明是在美国国立卫生研究院所授予的授权号NIA 1 RF1 AG057457-01和NIH/NIDDK P30DK116074下由政府支持作出的。政府对本发明拥有特定权利。

相关申请的交叉引用

本申请要求2018年2月27日提交的美国临时申请号62/635,946的权益，该申请的全部内容通过引用的方式被合并在本文中。

技术领域

本发明涉及与预测网络建模或者自顶向下和自底向上预测网络建模有关的计算机实施的系统和方法。更具体来说，本发明涉及用于计算机网络、生物科学和药物靶标发现中的应用的预测网络建模。

背景技术

推导出变量之间的因果性的问题，特别是从观测数据恢复因果网络，是一项特别具有挑战性的任务。贝叶斯网络是旨在从观测恢复一个变量集合之间的条件独立性的一种现有技术方法。但是众所周知的是，贝叶斯网络无法区分等价类中的因果结构，这是因为这些结构由具有相同的联合概率和条件独立性的多个成对变量编码链构成。

观测数据中的变量与系统噪声之间的非线性性质提供了推导出因果性方向的一个重要线索。因果性的方向是作为响应于数据中的波动产生最小残差的假设而推导出的。

发明内容

因此，本发明的一个目的是公开一种整体框架，从而使用一种集成的基于定性知识的贝叶斯推导(QKBI)(自底向上)方法将贝叶斯网络学习(自顶向下)算法与近来开发的因果性推导集成在一起。在之前的工作中已经表明，在给出生物网络拓扑的正确因果性配置的情况下，QKBI可以产生关于网络中的变量的边缘概率的准确预测。因此，自底向上QKBI方法将通过预测当前因果结构上的边缘概率并且将其与观测数据进行比较来直接评估关于因果性的假设。由于预测边缘概率对于网络中的因果性方向是敏感的，因此所提出的方法增强了贝叶斯网络学习，从而不仅推导出条件独立性而且还推导出因果性。

所提出的方法利用变换信号空间中的逐段带约束线性回归模型作为针对完全非参数模型的良好近似，后者也被称作Bernard非参数模型。本文中所公开的系统和方法的其中一个更加突出的优点在于，与复杂的非参数模型不同，带约束线性回归提供了模型中的因果假设的直观并且有效的表示，并且使得我们能够导出关于因果性模型拟合的残差的闭型表达式。第二，对于连续信号空间内的因果性推导问题，之前导出的条件概率分布上的约束与贝叶斯框架中的自底向上概率推导是调和的。这就允许我们把自底向上推导方法与现有技术贝叶斯网络结构学习(自顶向下)集成在一起，并且把成对等价结构类中的因果性推导嵌入到因果网络结构学习问题中。第三，与其中导出单个评估值并且将其用来估计因果性假设的先前所知的优化技术不同，本文中所公开的集成方法利用不把任何单个“最佳”拟合视为理所当然的完全贝叶斯概念来推导因果性的方向和函数。从带约束子空间内的一个模型集合计算拟合残差。与针对观测数据的其他拟合方法不同，本文中的系统和方法依赖于模型参数上的一个约束集合通过从模型参数空间内的受限子区段采样每一个可能的参数来实施拟合，并且在每一个模型样本中生成预测(边缘概率推导)。通过来自所有带约束模型的平均预测与观测数据之间的距离来测量最终拟合性。通过采用约束使得所述系统和方法不再严重依赖于观测数据，并且特别当数据稀疏时防止针对数据的过拟合。

胰岛素抗性(IR)是发展出2型糖尿病(T2D)的先兆，并且会增加心血管疾病的风险。虽然全基因组关联研究(GWAS)已经揭示了与T2D相关联的新的基因位点(loci)，但是其对于解释导致胰岛素敏感性降低的机制的贡献是非常有限的。因此，必须有新的方法来探索胰岛素抗性的基因架构。为此目的，跨越人体内的胰岛素敏感性光谱生成一个iPSC库。对于从100个个体导出的310个诱导性多能干细胞(iPSC)克隆体的基于RNA定序的分析允许我们识别出差异表达基因以及胰岛素抗性和敏感iPSC线之间的路径。对于共表达架构的分析揭示了几个胰岛素敏感性相关基因子网络，并且预测网络建模识别出调控这些共表达模块的关键驱动因素基因的集合。人类脂肪细胞和骨骼肌细胞(SKMC)中的功能验证证实了关键驱动因素候选基因对于胰岛素响应性的相关性。

附图说明

通过结合附图参考后面的详细描述将会更好地理解并且很容易获得关于本发明及其许多伴随优点的更加全面的认识，其中：

图1是本发明的一个示例性实施例，其中示出了在本发明的实现方式中所使用的硬件；

图2A示出了一个示例性贝叶斯网络，其中每一个变量在给定其亲代的情况下与其非后代条件独立；来自相同的等价结构类的三个因果结构变得不可由贝叶斯网络区分。

图2B示出了作为不可区分等价结构的结果的典型贝叶斯网络的示例性部分有向图；

图2C示例性地示出了可以通过预测网络学习方法进一步改进贝叶斯网络结构，这是因为在继续进行预测网络学习之后，来自贝叶斯网络学习的网络分数得到进一步改进；

图2D示例性地示出了使用本发明的算法在精确率和召回率方面对于标准贝叶斯技术的改进；

图3是说明本发明如何预测地达到因果推导的示例性逻辑流程图；

图4示出了与阿尔茨海默病的药物靶标有关的算法的示例性应用中的数据分析管线及之后的软件；

图5示出了来自与对应于RGS4的Abeta40和Abeta42结果有关的算法的示例性应用的结果；

图6示出了来自与对应于PDHB的Aβ结果有关的算法的示例性应用的结果；

图7示出了详述综合网络建模分析管线和功能分子验证的各个单独步骤的流程图：样本收集，数据生成，数据规范化，差异表达分析，共表达网络，预测网络和关键驱动因素分析，KD的优先级排序，以及分子和功能验证；

图8A示出了对于全部样本(AS)的差异表达分析，其中示出了具有多次测试校正p值<0.5的1338个差异表达基因。色标表明规范化残差表达水平(蓝色：低表达，橙色：高表达)，并且各列是根据样本/捐赠者来安排的(紫色：胰岛素抗性，青绿色：胰岛素敏感)；

图8B示出了对于患者平均(ApP)的差异表达分析：图中示出了前500个差异表达基因识别的ApP。色标表明规范化残差表达水平(蓝色：低表达，橙色：高表达)，并且各列是根据样本/捐赠者来安排的(紫色：胰岛素抗性，青绿色：胰岛素敏感)；

图9A示出了AS中的胰岛素抗性(IR)的拓扑重叠矩阵(TOM)。拓扑重叠矩阵(TOM)表明基因与其在AS和ApP中的IR和胰岛素敏感(IS)共表达网络的相应的通路富集分析(PEA)之间的共表达，并且每一个TOM中的色条表明不同的共表达模块。图中仅示出了包括在共表达模块中的基因。包括在灰色共表达模块(其中包含不可被指派到任何其他共表达模块的基因)中的基因未被示出。

图9B示出了AS中的IS的TOM。TOM表明基因与其在AS和ApP中的IR和IS共表达网络的相应PEA之间的共表达，并且每一个TOM中的色条表明不同的共表达模块。图中仅示出了包括在共表达模块中的基因。包括在灰色共表达模块(其中包含不可被指派到任何其他共表达模块的基因)中的基因未被示出。

图9C示出了ApP中的IR的TOM。TOM表明基因与其在AS和ApP中的IR和IS共表达网络的相应PEA之间的共表达，并且每一个TOM中的色条表明不同的共表达模块。图中仅示出了包括在共表达模块中的基因。包括在灰色共表达模块(其中包含不可被指派到任何其他共表达模块的基因)中的基因未被示出。

图9D示出了ApP中的IS的TOM。TOM表明基因与其在AS和ApP中的IR和IS共表达网络的相应PEA之间的共表达，并且每一个TOM中的色条表明不同的共表达模块。图中仅示出了包括在共表达模块中的基因。包括在灰色共表达模块(其中包含不可被指派到任何其他共表达模块的基因)中的基因未被示出。

图9E示出了对应于图9A中示出的共表达网络的相应PEA结果；

图9F示出了对应于图9B中示出的共表达网络的相应PEA结果；

图9G示出了对应于图9C中示出的共表达网络的相应PEA结果；

图9H示出了对应于图9D中示出的共表达网络的相应PEA结果；

图10A示出了一个预测因果网络(PCN)，其中从左到右对于ApP中的IS、ApP中的IR、AS中的IS和AS中的IR的PCN阐明了关键驱动因素和受测关键驱动因素的子网络；

图10B示出了一个预测因果网络(PCN)，其中对于在所有网络中都被复制的关键驱动因素阐明了关键驱动因素和受测关键驱动因素的子网络；

图10C示出了一个预测因果网络(PCN)，其中对应于图10A中的PCN的顺序从左到右对于关键驱动因素的子网络阐明了关键驱动因素和受测关键驱动因素的子网络；

图11A对于IR特定(442个基因)、IS特定(367个基因)以及由于HMGCR抑制所导致的IR/IS重叠(343个基因)DE基因作为维恩图示出了IR对比IS iPSC中的HMGCR抑制之后的差异通路富集；

图11B对于图11A中所定义的各个群组示出了通路富集。对于每一个群组示出了前10个重要通路。

图12A在阿托伐他汀(Atorvastain)治疗iPSC线中通过RNA定序数据示出了预测网络验证(AS网络)，其中随着与HMGCR的距离(层数)增加，DE基因的百分比减小；

图12B在阿托伐他汀治疗iPSC线中通过RNA定序数据示出了预测网络验证(AS网络)，其中在位于HMGCR下游的基因当中，随着与HMGCR的距离增加，倍数变化减小；

图12C在阿托伐他汀治疗iPSC线中通过RNA定序数据示出了预测网络验证(AS网络)，其中来自阿托伐他汀实验的差异表达分析的p值沿着从HMGCR向下的阶梯减小。上方线道是AS IR网络，下方线道是AS IS网络；

图13A在ApP网络中的IR和IS网络全部二者中通过位于HMGCR下游的基因当中的RNA定序数据示出了预测网络验证；

图13B在ApP网络中的IR和IS网络全部二者中通过位于HMGCR下游的基因当中的RNA定序数据示出了另一项预测网络验证；

图14A示出了关键驱动因素基因的功能验证，其中示出了成熟人类脂肪细胞中的胰岛素介导葡萄糖摄取试验。关于无胰岛素(-)示出了倍数变化值。每一组研究对象示出了不同浓度的抑制剂对于HMGCR(阿托伐他汀)、FDPS(阿仑膦酸钠)和SQLE(特比萘芬)的效果；

图14B示出了关键驱动因素基因的功能验证，其中示出了成脂分化试验。通过500nM下的油红O(Oil-O-Red)发射的吸收测量来测量不同浓度的抑制剂对于SGBS脂肪生成的效果；

图14C示出了关键驱动因素基因的功能验证，其中在人类SGBS前脂肪细胞中实施生长试验。在抑制剂存在/不存在的情况下的12天试验之后实施结晶紫染色；

图14D示出了关键驱动因素基因的功能验证，其中在成熟人类SKMC中实施胰岛素介导葡萄糖摄取试验；以及

图14E示出了关键驱动因素基因的功能验证，其中示出了人类SKMC中的生长试验。结果代表均值±标准偏差，并且通过与没有抑制剂的胰岛素条件(图14A和图14D)或者基础分化(图14C)相比较的单向ANOVA*p<0.05、**p<0.01、***p<0.001来评估统计意义。

具体实施方式

在描述附图中示出的本发明的优选实施例时，为了清楚起见将采用特定术语。但是本发明不意图受限于所选择的特定术语，并且应当理解的是，每一个特定术语包括按照类似方式操作来实现类似目的的所有技术等效方案。出于说明的目的将描述本发明的几个优选实施例，但是应当理解的是，可以通过未在附图中具体示出的其他形式来具体实现本发明。

胰岛素抗性对于发展出代谢综合症和2型糖尿病(T2D)是必要的，并且对于心血管疾病是一个重大风险因素，从而结合起来构成现代的大流行病。虽然全基因组关联研究(GWAS)已经识别出与T2D有关性状相关联的大量基因位点，但是这些信号当中的大多数与胰岛β细胞功能和胰岛素分泌相关联而不是与胰岛素抗性

为了填补这个差距，目标是利用从此外还经历了GWAS基因分型的跨越胰岛素敏感性光谱的个体导出的一个大型诱导性多能干细胞(iPSC)库，这些iPSC线是iPSC转录变异性的完全被表征并且被证明的决定因素。举例来说，已经发现同期群中的个体对变异性贡献当中的最高者对于代谢功能被富集。

这些结果提示我们更加具体地分析与iPSC捐赠者的胰岛素敏感性状态相关联的基因表达模式和网络。对于例如具有多基因易感性的胰岛素抗性之类的复杂情况，集成了例如基因和转录组学数据之类的多尺度-组学数据的系统生物学和网络建模提供了用以在其中解释基因与功能变化或疾病状态之间的关联的有用情境。因此，通过重建分子网络可以导致对于疾病底层的通路的更加系统化和数据驱动的表征，并且从而得到一种用以识别治疗靶标并且对其进行优先级排序的更加全面的方法。共表达和因果/预测网络建模中的近来的进展允许我们采取这样的方法。本文中的系统和方法把来自高度表征研究对象的复杂疾病表型关联到伴随的分子网络，所述分子网络随后可以被用来揭示最终确定临床表型的连贯的功能分子子网络及其关键驱动因素基因。

总而言之，在胰岛素抗性(IR)和胰岛素敏感(IS)iPSC之间实施差异表达分析，建立共表达网络以便把数据系统地组织到对于胰岛素敏感性关联功能所富集的连贯模块中，并且实施关键驱动因素分析以便识别出控制和调控IR和IS网络的关键方面的基因。最后，通过人类脂肪细胞和骨骼肌细胞(SKMC)中的胰岛素响应性关联功能试验，以经验方式验证iPSC中的所构造的网络以及经过优先级排序的关键驱动因素。

图1是预测网络系统的硬件的一个示例性实施例。在示例性系统100中，一个或多个外围设备110通过网络130连接到一台或多台计算机120。外围设备110的实例包括时钟、智能电话、平板设备、智能手表之类的可穿戴设备以及本领域内已知的任何其他联网设备。网络130可以是例如因特网之类的广域网，或者是例如内联网之类的局域网。由于网络130，外围设备110和计算机120的物理位置对本发明的硬件和软件的功能没有影响。除非另行表明，否则设想外围设备110和计算机120可以处在相同或不同的物理位置。系统的硬件之间的通信可以通过多种已知的方式来实现，例如使用调制解调器或以太网适配器之类的网络连接组件。外围设备110和计算机120将都包括或附着到通信装备。通信被设想为通过例如HTTP之类的行业标准协议发生。

每一台计算机120由中央处理单元122、存储介质124、用户输入设备126和显示器128构成。可以使用的计算机的实例有：商业可用个人计算机，开源计算设备(例如Raspberry Pi)，商业可用服务器，以及商业可用便携式设备(例如智能电话、智能手表、平板设备)。在一个实施例中，系统的每一个外围设备110和每一台计算机120可以具有与安装在其上的系统相关的软件。在这样的实施例中，数据可以被本地存储在联网计算机120上，或者替换地被存储在可由任何联网计算机120通过网络130访问的一个或多个远程服务器140上。远程服务器140可以存储可由所公开的本发明使用的科学或其他数据库。在替换实施例中，软件作为应用运行在外围设备110上。

本发明的系统和方法在关于成对因果性的现有工作的基础上建立多维设定中的因果性推导的完全光谱。在网络中，从相关辨识出直接因果性是特别具有挑战性的。学习因果网络的挑战包括两个子任务：i)推导出多个变量之间的条件独立性；以及ii)推导出两个或多个变量之间的直接因果性。在成对因果性推导的情况下，已经有采用非参数模型和信息几何方法的优秀研究。在他们的方法中，通过测试原因的边缘分布与给定原因的响应的条件分布之间的独立性来推导因果性。独立性基于相对熵被定义为信息空间中的正交性。这种算法在X和Y确定性地相关的情况下工作得特别好。但是如果存在具有更接近参考的分布的加性噪声，在f的非线性与噪声的宽度相比较小的情况下，具有参考流形的基于熵的IGCI则会失败。

在给定关于因果性方向和函数的情况下，基于条件概率分布上的一个约束集合，本发明利用了在自底向上贝叶斯推导的框架中加密的直观线性回归处理。所述约束概括了可能的线性回归拟合的全体。以此为基础，可以表明通过带约束贝叶斯置信传播对于X和Y的边缘概率的线性计算构成了针对非线性信息几何度量的一种有效并且直观的近似，以便从复杂非线性函数推导出因果性方向。

贝叶斯网络

贝叶斯网络代表假设是范畴的一个随机变量集合v＝{x

本发明使用以下概念：X

贝叶斯网络结构学习

贝叶斯网络代表假设是范畴的一个随机变量集合v＝{x

在给定观测数据D的情况下，本发明的一个目的是重建贝叶斯网络的结构。将一个分数S(G)指派给图G，以便评估网络G对数据集合D的拟合性。该分数由下面的等式A1的后验概率给出。

其中P(D|G)是边缘数据似然性，P(G)是结构G的先验概率，并且P(D)是归一化常数。可以如下面的等式A2计算边缘数据似然性。

P(D|G)＝∫

其中θ是参数集合，ξ表明先验背景信息。假设数据集合D由N个独立的数据样本d

因此，等式A2可以被写成下面的等式A4。

为了以闭型求解等式A4，作出五个假设。

假设1多项分布：设

假设2参数独立性：在给定网络结构G的情况下，与每一个变量相关联的参数彼此独立，从而使得P(θ|G，ξ)分解成下式

由于变量Xi的亲代的每一个实例是独立的，因此P(θ

其中q

假设3参数模块性：在给定网络结构G1和G2的情况下，如果Xi在G2和G2中具有相同的亲代，则有

假设4狄利克雷先验：在给定网络结构G的情况下，P(θ

其中Γ(*)标示伽马函数，r

N′

其中N′是等价样本尺寸，P(X

假设5完备数据：数据集合是完备的。也就是说，D不包含缺失的值或隐藏变量。根据多项样本假设和完备数据的假设，P(D|θ，G)可以被因式分解成下式：

其中N

由于狄利克雷分布对于此域是共轭的，因此每一个参数的后验保持在共轭族中。积分等于：

其中N

自底向上因果性推导

在X与Y之间是非线性关系的情况下，对于存在加性噪声情况下的因果推导，定义原因X与后果Y之间的关系的函数(即Y＝f(X))中的非线性也可以被考虑。所述非线性提供关于底层因果模型的信息，因此允许识别出真实因果机制的更多方面。之前提供了一种用于函数因果建模的方法，其中集成了贝叶斯网络固有的概率推导能力，以便使用假设的亲代(因果)节点的观测数据来生成关于假设的子代(响应)节点的预测。通过在概率空间中定义评估子代节点的预测分布与其观测值的分布之间的匹配良好程度的距离量度，可以对一个等价类内的不同因果性配置进行评估，以便确定数据所支持的最佳的一个。这种建模方法的一个关键优点在于，其所允许的亲代节点到子代节点的后果传播可以大于来自亲代的路径长度1，从而使得有可能在一个节点链条中或者在更加复杂的网络结构中推导出因果性。

为了描述本发明的方法，过程开始于对与传统贝叶斯网络方法中所定义的任何给定的图结构相对应的一般后验概率进行格式变换。设X代表被表示成网络中的节点的变量的矢量；E是证据；D标示观测数据；G是将要推导的图网络结构；并且θ是模型参数的矢量。由此可以将后验概率写成P(G|D)＝P(D|G)P(G)/P(D)，其中在给定特定的图结构G的情况下，边缘概率P(D|G)可以被表达成参数上的积分：P(D|G)＝∫

其中

导出自底向上数据似然性分数

为了计算和优化数据似然性P(D|G，θ)，置信传播作为一个子例程被合并在自底向上因果推导规程中，以便在对于给定的因果结构/配置(G)为预测因子变量给出观测数据的情况下，预测所有响应变量的边缘概率。在此情况下，通过置信传播来计算给定G和参数θ的情况下的X的边缘概率。后面的公开内容解释了如何示例性地使用贝叶斯置信推导P(X|E，G，θ)来计算等式1中的数据似然性P(D|G，θ)。其中E和D分别代表亲代和子代节点上的观测数据。首先，为了避免混淆，在本文中使用X

P(D|G，θ)＝∫

在给定G和θ的情况下，软证据作为经过调整尺度的观测数据D进入P(X

P(D|G，θ)＝P(D|H(G，θ))＝P(D|P(X

内部概率描述原始连续变量X被映射到的概率空间内的二进制变量X

其中

在此情况下，条件概率分布θ＝{P(B|A)，P(B|A)}是从[0，1]上的均匀分布采样的参数。应当提到的是，这一置信传播可以被视为P(B)＝βP(A)+C形式的线性回归，其中

由于真实值与其置信概率之间的确切函数映射是未知的，因此采用一个非参数量度(即库尔贝克-莱布勒(KL)散度)把真实观测的分布与预测边缘概率的分布进行比较。如果子代节点的预测和观测分布匹配良好，则可以得出结论基于G和θ的预测良好地反映观测数据D，从而得到KL散度的更小值。为了强制KL散度充当真实概率度量，对于定义在D和H上的该函数作出对称性和归一化修改，从而使得k(D；H)＝1-exp[-(KL(D||H)+KL(H||D))/2]。可以通过内核上的任何归一化单调递减函数来定义等式3中的数据似然性函数。对于模型选择，其中用以表示模型的后验分数的S＝-log(K(D，H))与内核值负相关。为了优化模型，该分数被最大化，即等价于最小化KL散度。

在给定特定的因果假设G和θ的情况下，KL散度的计算涉及把实值观测数据D与所推导出的置信概率H进行比较。可以通过任意函数

导出在等式4中给出的置信推导的一种方式是将亲代和子代节点之间的关系表示成三次样条，从而可以良好地近似一般的非线性关系。但是在假设线性关系的情况下，针对样条的一种更加直接明了的替换方案将是把X

为了导出用于以逐段方式把置信推导等式拟合到数据的规程，首先将定义一些术语。设X标示(多个)预测因子/(多个)亲代节点的矢量(在这种成对因果性设定中，对于所考虑的任何给定的马尔可夫等价结构，X仅代表单个变量)，并且设Y代表响应变量。设D∈R

其中对于在[0，1]中均匀地采样的θ有P(θ|G)＝1/M。对应于第k个片元中的第l个节段的预测概率

为了简单起见，在后面的实验部分中，为了获得[0，1]中的L个节段，每一个亲代节点的范围可以被均匀地划分成L个节段。

自顶向下和自底向上预测网络学习算法

在非常高的层级，本发明所使用的预测网络学习算法可以被概括成(1)检查动作后的结构是否等价于动作前；以及(2)如果检查返回为真，则对于动作后和动作前结构计算自底向上分数，并且使用该分数来确定因果性。

自顶向下和自底向上打分

自顶向下和自底向上预测网络算法把所讨论的传统的自顶向下贝叶斯打分与新颖的自底向上预测分数集成在一起，以便推导出因果网络拓扑。具体来说，自底向上分数将被用来推导等价结构之间的因果性，自顶向下分数则被用来推导条件独立性。

为了阐明所述算法，定义下面的概念：对于自底向上预测分数，假设原始观测(连续)数据是D

边缘化数据似然性可以被写成：

P(D

D

由于离散化的和经过调整尺度的数据是相关的，因此引入隐藏变量I以使得联合数据似然性是可分解的。

其中I代表仅取决于当前网络结构(G)的二进制指标。

如果当前网络结构G属于任何等价结构类，则I＝1，也就是P(I＝1|G，θ)＝P(I＝1|G)＝1并且P(I＝0|G，θ)＝P(I＝0|G)＝0；如果当前网络结构G不属于任何等价结构类，则I＝0，也就是P(I＝1|G，θ)＝P(I＝1|G)＝0并且P(I＝0|G，θ)＝P(I＝0|G)＝1。我们进一步定义：

因此，集成数据似然性可以被写成：

其中

因此，等式9中的总体自顶向下和自底向上分数可以被写成：

等式12中的第一项代表贝叶斯网络中的(自顶向下)贝叶斯-狄利克雷分数，并且等式12中的第二项标示如前面展开的(自底向上)预测分数。最终的自顶向下和自底向上分数可以被写成：

其中第一项中的变量如等式0中定义，第二项中的变量如等式7中定义。

自顶向下和自底向上(“TDBU”)方法中的最终后验是原始自顶向下和自底向上分组的组合。与贝叶斯网络中一样，由于可能的结构空间的尺寸，优化该分数函数(等式13)是一个NP难度问题。但是任何启发式搜索方法都可以被应用来优化该后验函数，比如蒙特卡洛-马尔可夫链(“MCMC”)或爬山算法等等。图2C和2D分别示例性地示出了网络打分以及使用本发明的算法在精确率和召回率方面对于标准贝叶斯技术的改进。

实例1：预测网络分析将HSPA2识别为新颖阿尔茨海默病靶标

据信，基因和蛋白质表达中的改变对于晚发性阿尔茨海默病(LOAD)的发展是至关重要的。在本文中，蛋白质在两个分开的系列中被检查并且合并到网络中，并且在两个不同的细胞系中对输出进行评估。所述管线包括以下步骤：(1)预测表达数量性状基因位点(eQTL)；(2)确定差异表达；(3)分析转录和肽关系的网络；以及(4)验证两个分开的细胞系中的效果。我们在两个分开的大脑系列中实施所有分析以便验证效果。所述两个系列在第一集合中包括345个样本(177个对照例，168个病例；年龄范围65-105；58％女性；KRONOSII同期群)，并且在复制集合中包括409个样本(153个对照例，141个病例，115个MCI；年龄范围66-107；63％女性；RUSH同期群)。首要靶标是在两个数据集合中被识别为关键驱动因素的热休克蛋白家族A成员2(HSPA2)。在两个细胞系中对HSPA2进行了验证，其中过表达驱动了APP变异细胞中的ABeta40和ABeta42水平的进一步升高，以及经过修饰的神经胶质瘤细胞系中的Tau蛋白和磷酸化Tau蛋白的显著升高。这一工作进一步证明了研究基因和蛋白质表达中的改变对于理解晚发性疾病是至关重要的，并且进一步将HSPA2指名为LOAD过程的一个具体的关键调控因子。

图3公开了根据本发明的一个实施例的示例性逻辑流程，其中将前面所提到的等式和过程应用于针对阿尔茨海默病(AD)的生物化学药物靶标识别。本发明的软件的操作在系统的任一个计算机120或外围设备110处发生。在步骤302处，算法开始于从外部数据库收集基因学、基因组学和蛋白质组学数据。在步骤304处，软件收集共表达和/或差异表达(DE)基因签名以及/或者其他生物标记物。在步骤306处，软件查询外部文献和通路数据库，从而为预测网络收集附加的数据(可选)。在步骤308处，根据前面的公开内容，软件使用所收集的所有数据来形成播种通路。在步骤310处合并表观遗传数据，比如而不限于ENCODE、RoadMap，从而如步骤312所示产生特定于组织的多尺度网络。

在步骤314处，从特定于组织的多尺度网络创建先验/初始网络。总的来说，先验/初始网络是可选的，不需要先验/初始网络后面的步骤也可工作。在没有初始网络可用或者之前没有构造初始网络的情况下，本发明的算法将在步骤322处把步骤302处的基因组学数据、蛋白质组学数据与基因型数据直接合并。如果构造了初始网络，则使得该先验/初始网络经过任何现有的启发式搜索方法，比如而不限于步骤316处的爬山算法、MCMC、步骤318处的全局MCMC以及步骤320处的基于顺序的MCMC。在启发式搜索的每一个步骤处，本发明的算法以交互方式计算集成的自顶向下分数(步骤328处)和自底向上预测分数(示为步骤332)。基因型数据322、eSNP/pSNP因果先验数据324以及GE/蛋白质组学临床数据326在328处被合并到自顶向下逆向工程分数中，从而得到候选因果多尺度网络330。类似地，自底向上预测分数332被用来计算预测GE和临床简档334。自顶向下分数(328和330)和预测分数(332和334)都被周期性地传递到步骤316处的爬山算法、MCMC、步骤318处的全局MCMC和步骤320处的基于顺序的MCMC，以便随着软件的结构学习的进程对其进行更新。这一学习过程的最终输出是通过优化组合的自顶向下和自底向上分数而得到的集成自顶向下和自底向上预测网络模型。最终预测网络被用来确定AD关键驱动因素分析336和AD因果多尺度网络338，预测GE和临床简档334则被用来计算关于GE和临床性状的计算机(in-silico)预测340。后面将更加详细地描述结果。

样本：KRONOSII是已经被给出(Corneveaux等人，2010年)的数据的一个子集，并且包含来自阿尔茨海默病研究中心资助的美国脑库以及6个欧洲和英国脑库。KRONOSII是一个便利同期群，其在受到LOAD影响的样本中具有低次要病理(即路易体病)和高病理载荷，并且对于对照例具有低病理载荷。第二个集合(RUSH)包括来自由Rush大学医学中心(Chicago，IL)的研究者维持的两个大型前瞻性随访同期群：即Religious Orders Study(宗教团体研究)和Memory and Aging Project(记忆和衰老项目)。RUSH集合是基于流行病学的同期群，并且具有病理学和病理分期的更大混合。对于KRONOSII同期群所收集的所有数据集合，其中有168个受LOAD影响的样本和177个未受影响的样本。对于RUSHII同期群所收集的所有数据集合，其中有141个受LOAD影响的样本和153个未受影响的样本。KRONOSII同期群的平均年龄是81岁，59％的研究对象是女性。RUSH同期群的平均年龄时88岁，63％的研究对象是女性。组织切片是从额叶(82％的样本)和颞叶(18％的样本)皮质区取得的。

数据收集：根据制造商的协议在人类全基因组SNP 6.0阵列(Affymetrix,Inc.，Santa Clara，CA)上分析基因组DNA样本。Birdsuite(Korn等人，2008年)被用来从CEL文件调用SNP基因型。DNA质量控制管线类似于在Andersen等人的著作(Anderson等人，2010年)中所描述的DNA质量控制管线。cRNA被杂交到Illumina HumanRefseq-HT-12 v2Expression BeadChip(Illumina，San Diego，CA)。使用BeadStudio软件提取表达简档，减去背景，并且估算缺失的珠粒类型(bead type)。使用lumi(Du等人，2008年)和limma(Ritchie等人，2015年)实施针对RNA简档的规范化。使用装备有定制电喷雾离子化(ESI)接口的Exactive Orbitrap质谱仪(Thermo Scientific，San Jose，CA)实施MS/MS分析。使用精确质量和时间(AMT)标签方法(Zimmer等人，2006年)对肽实施识别和量化。使用Decon2LS进行峰值挑选并且用于确定同位素分布和电荷状态(Jaitly等人，2009年)。去同位素的谱信息被加载到VIPER中以便找到并且匹配特征到AMT标签数据库中的肽识别(Monroe等人，2007年)。来自所提取出的离子色谱的曲线下方的面积被用作肽丰度(peptide abundance)的度量。

数据分析：在图4中示出了数据分析管线。这是用以揭示LOAD风险靶标并且将其置于上游调控(等位基因信息)和下游输出(转录和肽)的情境中的多遍选择规程。目标是识别出一个高置信度靶标的最小集合以进行验证。在规范化之后在KRONOSII和RUSH中分开实施该管线，以便确保独立的复制。

差异表达(DE)：使用比较LOAD和证实病理的对照例的limma(Ritchie等人，2015年)实施DE分析。每一个数据集合(KRONOSII，RUSH)被独立运行。使用Benjamini-Hochberg校正(5％FDR)实施多次测试调节。结果被用来定义用于下游分析的播种集合。

表达数量性状基因位点(eQTL)：MatrixeQTL(Shabalin，2012年)被用来预测等位基因-转录关系。每一个数据集合(KRONOSII，RUSH)被独立运行。使用Benjamini-Hochberg校正(5％FDR)实施多次测试调节。结果被用来定义用于下游分析的播种集合。

网络分析：除了从文献、通路数据库(mSIGDB，GO)和Roadmap倡议(RoadmapEpigenomics等，2015年)导出的外部数据之外，实施作为来自两个数据集合(KRONOSII，RUSH)的输入基因组、转录组和蛋白质组简档发生的网络分析。目标是产生在LOAD中失调的主要生物学过程的输出列表，以及影响LOAD相关联过程的首要KD的较小列表。KRONOSII和RUSH被作为独立的数据集合对待，并且跨越集合比较效果以便确定复制的靶标。该管线包括以下规程(图4步骤3到6，深橙色正方形)：1、构造COEN以便识别与LOAD病理相关联的共调基因的集合，并且确定在每一个网络模块中被富集的通路(步骤4b)；2、确定与LOAD病理相关联的播种基因集合(步骤4a，模块选择(MS)和模块富集(ME))；3、构建多尺度CPN(步骤5)；以及4、确定调制CPN子网络的状态的KD(步骤6)。

因果预测网络(CPN)：虽然COEN允许基因-蛋白质关系的描述性表征，但是为了将网络数据排序到关系的分级结构中从而允许KD分析，因果关系预测是必要的。COEN仅反映关联关系，BN则推导出代表信息流的方向的有向边。BN分析可以捕获非线性和组合交互。对于标准BN分析的一个限制在于，BN内的子结构有时是矛盾的，从而导致许多具有低置信度的有向边。为了解决这一固有的限制，开发了一种新颖的CPN方法，该方法把自顶向下BN与考虑到已知的因果关系的自底向上因果推导集成在一起，从而打破矛盾因果结构之间的对称性并且从而导致边方向中的更高置信度。

网络构建的复杂度是所考虑的节点数目和样本尺寸的函数。网络构造中的所有肽都被使用；但是鉴于被用来查询基因表达水平的探针的较大数目，在不损失重要LOAD基因和通路信息的情况下，将在CPN重建中使用的探针数目被减少。构建仅有基因的COEN并且识别出对于DE基因被富集的那些模块，并且随后将CPN构造限制到连贯的LOAD聚焦的基因集合的该子集。

搜索聚焦在与LOAD相关联的网络状态的关键驱动因素的识别上，因此仅使用了LOAD数据集合。用于KRONOSII和RUSH LOAD数据集合全部二者的播种基因集合包括对于DE转录靶标被富集的模块；因此针对LOAD病理的相关性的通路被选择。这些集合被扩展来包括多于仅仅DE转录的内容，这是通过包括来自基于文献的脑特定网络的先验。鉴于所测量的肽的适度数目及其把网络组成部分关联到包含DE转录的通路的潜力，所有的肽都被使用在网络模型中。转录数据被精简到最关键的靶标(模块选择)，并且随后通过包括来自专业数据库中的相同通路的附加靶标而被扩展(模块富集)。为了确保健壮(robust)的复制，KRONOSII和RUSH作为分开的集合被管线处理。

关键驱动因素分析(KDA)：在实施CPN分析之后，使用KDA算法对所得到的预测网络模型进行检查。KD是对于其他靶标的调控状态具有显著影响的靶标。对于KRONOSII和RUSH分开预测KD并且确定重叠。对之应用KDA的LOAD相关联的子网络是通过把三个不同的数据集合投影到网络上而生成的。首先，在KRONOSII与RUSH之间重叠的所有DE转录针对仅有转录和转录-肽因果网络全部二者被投影。第二，DE产物在对应于LOAD的子网络中分解，并且来自COEN的控制模块被投影，以便确定控制效果如何作用在LOAD预测结构中以及富集DE数据集合。最后，整个肽集合被投影到转录-肽因果网络上。

数据验证：在八个靶标当中，由于构造尺寸和成本，一个靶标(ST18)未被跟踪。其他七个构造在HEK293和H4线中被测试。一个构造(CP110)由于构造尺寸而给出低转导效率，并且未被跟踪(数据未示出)。在其他靶标当中，三个靶标(HSPA2，GNA12，COMT)在至少一个LOAD同期群中被过表达，并且两个靶标在LOAD中被欠表达(PDHB和RGS4)。这些构造随着LOAD状态被复制，过表达HSPA2、GNA12、COMT，并且击倒PDHB和RGS4。CCT5未被差异表达，但是作为复制的关键驱动因素肽被跟踪。CCT5作为第一遍复制被过表达。目标是在转导之后的多个时间点(48、72和96小时)获得关于Aβ和/或Tau蛋白的改变的一致度量。

对于RGS4，在所测量的所有时间点都存在Abeta40的显著下调，但是对于ABeta42没有效果(图5)。Aβ水平的这一下降在脑组织中看到的效果相反。在脑中具有更少的RGS4是有毒的，因此对于更少的RGS4，Aβ应当增加。或者，RGS4的更低水平可能是反映出晚期保护性补偿机制，在这种情境中所述结果是合理的。对应于PDHB的Aβ结果在大多数情况下无关紧要，仅有一个时间点在ABeta40中表现出差异(图6)。对于Tau蛋白，RGS4没有显著结果，在数据中也没有任何趋势(图6)。对于PDHB，在所测量的三个时间点当中的两个时间点处，在总的Tau蛋白和磷酸化Tau蛋白中没有显著改变。这些效果与脑组织数据一致，由于在LOAD脑中PDHB的表达更少，因此Tau蛋白应当随着击倒而增加。

实例2：胰岛素敏感性测量和IPSC生成

本发明的系统和方法还适用于人类诱导性多能干细胞中的预测网络建模，以便识别出对应于胰岛素响应性的关键驱动因素基因。

研究中的个体具有伴随的全基因组基因分型和胰岛素敏感性的黄金标准测量(也就是从胰岛素抑制测试导出的稳态血糖SSPG)。其他生物计量参数包括年龄、身体质量指数、性别和种族/民族(表1)。从每一个个体生成三到七条iPSC线，在IR和IS细胞之间没有重编程效率的明显差异。之前描述了对应于iPSC生成和质量控制的完整管线。简而言之，生成对应于来自101个个体的317条iPSC线的RNA定序数据，并且在质量控制之后，对来自100个个体的310个样本的RNA定序数据进行分析，其中48个个体是IS(149个样本)，52个个体是IR(161个样本)。基于之前的出版物，用以在IR或IS状态之间作出区分的SSPG截止值被设定在140mg/dl。平均SSPG值对于IS组是84mg/dl，并且对于IR组时210mg/dl。最后，对全部两组中的样本进行年龄和身体质量指数(BMI)匹配以避免可能的偏差(IS对比IR组中的平均年龄分别是57.7岁和59.5岁，并且IS对比IR组中的平均BMI是28.5和30)。

iPSC线已被证明概括了许多孟德尔疾病，包括由于胰岛素受体中的严重变异所导致的胰岛素抗性。但是iPSC概括对于胰岛素敏感性/抗性的常见多基因易感性的基因架构的程度是未知的。由于用以比较IR和IS状态的整体方法是基础的并且被扩展到整个转录网络(图7)差异表达分析被用来印证关于iPSC保持与个体捐赠者的胰岛素敏感性状态相关联的基因架构的组成部分的假设。为此，首先检查IR和IS iPSC之间的100个最显著的DE基因以用于KEGG通路的富集。在前10个最为富集的通路当中有糖代谢、胰岛素信号、糖酵解/糖异生以及II型糖尿病(表2)，这显示iPSC本身至少部分地反映出个体捐赠者的胰岛素敏感性状态和分子调控机制。

图7所示出的流程图示例性地描绘出本发明的一个示例性实施例的逻辑流程。在本发明的某些实施例中，该逻辑流程可以被实施在软件中。在步骤702处，从目标群组抽取血液，并且收集包括年龄、身体质量指数、性别和种族/民族在内的生物计量参数。在步骤704“iPSC重编程”处，即如本文中更加详细地解释的那样对iPSC进行重编陈。在步骤706“RNA定序”处，实施RNA定序。在步骤708“规范化和调节”处，对样本进行规范化和调节。

由于每个患者具有多个克隆体，因此分析的有效样本尺寸是患者的数目，而不是样本的数目。但是由于胰岛素抗性状态是个体水平的特性，因此无法在不移除感兴趣的信号的情况下针对捐赠者进行调节。因此，按照两种方式来分析数据。第一，在对于每一个个体不考虑多个克隆体的情况下使用所有样本(被称作AS分析)。第二，将对应于每一个个体的所有克隆体的残差表达求平均(被称作每患者平均ApP分析)。虽然可以使用一般化估计等式、混合模型或者类似的技术在考虑到来自相同个体的样本的相关的情况下计算残差表达，但是将不会解决对于每个个体具有多个克隆体的数据构造共表达和预测网络的问题。

这一分析包括步骤710“残差表达”，步骤712“IR和IS共表达网络”，步骤714“IR对比IS DE”，步骤716“预测网络”，步骤718“iPSC先验网络”，以及步骤720“KDA”。步骤710涉及对于基因的残差表达的分析。该分析继续到步骤714，其中实施AS和ApP DE，正如在后面的“差异表达分析”章节中所解释的那样，以及/或者步骤712，其中实施共表达分析，正如在后面的“共表达网络分析”章节中所解释的那样。在步骤716“预测网络”处，实施预测网络分析，正如在后面的“预测网络”章节中所解释的那样。在步骤718“iPSC先验网络”中，预测网络分析合并在先验样本集合中实施的先验网络分析。在步骤720“KDA”处，实施关键驱动因素分析(KDA)并且对样本进行排序，正如在后面的“关键驱动因素分析(KDA)和排序”章节中进一步解释的那样。最后，使用该数据，在步骤722“候选靶标”处通过算法识别出靶标基因候选。

在AS分析中识别出1338个基因在IR和IS样本之间差异表达基因(FDR调节p值<0.05)，而在ApP分析中则没有识别出显著差异表达的基因。但是对于来自全部两项分析的500个最为差异表达的基因的测试统计量的排序的比较显示出一致的结果，其中对于AS和ApP DE分析具有非常类似的排序(Spearman相关系数＝0.66，AS和ApP中的中值排序＝308，配对Wilcoxon签名排序测试P＝0.90)(图8A-B)。这一结果显示在ApP分析中缺少统计显著DE基因是由于减小的样本尺寸。因此我们决定使用AS和ApP方法全部二者，从而在所有分析中一致使用相同的残差表达水平。

根据本发明训练了四个共表达网络，如图9A-D中所示：对于AS和ApP调节的表达残差当中的每一个，为IR样本构建一个网络并且为IS样本构建另一个网络。共表达网络识别跨越样本具有高度相关的表达模式的基因群组(模块)，从而表明其涉及在类似的生物学过程中。使用基因集合富集分析(GSEA)，并且来自分子签名数据库(MSigDB)的基因本体(GO、C5生物学过程，v5.1)基因集合被用来测试每一个共表达网络的模块，以用于胰岛素或代谢相关通路中的富集。具体来说，在AS IR网络中识别出5个相关模块(对应于总共1565个基因)，在AS IS网络中识别出3个模块(430个基因)，在ApP IR网络中识别出4个模块(1689个基因)，并且在ApP IS网络中识别出5个模块(2791个基因)。这些模块在图9E-H中示出。

用于每一个预测网络的播种基因是通过扩张来自相应的共表达网络的所选模块中的基因集合而获得的，这是通过把从公共基因和蛋白质交互数据库(ConsensusPathDB(CPDB)和MetaCore(来自Thomson Reuters的v6.24))导出的先验细胞类型特定网络中的在k＝3或更少步内连接到任何所选模块基因的所有基因都包括在内。前面的播种基因选择过程确保影响胰岛素敏感性的基因被包括在内，同时通过排除不相关的基因以训练预测网络来精简特征空间。最终的播种基因集合分别由7250(AS IR)、3797(AS IS)、8183(ApP IR)和9712(ApP IS)个基因构成。这一最终播种基因列表被使用在自顶向下和自底向上预测网络管线中(参见在线方法)以便构建因果网络模型。

在四个预测网络上实施KDA，并且一共识别出AS预测网络中的281个关键驱动因素(KD)和ApP网络中的259个关键驱动因素。全部两个集合有45个共同的关键驱动因素基因。KDA要求指定一个起始基因集合，并且KDA被运行多次，其中对于关联到胰岛素敏感性的通路被富集的每一个共表达模块中的基因以及来自AS和ApP分析的DE基因(对应于AS的5％FDR和对应于ApP的前500个DE基因)被运行一次。这意味着给定的基因可以对于多于一个基因集合(每一个基因集合代表不同的通路)被识别为KD，并且基因越常被识别出来，在表型中蕴涵该基因的证据就越强。识别出在AS和ApP网络中3次或更多次作为KD的9个基因(BNIP3，CARS，IDH1，NDUFB1，HMGCR，HPN，FDPS，SLC27A1，TMEM54)(图10A-10B)。

对于这前9个KD，在四个网络中检车围绕其中每一个KD的局部子网络。具体来说，计算2个分数(表3)：从每一个关键驱动因素到AS网络中的显著差异表达基因和ApP网络中的前500个最为差异表达的基因的逆路径长度的总和(DE邻近度分数)，第二个分数是从每一个KD到其下游的其他KD的逆路径长度与到其上游的其他KD的逆路径长度的总和之间的差(KD优势度分数)。DE邻近度分数越高，从该KD到DE基因的路径越多，并且/或者这些路径越短；KD优势度分数越高，更直接处于该KD下游的其他KD就越多，并且/或者直接处于其上游的其他KD就越少。

最后，实施定向文献搜索，其中把不同的KD与术语胰岛素、葡萄糖和糖尿病相组合。在后面的表3中总结了将KD关联到胰岛素敏感性的相关出版物的结果。简而言之，BNIP3和SLC27A1都与胰岛素抗性表型强烈相关联。溶质载体家族27成员1(SLC27A1，也被称作FATP1)是在饮食诱导肥胖中所涉及的胰岛素敏感脂肪酸转运体并且已被关联到骨骼肌中的IR，并且BCL2/腺病毒E1B 19kDa交互蛋白3(BNIP3)对于脂肪细胞重构期间的线粒体生物能学是至关重要的，其调控肝脏中的线粒体功能和脂类代谢，并且与PPARγ相结合把线粒体融合-分裂平衡耦合到系统性胰岛素敏感性

上面的表3示出了对应于排序靠前的关键驱动因素的概要统计量和参考。网络出现：跨越对应于AS和ApP的IR和IS网络的出现次数。DE邻近度：从每一个关键驱动因素到AS网络中的显著差异表达基因和ApP网络中的前500个最为差异表达的基因的逆路径长度的总和。KD优势度：从每一个KD到其下游的其他KD的逆路径长度和到其上游的其他KD的逆路径长度的总和之间的差。KD：关键驱动因素。

根据本发明，使用来自iPSC细胞线中的HMGCR抑制实验的DE基因签名来验证因果IR/IS网络和关键驱动因素分析，其中所述iPSC细胞线是从IR(n＝6)和IS个体(n＝6)全部二者导出的。对于该实验中的每一条iPSC线(表4)，在存在或不存在阿托伐他汀的情况下生成RNA定序数据，其中阿托伐他汀是在患有高胆固醇血症的患者中广泛使用的HMGCR抑制剂(他汀)。之前的努力通过类似的方法验证了预测网络。

经过阿托伐他汀治疗与未经治疗的样本的比较得到包括3205个DE基因的一个列表，其中表明了对于他汀作用通路和胆固醇生物合成通路以及其他相关通路的最高富集，显示出HMGCR抑制触发了转录补偿响应以便平衡胆固醇通路中间体的减少。

接下来将比较IR对比IS iPSC对于阿托伐他汀治疗的具体响应。在独立地考虑IR或IS群组时，DE基因的数目被减少到785个(IR)和711个(IS)。如图11A中所示，IR与IS样本之间的DE基因仅仅部分地重叠(分别是343/785或343/711)，从而显示出基于捐赠者的IR/IS状态对于阿托伐他汀治疗的差异响应。此外，对于442个IR特定和367个IS特定的DE基因的通路富集分析证明了IR和IS iPSC之间的经过富集的通路中的显著差异(图11B)，这可能与T2D在处于他汀治疗下的胰岛素抗性患者中的不成比例的发生率有关。虽然阿托伐他汀治疗把代谢通路的扰动翻译成可测量的转录改变并且给出了关于HMGCR功能及其与胰岛素敏感性的关联的线索，但是需要新颖的附加分析来验证预测网络。

在本发明的某些实施例中，网络结构的验证优选地涉及多个步骤：(1)计算网络中的HMGCR的每一个下游层中的在HMGCR抑制实验中发生基因表达的显著变动(FDR＜0.05)的基因的百分比。发现对于IR和IS网络全部二者，HMGCR下游的第一层中的超过80％的基因是DE基因，并且这一百分比随着在网络中与HMGCR的距离增大而减小(图12A)；以及(2)把DE基因倍数改变(log FC)(图12B)和相关联的显著性(-logFDR)(图12C)与AS网络拓扑进行比较。在IR和IS网络全部二者中位于HMGCR下游的基因当中，倍数改变及其显著性随着与HMGCR的距离增大而减小(图12B，图12C)。在ApP网络中观察到相同的模式(图13A-13B)。最后，由于HMGCR抑制，所有四个预测网络(AS IR、AS IS、ApP IR和ApP IS)对于DE基因表现出HMGCR下游的显著富集。

因此，结果证实DE基因的百分比、DE基因倍数改变以及相关联的显著性随着与受扰动靶标的距离(离开的层数/步数)增大而减小，并且DE基因与非DE基因相比在HMGCR的下游各步中被富集。合在一起，HMGCR抑制实验验证了预测网络及其拓扑结构。

在本发明的其他实施例中，所述系统和方法在与胰岛素敏感性相关联的过程中对经过优先级排序的KD(HMGCR、FDPS和SQLE)实施验证，特别是与例如脂肪细胞和SKMC之类的相关代谢细胞类型中的胰岛素介导葡萄糖摄取相关联的过程。为此，Simpson-Golabi-Behmel综合症(SGBS)人类前脂肪细胞线和人类SKMC线HMCL-7304分别被分化到最终分化脂肪细胞和肌管。在某些实施例中，验证努力聚焦在HMGCR、FDPS和SQLE上，因为这三个关键驱动因素参与在相同的代谢通路(胆固醇生物合成)中，此外他汀(HMGCR抑制剂)正处在关于其对胰岛素抗性和II型糖尿病风险的效果的激烈讨论的中心。

为了在功能上抑制所述三个候选基因，应用已被明确描述并且广泛使用的化学抑制剂：阿托伐他汀(针对HMGCR)、阿仑膦酸钠(用于FDPS)和特比萘芬(用于SQLE)。如图14A中所示，所有三种抑制剂都减少人类脂肪细胞中的胰岛素介导葡萄糖摄取。但是仅有HMGCR抑制剂表现出前脂肪细胞生长和分化效率的可检测的降低(图14B，图14C)。沿着相同的线，阿托伐他汀同时抑制SKMC中的胰岛素介导葡萄糖摄取和细胞增殖，阿仑膦酸钠和特比萘芬则仅对葡萄糖摄取表现出显著效果(图14D，图14E)。结果显示出与关键驱动因素分析相组合的数据驱动共表达和预测网络是用于发现在IR中所涉及的新颖基因的强大工具。

虽然GWAS研究是针对T2D和胰岛素抗性相关联的血糖性状，但是在识别对于胰岛素抗性风险有贡献的新基因方面的成功是相当有限的。群组之前已证明，基于黄金标准测量的GWAS研究允许发现与胰岛素抗性

考虑到这一目标，利用反映出人类群体中的广谱胰岛素响应的胰岛素敏感性的准确测量生成一个iPSC库。虽然样本尺寸受到限制(对于一共近似300条iPSC线有52个IR对比48个IS个体)，差异基因表达分析突出了针对与胰岛素和葡萄糖代谢直接相关联的分子通路的富集，从而显示出iPSC线确实反映了从中导出所述iPSC线的个体的胰岛素抗性状态。

为了克服样本尺寸限制并且开发出与IR相关联的基因网络的更加整体性的视图，在某些实施例中对于IR和IS iPSC线全部二者分开构建共表达网络。此外，对于网络构造，实施一种双重方法，其中使用对应于所有样本的基因表达值(AS)(用以增强能力)或者每患者的平均值(ApP)(用以增强严密性)。所构造的网络突出了针对与胰岛素敏感性关联过程密切相关联的细胞功能(比如呼吸电子转运链、糖酵解、胆固醇和类固醇生物合成以及葡萄糖代谢)富集的共表达模块，以及从DE分析出现的通路富集签名。此外还实施预测网络和关键驱动因素分析来研究控制前面提到的模块和功能并且从而最有可能涉及在胰岛素抗性的病原学中的中心基因节点。

为了更好地对关键驱动因素列表进行界定和排序，在某些实施例中，只有在AS和ApP方法全部二者中都被视为KD的才被定义(45个基因)，随后考虑在4个所构造的网络中的总的出现次数(AS IR、AS IS、ApP IR和ApP IS)，从而给出前9个关键驱动因素。如表3中所突出的那样，IDH1、BNIP3和SLC27A1(FATP-1)被示出为参与在与胰岛素敏感性相关联的功能中，或者与胰岛素抗性或2型糖尿病直接相关联。在其余的所选关键驱动因素当中，代表给定的KD到DE基因和到其他KD的连接的具有最高DE路径和KD路径值的前两个关键驱动因素是协调参与在胆固醇生物合成通路中的法尼基焦磷酸合酶(FDPS)和3-羟基-3-甲基戊二酰辅酶A还原酶(HMGCR)。此外，参与在该通路中的另一个基因SQLE处于在全部两种方法之间共享的45个KD当中。对于他汀(HMGCR抑制剂)的临床试验的元分析显示出影响胰岛素抗性个体的T2D发生率以不成比例的方式增加。此外，降低LDL-C的HMGCR中的等位基因带来发展出T2D的更高风险，从而导致关于他汀会影响胰岛素敏感性或胰岛素分泌的推测，尽管仍然未能很好地理解导致这样的T2D风险升高的确切的细胞和分子机制。

预测网络不仅示出了相同通路中的共调基因，而且还表明了给定基因的因果性上游和下游。已经有了成功的努力来验证预测网络，因此在本发明的某些实施例中，示出捕获了预测网络中的关键驱动因素的下游效应因子基因是有用的。通过HMGCR抑制的经验网络验证表明对于DE基因的富集和HMGCR的下游附近的对数倍数改变都作用来验证网络的总体结构。

类似于鼠类前脂肪细胞线中的之前的结果，人类前脂肪(SGBS)细胞中的功能验证实验表明HMGCR抑制减少了前脂肪细胞增值以及成熟脂肪细胞中的分化和胰岛素介导葡萄糖摄取。此外，在受到阿托伐他汀治疗的人类SKMC(HMCL-7304)中，增值和胰岛素介导葡萄糖摄取都受到影响。FDPS和SQLE抑制在减少人类脂肪细胞和肌管中的胰岛素介导葡萄糖摄取方面也具有显著效果，同时不会影响增值或分化。这些结果显示出对于脂肪细胞和SKMC的增值和分化的效果与对于胰岛素介导葡萄糖摄取的效果是独立的。此外，与HMGCR抑制相比，SQLE抑制在胰岛素介导葡萄糖摄取上发挥了可比的效果，这显示出底层机制可以通过细胞内或膜结合胆固醇水平的去调控而被介导。

总而言之可以得出结论：(1)iPSC保留捐赠者特定签名；(2)IR与IS iPSC之间的差异基因表达分析显示出对于与胰岛素敏感性相关联的通路的富集；(3)与IS细胞相比，IRiPSC对HMGCR抑制具有差异响应；(4)与关键驱动因素分析相组合的共表达和预测网络揭示出参与在IR中的健壮候选。合在一起，这些结果表明根据本发明的iPSC提供了用于IR和相关联的心血管疾病的研究的一种新颖并且精密的模型，特别当相关的代谢(脂肪细胞或SKMC)和血管(内皮)细胞类型是通过胰岛素敏感性的准确测量从iPSC库生成时尤其是如此。

患者招募包括血液采样，并且通过根据Knowles等人的修饰胰岛素抑制测试来实施胰岛素敏感性测量。

STAR v2.4.0gl被用来把RNA定序读数对准到人类基因组版本GRCh37。使用featureCounts v1.4.4，唯一映射读数重叠基因被计数为由ENSEMBL v70注释。

除非另行表明，否则在该示例性实施例中使用R v.3.0.3来实施统计分析和数据处理步骤。本领域技术人员将很容易理解的是，在本发明的其他实施例中可以用任何类似软件来替代。

在实例2中所描述的所有分析中使用了每百万Log2计数(CPM)和TMM规范化(如在edgeR中实施)。正如本文中所讨论的那样，在至少30％的实验中，仅对于具有超过1的CPM的基因保留规范化计数。作为基因计数的最终规范化，在本发明的某些实施例中使用了来自limma R包的voom函数。

所有RNA定序数据分析都是在对于技术效果(定序批次和RNA准备工具套装、重编程来源细胞)和患者协变量(性别、民族、年龄、BMI)校正的表达残差上实施的。使用variancePartition R库作为随机效应来调节批次和RNA准备工具套装，同时使用limma包作为固定效应来调节重编程来源细胞和患者特性。由于每个患者具有多个克隆体，因此计算两个表达残差集合：一个(AS)使用来自对应于每一个患者的所有克隆体的所有样本，一个(ApP)对于每个患者的残差求平均。

类似于前面所描述的那些处理和分析步骤，根据本发明的eQTL分析是通过对基因型数据进行过滤而实施的，以便移除具有超过5％缺失条目的标记物，低于1％的次要等位基因频率，以及<10

遵循标准实践，只有欧洲血统的个体被包括在eQTL分析中，以避免由于血统与基因表达之间的相关而导致的假阳性。基于全基因组基因型数据的主分量分析对于eQTL识别出81个欧洲血统的个体。eQTL分析是由MatrixEQTL v2.2.1使用前5个基因型主分量作为协变量来实施的。使用PEER v1.0在基因表达数据中识别出隐变量。通过移除前20个PEER分量来计算表达残差。由于从每一个个体实验了多条iPSC线，因此对于给定的个体和基因，对应于每一个个体的表达值被总结为来自所述多条线的平均表达残差值。对应于每一个个体的平均值随后对于每一个基因被分位数规范化。Cis-eQTL分析考虑每一个基因的转录状态位点(transcription state site)的1Mb内的标记物。错误发现率是遵循Benjamini-Hochberg计算的。

对于初始差异表达分析，在对准和特征计数之后，在针对技术(定序批次、RNA提取方法；建模为随机效应)和生物学/人口(重编程来源细胞、性别、民族、年龄、BMI；建模为固定效应)协变量全部二者进行调节之后，对RNA定序计数进行规范化并且计算残差表达。所述数据集合随后被分离成个体的IS和IR群组，并且对于每一个群组中的患者ID使用线性模型分开进行调节，但是随后把来自模型的估计截距加回到残差表达值。随后将两个残差集合相组合，并且实施DE分析。由于截距是估计参数，即使真实的人口值在两个群组中是完全相同的，估计仍将会有略微不同的数值，这是因为数据集合是有限的，这意味着在该DE分析中，所有基因都将是显著的。但是假设是由于胰岛素敏感性状态所导致的平均表达中的差异比由于截距估计规程所导致的随机差异占优势。

对于被用于共表达和预测网络分析的表达残差，采用了2个分析流：一个在不针对患者ID进行调节的情况下使用了所有样本(AS)，一个对每患者残差求平均(ApP)。使用如在来自R语言limma包v3.18.13的lmFit函数中所实施的线性模型确定了差异表达基因。使用FDR调节p值上的0.05截止值对统计显著性进行了评估。

在示例性实施例中，使用coexpp R包构造了共表达网络，其中coexpp R包对于大量基因提供了针对WGCNA R包(v1.14-1，在这里与R v3.0.3一起使用)的优化工作流程。用于预测网络的播种基因(seeding genes)(具体来说是针对pathfinder的输入)被选择成对于胰岛素抗性相关性状(仅有生物学过程)是相关的GO项目被统计富集(FDR<0.05)的共表达网络模块中的基因。

在本发明中开发了自顶向下和自底向上预测网络管线以便构建因果预测网络模型，所述模型利用自底向上置信传播引擎作为子例程来推导因果性。传统的(自顶向下)贝叶斯网络无法捕获相反的因果性，因为这种方法无法在等价因果性之间作出区分。自底向上方法利用生物化学反应的非线性来推导分子交互的因果性，也就是说沿着真因果方向比假因果方向更好地拟合到数据，从而打断统计等价性。通过把新颖的自底向上因果性推导方向与(自顶向下)贝叶斯网络集成在一起，所述集成的自顶向下和自底向上预测网络平台将得到在等价结构之间解析出因果性的完整因果网络。通过集成多尺度“组学”数据(基因型和转录)来构造多尺度网络模型，这一管线集成了BN的优点。基因型数据作为cis-eQTL基因被合并在模型中，并且被约束到顶部节点(没有其他亲代)。我们使用来自共表达网络的所选模块中的基因作为用于预测网络建模的播种基因。

在本发明的该特定实施例中，使用R包KDA的经过修改的版本(https://github.com/zhukuixi/KDA)来实施关键驱动因素分析。第一，通过KDA定义一个背景子网络，这是通过在网络中的靶标基因列表中围绕每一个节点寻找一个K步上游邻域。第二，从该子网络中的每一个节点开始，KDA对于靶标基因列表评估下游邻域的富集(对于从1到K的每一个步长)。使用了K＝6。从AS和ApP网络中取得所识别出的关键驱动因素的重叠，并且进一步排序出前9个KD(在表3中示例性地描述)。

来自6个IS和6个IR个体的iPSC被保持在无饲养层条件下，其中在5％基质胶涂覆的6阱TC板上使用补充有1mM左旋谷酰胺、1mM青链霉素、0.1μg/ml二性霉素B的mTesr1(StemCell technologies,Inc)。为了进行传代，用PBS清洗一次细胞，并且用预热的1mM EDTA(Sigma)进行处理，在37度下培养1-5分钟，并且在新鲜的mTesr1介质2μM Thiazovivin(Millipore)中重新悬浮。在Thiazovivin中培养12小时之后，每天用新鲜mTesr1更换介质。细胞被生长到～90-100％融汇，用PBS清洗一次，并且用DMSO(D)或1μM阿托伐他汀(A)(Selleckhem)处理12小时。在PBS中清洗一次细胞，并且使用PureLInk RNA迷你工具套装(Thermo Fisher Scientific)收获以进行RNA提取。使用Nanodrop(Thermo Scientific)对总RNA进行量化。从进一步的处理中排除具有A260/280比<1.8或>2.3的RNA样本，并且使用Illumina HiSeq 2500系统对RNA进行定序。

Simpson-Golabi-Behmel综合症(SGBS)细胞(人类前脂肪细胞)由MartinWabitsch博士(Ulm University，Ulm，德国)提供。SGBS细胞在补充有10％FBS、33uM生物素和17uM泛酸的DMEM/F12中被培养。SKMC线HMCL-7304细胞由Institute of Child Health(ICH)(University College London)提供。细胞在SKMC生长介质(PromoCell)中被培养。iPSC如前面所描述的那样被生成和培养。

对于成脂分化，SGBS细胞被生长到融汇并且经过一个两步骤分化过程。首先将细胞暴露于介质3天，所述介质由补充有0.01mg/mL的转铁蛋白、20uM的胰岛素、100nM皮质醇、0.2nM的3,3′,5-三碘代-L-甲状腺原氨酸、25nM地塞米松、250uM的3-异丁基-1-甲基黄嘌呤和2uM罗格列酮的DMEM/F12构成。随后，将细胞暴露于补充有0.01mg/mL的转铁蛋白、20uM的胰岛素、100nM皮质醇和0.2nM的3,3′,5-三碘代-L-甲状腺原氨酸的DMEM/F12附加的12天。在实施葡萄糖摄取之前，在SKMC分化介质(PromoCell)存在的情况下将HMCL-7304分化4-5天。

为了量化对于成脂分化的效果，在分化的第0天以不同的浓度(10nM、100nM、1uM、10uM)添加用于HMGCR(阿托伐他汀)、FDPS(阿仑膦酸钠)和SQLE(特比萘芬)的化学抑制剂。利用油红O实施分化量化。在室温下用10％福尔马林将分化的SGBS细胞固定10分钟。在用60％异丙醇清洗之后，在室温下将样本在油红O(Sigma-Aldrich)中培养10分钟。用100％异丙醇将油红O洗脱10分钟。随后将溶液转移到一块96阱板中，并且在500nm下测量吸光度。

以不同浓度(10nM、100nM、1uM、10uM)的用于HMGCR(阿托伐他汀)、FDPS(阿仑膦酸钠)和SQLE(特比萘芬)的化学抑制剂对分化的SGBS或HMCL-7304细胞进行24小时预处理。在预培养之后，在37度下进行30分钟的100nM胰岛素刺激之前，将细胞饿2小时。在刺激之后，在室温下用Krebs-Ringer碳酸氢钠-HEPES(KRBH)缓冲液(130mM的NaCl，5mM的KCl，1.3mM的CaCl

利用闪烁计数器(型号ID：Beckman LS6500)确定放射性计数。对多余样本进行BCA实验以进行蛋白质量化和放射性计数的规范化。所有样本都被表示成与无刺激(无胰岛素)条件相比的倍数改变。

SGBS或HMCL-7304细胞在50或100细胞/cm2下被装盘在12阱盘中，并且在没有或存在10nM、100nM、1uM、10uM的阿托伐他汀、特比萘芬或阿仑膦酸钠的情况下被生长12到14天。在进行处理之后，细胞在冷甲醇中被固定15分钟，并且用结晶紫着色10分钟。用水清洗多余的染料，并且随后理解拍摄照片。

本发明的系统和方法可以通过允许访问来自电子信息来源的数据的计算机软件来实施。根据本发明的软件和信息可以处在单个独立计算机内，或者可以处在联网到一组其他计算机或其他电子设备的中央计算机中。信息可以被存储在计算机硬盘驱动器、CD-ROM盘或者任何其他适当的数据存储设备上。

前面的描述和附图应当被视为仅仅说明本发明的原理。本发明不意图受限于优选实施例，而是可以通过本领域技术人员将会想到的多种方式来实施。本领域技术人员将很容易想到本发明的许多应用。因此，不希望将本发明限制到所公开的具体实例或者所示出和描述的确切构造和操作。相反，可以采取落在本发明的范围内的所有适当的修改和等效方案。