CN105718999A - 一种启发式代谢共表达网络的构建方法及系统 - Google Patents

Abstract

本发明公开一种启发式代谢共表达网络的构建方法及系统。本发明根据最大依赖准则将多个代谢物特征的互信息作为适应度函数值，并使用启发式计算智能多模优化算法对最优特征子集进行寻优搜索。通过多次运行这一优化过程，将各次结果进行组合学习，构造共表达网络结构。最后，通过概率模型计算切分阈值，获得准确、稳定的代谢共表达网络。

Description

一种启发式代谢共表达网络的构建方法及系统

技术领域

本发明涉及代谢组学网络领域，尤其涉及一种启发式代谢共表达网络的构建方法及系统。

背景技术

代谢物是生物体内完成代谢过程的小分子有机化合物总称，包含了丰富的生理状态信息。代谢组学基于对代谢物的整体系统性研究，可有效揭示生理现象背后的真实机理，并更为全面地展示生命体的动态状态。因此获得了越来越多的重视，被广泛应用于诸多科研与实用领域中。而另一方面，传统机器学习方法往往难以应对代谢组学特征高维度、小样本、高噪声的数据特点。使用创新的网络结构描述代谢物间的相互关联，并以此进行准确、稳定的分析，是代谢组学未来发展的重要方向。

已有的代谢组学网络描述方法主要包括两类：

其一是全基因代谢网络重构方法。其以基因表达信息为基础，通过获取其可能产生的蛋白质列表，搜索EC(EnzymeCommissionNumber)数据库得到对应的蛋白酶，并根据代谢途径(Pathway)数据库取得所有可能的化学反应，使用连接算法组合为包含高假阳性的代谢网络草图。而后根据在特定条件下的实验表达信息，对草图进行修正与剪裁，最终获得较为准确的网络结构。

其二是代谢网络的共表达构建方法(MetabolicCo-expressionNetwork)。直接评估不同代谢物特征在各实验条件下的表达差异，通过计算相关性参数(CorrelationCoefficient)形成权值矩阵。而后人为设定或使用适应性算法确定切分阈值，对矩阵进行简化，最终映射为网络结构。

一般认为，代谢共表达网络可更为有效地描述未知的生理关联信息，且对先验知识要求较低，更适合于非针对性代谢组学研究，是发掘与分析代谢组学新知识的有力工具。但在生物数据中，其相关性参数的计算往往有着较大误差，且人为设定的切分阈值缺乏理论依据，导致最终结果难以令人满意。针对这一问题，近年来提出了基于特征选择的共表达网络构建方法，获得了学术界的广泛重视。

现有的全基因代谢网络重构方法，其缺点在于：

第一，其包含了现有数据库中所有可能的代谢反应，具有极高的假阳性。尽管实验数据可部分消除此类网络连接，但若要对其进行准确修正，所需的样本量过高，成本较高。

第二，其严重依赖于现有的基因表达、酶催化及代谢途径等先验知识。而此类知识，特别是代谢组学相关的数据库仍有着大量信息缺失。导致所构建的网络具有高假阴性。此外，其网络完全基于现有知识，难以用于新生物信息的发掘。

现有的代谢共表达网络构建方法，其缺点在于：

第一，基于相关性参数的方法，如Pearson相关系数、Spearman相关系数等，其参数计算所需样本量较高，在生物实验中往往难以满足。导致所估计的相关程度有所偏差，构造网络的鲁棒性不佳。且人为设定切分阈值并无理论支撑，容易再次引入误差，影响分析结果。

第二，现有算法仅能估计两两特征(PairwiseFeatures)间的关联信息。而在真实生命体中，多个代谢物往往会相互连接形成功能模块，整体对生理过程进行调控。现有方法并不能有效描述这一特点。

第三，现有基于特征选择的网络构建方法一般使用确定性搜索方法，对于相同数据集仅能获得唯一的特征子集。而对于高维代谢组学数据，此类解往往不是最优的。此类方法也无法通过多次运行程序搜索更佳的结果。

因此，现有技术还有待于改进和发展。

发明内容

鉴于上述现有技术的不足，本发明的目的在于提供一种启发式代谢共表达网络的构建方法及系统，旨在解决现有的构建方法准确性低、稳定性差、成本高等问题。

本发明的技术方案如下：

一种启发式代谢共表达网络的构建方法，其中，包括步骤：

A、对原始的代谢特征数据集F^*进行正规化预处理，使其中所有M个代谢特征矢量在每个维度上都具有0均值及单位方差:

$F_m=\frac{F_m^{*}-{\mathrm{\μ}}_m}{{\mathrm{\δ}}_m},F_m^{*}\∈F^{*};$

F＝{F_m；m＝1,2,…,M}为预处理后的代谢特征数据集，μ_m、δ_m分别为第m个原始的代谢特征矢量F^* _m的均值与方差；

B、设定特征选择总运行次数为K，初始化运行计数器k＝1；

C、构造多模优化的进化种群ps，将其所包含的每个寻优个体X_i∈ps初始化为范围R＝[0,1]内均匀分布的M维随机矢量；

D、设定算法迭代总次数为G，初始化迭代计数器g＝1；

E、计算进化种群ps中每个寻优个体的共享适应度函数值；

F、在计算所有寻优个体的共享适应度函数值后，使用启发式计算智能算法优化进化种群ps；

G、更新迭代计数器g＝g+1，若g<G，则返回步骤E；否则本次优化结束，进入步骤H；

H、对于优化后进化种群ps中的每个寻优个体X_i，将其映射为选择矢量S_i；

I、构造对称的共表达权值矩阵W_k＝{w_p,q}_M×M，其中对角线元素w_p,p为所有S_i中代谢特征矢量F_p被选中的次数，p∈M：

$w_{p,p}=\underset{i\&Element;\left|ps\right|}{\&Sigma;}{s}_p\&Element;S_i$

其余元素w_p,q则为S_i中代谢特征矢量F_p与F_q被同时选中的次数,p,q∈M,p≠q：

w_p,q＝∑_i∈|ps|s_p∩s_q；s_p,s_q∈S_i；

J、更新运行计数器k＝k+1，若k<K，则返回步骤C，否则特征选择完成，进入步骤K；

