1、引言

癌症的发生从微观方面来讲是一种异常细胞不受控制地分裂，并可通过血液和淋巴系统侵入附近组织和身体其他部位的疾病。在过去的十年中，一些大规模的癌症基因组学项目已经发表了来自数千名癌症患者的基因组、表观基因组、转录组和蛋白质组的数据。这些项目包括癌症基因组图谱(TCGA)、国际癌症基因组联盟(ICGC)和全基因组泛癌症分析(PCAWG)。癌症微阵列技术在基因组研究、计算生物学、统计学和机器学习等学科中开辟了广泛的多学科研究领域。在微阵列癌症数据领域开展研究，对于癌症患者的诊断、癌症亚型的识别和区分具有重要意义[1]。

通过高通量测序技术[2]获得的基因表达的数据结构复杂，信息冗余度高，基因之间的相关性很强。传统的生物学研究方法难以有效地处理基因表达数据，所以从基因表达数据中挖掘基因之间的相互作用关系，构建基因调控网络，成为生物信息学中重要的研究课题之一。Chen等[3]运用皮尔逊相关系数的方法建立了基因调控网络，该方法需要不断尝试设置最优阈值，以此来保留调控网络中相关性强的基因间关系。加权基因共表达网络分析(Weighted Geneco-expression Network Analysis, WGCNA)方法基于表达模式类似的分子可能参与特定生物学功能的理论，最初由 Zhang 和Horvath[5]提出。它被广泛研究并用于预测新的基因功能[4],发现新的疾病生物标记物，以及检测癌症中的遗传变异。

近年来，人们提出了不同的计算方法来检测癌症亚型。这些方法通常建立在特征工程的基础上，对患者进行聚类或分类。由于生物数据高维且样本量小，早期的方法在一定程度上减少了样本的特征，并利用这些特征来聚类癌症亚型。Li等[6]将L1正则化惩罚添加到目标函数中，该方法在处理噪声和异常值方面是非常有效的。Guo等[7]提出了两步L1正则化方法对微阵列数据进行分类，该文定义了一种新的L1正则化特征选择方法，以去除不相关和冗余的特征并选择重要的特征。在我们的工作中，为了处理小样本的维度灾难，引入了L1正则化的特征选择，然后使用深度学习技术对多类微阵列癌症数据进行分类。

随着深度学习技术的发展，一些诸如卷积神经网络(CNNs)[8,9,10]、堆栈自编码器(SAE)[11]、生成对抗模型(GANs)、卷积自编码器(CAE)、变分自编码器(VAE)[12]等的深度学习方法已经应用于癌症研究领域[13,14]。基于分类的方法通常是在一些已知亚型标签的癌症样本上训练一个模型，然后使用该模型预测新的癌症样本的亚型。根据特征提取和分类的差异，基于分类的方法分为两类：两阶段方法和端到端方法。两阶段方法需要先对基因表达数据进行特征选择，然后再进行分类预测。端到端的方法可以同时进行选择和分类。DeepType[15]将癌症样本数据输入多层神经网络用以降低数据维数，同时将监督分类与非监督聚类相结合，用聚类结构确定癌症相关数据的表示形式。Dai等[16]设计的端到端的深度学习方法ERGCN,使用癌症样本相似性网络和残差图卷积网络来对癌症亚型进行分类，得到了较优的分类结果。本文使用两阶段的分类方法，将特征选择与端到端的分类方法相结合，试图进一步提升癌症预测的准确性。目前流行的深度神经网络在其架构中仅考虑基因之间的交互和关联，不考虑患者之间的交互和关联。考虑患者之间的关系是有益的，因为它有助于一起分析和研究相似的患者队列。通过将患者视为节点，将他们的相互作用关系视为边，图表提供了一种自然的方式来表示患者之间的相互作用。

为了解决这一问题，我们提出一种新的预测模型：基于加权基因相似网络和样本相似网络融合算法的癌症亚型预测模型(WGCSS)。在数据预处理中，为了减低基因数据的维度，采用L1正则进行特征选择；在构建基因相似网络时，模型采用自适应选择策略进行阈值选择；在构建样本相似网络后，本文融合了这两个相似网络中的信息；在本文中使用图卷积网络模型进行预测，为了避免训练过程中的过度平滑，在模型中构建了残差层。实验结果表明本文模型优于现有方法。

2、提出的方法

