CN111009292B - 基于单样本sKLD指标检测复杂生物系统相变临界点的方法 - Google Patents

Abstract

本发明公开了一种基于单样本sKLD指标检测复杂生物系统相变临界点的方法，通过从高通量的数据中挖掘丰富的动态信息并利用正常状态和临界状态之间的不同动态特性，量化单一样本对参照组样本分布的干扰，从而确定临界状态或相变的早期预警信号。为了验证算法的有效性，本发明将该算法应用于一个基于模拟数据集的调控网络和六个真实的数据集。这六个真实数据集分别是：小鼠实验中通过吸入羰基氯引发肺水肿生成的基因表达谱的数据集、5个TCGA数据库的癌症数据集(肺鳞癌、肺腺癌、胃腺癌、甲状腺癌、结肠癌)。

Description

基于单样本sKLD指标检测复杂生物系统相变临界点的方法

技术领域

本发明涉及生物系统相变临界点检查的技术领域，具体涉及一种基于单样本sKLD指标检测复杂生物系统相变临界点的方法。

背景技术

生物系统的发展进程通常是一个非线性过程，具有三个阶段，即正常状态、临界状态和疾病状态，其中临界状态是正常状态进入疾病状态的临界点。传统的生物标记物旨在通过利用观察到的分子的差异表达信息来识别疾病状态，但是由于正常和临界状态之间通常没有显著差异，所以可能无法检测到复杂生物系统相变的临界。因此，对临界状态进行信号预警是一个挑战，这实际上意味着对复杂生物系统相变临界点进行预测。

下面介绍计算方法的理论推导如下：

在临界相变之前和附近不同的动态特性：

复杂生物系统相变临界点进展的动力学可以通过以下非线性离散时间动态系统来表示：Z(t)＝f(Z(t-1)；P),这里Z(t)＝ (z₁(t),z₂(t),…,z_n(t))是n维状态向量或在时间t＝1,2,…的变量，P＝ (P₁,…,P_s)是一个代表缓慢变化因素的参数向量或者驱动因素，例如，遗传因素(SNP,CNV等)，表观遗传因素(methylation，acetylation等)或者环境因素。f：Rⁿ×R^s×Rⁿ是一个非线性函数。对于这样一个非线性系统，该系统在

处将经历一个相变或者是一种当参数P达到阈值P_c时来自稳定平衡的分叉(Gilmore,1993)。

对于z附近的系统(1)，在P到达P_c之前，系统应该保持稳定的平衡

因而所有的特征值的绝对值都在(0,1)内。使系统状态发生移位的参数值P_c称为一个分岔参数值或一个临界值，而在这种分歧之前的状态被称为前疾病状态。为了从理论和数学上描述生物系统的动力学，通常将其演化建模为一个含时非线性动力学系统，在该系统中，突变被视为分叉点处的时期转移。在小噪声的理想情况下，当一个复杂系统接近临界点时，在所有观测变量中，存在一个被定义为动态网络生物标记物生物分子的优势群，它基于观测数据满足以下三个条件(Chen et al.,2012；Liu et al.,2012,2013a, 2014b)

1.动态网络生物标记物组中任何一对分子之间的相关性迅速增加；

2.动态网络生物标记物组中任何的一个分子与任何其他非动态网络生物标记物分子之间的相关性迅速降低；

3.动态网络生物标记物组中任何一个分子的标准差急剧增加。

因此，从上述三个性质可以看出，系统的临界转变实际上可以由“分布转变”来表示，即对于动态网络生物标记物，当系统接近临界转变点时，它们的分布会发生显著变化。因此，通过探索这些变量的分布，我们可以预测即将到来的定性时期转变。另一方面，一个高通量数据样本使研究者能够同时分析成千上万种生物分子的表达。这样一个高维的样本实际上富含了累积效应的动态信息，比如相关生物系统在长期发展中的生物分子的相互作用。

