CN110473591B - 基于量子计算的基因网络功能模块挖掘及分析方法 - Google Patents

Abstract

本发明公开了一种基于量子计算的基因网络功能模块挖掘及分析方法。本发明从生物复杂网络的角度，利用重整化方法研究生物由正常到疾病再到治疗时，相应复杂网络重整化特征的变化以及新属性、新结构的产生，探讨生物复杂系统疾病发生、发展和治疗过程中涌现产生的机制。本发明利用量子计算的退火算法代替了传统的经典模拟退火算法，量子退火中的势能函数包含了经典模拟退火中的评价函数，又引入横向场使其包含量子动能函数，这样可以产生隧穿效应，平行改变所有状态幅度的量子力学概率，提高了计算效率。

Description

基于量子计算的基因网络功能模块挖掘及分析方法

技术领域

本发明涉及生命科学技术领域，具体涉及一种基于量子计算的基因网络功能模块挖掘及分析方法。

背景技术

在生命科学领域，有多种多样的生物复杂网络，有反映转录因子和基因间调控关系的基因调控网络，反映基因相互间通过共同表达来控制生命体微观活动和多样性的基因共表达网络等。这些基因网络由基因构成功能模块行使功能，对于基因网络中的模块进行分析有助于预测和指派未知基因的功能。基因功能模块的挖掘与分析，有较多的优良算法，比如谱聚类方法、模块度(Modularity)的优化和扩展方法、图论方法、随机游走方法等等。这些方法虽然各具优点，但是不管在基因网络，还是在算法方面都存在局限或缺陷。如何高效和准确地挖掘基因网络的功能模块，并进行研究分析，是生物信息学研究的重要的问题。

传统基因网络功能模块挖掘算法主要采用谱聚类、经典模拟退火、遗传算法等，时间复杂度都较高。比如经典模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解；遗传算法是一种通过模拟达尔文自然选择和遗传学机理的生物进化过程搜索最优解的方法。经典模拟退火算法陷入局部最优势阱后，要到达全局最优需要翻越能量势垒，具体的实现是按随机方式，以很小的概率作选择以跳出势垒，摆脱局部最优，然后再搜索到全局最优。遗传算法也有类似情况，它摆脱局部最优解，一般采取变异方式，基本是随机选择。所以，它们的共同点是能摆脱局部最优，但是会耗费大量时间，性能不高。与经典模拟退火算法和遗传算法相比，量子退火算法则不同，它利用量子隧穿效应(指在量子跃迁过程中，由于其具有波动性，在一定条件下，能够直接穿透比它能量更高的势垒现象)，由局部最优解穿过势垒到达全局最优解。

在模拟退火中，温度确定从单个当前状态移动到更高“能量”状态的概率；在量子退火中，横向场的强度决定了平行改变所有状态的幅度的量子力学概率；因此，量子退火可能比经典模拟退火算法和经典遗传算法有更高的性能(甚至能达到指数级的差别)。在实验和理论上已经证明，很多情况下，量子退火确实优于模拟退火，特别是在围绕局部极小值非常高但很薄的势垒情况下。

基因网络往往具有多层次结构，层次结构不仅表现为层次模块结构，同时也表现为模块内部元素之间的层次结构，它最主要特点是网络成份的异质性和系统多层次的复杂性，基因调控网络和基因共表达网络都具有典型的分层结构。从系统生物学角度出发，基因间共表达或调控关系实际上对应于基因信息在不同尺度空间上的复杂关联过程。从多层次多尺度去认识基因信息在不同尺度上的关联，探索尺度到尺度之间的转换和推演，有助于认识和理解基因间以及基因与生命现象表现间的复杂关联过程。

当生物复杂系统由正常变为疾病或由疾病治疗后逐步恢复正常的过程中，会在临界点附近出现自组织涌现，即其对应的复杂网络比如基因网络结构特征指数会变化，然后由无序进入有序状态。重整化作为一个研究复杂系统多层次多尺度现象的重要理论，是量子场论的方法，提供了不断粗粒化的途径，也是不断变换标度的重标度过程。复杂网络在不同尺度上，特别在临界点附近所表现出的自相似性恰好能用重整化的方法来描述，其自相似性具有标度不变性，重整化能够发现其结构特征指数。

