阅读说明：本技术 基于多网络的生物标志物识别方法 (Multi-network-based biomarker identification method ) 是由 李幸一 李敏 项炬 王建新 于 2019-09-18 设计创作，主要内容包括：本发明公开了一种基于多网络的生物标志物识别方法；考虑到样本异质性的影响,本发明方法首先对基因表达谱数据进行标准化处理并对样本进行主成分分析,利用前两个主成分对样本通过高斯混合模型聚类；对于每一类样本,构建一个基于多网络的网络传播模型对网络中的所有基因进行排序,用于初步筛选重要基因；为了获得具有最大区分能力、最小冗余的生物标志物,通过接收者操作特征曲线的线下面积(AUC)优化模型在上一步得到的重要特征中进一步对基因打分排序,获得生物标志物。本发明的方法充分利用了多源生物网络信息,可以有效识别具有最大分类能力、最小冗余和生物可解释性的生物标志物,用于异质的复杂疾病分析。(The invention discloses a biomarker identification method based on multiple networks; considering the influence of sample heterogeneity, the method firstly carries out standardization processing on gene expression profile data and carries out principal component analysis on the sample, and the first two principal components are utilized to cluster the sample through a Gaussian mixture model; for each type of samples, constructing a network propagation model based on multiple networks to sequence all genes in the network for primarily screening important genes; in order to obtain biomarkers with maximum discriminatory power and minimum redundancy, the biomarkers were obtained by further ranking the genes in the important features obtained in the previous step by the receiver operating the area under the line (AUC) optimization model of the characteristic curve. The method of the invention fully utilizes the multi-source biological network information, can effectively identify the biomarker with maximum classification capability, minimum redundancy and biological interpretability, and is used for heterogeneous complex disease analysis.)

基于多网络的生物标志物识别方法

技术领域

本发明涉及生物信息学领域，特别是一种基于多网络的生物标志物识别方法。

背景技术

复杂疾病异质性强，易受环境因素影响，这给复杂疾病的诊断和治疗带来困难。因此，异质复杂疾病分析已成为现代医学研究的重点之一。生物标志物是客观测量和评价正常生物过程、病理过程或药物干预反应的指标，也是机体受到损伤时的重要预警指标。从多种生物数据中挖掘有效的生物标志物是解决复杂疾病的关键。

随着系统生物学的深入研究和高通量技术的迅速发展，大量的生物相互作用网络被获得，如蛋白质相互作用网络、代谢网络、调控网络、信号网络、蛋白质复合物等。网络信息反映了生命系统各组成元件之间复杂的相互作用。基于生物网络进行复杂疾病分析可以帮助研究人员了解疾病发生的规律并从网络层面揭示特定的疾病发病过程。近年来，网络传播模型已成功地用于识别与特定疾病相关的生物标志物。例如，Sun等人^[1]提出一个改进的随机游走的MarkRank算法，用于在蛋白质相互作用网络上挖掘协同合作的网络生物标志物。Winter等人^[2]提出一个与佩奇排序相似的NetRank算法从而提高了生物标志物预测的准确性和可再生性。然而，网络传播模型很大程度上取决于网络数据的可靠性和完整性，基于单一生物网络可能由于网络的不完全性导致结果不准确。

另外一些方法将网络传播模型应用于多信息网络。例如，Choi等^[3]利用改进的佩奇排序从包含通路、蛋白质相互作用网络、基因共表达网络和基因本体等信息的功能相互作用网络中的异质生物标志物。但是，直接将多种不同生物相互作用信息整合为一个网络，可能存在不同类型的数据相互干扰、无法很好地反映每种类型网络自身的特征和拓扑结构等问题。

基于目前研究的不足，有必要提供一种基于多网络的生物标志物识别方法。

发明内容

本发明所要解决的技术问题是，针对现有技术不足，提供一种基于多网络的生物标志物识别方法，以提高疾病样本分类的准确性。

为解决上述技术问题，本发明所采用的技术方案是：一种基于多网络的生物标志物识别方法，包括以下步骤：