K、对每次运行所获得的共表达权值矩阵进行平均，并计算其对应的概率，得到最终的共表达权值矩阵为Ω＝{ω_p,q}_M×M，其中|ps|为进化种群ps中的寻优个体总数：

${\mathrm{\&omega;}}_{p,q}=\frac{1}{K\left|ps\right|}{\mathrm{\&Sigma;}}_{k\&Element;K}{w}_{p,q}\&Element;W_k;$

L、将每次特征选择中最终输出的每个S_i视作优化算法对于代谢特征数据集空间的一次采样，其s_m∈S_i服从概率p_m的伯努利分布，则w_p,p即为服从B(|ps|,p_m)二次分布的随机变量；

M、将最终的共表达权值矩阵视作组合学习投票方法的稳态结果；

N、使用最终的共表达权值矩阵中的对角线元素ω_p,p作为节点p的重要性权重，其余任意ω_p,q,p≠q作为节点F_p与F_q之间的连接权重，构建全连通加权网络G，而后移除其中权重小于阈值ω_t的节点与边，形成该原始的代谢特征数据集F^*的代谢共表达网络；

O、输出所述代谢共表达网络作为结果。

所述的启发式代谢共表达网络的构建方法，其中，所述步骤E具体包括：

E1、设若输入个体为X_i＝{x_m；m＝1,2,...,M}，其各维度上均为R范围内的实数值，将其二值化为离散的选择矢量S_i＝{s_m；m＝1,2,...,M}：