发明内容

针对现有技术中的上述不足，本发明提供的一种基于量子计算的基因网络功能模块挖掘及分析方法解决了如何高效挖掘基因网络中功能模块并分析的问题。

为了达到上述发明目的，本发明采用的技术方案为：一种基于量子计算的基因网络功能模块挖掘及分析方法，包括以下步骤：

S1、通过小鼠样本建立基因共表达网络；

S2、挖掘基因共表达网络中的功能模块；

S3、通过多层次多尺度方法对基因网络进行分析。

进一步地：所述步骤S1的具体步骤为：

S11、由小鼠正常组、便秘组、便秘自发肠癌组和肠癌云南参治疗组样本通过转录组高通量方法进行测序，获取转录组基因表达数据；

S12、计算转录组基因表达数据中任意两个基因之间的表达相关性，并根据表达相关性建立基因共表达网络。

进一步地：所述步骤S12中表达相关性的计算公式为：

上式中，r_xy为基因x与基因y之间的表达相关性，k为样本计数，m为基因表达数据的个数，x_k为不同样本同一基因的表达数量，

为同一基因表达平均数量，y_k为不同样本另一种基因的表达数量，

为此种基因表达平均数量，其中，

进一步地：所述步骤S12中建立基因共表达网络的具体方法为：

将表达相关性转换为：

上式中，n为计算基因x与基因y的表达相关性所使用的数据点数，r₁为转换后的表达相关性；

当r₁大于预设的p-value对应的t分布值，则基因x与基因y显著表达相关，将基因共表达网络中连接基因x与基因y之间的边被赋予表达相关性r_xy作为其权重，否则连接基因x与基因y之间的边的权重为0。

进一步地：所述步骤S2的具体步骤为：

S21、计算基因共表达网络的模块性度量参数；

S22、由模块性度量参数构建挖掘功能模块的目标函数；

S23、利用量子退火算法和目标函数挖掘功能模块。

进一步地：所述步骤S21中模块性度量参数的计算公式为：

上式中，Q为模块性度量参数，A_vw为基因网络对应的邻接矩阵A中v结点和w结点间的带权连接值，k_v和k_w分别为v结点和w结点的度数，δ(C_v,C_w)为通过dalta函数判断v结点和w结点是否属于同个模块的取值，相同时取值为1，否则取值为0。

进一步地：所述步骤S23中量子退火算法的具体步骤为：

S231、设定初始动能值E₀和执行次数N；

S232、根据量子退火算法的评价函数采用蒙特卡洛方法搜索最优解；

S233、根据Γ(t)＝Γ₀β^t(0<β<1)策略和最优解逐步减小初始动能值E₀或Γ(t)，当达到执行次数N或达到了基态能量，进入步骤S234，否则返回步骤S232；

Γ(t)为引起不同状态之间跃迁的横向场，Γ₀为算法赋予的初始能量，β为衰减系数，t为时间；

S234、得到基态能量或执行次数N的最小能量，作为量子化的最优解；

S235、将量子化最优解转化为经典最优解，作为挖掘出的功能模块。

进一步地：所述步骤S232中量子退火算法的评价函数H_q为：

H_q＝H_pot+H_kin

上式中，H_pot为势能，对应经典模拟退火算法中的评价函数，H_kin为动能，是使得系统产生隧道效应跳出局部最优的能量。

进一步地：所述步骤S3的具体步骤为：

S31、利用重整化方法计算标度变化下的特征指数；

S32、利用重整化特征指数对基因网络的涌现机制进行分析。

进一步地：所述步骤S31中的重整化方法为基于贪婪着色的盒子覆盖法，其具体步骤为：

以任意结点作为开始结点，遍历基因网络图中的未被着色的每个结点，如果一个结点的邻接点都未用颜色a着色，则这个结点可以用颜色a着色，当没有结点能以颜色a着色时，选择颜色b和一个未被着色的结点作为开始结点，用颜色b为最多的结点着色，如果还有未着色的结点，选择颜色c为最多的结点着色，直到所有结点都着上色，统计使用颜色的数量，即为盒子的数量N_b，通过盒子的数量N_b即可计算特征指数d_b，计算公式为：