首先，采用L1正则进行特征选择，它通过在线性回归的最小二乘误差中加入惩罚项来应用收缩策略，其将不相关的特征赋零系数，仅考虑非零系数变量。该方法解决了模型的过拟合问题，减小了预测误差，简化了模型的复杂度，提升了计算稳定性。然后，基于WGCNA计算癌症患者基因表达之间的相似度，构建一个以基因为节点的加权基因相似网络。样本相似网络是在样本空间中基于皮尔逊相关系数构建的，本文的基因图和样本图都是提前构造的，这可以大大降低计算的复杂度，保持样本之间关系的稳定性。最后，采用残差图卷积神经网络算法进行分层消息传递，通过网络扩散节点的特征信息，学习每个节点的特征表示。残差层的构建主要是为了防止在两个空间聚合信息时的过度平滑问题。WGCSS根据患者的特征表征来预测癌症亚型，该模型由两层卷积层、两层残差层和三层线性层构成。图1给出了WGCSS模型概述。

图1 WGCSS模型流程图

2.1 基于L1正则的特征选择

L1正则化特征选择是微阵列数据分析中引入的一种新的特征选择方法。基于L1的特征选择使用LSVM拟合数据，并返回将数据划分为类别的最佳拟合超平面。它利用局部最优解去除系数为零的特征。在本研究中，使用L1正则化支持向量机算法进行特征选择。

假设一个数据集S有n个实例，如式(1)所示：

其中，xi是第i个实例，具有n个特征和一个类标签yi。xi的表达式如式(2)所示：

xi={xi1,xi2,…,xin} (2)

其中，xij是实例xi中第j个特征的值。

对于二分类问题，支持向量机SVM的原理要找到最优分离超平面h×x=b,支持向量到超平面的距离为

其中ω是权重向量，b是偏置项，其表达式如下：

但须符合以下条件：

Bradley和Mangasarian[17]提出了L1-SVM算法，根据所得到的稀疏解进行特征选择，如式(5)所示：

其中，εi为松弛变量，α>0为错误惩罚参数。但须符合以下条件：

L1正则支持向量机的优化问题如式(7)所示：

该算法利用控制参数α来控制数据的稀疏性。稀疏性允许矩阵中的少数特征具有较大的非零系数值。本文中设定α=0.5。

2.2 加权基因相似网络

我们引用WGCNA中的相似度矩阵概念，以构建加权基因相似网络。首先，计算所有基因的相似度共表达矩阵sij,其中sij=|r(xi,xj)|为节点i和节点j的基因表达谱之间相关系数的绝对值。本文采用距离相关法计算相关系数。因为距离相关系数总是正的，所以定义了一个无符号网络，其中正相关和负相关是同等的。然后，通过设置aij=sijβ,以β为软阈值幂，将相似度共表达式矩阵转化为邻接矩阵aij,然后计算TOM拓扑重叠矩阵(使用R中的pick Soft Threshold 函数，加权是指对相关性值进行幂次运算，幂次的值即为软阈值)。

2.3 样本相似网络

我们根据患者的相似性构建了一个患者网络。通过患者基因表达谱的Pearson相关系数(PCC)计算患者间的相似性(见式(8))。

Pearson相关系数(r(x,y))可以衡量两个患者之间的线性相关程度，输出范围为-1到1,当r的绝对值接近1时，两个患者呈正相关或负相关，否则两个患者无相关性。因此，我们认为两个患者可由一条边连接，如果两个样本之间的Pearson系数绝对值大于阈值θ,则将邻接矩阵中对应的值设为1,否则两个患者之间没有边连接，对应的邻接矩阵值为零。在本文中，θ取值为0.8。

2.4 残差图卷积网络

在本文的模型中，构造了两个图：G(G)=(V(G),E(G))和G(s)=(V(s),E(s))。G(G)为基因相互作用图，其中V(G)为基因集合，E(G)为表示基因间相互作用的边集合。与G(G)不同，G(s)是样本相似度图，其中V(s)是样本的集合，E(G)是表示样本之间相似度的边的集合。A和N分别为G(G)和G(s)的邻接矩阵，Dg和Ds分别为A和N的度矩阵，X为n×D的网络节点初始属性矩阵。X为n×D的网络节点初始属性矩阵。本章的GCN模型需要3个输入：两个存储节点连接的邻接矩阵和一个节点初始属性矩阵。在特征空间中的GCN定义如下：

在基于GCN的方法中，过度平滑是一个常见的挑战，在两个空间中聚合信息时会非常严重。因为要先在特征空间(垂直)中聚合，然后在样本空间(水平)中聚合，每个元素将聚合两次。为了防止模型过度平滑，本章同样也在GCN中增加了一层残差数据。为了保证输入特征的维度与一层GCN的节点特征维度一致，我们将初始输入特征通过一个独立的线性层直接连接到GCN层的输出。可以写成如下形式：