$s_m=\left\{\begin{array}{cc}1,& \begin{array}{cc}if& x_m>0.5\end{array}\\ 0,& otherwise\end{array}\right.,s_m\&Element;S_i;$

E2、对于S_i中任意第m个选择值s_m，若其值为1，则选择对应代谢特征矢量F_m包含于所构造的特征子集F_S中，否则F_m将不会被选中；

F_S＝{F_m；m＝1,2,…,M,s_m＝1}；

E3、计算F_S中多变量的近似互信息值作为原始适应度函数值；

E4、定义稀疏适应度函数值为矢量X_i的1范数：

f_spr.(X_i)＝‖X_i‖₁；

E5、计算当前个体X_i的总体适应度函数值为：

f(X_i)＝f_raw(X_i)+λf_spr.(X_i)；

其中λ为拉格朗日乘数；

E6、若每个寻优个体的总体适应度函数值都已计算，则转至步骤E7，否则转至步骤E1；

E7、计算各寻优个体的共享适应度函数值：

$f_{share}\left(X_i\right)=f\left(X_i\right)(1+{\mathrm{\&Sigma;}}_{X_j\&Element;ps,\left|\right|X_i-X_j|{|}_2<r,j\&NotEqual;i}{\left(1-\frac{\left|\right|X_i-X_j|{|}_2}{r}\right)}^{\&Element;}),X_i\&Element;ps;$

其中r为聚集半径，∈为驱散因子。

所述的启发式代谢共表达网络的构建方法，其中，所述步骤E3具体包括：

E31、设若C为F的N个样本所对应的类标矢量，则F_S的互信息计算为：

I(F_S；C)＝H(F_S)-H(F_s|C)＝H(F_S)-∑_c∈Cp(c)H(F_s|c)；

其中p(c)表示类标c的出现概率，H()表示变量的熵；

E32、使用F_S中的N个样本作为节点，其相互欧氏距离作为边的权值，构造最小生成树，则L_γ(F_S)为此最小生成树的边权值总和：

$L_{\mathrm{\&gamma;}}\left(F_S\right)={\mathrm{\&Sigma;}}_{e_{i,j}\&Element;MST\left(F_S\right)}\left|\right|e_{i,j}|{|}^{\mathrm{\&gamma;}};$

其中γ为接近于0的正值常数；

E33、计算F_S的多变量互信息为：

I_appx.(F_S；C)＝L_γ(F_S)-∑_c∈Cp(c)L_γ(F_S|c)；

从而原始适应度函数值定义为：

f_raw(X_i)＝-I_appx.(F_S；C)。

一种启发式代谢共表达网络的构建系统，其中，包括：

正规化模块，用于对原始的代谢特征数据集F^*进行正规化预处理，使其中所有M个代谢特征矢量在每个维度上都具有0均值及单位方差:

$F_m=\frac{F_m^{*}-{\mathrm{\&mu;}}_m}{{\mathrm{\&delta;}}_m},F_m^{*}\&Element;F^{*};$

F＝{F_m；m＝1,2,…,M}为预处理后的代谢特征数据集，μ_m、δ_m分别为第m个原始的代谢特征矢量F^* _m的均值与方差；

运行计数器初始化模块，用于设定特征选择总运行次数为K，初始化运行计数器k＝1；

进化种群构造模块，用于构造多模优化的进化种群ps，将其所包含的每个寻优个体X_i∈ps初始化为范围R＝[0,1]内均匀分布的M维随机矢量；

迭代计数器初始化模块，用于设定算法迭代总次数为G，初始化迭代计数器g＝1；

适应度函数值计算模块，用于计算进化种群ps中每个寻优个体的共享适应度函数值；

种群优化模块，用于在计算所有寻优个体的共享适应度函数值后，使用启发式计算智能算法优化进化种群ps；

迭代计数器更新模块，用于更新迭代计数器g＝g+1，若g<G，则返回适应度函数值计算模块；否则本次优化结束，进入映射模块；

映射模块，用于对于优化后进化种群ps中的每个寻优个体X_i，将其映射为选择矢量S_i；

共表达权值矩阵构造模块，用于构造对称的共表达权值矩阵W_k＝{w_p,q}_M×M，其中对角线元素w_p,p为所有S_i中代谢特征矢量F_p被选中的次数，p∈M：

$w_{p,p}=\underset{i\&Element;\left|ps\right|}{\&Sigma;}{s}_p\&Element;S_i$

其余元素w_p,q则为S_i中代谢特征矢量F_p与F_q被同时选中的次数,p,q∈M,p≠q：

w_p,q＝∑_i∈|ps|s_p∩s_q；s_p,s_q∈S_i；

运行计数器更新模块，用于更新运行计数器k＝k+1，若k<K，则返回进化种群构造模块，否则特征选择完成，进入平均模块；

平均模块，用于对每次运行所获得的共表达权值矩阵进行平均，并计算其对应的概率，得到最终的共表达权值矩阵为Ω＝{ω_p,q}_M×M，其中|ps|为进化种群ps中的寻优个体总数：

${\mathrm{\&omega;}}_{p,q}=\frac{1}{K\left|ps\right|}{\mathrm{\&Sigma;}}_{k\&Element;K}{w}_{p,q}\&Element;W_k;$

采样模块，用于将每次特征选择中最终输出的每个S_i视作优化算法对于代谢数据集空间的一次采样，其s_m∈S_i服从概率p_m的伯努利分布，则w_p,p即为服从B(|ps|,p_m)二次分布的随机变量；

稳态结果输出模块，用于将最终的共表达权值矩阵视作组合学习投票方法的稳态结果；

代谢共表达网络计算模块，用于使用最终的共表达权值矩阵中的对角线元素ω_p,p作为节点p的重要性权重，其余任意ω_p,q,p≠q作为节点F_p与F_q之间的连接权重，构建全连通加权网络G，而后移除其中权重小于阈值ω_t的节点与边，形成该原始的代谢特征数据集F^*的代谢共表达网络；

代谢共表达网络输出模块，用于输出所述代谢共表达网络作为结果。

所述的启发式代谢共表达网络的构建系统，其中，所述适应度函数值计算模块具体包括：

二值化单元，用于设若输入个体为X_i＝{x_m；m＝1,2,...,M}，其各维度上均为R范围内的实数值，将其二值化为离散的选择矢量S_i＝{s_m；m＝1,2,...,M}：

$s_m=\left\{\begin{array}{cc}1,& \begin{array}{cc}if& x_m>0.5\end{array}\\ 0,& otherwise\end{array}\right.,s_m\&Element;S_i;$

选择单元，用于对于S_i中任意第m个选择值s_m，若其值为1，则选择对应代谢特征矢量F_m包含于所构造的特征子集F_S中，否则F_m将不会被选中；

F_S＝{F_m；m＝1,2,…,M,s_m＝1}；

原始适应度函数值计算单元，用于计算F_S中多变量的近似互信息值作为原始适应度函数值；

定义单元，用于定义稀疏适应度函数值为矢量X_i的1范数：

f_spr.(X_i)＝‖X_i‖₁；

总体适应度函数值计算单元，用于计算当前个体X_i的总体适应度函数值为：

f(X_i)＝f_raw(X_i)+λf_spr.(X_i)

其中λ为拉格朗日乘数；

判断单元，用于若每个寻优个体的总体适应度函数值都已计算，则转至共享适应度函数值计算单元，否则转至二值化单元；

共享适应度函数值计算单元，用于计算各寻优个体的共享适应度函数值：

$f_{share}\left(X_i\right)=f\left(X_i\right)(1+\underset{X_j\&Element;ps,\left|\right|X_i-X_j|{|}_2<r,j\&NotEqual;i}{\&Sigma;}{\left(1-\frac{\left|\right|X_i-X_j|{|}_2}{r}\right)}^{\&Element;}),X_i\&Element;ps$

其中r为聚集半径，∈为驱散因子。

所述的启发式代谢共表达网络的构建系统，其中，所述原始适应度函数值计算单元具体包括：

互信息计算子单元，用于设若C为F的N个样本所对应的类标矢量，则F_S的互信息计算为：

$I(F_S;C)=H\left(F_S\right)-H\left(F_s\right|C)=H\left(F_S\right)-\underset{c\&Element;C}{\&Sigma;}p\left(c\right)H\left(F_s\right|c)$

其中p(c)表示类标c的出现概率，H()表示变量的熵；

边权值计算子单元，用于使用F_S中的N个样本作为节点，其相互欧氏距离作为边的权值，构造最小生成树，则L_γ(F_S)为此最小生成树的边权值总和：

$L_{\mathrm{\&gamma;}}\left(F_S\right)=\underset{e_{i,j}\&Element;MST\left(F_S\right)}{\&Sigma;}\left|\right|e_{i,j}|{|}^{\mathrm{\&gamma;}}$

其中γ为接近于0的正值常数；

函数值计算子单元，用于计算F_S的多变量互信息为：

$I_{appx.}(F_S;C)=L_{\mathrm{\&gamma;}}\left(F_S\right)-\underset{c\&Element;C}{\&Sigma;}p\left(c\right)L_{\mathrm{\&gamma;}}\left(F_S\right|c)$

从而原始适应度函数值定义为：

f_raw(X_i)＝-I_appx.(F_S；C)。

有益效果：本发明根据最大依赖准则(Max-DependencyCriteria)将多个代谢物特征的互信息(MultivariateMutualInformation)作为适应度函数值，并使用启发式计算智能(ComputationalIntelligence)多模优化算法(MultimodalOptimization)对最优特征子集进行寻优搜索。通过多次运行这一优化过程，将各次结果进行组合学习，构造共表达网络结构。最后，通过概率模型计算切分阈值，获得准确、稳定的代谢共表达网络。

附图说明

图1为本发明一种启发式代谢共表达网络的构建方法较佳实施例的流程图。

图2为本发明中使用F_S中的样本点构造最小生成树的具体流程图。

图3为本发明中使用阈值切分构造代谢共表达网络的具体流程图。

具体实施方式

本发明提供一种启发式代谢共表达网络的构建方法及系统，为使本发明的目的、技术方案及效果更加清楚、明确，以下对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

请参阅图1，图1为本发明一种启发式代谢共表达网络的构建方法较佳实施例的流程图，如图所示，其包括步骤：

1)、对原始的代谢特征数据集F^*进行正规化预处理，使其中所有M个代谢特征矢量在每个维度上都具有0均值及单位方差:

$F_m=\frac{F_m^{*}-{\mathrm{\&mu;}}_m}{{\mathrm{\&delta;}}_m},F_m^{*}\&Element;F^{*}$

F＝{F_m；m＝1,2,…,M}为预处理后的代谢特征数据集，μ_m、δ_m分别为第m个原始的代谢特征矢量F^* _m的均值与方差；

2)、设定特征选择总运行次数为K，初始化运行计数器k＝1；