上式中，l_b为盒子的大小。

本发明的有益效果为：

(1)本发明从生物复杂网络的角度，利用重整化方法研究生物由正常到疾病再到治疗时，相应复杂网络重整化特征的变化以及新属性、新结构的产生，探讨生物复杂系统疾病发生、发展和治疗过程中涌现产生的机制。

(2)本发明利用量子计算的退火算法代替了传统的经典模拟退火算法，量子退火中的势能函数包含了经典模拟退火中的评价函数，又引入横向场使其包含量子动能函数，这样可以产生隧穿效应，平行改变所有状态幅度的量子力学概率，提高了计算效率。

附图说明

图1为本发明流程图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，一种基于量子计算的基因网络功能模块挖掘及分析方法，包括以下步骤：

S1、通过小鼠样本建立基因共表达网络；具体步骤为：

S11、由小鼠正常组、便秘组、便秘自发肠癌组和肠癌云南参治疗组样本通过转录组高通量方法进行测序，获取转录组原始数据和转录组基因表达数据。

基因共表达网络(Gene Co-expression Network)是用来展现和分析基因间相互作用关系的一种手段，也是基于基因间表达数据的相似性而构建的调控网络图。它可以根据基因表达信号值的动态变化，计算基因间的共表达关系，建立基因转录调控模型，得到基因间的表达调控关系及调控方向，通过基因表达的相似性可分析基因产物可能的相互作用关系，从而理清基因间相互作用脉络及寻找核心基因。Pearson相关系数(PearsonCorrelation Coefficient)是用来衡量两个数据集合间的相关性的，相关系数的绝对值越大，相关性越强：相关系数越接近于1(正相关)或-1(负相关)，相关度越强，相关系数越接近于0，相关度越弱。这里，应用Pearson相关系数构建基因共表达网络。假设共有m个不同的基因表达数据，为了得到相应的共表达网络，需要判断转录组基因表达数据中任意两个基因之间是否显著表达相关。

S12、计算转录组基因表达数据中任意两个基因之间的表达相关性，并根据表达相关性建立基因共表达网络。表达相关性的计算公式为：

上式中，r_xy为基因x与基因y之间的表达相关性，k为样本计数，m为基因表达数据的个数，x_k为不同样本同一基因的表达数量，

为同一基因表达平均数量，y_k为不同样本另一种基因的表达数量，

为此种基因表达平均数量，其中，

建立基因共表达网络的具体方法为：

将表达相关性转换为：

上式中，n为计算基因x与基因y的表达相关性所使用的数据点数，r₁为转换后的表达相关性；

当r₁大于预设的p-value对应的t分布值，则基因x与基因y显著表达相关，将基因共表达网络中连接基因x与基因y之间的边被赋予表达相关性r_xy作为其权重，否则连接基因x与基因y之间的边的权重为0。

S2、挖掘基因共表达网络中的功能模块；具体步骤为：

模块的基本要求是其内部高内聚，外部低耦合，按照模块度的定义，需要计算内聚性和耦合性度量。Newman模块度Q基本满足综合度量内聚性和耦合性要求，

S21、计算基因共表达网络的模块性度量参数；模块性度量参数的计算公式为：

上式中，Q为模块性度量参数，A_vw为基因网络对应的邻接矩阵A中v结点和w结点间的带权连接值，k_v和k_w分别为v结点和w结点的度数，δ(C_v,C_w)为通过dalta函数判断v结点和w结点是否属于同个模块的取值，相同时取值为1，否则取值为0。

S22、由模块性度量参数构建挖掘功能模块的目标函数；

模块性度量参数Q可以用来衡量模块划分优劣，所以，它也可以作为基因网络挖掘功能模块的目标函数。引入模块性度量作为目标函数后，基因网络模块挖掘方法就成了优化方法。以模块度参数为优化的目标函数，每次合并使得Q增加最大的两个结点，直到得到的R减小时，停止合并。但是，这样优化模块度Q实际是一个NP问题，只能进行近似求解，因此，可以引入基于贪婪策略的算法，即首先按照期望最大的原则，选择尽可能多地影响模块性度量的结点，从中选择两个使目标函数Q增加最大的结点进行合并，这样就有了基因网络功能模块的最优近似挖掘方法。