1)对基因表达谱数据的每一个样本进行Z分数标准化，并进行主成分分析，取前两个主成分用于后续分析；

2)根据前两个主成分，对样本通过高斯混合模型聚类；

3)对于每一类，我们构建了一个基于多网络的网络传播模型对网络中的所有基因进行排序，用于初步筛选重要基因；

4)通过AUC优化模型在上一步得到的重要特征中进一步对基因打分排序，获得生物标志物。

通过Z分数对基因表达数据进行标准化的表达式为：x表示每一个样本的基因原始表达值；μ表示每一个样本的所有基因原始表达数据的均值；σ为每一个样本的所有基因原始表达数据的标准差。

步骤1)的具体实现过程包括：

1)求标准化后的基因表达数据中特征的协方差矩阵；

2)求协方差矩阵的特征值及对应的特征向量；

3)将特征值按照从大到小的顺序排序，选择其中最大的两个，然后将其对应的两个特征向量分别作为列向量，组成特征向量矩阵；

4)将样本点投影到选取的特征向量上，设标准化后的基因表达数据是样本数为m，基因数为n的矩阵a，协方差矩阵为n*n维，选取的两个特征向量组成的矩阵为b，那么投影后的数据为：Y＝a×b，从而将原始样本的n维特征变成了两维。

步骤2)中，对于每个样本x，高斯混合模型的表达式为：

其中k表示聚类的个数，

表示高斯混合模型中的第k个分量的高斯分布；π_k是混合系数。

步骤3)中，首先对节点通过t统计量打分，作为节点的初始状态，表达式为：

为两类样本的平均数，

为两类样本的方差，n₁、n₂为两类样本的容量。一个网络可以表示为一个邻接矩阵为A的图：G＝(V,E)，其中，A_ij表示节点v_i和v_j之间存在相互作用边。V＝{v₁,v₂,...,v_n}为网络中的节点，E＝{(v_i,v_j),i,j＝1,...,n,A_ij≠0}表示网络中的边。对于一个由S个网络构成的多网络，其邻接矩阵可以表示为：

其中，I表示n×n的单位矩阵，A_s表示第s层网络的邻接矩阵，μ可以量化节点在层内或在层间的概率。基于多网络的网络传播模型的表达式为：F_t+1＝αW'F_t+(1-α)T，其中，W'为多网络的邻接矩阵经过行标准化后的矩阵，α可以量化初始状态的相对重要性。通过基于多网络的网络传播模型对多网络中的节点打分排序后，每个节点可以得到S个重要性度量值，通过几何平均求得节点的打分值，根据该打分值计算该节点在该簇中的排序值。最后通过计算每个节点在所有样本簇的平均排序值，作为该节点最终的排序值。

步骤4)中，AUC优化模型的表达式可以通过一个线性规划方程来表示：

该线性规划方程服从于：

其中X₁,X₂,...,X_p为正样本，Y₁,Y₂,...,Y_q为负样本，每个样本都由一个k维特征向量表示。ξ_ij+、ξ_ij-分别为σ^T(X_i-Y_j)的正、负部分。σ是一个k维向量，表示为：σ＝{σ_k,k＝1,..,K}；θ是一个K维向量，表示这些K个特征被选为标志物的概率；当给出正负样本的中心指标，即

和

ω可以表示为：ω＝β₊-β_-。λ设为0.00001。

与现有技术相比，本发明所具有的有益效果为：本发明为了减小样本异质性的影响，首先对样本的基因表达数据进行标准化和主成分分析，取前两个主成分、用高斯混合模型对样本进行聚类；对于每一类，我们构建了一个基于多网络的网络传播模型对网络中的所有基因进行排序，用于初步筛选重要基因；为了获得具有最大区分能力、最小冗余的生物标志物，我们通过AUC优化模型在上一步得到的重要特征中进一步对基因打分排序，获得生物标志物。本发明的方法可以有效识别具有最大分类能力和生物可解释性的生物标志物，用于异质的复杂疾病分析。