3)、构造多模优化的进化种群ps，将其所包含的每个寻优个体X_i∈ps初始化为范围R＝[0,1]内均匀分布的M维随机矢量；

4)、设定算法迭代总次数为G，初始化迭代计数器g＝1；

5)、计算进化种群ps中每个寻优个体的共享适应度函数值；

6)、在计算所有寻优个体的共享适应度函数值后，使用启发式计算智能算法优化进化种群ps；

7)、更新迭代计数器g＝g+1，若g<G，则返回5)；否则本次优化结束，进入步骤8)；

8)、对于优化后进化种群ps中的每个寻优个体X_i，将其映射为选择矢量S_i；

9)、构造对称的共表达权值矩阵W_k＝{w_p,q}_M×M，其中对角线元素w_p,p为所有S_i中代谢特征矢量F_p被选中的次数，p∈M：

$w_{p,p}=\underset{i\&Element;\left|ps\right|}{\&Sigma;}{s}_p\&Element;S_i$

其余元素w_p,q则为S_i中代谢特征矢量F_p与F_q被同时选中的次数,p,q∈M,p≠q：

w_p,q＝∑_i∈|ps|s_p∩s_q；s_p,s_q∈S_i；

10)、更新运行计数器k＝k+1，若k<K，则返回步骤3)，否则特征选择完成，进入步骤11)；

11)、对每次运行所获得的共表达权值矩阵进行平均，并计算其对应的概率，得到最终的共表达权值矩阵为Ω＝{ω_p,q}_M×M，其中|ps|为进化种群ps中的寻优个体总数：

${\mathrm{\&omega;}}_{p,q}=\frac{1}{K\left|ps\right|}{\mathrm{\&Sigma;}}_{k\&Element;K}{w}_{p,q}\&Element;W_k;$

12)、将每次特征选择中最终输出的每个S_i视作优化算法对于代谢数据集空间的一次采样，其s_m∈S_i服从概率p_m的伯努利分布，则w_p,p即为服从B(|ps|,p_m)二次分布的随机变量；

13)、将最终的共表达权值矩阵视作组合学习投票方法的稳态结果；

14)、使用最终的共表达权值矩阵中的对角线元素ω_p,p作为节点p的重要性权重，其余任意ω_p,q,p≠q作为节点F_p与F_q之间的连接权重，构建全连通加权网络G，而后移除其中权重小于阈值ω_t的节点与边，形成该原始的代谢特征数据集F^*的代谢共表达网络；

15)、输出所述代谢共表达网络作为结果。

具体来说，在步骤1)中，在进行特征选择前，对原始的代谢特征数据集F^*进行正规化预处理(Standardization)，使其所有M个代谢特征矢量F^* _m在每个维度上都具有0均值及单位方差:

$F_m=\frac{F_m^{*}-{\mathrm{\&mu;}}_m}{{\mathrm{\&delta;}}_m},F_m^{*}\&Element;F^{*}$

其中F＝{F_m；m＝1,2,…,M}为预处理后的代谢特征数据集，μ_m、δ_m分别为第m个原始的代谢特征矢量F^* _m的均值与方差。

在步骤2)中，设定特征选择总运行次数为K，初始化运行计数器k＝1；

在步骤3)中，构造多模优化的进化种群ps，所包含的每个寻优个体Xi∈ps初始化为范围R＝[0,1]内均匀分布的M维随机矢量；

在步骤4)中，开始特征选择子集的优化设计。设定算法迭代总次数为G，初始化迭代计数器g＝1。

在步骤5)中，计算进化种群ps中每个寻优个体的共享适应度函数值。

所述步骤5)具体包括：

a、设若输入个体(即输入的寻优个体)为X_i＝{x_m；m＝1,2,...,M}，其各维度上均为R范围内的实数值，将其二值化为离散的选择矢量S_i＝{s_m；m＝1,2,...,M}：oherwise即指除x_m＞0.5之外的其他情况。

$s_m=\left\{\begin{array}{cc}1,& \begin{array}{cc}if& x_m>0.5\end{array}\\ 0,& otherwise\end{array}\right.,s_m\&Element;S_i;$

b、对于S_i中任意第m个选择值s_m，若其值为1，则选择对应代谢特征矢量F_m包含于所构造的特征子集F_S中，否则F_m将不会被选中；

F_S＝{F_m；m＝1,2,…,M,s_m＝1}；

c、计算F_S中多变量的近似互信息值作为原始适应度函数值；

d、定义稀疏适应度函数值为矢量X_i的1范数：

f_spr.(X_i)＝‖X_i‖₁；

其引入可使算法选择最为重要的核心代谢物特征。

e、计算当前个体X_i的总体适应度函数值为：

f(X_i)＝f_raw(X_i)+λf_spr.(X_i)

其中λ为拉格朗日乘数；

f、若每个寻优个体的总体适应度函数值都已计算，则转至步骤5).g)，否则转至步骤5).a)；

g、使用FitnessSharing(适应度函数值共享)方法计算各寻优个体的共享适应度函数值：

$f_{share}\left(X_i\right)=f\left(X_i\right)(1+\underset{X_j\&Element;ps,\left|\right|X_i-X_j|{|}_2<r,j\&NotEqual;i}{\&Sigma;}{\left(1-\frac{\left|\right|X_i-X_j|{|}_2}{r}\right)}^{\&Element;}),X_i\&Element;ps$

其中r为聚集半径，∈为驱散因子。此方法可使搜索算法进行多模优化，获得特征空间中所有的全局或局部最优解(亦即特征选择子集)。

所述步骤c具体包括：

i、设若C为F的N个样本所对应的类标矢量，则F_S的互信息计算为：

$I(F_S;C)=H\left(F_S\right)-H\left(F_s\right|C)=H\left(F_S\right)-\underset{c\&Element;C}{\&Sigma;}p\left(c\right)H\left(F_S\right|c)$

其中p(c)表示类标c的出现概率，其值可根据数据集样本进行估算，H()表示变量的熵；可使用Rényi’sα-Entropy取得：