S23、利用量子退火算法和目标函数挖掘功能模块。

量子退火算法模型由量子势能和量子动能两部分构成，前者部分将优化问题映射成量子系统，优化的目标函数映射成势场；后者部分引入量子波动可控的穿透场(横向场)。其量子系统的演化描述可以用含时薛定谔方程：

H(t)＝H_pot+H_kin()

上式中，H_pot为势能项，H_kin(t)为含时动能项，其值越大，量子波动也较大。因为量子波动有穿透性，会产生隧穿效应，所以量子退火采用量子波动重构退火算法实现最优化，可以使较低能量的量子穿过较高能量的势垒或摆脱较高能量的势阱约束。量子退火在经典物理系统中引入横向场(Tunneling Field)，开始时横向场能量较大，粒子有足够能量隧穿，因而能够搜索全部空间。然后，以某种方式减小横向场能量，直到系统最终停留的能量基态，即能量最低状态。横向场能量对应的就是上面公式中的含时动能项。量子退火算法改造了经典退火算法，引入动能项所产生的量子隧穿效应摆脱局部最优，以较小能量值实现了最优化策略。这里，要用计算机模拟动能项变化求解含时薛定谔方程，因此采用量子蒙特卡洛方法来模拟量子退火的随机过程。

这里，可以将问题中的变量视为经典自由度，并且将成本函数视为势能函数(经典哈密顿量)，在哈密顿量中引入隧道横向场的作用作为动力学部分构造量子哈密尔顿量进行模拟。

量子退火算法的具体步骤为：

S231、设定初始动能值E₀和执行次数N；

S232、根据量子退火算法的评价函数采用蒙特卡洛方法搜索最优解；

评价函数H_q为：

H_q＝H_pot+H_kin

上式中，H_pot为势能，对应经典模拟退火算法中的评价函数，H_kin为动能，是使得系统产生隧道效应跳出局部最优的能量。

S233、根据Γ(t)＝Γ₀β^t(0<β<1)策略和最优解逐步减小初始动能值E₀或Γ(t)，当达到执行次数N或达到了基态能量，进入步骤S234，否则返回步骤S232；

Γ(t)为引起不同状态之间跃迁的横向场，Γ₀为算法赋予的初始能量，β为衰减系数，t为时间；

S234、得到基态能量或执行次数N的最小能量，作为量子化的最优解；

S235、将量子化最优解转化为经典最优解，作为挖掘出的功能模块。

体系的Hamilton量可以表示为：

式中，H₀为体系在无外力作用时的Hamilton量；

为外力作用的结果，设定为横向场，表示为一个引起必要的量子跃迁的适当动能；Γ(t)表示表示引起不同状态之间跃迁的横向场，起控制作用；

表示第i个粒子在x轴的自旋泡利表象。令

它控制着量子跃迁的概率，保证在保留搜索到的模型参数时有较大的概率，从而能发挥量子的隧穿效应。为了使目标函数达到全局最小，需要缓慢地减小K值。其中，Γ(t)基本只与β^t相关，比较容易确定，但是

关于粒子间的相互作用，较为复杂，也很重要，这里，需要对伊辛(Ising)模型进行计算机模拟，然后蒙特卡洛方法确定它们的数值。

S3、通过多层次多尺度方法对基因网络进行分析。具体步骤为：

S31、利用重整化方法计算标度变化下的特征指数；

基因网络拥有大量自由度，对其进行研究时，需要寻找真正重要的自由度。如果利用统计学习方法，一般采用降维方法来对系统的自由度进行约简，构造新的主要变量。当然，也可以用重整化(Renormalization)的思路来求解这一问题。降维方法比如主成分分析方法就是寻找被研究复杂系统的主成分方向，而重整化方法关注的是复杂系统的自相似特征。基因网络基本属于无标度网络，它们的结点的度符合幂律分布，也有部分相似于整体的自相似特征。重整化就是想要提取隐藏在这些自相似以及多层次多尺度中的特征，希望挖掘出系统在不同的层次和尺度之间的关系。