附图说明

图1本发明方法的流程图；

图2为本发明方法和两个生物标志物识别方法：NetRank和MarkRank,以及三个基准方法：带重启的随机游走(RWR),标签传播(LP),主成分分析(PCA)的AUC值在六个数据集的比较图；(1)BRCA(2)NKI(3)GSE3365(4)GSE36807(5)GSE35713(6)GSE9006；

图3为本发明方法和本发明方法基于单网络的比较图。

图4为本发明方法和本发明方法基于融合网络的比较图。

图5为本发明方法和NetRank,MarkRank方法的已知致病基因和差异表达基因富集性分析的比较图。

具体实施方式

一、基因表达数据的预处理

读入基因表达数据文件，通过Z分数对基因表达数据进行标准化：

x表示每一个样本的基因原始表达值；μ表示每一个样本的所有基因原始表达数据的均值；σ为每一个样本的所有基因原始表达数据的标准差。

二、基因表达数据的主成分分析

基于标准化后的基因表达数据，通过主成分分析获得基因表达矩阵前两个主成分的具体过程为：

1)求标准化后的基因表达数据中特征的协方差矩阵；

2)求协方差矩阵的特征值及对应的特征向量；

3)将特征值按照从大到小的顺序排序，选择其中最大的两个，然后将其对应的两个特征向量分别作为列向量组成特征向量矩阵。

4)将样本点投影到选取的特征向量上。假设标准化后的基因表达数据是样本数为m，基因数为n的矩阵a(m*n维)，协方差矩阵为n*n维，选取的两个特征向量组成的矩阵为b(n*2维)，那么投影后的数据为：Y＝a×b，从而将原始样本的n维特征变成了两维。

三、高斯混合模型对样本聚类

设有随机变量X，高斯混合模型可以表示为：

表示混合模型中的第k个分量的高斯分布；π_k是混合系数，代表每个簇被选为随机变量所在族的概率。

四、对网络节点打分排序

对于聚类后的每一类样本，通过基于多网络的网络传播模型对网络节点进行打分排序。具体过程为：

1)首先对节点通过t统计量打分，作为节点的初始状态，表达式为：

为两类样本的平均数，

为两类样本的方差，n₁、n₂为两类样本的容量。

2)一个网络可以表示为一个邻接矩阵为A的图：G＝(V,E)，其中，A_ij表示节点v_i和v_j之间存在相互作用边。V＝{v₁,v₂,...,v_n}为网络中的节点，E＝{(v_i,v_j),i,j＝1,...,n,A_ij≠0}表示网络中的边。对于一个由S个网络构成的多网络，其邻接矩阵可以表示为：

其中，I表示n×n的单位矩阵，A_s表示第s层网络的邻接矩阵，μ可以量化节点在层内或在层间的概率。基于多网络的网络传播模型的表达式为：

F_t+1＝αW'F_t+(1-α)T (4)

其中，W'为多网络的邻接矩阵经过行标准化后的矩阵，α可以量化初始状态的相对重要性，α取值为0.7。

3)通过基于多网络的网络传播模型对多网络中的节点打分排序后，每个节点可以得到S个重要性度量值，通过几何平均求得节点的打分值，根据该打分值计算该节点在该簇中的排序值。

4)通过计算每个节点在所有样本簇的平均排序值，作为该节点最终的排序值。取前1％的基因作为初步筛选后的重要基因。

五、AUC优化模型获得生物标志物

根据上一步得到的重要基因，通过AUC优化模型进一步优化获得生物标志物。AUC优化模型的表达式可以通过一个线性规划方程来表示：

该线性规划方程服从于：

其中X₁,X₂,...,X_p为正样本，Y₁,Y₂,...,Y_q为负样本，每个样本都由一个k维特征向量表示。ξ_ij+、ξ_ij-分别为σ^T(X_i-Y_j)的正、负部分。σ是一个k维向量，表示为：σ＝{σ_k,k＝1,..,K}；θ是一个K维向量，表示这些K个特征被选为标志物的概率；当给出正负样本的中心指标，即