$H\left(F_S\right)=\frac{1}{1-\mathrm{\&alpha;}}\&lsqb;log\frac{L_{\mathrm{\&gamma;}}\left(F_S\right)}{N^{\mathrm{\&alpha;}}}-log\mathrm{\&beta;}\&rsqb;$

其中α为趋近于1的常数，β为与概率分布无关的偏差修正值，因此有：

H(F_S)∝L_γ(F_S)

呈正相关关系。

ii、使用F_S中的N个样本作为节点(Vertices)，其相互欧氏距离作为边(Edges)的权值，构造最小生成树(MinimumSpanningTree,MST)，则L_γ(F_S)为此MST的边权值总和：

$L_{\mathrm{\&gamma;}}\left(F_S\right)=\underset{e_{i,j}\&Element;MST\left(F_S\right)}{\&Sigma;}\left|\right|e_{i,j}|{|}^{\mathrm{\&gamma;}}$

其中γ为接近于0的正值常数；常用的MST构造算法包括Prim算法等。

如图2所示，F_S＝{pt₁＝(9,3),pt₂＝(3,5),pt₃＝(7,7),pt₄＝(5,10),pt₅＝(10,12)}，F_S由5个样本点构成，则其MST有：

e_1,3＝‖pt₁-pt₃‖＝4.47；

e_2,3＝‖pt₂-pt₃‖＝4.47；

e_3,5＝‖pt₃-pt₅‖＝4.47；

e_3,4＝‖pt₃-pt₄‖＝4.47；

L₁(F_S)＝4.47+4.47+5.83+3.60＝18.37。

iii、计算F_S的多变量互信息为：

$I_{appx.}(F_S;C)=L_{\mathrm{\&gamma;}}\left(F_S\right)-\underset{c\&Element;C}{\&Sigma;}p\left(c\right)L_{\mathrm{\&gamma;}}\left(F_S\right|C)$

其值越大则表示该代谢特征子集与目标生理状态的关联更为显著，从而原始适应度函数值定义为：

f_raw(X_i)＝-I_appx.(F_S；C)。

在步骤6)中，在计算所有寻优个体的共享适应度函数值后，使用启发式计算智能算法优化进化种群ps；常用方法为差分进化(DifferentialEvolution,DE)与Memetic算法(MemeticAlgorithm,MA)等。

在步骤7)中，更新迭代计数器g＝g+1。若g<G，则返回5)。否则本次优化结束，转至8)。

在步骤8)中，对于优化后ps中的每个寻优个体Xi，使用5).a)中的方法映射为选择矢量S_i。

在步骤9)中，构造对称的共表达权值矩阵W_k＝{w_p,q}_M×M，其中对角线元素w_p,p,p∈M为所有S_i中代谢特征矢量F_p被选中的次数：

$w_{p,p}=\underset{i\&Element;\left|ps\right|}{\&Sigma;}{s}_p\&Element;S_i$

其余元素w_p,q,p,q∈M,p≠q则为S_i中特征F_p与F_q被同时选中的次数：

$w_{p,q}=\underset{i\&Element;\left|ps\right|}{\&Sigma;}{s}_P\&cap;s_q;s_p,s_q\&Element;S_i$

在步骤10)中，更新运行计数器k＝k+1。若k<K，则返回3)。否则特征选择完成，转至9)。

在步骤11)中，对每次运行所获得的共表达权值矩阵W_k进行平均，并计算其对应的概率，得到最终的共表达权值矩阵为Ω＝{ω_p,q}_M×M，其中|ps|为进化种群中的寻优个体总数：

${\mathrm{\&omega;}}_{p,q}=\frac{1}{K\left|ps\right|}\underset{k\&Element;K}{\&Sigma;}{w}_{p,q}\&Element;W_k$

在步骤12)中，设若在每次特征选择中，将最终输出的每个S_i视作优化算法对于代谢数据集空间的一次采样，其s_m∈S_i服从概率p_m的伯努利分布。则w_p,p即为服从B(|ps|,p_m)二次分布的随机变量。则在种群规模|ps|设置为:

的条件下，其可视作服从均值为μ＝|ps|p_m，方差为σ＝|ps|p_m(1-p_m)的正态分布N(μ,σ)。从而总运行次数K可通过以下公式求得：

其中z^*为置信值，ε为均值的最大误差范围。

例如，假设F_m的选择概率p_m∈[0.05,0.95]，则在每次特征选择使用|ps|＝100个寻优个体的条件下，重复运行K＝6次，即可在98％的置信区间内(z^*＝2.33)，保证ω_p,p的均值误差不超过ε＝5％。

在步骤13)中，在该置信值下，可将最终的共表达权值矩阵Ω视作组合学习BootstrapAggregating投票方法的稳态结果，例如可设置切分阈值为ω_t＝0.5。

在步骤14)中，如图3所示，使用最终的共表达权值矩阵Ω中的对角线元素ω_p,p作为节点p(代谢物特征F_p)的重要性权重，其余任意ω_p,q,p≠q作为节点F_p、F_q间的连接权重，构建全连通加权网络G。而后移除其中权重小于阈值ω_t的节点与边，形成该代谢特征数据集F^*的代谢共表达网络。

在步骤15)中，输出所述代谢共表达网络作为结果。

基于上述方法，本发明还提供一种启发式代谢共表达网络的构建系统，其包括：

正规化模块，用于对原始的代谢特征数据集F^*进行正规化预处理，使其中所有M个代谢特征矢量在每个维度上都具有0均值及单位方差:

$F_m=\frac{F_m^{*}-{\mathrm{\&mu;}}_m}{{\mathrm{\&delta;}}_m},F_m^{*}\&Element;F^{*};$

F＝{F_m；m＝1,2,…,M}为预处理后的代谢特征数据集，μ_m、δ_m分别为第m个原始的代谢特征矢量F^* _m的均值与方差；

运行计数器初始化模块，用于设定特征选择总运行次数为K，初始化运行计数器k＝1；

进化种群构造模块，用于构造多模优化的进化种群ps，将其所包含的每个寻优个体X_i∈ps初始化为范围R＝[0,1]内均匀分布的M维随机矢量；