重整化方法为基于贪婪着色的盒子覆盖法，其具体步骤为：

以任意结点作为开始结点，遍历图中的未被着色的每个结点，如果一个结点的邻接点都未用颜色a着色，则这个结点可以用颜色a着色，当没有结点能以颜色a着色时，选择颜色b和一个未被着色的结点作为开始结点，用颜色b为最多的结点着色，如果还有未着色的结点，选择颜色c为最多的结点着色，直到所有结点都着上色，统计使用颜色的数量，即为盒子的数量N_b，通过盒子的数量N_b即可计算特征指数d_b，计算公式为：

上式中，l_b为盒子的大小。

重整化过程通过对网络进行不同盒子大小l_b的着色来实现不同尺度上的缩放，而指数d_b即是重整化下的生物网络的尺度不变性的特征。

幂律函数指数d_b的尺度不变性表明不同尺度下的重整化网络的度分布仍然服从幂律分布。因为基因网络已经划分成多个模块，所以每个模块实际构成了它的子网，则每个模块可以计算得出自己子网的幂律函数特征指数d_b。然后，根据特征d_b的差异，利用支持向量机SVM进行模块的聚类，得出各自模块的构造相似子网和差异子网。接着，将每个模块看作一个结点，在更高的层次上构建网络，继续采用上述方法挖掘模块和聚类，可以得出不同层次上模块划分和聚类。这样，可以在不同层次上发现特征相似的不同模块，进而可以在功能比如代谢通路或信号通路上寻找它们间的联系。

S32、利用重整化特征指数对基因网络的涌现机制进行分析。

由大量子系统组成的系统的可测宏观量在每一时刻的实际测度相对平均值或多或少有些偏差，这些偏差就叫涨落，涨落是偶然的、杂乱无章的、随机的。在正常情况下，涨落相对于平均值是很小的，即使偶尔有大的涨落也会很快耗散掉，系统基本回到平均值附近，这些涨落不会对宏观的实际测量产生影响，因而可以被忽略掉。然而，如果持续性涨落很大，越过了临界点，情况就大不相同了，这时涨落不可能自生自灭，而是在不稳定的情况下被系统放大，最后促使系统达到新的宏观态。

复杂系统由大量相互作用的单元组成，其活动呈现非线性，往往形成具备多层级的复杂组织。复杂系统平衡态改变后生成新的平衡态时，会伴随结构或属性等的变化，一个比较明显的现象是涌现，即系统可能由一种自相似结构自组织变化为另一种自相似结构，这时，重整化特征指数(幂律指数)会发生差异性变化。生物复杂系统由正常态转变为疾病态，再转变为治疗态，实际也会在临界点附近出现自组织涌现，伴随着基因共表达网络中幂律指数d_b的显著变化，这也就是它们重整化特征的变化。

利用小鼠建模，正常状态、便秘状态、疾病状态、治疗状态是小鼠不同生理阶段的不同平衡状态，这些平衡态都会存在系统涨落。只有涨落较大，持续时间长，持续越过临界点后，才会在新平衡点附近达到新的平衡。当新平衡态形成时，会在临界点附近出现自组织涌现，形成新的自相似结构，基因网络功能模块结构特征也会发生显著变化。小鼠模型由正常态到便秘态，正常态到疾病态，疾病态到治疗态都可以发现基因共表达网络结构的差异，得到不同的幂律指数d_b，也就可以区别正常、疾病和治疗它们之间重整化特征的不同。这样，就能区别基因共表达网络中重整化结构特征指数d_b在产生涌现时的差别，也就可以从系统科学复杂网络角度探讨疾病发生、发展和治疗过程中涌现产生的机制。

Claims (6)

1.一种基于量子计算的基因网络功能模块挖掘及分析方法，其特征在于，包括以下步骤：

S1、通过小鼠样本建立基因共表达网络；

所述步骤S1的具体步骤为：