和ω可以表示为：ω＝β₊-β_-。λ设为0.00001。

通过AUC优化模型，基因可以根据概率排序，一组具有高概率的基因被选为生物标志物。

为了验证本方法的有效性，我们从疾病诊断、癌症预后、相似疾病分类三个方面对异质性复杂疾病进行了综合分析。我们在六个疾病的真实数据集上进行了测试，其中包含两个包含一型糖尿病患者和正常人样本的数据集，来源于GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)，分别为GSE35713和GSE9006；两个来源于GEO数据库的包含局限性肠炎和溃疡性结肠炎这两个相似疾病的数据集，分别为GSE3365和GSE36807；一个来源于TCGA数据库(https://portal.gdc.cancer.gov/projects)的乳腺癌高通量测序(HTS)数据BRCA和一个来源于Van De Vijver等人发表在New England Journal of Medicine上的用于乳腺癌病人生存分析的公开数据集NKI。如果患者存活10年以上，样本被标记为预后良好(对于BRCA，患者存活5年以上标记为预后良好)，如果患者生存时间不超过5年，则标记为预后不良。

Menche等人发表在Science上的六种生物相互作用公开数据集：(1)二元相互作用的酵母双杂交数据；(2)文献验证的低通量实验相互作用数据；(3)调控相互作用数据；(4)代谢相互作用数据；(5)蛋白质复合物；(6)信号网络相互作用数据。

为了评价本方法分类的准确性及生物可解释性，进行以下三项分析：

(1)样本分类的准确性

对于每一个数据集，首先通过五折交叉验证将样本分为训练集和测试集，在训练集上，基于本发明方法与NetRank,MarkRank,RWR,LP,PCA中的每一种方法提取到的特征，通过随机森林分类器在测试集上评价方法的准确性；为了获得稳定的分类结果，重复进行100次。采用AUC指标评价分类结果，AUC值的实验结果分别如图2所示。从图2可以看出，我们的方法得到的AUC值除了在GSE35713数据集上略低于MarkRank，其余数据集上均优于其他的方法。由此可见我们提出的方法具有很好的分类能力和稳定性。

(2)基于多网络与基于单一网络和融合网络方法的比较

为了进一步验证多网络的有效性，我们将本发明方法框架应用于每个单一网络，这意味着在网络传播模型中，只基于单一生物网络。结果如图3所示。在所有六个数据集上，任何单一网络的效果都不如多网络，这证明了我们将不同生物学意义的网络组合在一起的有效性。

同时，我们比较了基于多网络和基于融合网络的性能。融合网络依然包含六种类型的生物相互作用。在六个数据集上，基于多网络和基于融合网络的AUC结果如图4所示。结果表明，基于多网络的性能明显优于基于融合网络的性能，这证实了我们的假设，即不同类型的网络直接融合成一个网络可能导致网络信息的丢失，从而影响生物标志物的识别。

(3)预后生物标志物的生物可解释性

为了检验方法获得的生物标志物的生物可解释性，我们分析了获得的生物标志物对于已知致病基因和差异表达基因的富集程度。对于基因表达数据中的每一个基因，采用t检验获得差异表达基因(P值小于0.01)。通过超几何检验被计算已知致病基因和差异表达基因在生物标志物中富集程度的P值：

其中，N为所有基因的数量，M为所有基因中的已知致病基因和差异表达基因的数量，n为生物标志物的数量，m为在生物标志物中的已知致病基因和差异表达基因的数量。P值越小，说明已知致病基因和差异表达基因在生物标志物中的富集程度越高。基于六数据集得到的-log₁₀P的结果如图5所示，从图5可以看出，本发明方法的-log₁₀P总体上比其他方法富集程度更高，说明已知致病基因和差异表达基因在本发明方法获得的生物标志物中的显著富集，即具有良好的生物可解释性。