迭代计数器初始化模块，用于设定算法迭代总次数为G，初始化迭代计数器g＝1；

适应度函数值计算模块，用于计算进化种群ps中每个寻优个体的共享适应度函数值；

种群优化模块，用于在计算所有寻优个体的共享适应度函数值后，使用启发式计算智能算法优化进化种群ps；

迭代计数器更新模块，用于更新迭代计数器g＝g+1，若g<G，则返回适应度函数值计算模块；否则本次优化结束，进入映射模块；

映射模块，用于对于优化后进化种群ps中的每个寻优个体X_i，将其映射为选择矢量S_i；

共表达权值矩阵构造模块，用于构造对称的共表达权值矩阵W_k＝{w_p,q}_M×M，其中对角线元素w_p,p为所有S_i中代谢特征矢量F_p被选中的次数，p∈M：

$w_{p,p}=\underset{i\&Element;\left|ps\right|}{\&Sigma;}{s}_p\&Element;S_i$

其余元素w_p,q则为S_i中代谢特征矢量F_p与F_q被同时选中的次数,p,q∈M,p≠q：

w_p,q＝∑_i∈|ps|s_p∩s_q；s_p,s_q∈S_i；

运行计数器更新模块，用于更新运行计数器k＝k+1，若k<K，则返回进化种群构造模块，否则特征选择完成，进入平均模块；

平均模块，用于对每次运行所获得的共表达权值矩阵进行平均，并计算其对应的概率，得到最终的共表达权值矩阵为Ω＝{ω_p,q}_M×M，其中|ps|为进化种群ps中的寻优个体总数：

${\mathrm{\&omega;}}_{p,q}=\frac{1}{K\left|ps\right|}{\mathrm{\&Sigma;}}_{k\&Element;K}{w}_{p,q}\&Element;W_k;$

采样模块，用于将每次特征选择中最终输出的每个S_i视作优化算法对于代谢数据集空间的一次采样，其s_m∈S_i服从概率p_m的伯努利分布，则w_p,p即为服从B(|ps|,p_m)二次分布的随机变量；

稳态结果输出模块，用于将最终的共表达权值矩阵视作组合学习投票方法的稳态结果；

代谢共表达网络计算模块，用于使用最终的共表达权值矩阵中的对角线元素ω_p,p作为节点p的重要性权重，其余任意ω_p,q,p≠q作为节点F_p与F_q之间的连接权重，构建全连通加权网络G，而后移除其中权重小于阈值ω_t的节点与边，形成该原始的代谢特征数据集F^*的代谢共表达网络；

代谢共表达网络输出模块，用于输出所述代谢共表达网络作为结果。

其中，所述适应度函数值计算模块具体包括：

二值化单元，用于设若输入个体为X_i＝{x_m；m＝1,2,...,M}，其各维度上均为R范围内的实数值，将其二值化为离散的选择矢量S_i＝{s_m；m＝1,2,...,M}：