大多数生物分子通过与功能模块或模块之间的其它生物分子的相互作用来执行其功能。这种模块间和模块内的互连性表明，特定遗传异常的影响不仅影响携带它的基因产物的活性，而且可以沿着由生物分子组成的网络的链接延伸，改变其它基因产物的活性。因此，了解生物分子的相互作用网络环境对于确定影响生物分子的缺陷的表型至关重要。

发明内容

本发明的目的是通过利用正常状态和临界状态之间的不同动态特性，提出一种基于单样本sKLD指标检测复杂生物系统相变临界点的方法。同时，高通量技术的快速发展为计算分析提供了新的视角，即使只有一个样本也可以识别临界状态。实际上，基于高通量数据样本，可以同时测量数千个生物分子的表达，这种全生物分子组范围内的高维观测不仅提供了一个生物系统的全局视图，而且还展示了其长期动态的累积效应。基于这一点，本研究开发了一种数据驱动的计算方法，实现了从高通量的数据中挖掘丰富的动态信息以实现基于单样本识别复杂生物系统相变临界点的目的。

本发明的目的可以通过采取如下技术方案达到：

一种基于单样本的sKLD指标检测复杂生物系统相变临界点的方法，给定至少6个参照样本(来自正常时期的样本)，仅使用一个实验样本来识别临界状态，所述的方法包括如下步骤：

S1、准备一组参照样本。以正常时期提取的样本作为参照样本，也就是把处于正常时期的相对健康个体作为背景。对于数值模拟，将初始时间点采集的样本被视为参照样本。对于真实的数据集，选择来自正常组织的样本作为参照样本；

S2、根据参照样本，拟合生物分子的分布。具体来说，对于生物分子 g_i，基于参照样本{s1，s2，…，sk}中的表达水平拟合高斯分布。然后，得到一个k维向量(area(D_gi(S₁)),area(D_gi(S₂)),…,area(D_gi(S_k)))，其中， area(D_gi(S_k))代表第k个样本中的生物分子gi被高斯分布所确定的累积面积；

S3、按照如下公式构造参照分布P

其中，area(D_gi(S_k))代表第k个样本中的生物分子gi被高斯分布所确定的累积面积，对于分布P，有

S4、对于单个实验样本，基于这单个实验样本构造扰动分布Q，如下所示

其中，area(D_gi(S_case))代表的是实验样本的生物分子gi被对应的高斯分布所确定的累积面积，对于分布Q,有

S5、根据公式(A3)计算sKLD，sKLD用于评估参照分布P和扰动分布Q之间的差异

(A3)这里P和Q代表是离散概率分布，分布P和Q分别满足

其中P(k)和Q(k)分别是参照分布和扰动分布的第k个生物分子的概率值。

根据动态网络生物标记物理论，当复杂生物系统接近临界时期时，动态网络生物标记物生物分子表现出明显的集体波动行为，这导致了动态网络生物标记物生物分子处于前疾病时期的分布与处于正常时期的明显不同。因此，一组参照/正常样本的背景分布与一个新样本的扰动分布显著不同，导致公式(A3)中sKLD明显增加。因此，sKLD可以提供临界转变的预警信号。从上面的算法可以看出，该方法是数据驱动的，因此不需要模型。

进一步地，该方法需要至少6个正常样本作为参照样本。

进一步地，利用公式(A3)(sKLD指标)评估参照分布P和扰动分布 Q之间的差异，sKLD指标突然迅速增加表示即将到来的临界转变，或者复杂生物系统相变临界点的出现。

本发明相对于现有技术具有如下的优点及效果：

本发明给出了一个基于单样本的sKLD指标的计算方法，用来识别即将到来的临界转变，它被真实数据集证明是有效的。值得注意的是，本发明的目的是检测临界状态产生的早期预警信号，而不是找到发生质变的疾病状态的迹象。本发明的创新之处有：