H(P)=H(1)+tanh(linear(xT)) (10)

在样本空间中的GCN定义如下：

同样，在样本空间的GCN中也增加了一层残差数据，如下所示：

H(q)=H(2)+tanh(linear(xT)) (12)

使用交叉熵损失函数来量化癌症亚型预测损失：

给定训练集T中的患者k,Ykf为符号函数(0或1),Zdf为癌症患者k属于f类的预测概率。我们将癌症亚型预测损失最小化，以优化WGCSS。

3、实验和结果分析

3.1 实验数据集

为了验证WGCSS方法的有效性，使用了来自TCGA(the Cancer Genome Atlas)的乳腺浸润性癌(BRCA)、多形性胶质母细胞瘤(GBM)和肺癌(LUNG)的肿瘤数据集，在Wang[18]的补充文件中下载了3种癌症类型的基因表达数据和生存信息。使用R包TCGAbiolinks从TCGA中检索癌症亚型信息。通过匹配基因表达数据的样本ID,将其样本基因表达数据与他们的癌症亚型结合起来。最后整理出102例乳腺浸润性癌患者、213例多形性胶质母细胞瘤患者和85例肺癌患者都有4种癌症亚型，详细信息如表1所列。

表1 数据集说明

3.2 评价指标

外部评价指标：通过将预测分类结果与实际分类结果进行比较来评价算法的有效性。每个性能指标都有其优缺点。为了缓解这种限制，使用了诸如准确度、召回率、精度、F1分数、混淆矩阵和宏观平均ROC等性能度量进行比较。这些指标的公式如下：

其中，TP为模型预测的阳性样本为阳性，TN为模型预测的阴性样本为阴性，FP为模型预测的阴性样本为阳性，FN为模型预测的阳性样本为阴性。

我们还使用了一个内部评价指标——Davies-Boulding指数：

其中，n为聚类个数，avg(cj)为第i类样本到其聚类质心的平均距离，dis(ci,cj)为ci类中心到cj类中心的距离。DBI的下限为0,且DBI值越小，聚类效果越好。

3.3 实验设置

在特征选择阶段，设置系数性参数α=0.5。使用Adam优化器函数，学习率设置为0.001。为了避免“选择偏差”,在所有微阵列数据集的实验中使用了5倍交叉验证(CV)。每个数据集的样本被随机分成5个大小相等的子样本，其中4个子样本用于训练，其余一个子样本用于测试。该过程重复5次，5个子样本中的每一个都用作测试数据。所有比较方法都在由5倍CV随机生成的同一数据集上实现。在构建基因相似网络时，本文要求k与p(k)的相关性达到0.85时的power作为β值。根据无标度拟合指数平均连接度分析BRCA的β值为2,GBM的β值为6,LUNG系统没有给出合适的β值，我们手动调节软阈值并给出它的β值为3。构建样本相似网络时的θ取值为0.8。

3.4 实验结果与分析

本文对3种标准的多类微阵列癌症数据集进行了实验，结果如表2所列(将5倍交叉验证测试集的平均结果作为评价指标)。表2显示了本文的模型在BRCA,GBM和LUNG数据集上关于准确度、召回率、精度、F1分数和DBI的性能。经过特征选择，3个数据集的特征都大规模减少。所提出的方法在乳腺浸润性癌和肺癌数据集上显示出完美的分类性能评分1.00,在多形性胶质母细胞瘤数据集上的分类准确度为0.97。文中方法在乳腺浸润性癌和肺癌数据集实现了100%的精度、召回率和F1度量；在多形性胶质母细胞瘤数据集中精确度、f度量和召回率都达到了97%。DBI指数在3个数据集上都较低，证明本章的模型有较优的聚类效果。

表2 WGCSS模型分类结果

通过对3种标准的多类微阵列癌症数据集进行实验，我们发现本文提出的算法可以提高模型的分类性能和鲁棒性。该算法在多个数据集上取得了较好的分类效果，表明其在微阵列癌症分类领域具有广阔的应用前景。

在分类问题中，ROC曲线被广泛用于检验模型的性能。3类数据集的五倍交叉平均ROC曲线如图2所示。本文模型在多形性胶质母细胞瘤数据集上所得的AUC评分为99.15%,在乳腺浸润性癌和肺癌数据集上也都实现了完美的分类精度，因此这两类数据集的AUC评分为100%。每次交叉验证的结果都是较优的，证明本文的模型也是较为优越的。

图2 五倍交叉验证后的平均ROC曲线

3.5 对比分析