$s_m=\left\{\begin{array}{cc}1,& \begin{array}{cc}if& x_m>0.5\end{array}\\ 0,& otherwise\end{array}\right.,s_m\&Element;S_i;$

选择单元，用于对于S_i中任意第m个选择值s_m，若其值为1，则选择对应代谢特征矢量F_m包含于所构造的特征子集F_S中，否则F_m将不会被选中；

F_S＝{F_m；m＝1,2,…,M,s_m＝1}；

原始适应度函数值计算单元，用于计算F_S中多变量的近似互信息值作为原始适应度函数值；

定义单元，用于定义稀疏适应度函数值为矢量X_i的1范数：

f_spr.(X_i)＝‖X_i‖₁；

总体适应度函数值计算单元，用于计算当前个体X_i的总体适应度函数值为：

f(X_i)＝f_raw(X_i)+λf_spr.(X_i)

其中λ为拉格朗日乘数；

判断单元，用于若每个寻优个体的总体适应度函数值都已计算，则转至共享适应度函数值计算单元，否则转至二值化单元；

共享适应度函数值计算单元，用于计算各寻优个体的共享适应度函数值：

$f_{share}\left(X_i\right)=f\left(X_i\right)(1+\underset{X_j\&Element;ps,\left|\right|X_i-X_j|{|}_2<r,j\&NotEqual;i}{\&Sigma;}{\left(1-\frac{\left|\right|X_i-X_j|{|}_2}{r}\right)}^{\&Element;}),X_i\&Element;ps$

其中r为聚集半径，∈为驱散因子。

其中，所述原始适应度函数值计算单元具体包括：

互信息计算子单元，用于设若C为F的N个样本所对应的类标矢量，则F_S的互信息计算为：

$I(F_S;C)=H\left(F_S\right)-H\left(F_s\right|C)=H\left(F_S\right)-\underset{c\&Element;C}{\&Sigma;}p\left(c\right)H\left(F_s\right|c)$

其中p(c)表示类标c的出现概率，H()表示变量的熵；

边权值计算子单元，用于使用F_S中的N个样本作为节点，其相互欧氏距离作为边的权值，构造最小生成树，则L_γ(F_S)为此最小生成树的边权值总和：

$L_{\mathrm{\&gamma;}}\left(F_S\right)=\underset{e_{i,j}\&Element;MST\left(F_S\right)}{\&Sigma;}\left|\right|e_{i,j}|{|}^{\mathrm{\&gamma;}}$

其中γ为接近于0的正值常数；

函数值计算子单元，用于计算F_S的多变量互信息为：

$I_{appx.}(F_S;C)=L_{\mathrm{\&gamma;}}\left(F_S\right)-\underset{c\&Element;C}{\&Sigma;}p\left(c\right)L_{\mathrm{\&gamma;}}\left(F_S\right|c)$

从而原始适应度函数值定义为：

f_raw(X_i)＝-I_appx.(F_S；C)。

应当理解的是，本发明的应用不限于上述的举例，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (6)

1.一种启发式代谢共表达网络的构建方法，其特征在于，包括步骤：

A、对原始的代谢特征数据集F^*进行正规化预处理，使其中所有M个代谢特征矢量在每个维度上都具有0均值及单位方差:

$F_m=\frac{F_m^{*}-{\mathrm{\&mu;}}_m}{{\mathrm{\&delta;}}_m},F_m^{*}\&Element;F^{*};$

F＝{F_m；m＝1,2,…,M}为预处理后的代谢特征数据集，μ_m、δ_m分别为第m个原始的代谢特征矢量F^* _m的均值与方差；

B、设定特征选择总运行次数为K，初始化运行计数器k＝1；

C、构造多模优化的进化种群ps，将其所包含的每个寻优个体X_i∈ps初始化为范围R＝[0,1]内均匀分布的M维随机矢量；

D、设定算法迭代总次数为G，初始化迭代计数器g＝1；

E、计算进化种群ps中每个寻优个体的共享适应度函数值；

F、在计算所有寻优个体的共享适应度函数值后，使用启发式计算智能算法优化进化种群ps；

G、更新迭代计数器g＝g+1，若g<G，则返回步骤E；否则本次优化结束，进入步骤H；

H、对于优化后进化种群ps中的每个寻优个体X_i，将其映射为选择矢量S_i；

I、构造对称的共表达权值矩阵W_k＝{w_p,q}_M×M，其中对角线元素w_p,p为所有S_i中代谢特征矢量F_p被选中的次数，p∈M：

$w_{p,p}=\underset{i\&Element;\left|ps\right|}{\&Sigma;}{s}_p\&Element;S_i$

其余元素w_p,q则为S_i中代谢特征矢量F_p与F_q被同时选中的次数,p,q∈M,p≠q：

w_p，q＝Σ_i∈|ps|s_p∩s_q；s_p，s_q∈S_i；

J、更新运行计数器k＝k+1，若k<K，则返回步骤C，否则特征选择完成，进入步骤K；

K、对每次运行所获得的共表达权值矩阵进行平均，并计算其对应的概率，得到最终的共表达权值矩阵为Ω＝{ω_p,q}_M×M，其中|ps|为进化种群ps中的寻优个体总数：

${\mathrm{\&omega;}}_{p,q}=\frac{1}{K\left|ps\right|}{\mathrm{\&Sigma;}}_{k\&Element;K}{w}_{p,q}\&Element;W_k;$

L、将每次特征选择中最终输出的每个S_i视作优化算法对于代谢特征数据集空间的一次采样，其s_m∈S_i服从概率p_m的伯努利分布，则w_p,p即为服从B(|ps|,p_m)二次分布的随机变量；

M、将最终的共表达权值矩阵视作组合学习投票方法的稳态结果；

N、使用最终的共表达权值矩阵中的对角线元素ω_p,p作为节点p的重要性权重，其余任意ω_p,q,p≠q作为节点F_p与F_q之间的连接权重，构建全连通加权网络G，而后移除其中权重小于阈值ω_t的节点与边，形成该原始的代谢特征数据集F^*的代谢共表达网络；

O、输出所述代谢共表达网络作为结果。

2.根据权利要求1所述的启发式代谢共表达网络的构建方法，其特征在于，所述步骤E具体包括：

E1、设若输入个体为X_i＝{x_m；m＝1,2,...,M}，其各维度上均为R范围内的实数值，将其二值化为离散的选择矢量S_i＝{s_m；m＝1,2,...,M}：

$s_m=\left\{\begin{array}{cc}1,& if{x}_m>0.5\\ 0,& otherwise\end{array}\right.,s_m\&Element;S_i;$

E2、对于S_i中任意第m个选择值s_m，若其值为1，则选择对应代谢特征矢量F_m包含于所构造的特征子集F_S中，否则F_m将不会被选中；

F_s＝{F_m；m＝1,2,…,M,s_m＝1}；

E3、计算F_S中多变量的近似互信息值作为原始适应度函数值；

E4、定义稀疏适应度函数值为矢量X_i的1范数：

f_spr.(X_i)＝‖X_i‖₁；

E5、计算当前个体X_i的总体适应度函数值为：

f(X_i)＝f_raw(X_i)+λf_spr.(X_i)；

其中λ为拉格朗日乘数；

E6、若每个寻优个体的总体适应度函数值都已计算，则转至步骤E7，否则转至步骤E1；

E7、计算各寻优个体的共享适应度函数值：

$f_{share}\left(X_i\right)=f\left(X_i\right)(1+{\mathrm{\&Sigma;}}_{X_j\&Element;ps,\left|\right|X_i-X_j|{|}_2<r,j\&NotEqual;i}{\left(1-\frac{\left|\right|X_i-X_j|{|}_2}{r}\right)}^{\&Element;}),X_i\&Element;ps;$

其中r为聚集半径，ε为驱散因子。

3.根据权利要求2所述的启发式代谢共表达网络的构建方法，其特征在于，所述步骤E3具体包括：

E31、设若C为F的N个样本所对应的类标矢量，则F_S的互信息计算为：

I(F_s；C)＝H(F_s)-H(F_s|C)＝H(F_s)-∑_c∈cp(c)H(F_s|c)；

其中p(c)表示类标c的出现概率，H()表示变量的熵；

E32、使用F_S中的N个样本作为节点，其相互欧氏距离作为边的权值，构造最小生成树，则L_γ(F_S)为此最小生成树的边权值总和：

L_γ(F_s)＝∑_{ei,k∈MST(Fs)}||e_i,j||^I；

其中γ为接近于0的正值常数；

E33、计算F_S的多变量互信息为：

I_appx.(F_s；C)＝L_γ(F_s)-∑_c∈cp(c)L_γ(F_s|c)；

从而原始适应度函数值定义为：

f_raw(X_i)＝-I_appX.(F_s；C)。

4.一种启发式代谢共表达网络的构建系统，其特征在于，包括：

正规化模块，用于对原始的代谢特征数据集F^*进行正规化预处理，使其中所有M个代谢特征矢量在每个维度上都具有0均值及单位方差:

$F_m=\frac{F_m^{*}-{\mathrm{\&mu;}}_m}{{\mathrm{\&delta;}}_m},F_m^{*}\&Element;F^{*};$

F＝{F_m；m＝1,2,…,M}为预处理后的代谢特征数据集，μ_m、w_m分别为第m个原始的代谢特征矢量F^* _m的均值与方差；

运行计数器初始化模块，用于设定特征选择总运行次数为K，初始化运行计数器k＝1；

进化种群构造模块，用于构造多模优化的进化种群ps，将其所包含的每个寻优个体X_i∈ps初始化为范围R＝[0,1]内均匀分布的M维随机矢量；

迭代计数器初始化模块，用于设定算法迭代总次数为G，初始化迭代计数器g＝1；

适应度函数值计算模块，用于计算进化种群ps中每个寻优个体的共享适应度函数值；

种群优化模块，用于在计算所有寻优个体的共享适应度函数值后，使用启发式计算智能算法优化进化种群ps；

迭代计数器更新模块，用于更新迭代计数器g＝g+1，若g<G，则返回适应度函数值计算模块；否则本次优化结束，进入映射模块；

映射模块，用于对于优化后进化种群ps中的每个寻优个体X_i，将其映射为选择矢量S_i；

共表达权值矩阵构造模块，用于构造对称的共表达权值矩阵W_k＝{w_p,q}_M×M，其中对角线元素w_p,p为所有S_i中代谢特征矢量F_p被选中的次数，p∈M：

$w_{p,p}=\underset{i\&Element;\left|ps\right|}{\&Sigma;}{s}_p\&Element;S_i$

其余元素w_p,q则为S_i中代谢特征矢量F_p与F_q被同时选中的次数,p,q∈M,p≠q：

w_p,q＝∑_i∈|ps|s_p∩s_q；s_p,s_q∈s_i；

运行计数器更新模块，用于更新运行计数器k＝k+1，若k<K，则返回进化种群构造模块，否则特征选择完成，进入平均模块；

平均模块，用于对每次运行所获得的共表达权值矩阵进行平均，并计算其对应的概率，得到最终的共表达权值矩阵为Ω＝{ω_p,q}_M×M，其中|ps|为进化种群ps中的寻优个体总数：

${\mathrm{\&omega;}}_{p,q}=\frac{1}{K\left|ps\right|}{\mathrm{\&Sigma;}}_{k\&Element;K}{w}_{p,q}\&Element;W_k;$

采样模块，用于将每次特征选择中最终输出的每个S_i视作优化算法对于代谢数据集空间的一次采样，其s_m∈S_i服从概率p_m的伯努利分布，则w_p,p即为服从B(|ps|,p_m)二次分布的随机变量；

稳态结果输出模块，用于将最终的共表达权值矩阵视作组合学习投票方法的稳态结果；

代谢共表达网络计算模块，用于使用最终的共表达权值矩阵中的对角线元素ω_p,p作为节点p的重要性权重，其余任意ω_p,q,p≠q作为节点F_p与F_q之间的连接权重，构建全连通加权网络G，而后移除其中权重小于阈值ω_t的节点与边，形成该原始的代谢特征数据集F^*的代谢共表达网络；

代谢共表达网络输出模块，用于输出所述代谢共表达网络作为结果。

5.根据权利要求4所述的启发式代谢共表达网络的构建系统，其特征在于，所述适应度函数值计算模块具体包括：

二值化单元，用于设若输入个体为X_i＝{x_m；m＝1,2,...,M}，其各维度上均为R范围内的实数值，将其二值化为离散的选择矢量S_i＝{s_m；m＝1,2,...,M}：

$s_m=\left\{\begin{array}{cc}1,& if{x}_m>0.5\\ 0,& otherwise\end{array}\right.,s_m\&Element;S_i;$

选择单元，用于对于S_i中任意第m个选择值s_m，若其值为1，则选择对应代谢特征矢量F_m包含于所构造的特征子集F_S中，否则F_m将不会被选中；

F_s＝{F_m；m＝1,2,…,M,s_m＝1}；

原始适应度函数值计算单元，用于计算F_S中多变量的近似互信息值作为原始适应度函数值；

定义单元，用于定义稀疏适应度函数值为矢量X_i的1范数：

f_spr.(X_i)＝‖X_i‖₁；

总体适应度函数值计算单元，用于计算当前个体X_i的总体适应度函数值为：

f(X_i)＝f_raw(X_i)+λf_spr.(X_i)

其中λ为拉格朗日乘数；

判断单元，用于若每个寻优个体的总体适应度函数值都已计算，则转至共享适应度函数值计算单元，否则转至二值化单元；

共享适应度函数值计算单元，用于计算各寻优个体的共享适应度函数值：

$f_{share}\left(X_i\right)=f\left(X_i\right)(1+\underset{X_j\&Element;ps,\left|\right|X_i-X_j|{|}_2<r,j\&NotEqual;i}{\mathrm{\&Sigma;}}{\left(1-\frac{\left|\right|X_i-X_j|{|}_2}{r}\right)}^{\&Element;}),X_i\&Element;ps$

其中r为聚集半径，ε为驱散因子。

6.根据权利要求5所述的启发式代谢共表达网络的构建系统，其特征在于，所述原始适应度函数值计算单元具体包括：

互信息计算子单元，用于设若C为F的N个样本所对应的类标矢量，则F_S的互信息计算为：

$I(F_S;C)=H\left(F_S\right)-H\left(F_S\right|C)=H\left(F_S\right)-\underset{c\&Element;C}{\&Sigma;}P\left(c\right)H\left(F_S\right|c)$

其中p(c)表示类标c的出现概率，H()表示变量的熵；

边权值计算子单元，用于使用F_S中的N个样本作为节点，其相互欧氏距离作为边的权值，构造最小生成树，则L_γ(F_S)为此最小生成树的边权值总和：

$L_{\mathrm{\&gamma;}}\left(F_S\right)=\underset{e_{i,j}\&Element;MST\left(F_S\right)}{\&Sigma;}\left|\right|e_{i,j}|{|}^{\mathrm{\&gamma;}}$

其中γ为接近于0的正值常数；

函数值计算子单元，用于计算F_S的多变量互信息为：

$I_{appx.}(F_S;C)=L_{\mathrm{\&gamma;}}\left(F_S\right)-\underset{c\&Element;C}{\&Sigma;}p\left(c\right)L_{\mathrm{\&gamma;}}\left(F_S\right|c)$

从而原始适应度函数值定义为：

f_raw(X_i)＝-I_appx.(F_s；C)。