1、传统的方法只能判断个体是处于正常状态还是疾病状态，而对于正常状态的极限状态也即临界转变关键期不能有效地察觉，本发明能够准确地反映出复杂生物系统发展过程中的临界状态的发生；

2、本发明方法是一种无模型方法，这意味着既没有特征选择，也没有模型/参数训练过程。因此，它不同于传统的机器学习或分类方法，在学习过程中产生一个鲁棒模型，需要大量的样本来避免过度拟合问题；

3、现有技术中单变量或少数变量受噪声影响太大，临界点信号不明显，而本发明方法可以克服；

4、本发明方法朝着检测复杂生物系统相变临界点进展中的关键转变的预警信号迈进了一步，它为识别生物系统进展中的临界转变的早期预警信号开辟了一条新的途径；

5、给定一些参照样本，sKLD指标基于单一样本也可以识别临界状态，即使只有一个样本也可以识生物系统的临界状态，在个体化的预防医学中具有巨大的潜力。

附图说明

图1是基于sKLD指标识别疾病前状态的流程示意图；

图2(A)是由模型控制八个节点的网络的示意图；

图2(B)是通过数值仿验证sKLD指标的验证图；

图2(C)是参照分布P和扰动分布Q的频率分布的分布图；

图3(A)是基于来自每个时间点的急性肺损伤数据集的sKLD指标的验证图；

图3(B)是急性肺损伤数据集的对照组和病例组的频率分布的分布图；

图3(C)是累积面积的变化最显著的前5％生物分子组成的网络的动力学演化示意图；

图4(A)是肺鳞癌的临界时期的识别示意图；

图4(B)是肺腺癌的临界时期的识别示意图；

图4(C)是胃腺癌的临界时期的识别示意图；

图4(D)是甲状腺癌的临界时期的识别示意图；

图4(E)是结肠癌的临界时期的识别示意图；

图4(F)是肺鳞癌的临界时期前后生存曲线的比较示意图；

图4(G)是肺腺癌的临界时期前后生存曲线的比较示意图；

图4(H)是胃腺癌的临界时期前后生存曲线的比较示意图；

图4(I)是甲状腺癌的临界时期前后生存曲线的比较示意图；

图4(J)是结肠癌的临界时期前后生存曲线的比较示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例一

基于sKLD指标在数值模拟中的验证

如图1所示，本发明所实施的案例公开了一种基于sKLD指标检测复杂生物系统相变前的临界状态的方法。根据图1公开的流程示意图，该实施案得到的结果如下：

1、识别八节点网络的疾病前状态：

使用八节点人工网络的模型(图2(A))来验证所提出的计算方法。该网络是一组八个生物分子的调控表示，由八个随机微分方程控制。

本实施例使用一个包含8个基因的调控网络(见图2(A))进行数值模拟，以使用SKLD检测生物系统发展过程的临界时期。这些类型的基因分子调控网络通常用于研究各种生物过程，包括影响基因调控活性的转录、翻译、扩散和易位过程。下面的8个微分方程代表了网络中8个基因的基因调控。在网络中，除降解率与相应基因浓度成线性比例外，基因调控以 michaelis-menten形式表示。

其中s是标量控制参数，

(i＝1，2，…，8)是均值和协方差为零的高斯噪声。z_i(t)(i＝1，2，…，8)表示mRNA-i的浓度。在式(1)中，

表示为mRNA-i的降解率的大小，

是微分方程式(1)的稳定平衡点。利用Euler 格式可以将微分方程(1)转化为差分方程。结果如下：

很容易注意到Z(k)是时间瞬间的向量。式(2)的雅可比矩阵可以定义为

并且

J＝e^Δt·m (3)

其中m是方程(2)线性化后的矩阵，从式(3)中，取Δt＝1，方程有8个不同的特征值，当_s→0时，方程的最大特征值0.66→1。因此，当s∈(0,1] 时，平衡点是稳定的，但是，当S_c＝0，系统会发生临界的变化。在理论模型方程(2)的基础上，收集了8个节点在各个参数条件下的时间序列数据。然后，模拟了平均sKLD曲线，如图2(B)所示。