S11、由小鼠正常组、便秘组、便秘自发肠癌组和肠癌云南参治疗组样本通过转录组高通量方法进行测序，获取转录组基因表达数据；

S12、计算转录组基因表达数据中任意两个基因之间的表达相关性，并根据表达相关性建立基因共表达网络；

S2、挖掘基因共表达网络中的功能模块；

所述步骤S2的具体步骤为：

S21、计算基因共表达网络的模块性度量参数；

S22、由模块性度量参数构建挖掘功能模块的目标函数；

S23、利用量子退火算法和目标函数挖掘功能模块；

S3、通过多层次多尺度方法对基因网络进行分析；

所述步骤S3的具体步骤为：

S31、利用重整化方法计算标度变化下的特征指数；

所述步骤S31中的重整化方法为基于贪婪着色的盒子覆盖法，其具体步骤为：

以任意结点作为开始结点，遍历基因网络图中的未被着色的每个结点，如果一个结点的邻接点都未用颜色a着色，则这个结点用颜色a着色，当没有结点能以颜色a着色时，选择颜色b和一个未被着色的结点作为开始结点，用颜色b为最多的结点着色，如果还有未着色的结点，选择颜色c为最多的结点着色，直到所有结点都着上色，统计使用颜色的数量，即为盒子的数量N_b，通过盒子的数量N_b计算特征指数d_b，计算公式为：

上式中，l_b为盒子的大小；

S32、利用重整化特征指数对基因网络的涌现机制进行分析。

2.根据权利要求1所述的基于量子计算的基因网络功能模块挖掘及分析方法，其特征在于，所述步骤S12中表达相关性的计算公式为：

上式中，r_xy为基因x与基因y之间的表达相关性，k为样本计数，m为基因表达数据的个数，x_k为不同样本同一基因的表达数量，

为同一基因表达平均数量，y_k为不同样本另一种基因的表达数量，

为此种基因表达平均数量，其中，

3.根据权利要求2所述的基于量子计算的基因网络功能模块挖掘及分析方法，其特征在于，所述步骤S12中建立基因共表达网络的具体方法为：

将表达相关性转换为：

上式中，n为计算基因x与基因y的表达相关性所使用的数据点数，r₁为转换后的表达相关性；

当r₁大于预设的p-value对应的t分布值，则基因x与基因y显著表达相关，将基因共表达网络中连接基因x与基因y之间的边的表达相关性r_xy作为其权重，否则连接基因x与基因y之间的边的权重为0。

4.根据权利要求1所述的基于量子计算的基因网络功能模块挖掘及分析方法，其特征在于，所述步骤S21中模块性度量参数的计算公式为：

上式中，Q为模块性度量参数，A_vw为基因网络对应的邻接矩阵A中v结点和w结点间的带权连接值，k_v和k_w分别为v结点和w结点的度数，δ(C_v,C_w)为通过dalta函数判断v结点和w结点是否属于同个模块的取值，相同时取值为1，否则取值为0。

5.根据权利要求1所述的基于量子计算的基因网络功能模块挖掘及分析方法，其特征在于，所述步骤S23中量子退火算法的具体步骤为：

S231、设定初始动能值E₀和执行次数N；

S232、根据量子退火算法的评价函数采用蒙特卡洛方法搜索最优解；

S233、根据Γ(t)＝Γ₀β^t策略和最优解逐步减小初始动能值E₀或Γ(t)，当达到执行次数N或达到了基态能量，进入步骤S234，否则返回步骤S232；其中0<β<1；

Γ(t)为引起不同状态之间跃迁的横向场，Γ₀为算法赋予的初始能量，β为衰减系数，t为时间；

S234、得到基态能量或执行次数N的最小能量，作为量子化的最优解；

S235、将量子化最优解转化为经典最优解，作为挖掘出的功能模块。

6.根据权利要求1所述的基于量子计算的基因网络功能模块挖掘及分析方法，其特征在于，所述步骤S232中量子退火算法的评价函数H_q为：

H_q＝H_pot+H_kin

上式中，H_pot为势能，对应经典模拟退火算法中的评价函数，H_kin为动能，是使得系统产生隧道效应跳出局部最优的能量。

Images