Dai等提出的ERGCN模型是基于PCC构建的样本相似网络，其只考虑了样本空间的信息而忽略了特征空间中的基因信息。我们将经过特征选择后的基因组数据输入该模型，其结果如表3所列(加粗字体为WGCSS模型的结果)。可以看到，在BRCA数据集上两个模型都达到了最优，在GBM数据集上的各项评价指标上，本文模型比ERGCN模型提高了2%左右；在LUNG数据集上的各项评价指标上，本文模型比ERGCN模型提高了1.5%左右。对于DBI指数，本文的模型在BRCA数据集上比ERGCN模型降低了6.76%,在GBM和LUNG数据集上各降低了约13.13%和8.92%,说明本文的模型具有更优的聚类效果。上述结果说明基因之间和样本之间的相互作用都包含了有价值的癌症样本分类信息。

表3 ERGCN模型和WGCSS模型的分类结果比较

在分类精度方面，将本文模型与现有的两阶段分类模型进行比较，如表4所列。可以看到，文中模型在BRCA数据集上比DFN Forest模型提高了19.14%,比SGL-SVM模型提高了14.35%,比CFN Forest模型提高了5.64%;在GBM数据集上比DFN Forest模型提高了10.74%,比SGL-SVM模型提高了14.63%,比CFN Forest模型提高了5.99%;在LUNG数据集上比CFN Forest模型提高了9.11%,比SGL-SVM模型提高了4.27%,比MOEDA模型提高了4%,比RNBC模型提高了11.5%。

图4 WGCSS模型和ERGCN模型结果的对比图

表4 不同两阶段分类模型精度的比较

图5 不同两阶段模型的精度比较

3.6 生存分析

为了进一步探讨所识别的亚型之间的关系，本节对WGCSS模型的结果进行了生存分析。从理论上讲，不同的癌症亚型应有不同的生存曲线。图6展示了文中模型在BRCA,GBM和LUNG数据集上的Kaplan-Meier生存曲线，横轴(x轴)表示以天为单位的时间，纵轴(y轴)表示生存的概率或生存人口的比例。图中曲线代表两组病人的生存曲线。曲线的垂直下降表示事件。曲线上的垂直刻度表示这个病人在这个时候被审查了。本文计算了不同亚型生存曲线上log-rank检验的p值，还在曲线上绘制了中位生存时间。这3类癌症参与者在时间为0时，生存概率是1.0,即所有的参与者都活着。在BRCA数据集上，第一种亚型的中位生存时间为1 563天，第二种亚型为3 418天；第三种亚型为NA;第四种亚型为2 227天。NA表示第三组大多数患者无法活过中位生存时间。在该数据集上第二种亚型的生存率明显优于其他3种亚型，且每个类别之间的中位生存时间差距都很大，表明在该数据集上的癌症亚型的生存曲线存在显著差异。对于GBM数据集，癌症亚型间的差异不是很明显。第一种亚型患者中位生存时间为394天；第二种亚型为455天；第三种亚型为440天；第四种亚型为362天。在LUNG数据集上，第一种亚型的中位生存时间为888天；第二种亚型为761天；第三种亚型为631天，第四类是1 456天。在该数据集上，第四种亚型的生存率明显优于其他3种亚型，且每个类别之间的中位生存时间差距也都很大，表明在该数据集上的癌症亚型的生存曲线存在显著差异。在BRCA、GBM和LUNG数据集生存曲线上，log-rank检验的p值分别为0.034,0.087,0.056,因此文中模型在3个癌症数据集上癌症亚型的生存曲线存在显著差异。

图6 不同亚型的生存时间

4、结束语

本文提出了一种基于加权基因相似网络和样本相似网络融合算法的癌症亚型预测模型(WGCSS)。该模型不仅考虑了基因之间的相互作用关系，同时也考虑了样本之间的相互作用关系，通过聚合两个空间中的基因数据信息，使得预测的结果更加准确。与仅考虑样本空间或特征空间的模型相比，该模型有更优的准确率和聚类效果，说明样本之间和基因之间的相互作用都包含了有价值的癌症样本分类信息。最后进行的生存分析也表明，在该模型，癌症亚型的生存曲线是有显著差异的。本文研究为癌症亚型分类预测提供了一种新的方法，也为精准医疗提供了新的可能。

基金资助:国家自然科学基金项目(11571009); 山西省应用基础研究项目(201901D111086); 山西省重点研发计划项目(202102020101004); 山西省回国留学人员科研资助项目(2022-074)~~;

文章来源:张晓茜,李东喜.基于相似网络融合算法的癌症亚型预测[J].计算机科学,2024,51(S1):639-645.