这样一个模型是以michaelis-menten形式表示的。这类调控网络通常用于研究遗传调控，包括转录和翻译过程，以及多稳定性和非线性生物过程。此外，michaelis-menten形式的分叉经常被用来模拟基因调控网络的时期转移。参数s在-0.5到0.2之间变化。在此模型的基础上，生成了数值模拟数据集。在图2(B)中可以看到，当系统接近一个特殊的参数值s＝0时，sKLD突然急剧上升，该参数值被设置为分岔点。换言之，在参数值s＝0 附近，sKLD的突然急剧上升表示参照分布P与由单个实验样本生成的扰动分布Q显著不同。为了更好地说明正常时期和临界时期之间的不同分布，在图2(C)中用一系列参数值演示了P和Q的频率分布的动态过程，其中每个频率分布是基于一万个模拟的统计图。从这些频率分布图可以看出，临界时期(s＝0)的扰动频率分布Q呈现两个峰值，即当网络系统处于临界时期时，一些节点的表达式以强烈的集体方式剧烈波动，从而产生的分布与正常时期的参照分布显著不同。sKLD能够准确地检测到这一临界现象，为识别即将到来的分岔点提供了一个定量的衡量。数值实验验证了基于sKLD指标在疾病前状态信号传递过程中的可靠性和准确性。

2、预测真实数据集的临界点：

本实施例将基于sKLD指标应用于六个真正的实验数据集，即小鼠实验中通过吸入羰基氯引发肺水肿生成的基因表达谱数据集(GSE2565)和5 个TCGA数据库的癌症数据集(肺鳞癌、肺腺癌、胃腺癌、甲状腺癌、结肠癌)。

实施例二

基于sKLD指标在光气吸入致急性肺损伤数据集中的应用

sKLD指标算法已应用于(GSE2565)数据集的微阵列数据，该数据集来源于光气诱导急性肺损伤的小鼠实验。在原实验中，实验样本的基因表达数据来源于光气暴露72小时的cd-1雄性小鼠肺组织，而对照样本的基因表达数据来源于暴露于空气中的cd-1雄性小鼠肺组织。在实验过程中，实验组和对照组共有9个采样点，分别为0、0.5、1、4、8、12、24、48和72小时，在每个采样点取6-8只小鼠的肺组织。将第一时间点(0小时)的样本作为参照组样本。如图3(A)所示，sKLD在1到4h之间突然增加，并在8h达到峰值，这表明在8h左右存在一个临界转变。为了证明结果的准确性，利用leave-one的方法生成了六个数据集。将sKLD算法分别应用于这些数据集，得到6个sKLD平均得分，并绘制为图3(A)中的曲线。可以看出，基于重新采样数据集的这些sKLD曲线都表示8h时是临界点。在图3(B)中，它显示了实验组和对照组分布的动态变化。显然，对于对照组，扰动分布几乎没有动态差异，而对于实验组，在第4个采样时间点 (8h)的扰动分布与其他采样时间点的扰动分布明显不同，导致sKLD的显著变化。在图3(C)中，展示了由实验样本累积面积的变化最显著的前 5％基因组成的网络的动态演化。显然，在8小时左右，网络结构发生了明显的变化，从网络层面角度预示着即将到来的临界转变。

简而言之，这些分析表明在吸入光气后的前8小时内，实验组的主要病理过程导致了支气管肺泡灌洗液蛋白水平升高，进而引发肺水肿，最终增加死亡率。严重的光气所致急性肺损伤大约出现在第8小时，暴露后持续至12小时，随着光气接触的继续，12小时后死亡率为50％-60％，24小时后死亡率为60％-70％。

实施例三

基于sKLD指标在5个肿瘤数据集的应用

为了进一步证明该方法的有效性，将其应用于5个肿瘤数据集：肺鳞状细胞癌、肺腺癌、胃腺癌、甲状腺癌、结肠癌，所有这些数据都来自于 TCGA癌基因图谱，由肿瘤和肿瘤邻近样本组成。根据TCGA相应的临床资料，将肿瘤分为不同的分期。肺鳞状细胞癌、肺腺癌和胃腺癌可以分为 7个时期，甲状腺癌、结肠癌可以分为4个时期。在所有5个数据集中，肿瘤邻近样本被用作正常/参照样本。然后根据sKLD算法计算每个单个肿瘤样本的sKLD。最后，取每个阶段的平均sKLD来确定肿瘤的临界时期。

sKLD指标成功地识别了所有五种癌症恶化前的关键阶段(图4(A) -图4(E))。为了验证所识别的临界时期，需对临界转化前和转换临界转化后的样本进行kaplan-meier(log-rank)生存分析进行比较(图4(F)-图 4(J))。临界转化前的样本的预后寿命通常比临界转化后的样本高。具体地说，对于肺鳞状细胞癌，从图4(F)可以看出，临界时期前样本(IA-IIA 期的样本)的存活时间比临界时期后样本(IIB-IV阶段的样本)的存活时间长得多，两组样本生存曲线之间存在显著差异(显著值p＝0.0034)。肺腺癌的ⅡB期前后样本的生存曲线有显著性差异(p＝3e-07，图4(G))，临界前样本(IA-IIB期的样本)的存活时间比临界时期后样本(IIIA-IV期的样本)长得多。同样地，对于胃腺癌、甲状腺癌和结肠癌，也可以得出相同的结论。这些结果表明所确定的临界阶段是准确的，并且与预后密切相关。

综上所述，本发明利用观测到的正常和临界状态状态分子间的分布差异关联信息，提出的基于单样本sKLD指标的计算方法，能够准确地反映生物系统临界状态的发生。这种方法与现有方法的区别在于，它为计算分析提供了新的视角，在只有单一病例样本的情况下，也可以识别复杂生物系统相变临界点。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (1)

1.一种基于单样本sKLD指标检测复杂生物系统相变临界点的方法，其特征在于，所述的方法包括如下步骤：

S1、准备参照样本，以正常时期提取的样本为参照样本，也就是把处于正常时期的相对健康个体作为背景，对于数值模拟，将初始时间点采集的样本被视为参照样本；对于真实的数据集，选择来自正常组织的样本作为参照样本；

S2、根据参照样本，拟合生物分子的分布，具体如下

对于生物分子g_i，基于参照样本{s1，s2，…，sk}中的表达水平拟合高斯分布；然后，得到一个k维向量(area(D_gi(S₁))，area(D_gi(S₂))，…，area(D_gi(S_k)))，其中，area(D_gi(S_k))代表第k个样本中的生物分子g_i被高斯分布所确定的累积面积；

S3、按照如下公式构造参照分布P

其中，area(D_gi(S_k))代表第k个样本中的生物分子g_i被对应的高斯分布所确定的累积面积，对于分布P，有

S4、对于单个实验样本，基于这单个实验样本构造扰动分布Q，如下所示

其中，area(D_gi(S_case))代表的是实验样本的生物分子gi被对应的高斯分布所确定的累积面积，对于分布Q,有

S5、根据公式(A3)计算sKLD，sKLD用于评估参照分布P和扰动分布Q之间的差异；

这里P和Q代表是离散概率分布，分布P和Q分别满足

其中P(k)和Q(k)分别是参照分布和扰动分布的第k个生物分子的概率值；

该检测复杂生物系统相变临界点的方法至少需要6个正常样本作为参照样本；

所述的sKLD指标突然迅速增加预示即将到来的临界转变，或者复杂生物系统相变临界点的出现。

Images