1. 引言
随着多种组学技术的快速发展,产生了大量的不同层面的组学数据,包括基因组、转录组、蛋白质组和代谢组等。癌症是一种高度异质性的基因疾病,涉及多种生物学过程的失调。单组学分析只能在单一水平上揭示肿瘤的分子信息。多组学整合分析可以综合不同组学的信息,能够全面系统地理解肿瘤的分子生物学机制,有助于癌症诊断和临床治疗。多组学整合分析逐渐应用于癌症生存分析的研究中[1] [2]。
生存分析(survival analysis)是研究事件发生前的持续时间数据的统计方法,在生物医学、工程学和经济学等领域都有广泛的应用[3]。Cox比例风险模型(Cox proportional-hazards, Cox-PH)是生存分析中使用最广泛的统计模型,由英国统计学家D.R.Cox于1972年首次提出[4]。Cox模型是以最终结局和生存时间为因变量的半参数回归模型,通过估计风险比来评估不同因素对生存时间的影响。随着Cox模型的推广,衍生出基于神经网络的非线性Cox模型,例如Cox-nnet,DeepSurv,CoxPASNet和Cox-Time等[5]-[8]。
深度学习技术具有强大的自动特征学习能力,在癌症基因组学方面具有巨大的潜力。研究者陆续提出了基于深度学习和多组学融合的生存分析模型[9]-[12]。为了全面平衡不同组学之间的相似性和差异性,准确描述复杂多样的多组学数据,研究者提出在多组学特征提取中引入对抗性学习。例如,Hai等人[13]提出了基于多组学融合和生成对抗网络的癌症亚型分类Subtype-GAN模型。模型中的生成器对四种组学数据进行低维表示,判别器利用共享层隐变量进行癌症亚型分类;Mondoi等人[14]提出了基于对抗自编码器的双阶段神经网络架构AFExNet,从高维单组学基因表达数据中进行特征提取,再应用各种分类器进行癌症的亚型分类。
受上述文献的启发,本文提出一种基于多组学融合和生成对抗网络的生存分析GAN-1DCCox模型。该模型将基于对抗自编码器的多组学特征提取与1D-CNNCox模型相结合进行生存分析。模型的创新之处:基于对抗性学习的多组学特征提取,充分利用了不同组学数据之间的互补性和差异性。在8种癌症数据集上进行了实验,结果表明GAN-1DCCox模型具有很好的多组学融合能力,提升了模型的预测性能。
2. 研究方法
2.1. Cox比例风险模型
生存分析常用于预测某种终点事件发生之前的等待时间,比如患者死亡、疾病复发、机器故障、客户流失等。生存分析的目标是充分利用删失数据所提供的不完全信息,对影响生存时间的因素进行分析。随访中未观察到终点事件发生的数据称为删失(censored)。删失数据的真实生存时间未知,只知道大于观察到的删失时间。常用三元组
表示生存数据,
表示协变量,
表示终点事件发生,
对应生存时间
;
表示删失,
对应删失时间
。
Cox比例风险模型假定风险比率为固定值,即协变量对生存概率的影响不随时间而变化:
(1)
其中
表示基线风险函数,
称为效应参数(偏回归系数),
称为风险评分或预后指数(Prognostic Index, PI)。Cox比例风险模型的偏似然函数定义为所有事件发生的联合概率,即:
(2)
其中
表示仍处于死亡风险的患者集合,即生存时间大于第i个患者生存时间的集合。通过最大化对数偏似然函数,对回归系数
进行估计,即:
(3)
2.2. 对抗自编码器
生成对抗网络(Generative Adversarial Networks, GAN)最早由Ian Goodfellow于2014年提出,以其优越的性能迅速成为深度学习领域的研究热点[15]。GAN包含两个神经网络:生成器
和判别器
。生成器的目标是生成尽量真实的虚假样本
,欺骗判别器将其标记为真实数据。判别器的目标是尽量区分虚假样本,输出[0, 1]的概率值(是否为真实数据的概率)。生成器试图使
接近于1,判别器试图使
接近于0。通过对抗博弈的方式训练判别器和生成器,其对抗损失函数为:
(4)
第一项是对真实数据的期望,第二项是对噪声数据的期望。上述minmax优化问题通过交替优化求解:首先固定生成器G,最大化损失函数优化判别器D;然后固定判别器D,最小化损失函数优化生成器G。
Makhzani等人在2015年提出了引入对抗性学习的对抗自编码器(Adversarial Autoencoder, AAE) [16]。AAE模型包含三个组件:编码器(Encoder)、解码器(Decoder)和判别器(Discriminator)。AAE的框架是将对抗网络添加在自编码器隐藏层的顶部,AAE的结构见图1。Encoder将真实数据x映射为隐编码z,Decoder使用隐编码z得到重构数据x'。Encoder在自编码器和对抗网络之间共享,Encoder同时是对抗网络的生成器,试图生成尽量真实的隐编码z。Discriminator试图将从先验分布采样的判别为真实,将编码器生成的隐编码z判别为虚假。AAE的训练过程分成两个阶段:一是自编码器的重构阶段,通过最小化重构损失更新Encoder和Decoder;二是对抗网络的正则化阶段,使用对抗训练方法更新判别器和生成器,通过最小化对抗损失函数交替更新Discriminator和Encoder,以提高编码器混淆判别器的能力。
Figure 1. The network architecture of an adversarial autoencoder
图1. 对抗自编码器的网络结构
2.3. 基于多组学融合和对抗自编码器的生存分析模型
本文提出了一种基于多组学融合和生成对抗网络的生存分析模型,简称GAN-1DCCox模型。该模型由两部分构成:基于对抗自编码器的多组学特征提取网络和基于一维卷积的1D-CNNCox生存分析模型。
2.3.1. 基于对抗自编码器的多组学特征提取
为了多组学数据的有效融合,本文设计了基于对抗自编码器的多组学特征提取网络。该网络由三个子网络构成:编码器、解码器和判别器,其结构如图2所示。编码器阶段进行多组学数据的低维嵌入。首先每种组学分别连接到组学独立层,输出的隐特征拼接后组成一个共享层。组学独立层是一层的全连接网络,假设
为第k种组学输入的第i个神经元,
为第k种组学输出的第j个神经元,则有
(5)
(6)
为了使训练过程不对特征尺度敏感,提高训练速度,在组学独立层之后添加了批归一化层。批归一化(batch normalization)是指对每个小批量训练数据进行标准化,即:
(7)
(8)
解码器具有与编码器对称的网络结构,将共享层隐特征通过组学独立层,重建原始的多组学输入。假设
表示共享层到第k种组学的重建输出,重构阶段的损失函数L1为:
(9)
在对抗网络的正则化阶段,通过对抗性学习训练编码器和判别器,对应的对抗损失函数L2为:
(10)
对抗自编码器网络的损失函数是重构损失函数和对抗损失函数的加权和:
(11)
最后,在模型训练完成后,取共享层的隐编码向量z作为下游生存分析模块的输入。
2.3.2. GAN-1DCCox生存分析模型
卷积神经网络(Convolutional Neural Network, CNN)利用卷积滤波器进行空间局部特征的自动提取,通过权重共享和局部连接减少了参数量,提高了模型的效率,已成为计算机视觉领域最流行的网络架构。虽然多数卷积网络采用2维卷积Conv2D,但Conv2D应用一维序列数据时存在局限性,1维卷积1D-CNN已成功应用于语音识别和机器翻译[17]。受Mostavi等人[18]提出的用于癌症分类的卷积网络模型的启发,Yin等人提出了生存分析的卷积网络模型:基于Conv2D的CNN-Cox模型和基于1D-CNN的1D-CNNCox模型,具体的网络架构详见文献[12]。1D-CNN具有浅层网络结构,相对Conv2D训练的计算成本较低,非常适合高维基因组数据的研究。因此,本文选取了1D-CNNCox模型搭建新模型的生存分析网络层。
本文提出的GAN-1DCCox模型由两个模块构成:基于对抗自编码器的特征提取模块和1D-CNNCox生存分析模块,具体架构见图2。特征提取模块是一个基于多组学数据的多输入多输出GAN网络,目的是充分利用不同组学之间的互补性,提取有效的关键基因。基于1D-CNNCox模型的生存分析模块进一步提取预后基因,提升模型的准确性和稳健性,使其在不同癌症数据集上均能取得优异的预测性能。
Figure 2. The network architecture of GAN-1DCCox model
图2. GAN-1DCCox模型的网络结构
3. 实验结果分析
3.1. 数据集和预处理
本文数据来源于癌症基因组图谱TCGA数据库中的8种癌症类型,包括BLCA (膀胱癌),BRCA (乳腺癌),CESC (宫颈鳞状细胞癌),KIRC (肾透明细胞癌),LIHC (肝细胞癌),SKCM (黑色素瘤),STAD (胃癌)和UCEC (子宫内膜癌)。利用UCSC Xena数据库下载TCGA数据集(https://xenabrowser.net/datapages/)。每个数据集包含样本的生存信息(生存结局和生存时间)和三种组学数据:RNA-seq,miRNA-seq和拷贝数变异(Copy Number Variation, CNV)。对于三种组学数据进行如下的预处理:1) 将不同组学以及生存信息的样本数据取交集,去除重复样本和筛选共同样本;2) 对于RNA-seq数据(60,483维),首先将表达值为0的基因剔除,其次基于中位数绝对偏差(Median Absolute Deviation, MAD)选取前18,390个基因,最后根据风险分组做进一步的基因筛选;3) 对于CNV数据(19,729维)和miRNA-seq数据(1881维),仅进行基于风险分组的基因筛选。
根据生存时间进行风险分组的过程如下:首先选取合适的生存时间阈值M,将OS.time ≤ M的死亡患者划为高风险组,OS.time ≤ M的患者划为低风险组,再利用sklearn.feature_selection库的SelectKBest函数进行基因选择。SelectKBest是一种filter-based特征选择方法,依赖统计指标对特征进行评分和排序。由于风险分组属于分类任务,使用基于方差分析的评分函数f_classif,选择得分最高的k个特征。对于RNA-seq,CNV和miRNA-seq,分别取k = 3000,3000,400。8种不同癌症数据集的样本量、删失率和风险分组情况等信息见表1。
Table 1. The sample information of 8 different cancer datasets
表1. 8种不同癌症数据集的样本信息
癌症类型 |
样本量 |
删失率 |
生存时间阈值M |
高风险组 |
低风险组 |
BLCA |
397 |
55.92% |
700 |
136 |
141 |
BRCA |
1032 |
86.3% |
1095 |
66 |
419 |
CESC |
279 |
75.63% |
3046 |
51 |
89 |
KIRC |
505 |
66.93% |
1095 |
107 |
277 |
LIHC |
359 |
64.62% |
1095 |
103 |
88 |
SKCM |
432 |
50.93% |
1095 |
106 |
222 |
STAD |
345 |
58.26% |
650 |
117 |
114 |
UCEC |
521 |
83.49% |
1095 |
67 |
199 |
3.2. 模型的有效性
为了评估GAN-1DCCox模型的有效性,我们进行了两种实验:第一种是消融实验,以验证GAN特征提取模块对多组学整合的有效性;第二种实验是将其与流行的生存分析模型进行性能比较,包括带ElasticNet惩罚的Cox模型(Cox-EN) [19],随机生存森林(RSF) [20],梯度提升机(GBM) [21],生存支持向量机(SSVM) [22],1D-CNNCox [12]。调用scikit-survival库的相应函数来实现前四种基准模型,利用Tensorflow和Keras框架实现1D-CNNCox和GAN-1DCCox模型。
C指数(concordance index, C-index)是生存分析中最常用的评价指标,用于衡量预测结果与实际生存结局的一致性。C指数越高模型的风险区分度越好。C指数定义为预测结果与实际结果一致的样本对在所有可比较的样本对中的比例。如果两个样本都经历过事件,或者生存时间较短的样本经历过事件,生存时间较长的另一个样本未经历事件,这样的样本对称为可比较的样本对。C指数的计算公式如下:
(12)
本文还计算了8种癌症数据集上的微平均C指数,其定义如下:
(13)
本文使用5折交叉验证评估模型的性能,使用网格搜索法对模型的超参数进行优化。GAN特征提取模块的超参数包括训练迭代轮次、网络层数和隐藏层节点数。生存分析模块的超参数包括一维卷积核尺寸、卷积滤波器个数和全连接层的节点数。一维卷积核的尺寸设为1 × 12,滤波器个数和全连接层节点数的搜索范围分别为[8, 16, 32, 64, 128]和[16, 32, 64, 128, 512]。GAN特征提取模块的超参数设置见表2。
Table 2. Hyperparameter settings of GAN feature extraction module
表2. GAN特征提取模块的超参数设置
网络层 |
编码器 |
解码器 + 判别器 |
1D-CNNCox |
Input |
RNA 3000 |
CNV 3000 |
miRNA 400 |
144 |
12 × 12 |
Hidden |
64 |
64 |
14 |
144 |
— |
Concat + Batch Norm |
142 |
144 |
— |
Output |
144 |
RNA 3000 |
CNV 3000 |
MiRNA 400 |
Outcome 1 |
C-index 1 |
Figure 3. The impact of feature extraction dimension on the performance of GAN-1DCCox
图3. 特征提取维数对GAN-1DCCox模型性能的影响
GAN特征提取模块中共享层隐特征的维数选为144,因为在特征提取维数对模型性能影响的实验中(40000/400/196/144/100/81/49),特征维数选为144时GAN-1DCCox模型的性能最优,即在8种癌症数据集上的微平均C指数最高,见图3。
图4(a)展示了基于GAN特征提取的144个特征,再经由1D-CNNCox、Cox-EN、GBM、RSF、SSVM五种生存分析模型,在8种癌症数据集上的C指数箱线图。可以看到本文所提模型GAN-1DCCox在所有数据集上的性能均显著优于其他基准模型。为了进一步检验GAN特征提取模块的有效性,将本文模型GAN-1DCCox与基于40,000个原始特征的1D-CNNCox、Cox-EN、GBM、RSF、SSVM五种基准模型在8种癌症数据集上的C指数进行比较,箱线图如图4(b)所示。可以发现,在8种癌症数据集上GAN-1DCCox的C指数值总是最高的。
Figure 4. The effectiveness of GAN feature extraction module
图4. GAN特征提取模块的有效性
我们还进行了消融实验,将本文GAN-1DCCox模型与GAN + SSVM,GAN + Cox-EN,GAN + RSF进行生存预测性能的比较,各个模型取得的C指数值见表3。实验结果表明,GAN-1DCCox的性能显著优于其他模型,证实了GAN特征提取的有效性。这揭示了GAN能够提取重要特征,一定程度上缓解多组学数据融合的过拟合问题,从而提高模型的预测性能和泛化能力。
Table 3. C-index values of all models on 8 cancer datasets
表3. 各个模型在8种癌症数据集上的C指数值
Model |
BLCA |
BRCA |
CESC |
KIRC |
LIHC |
SKCM |
STAD |
UCEC |
GAN-1DCCox |
0.683 |
0.721 |
0.738 |
0.700 |
0.716 |
0.677 |
0.619 |
0.725 |
1D-CNNCox |
0.562 |
0.566 |
0.591 |
0.660 |
0.631 |
0.581 |
0.529 |
0.633 |
GAN + GBM |
0.669 |
0.678 |
0.726 |
0.691 |
0.700 |
0.654 |
0.592 |
0.693 |
GAN + SSVM |
0.645 |
0.679 |
0.681 |
0.620 |
0.616 |
0.621 |
0.599 |
0.605 |
GAN + Cox-EN |
0.641 |
0.690 |
0.672 |
0.634 |
0.648 |
0.639 |
0.602 |
0.615 |
GAN+RSF |
0.668 |
0.674 |
0.717 |
0.694 |
0.695 |
0.633 |
0.603 |
0.705 |
GBM |
0.651 |
0.623 |
0.597 |
0.690 |
0.663 |
0.641 |
0.521 |
0.685 |
SSVM |
0.633 |
0.668 |
0.673 |
0.696 |
0.647 |
0.589 |
0.572 |
0.625 |
Cox-EN |
0.606 |
0.642 |
0.660 |
0.677 |
0.624 |
0.591 |
0.517 |
0.613 |
RSF |
0.596 |
0.607 |
0.594 |
0.679 |
0.640 |
0.583 |
0.535 |
0.660 |
本文进行了GAN-1DCCox模型在不同组学数据上的消融实验,见图5和表4。结果显示,在8种癌症数据集上基于3种组学数据的C指数均高于双组学和单组学。这表明3种组学的融合比单组学和双组学在癌症基因特征提取方面更具优势。
8种癌症数据集上的生存曲线及其对数秩检验结果,见图6。结果表明,高风险组的总体生存率明显低于低风险组,高低风险组的预后差异在6种癌症数据集上高度显著(对数秩检验p值低于0.005,KIRC,BRCA除外)。这证实GAN-1DCCox模型可以将不同癌症类型的患者进行高低风险组的准确分类。
Table 4. C-index values of GAN-1DCCox model based on integration of different omics
表4. 基于不同组学融合的GAN-1DCCox模型的C指数值
Omics |
BLCA |
BRCA |
CESC |
KIRC |
LIHC |
SKCM |
STAD |
UCEC |
Three |
0.683 |
0.721 |
0.738 |
0.700 |
0.716 |
0.677 |
0.619 |
0.725 |
RNA + CNV |
0.667 |
0.688 |
0.712 |
0.686 |
0.705 |
0.647 |
0.611 |
0.690 |
miRNA + CNV |
0.607 |
0.640 |
0.648 |
0.669 |
0.667 |
0.632 |
0.590 |
0.690 |
RNA + miRNA |
0.671 |
0.701 |
0.690 |
0.694 |
0.693 |
0.642 |
0.583 |
0.701 |
RNA |
0.665 |
0.684 |
0.713 |
0.691 |
0.679 |
0.650 |
0.609 |
0.697 |
miRNA |
0.599 |
0.660 |
0.648 |
0.667 |
0.673 |
0.605 |
0.555 |
0.652 |
CNV |
0.554 |
0.638 |
0.620 |
0.611 |
0.623 |
0.613 |
0.600 |
0.676 |
Figure 5. Performance comparison of GAN-1DCCox under different omics
图5. GAN-1DCCox模型在不同组学数据下的性能比较
Figure 6. Kaplan-Meier curves of GAN-1DCCox model on 8 cancer datasets
图6. GAN-1DCCox模型在8种癌症数据集上的Kaplan-Meier曲线
3.3. 预后基因的识别
识别具有生物学意义的基因子集,对于揭示癌症疾病的潜在机制至关重要。作为模型解释的示例,本文研究了乳腺癌BRCA数据集的预后相关基因。首先,进行基因集富集分析GSEA,筛选高风险组和低风险组之间的差异基因。然后进行蛋白质–蛋白质相互作用PPI网络分析,分析枢纽基因的生物学功能。
GSEA是一种评估基因子集在两种生物状态(高低风险组)之间表达是否显著差异的方法。BRCA数据集的基因富集图见图7(a),归一化富集分数ES值为1.24,显著性p值为0.032,每个表型前50个基因特征的热图见图7(b)。最终获得了747个富集分数小于−0.23或大于0.25的差异基因。再将747个差异基因导入STRING数据库构建PPI网络。当置信度阈值设置为0.7时,得到39个节点和32条边,见图7(c)。使用Cytoscape中的CytoHubba插件,识别PPI网络中的功能簇和高度连通的10个hub基因,包括CD40LG,FOXP3,FLT3LG,CCL22,CD1C,STAT5A,LTB,TNFRSF17,SLAMF1,CLEC10A,见图7(d)。
Figure 7. GSEA and PPI network analysis on BRCA datasets
图7. BRCA数据集的GSEA和PPI网络分析
枢纽hub基因是PPI网络中的核心基因,是疾病潜在的生物标志物。使用Gene Ontology分析10个枢纽基因富集的GO生物学过程,包括适应性免疫反应(adaptive immune response)、细胞因子介导的信号通路(cytokine-mediated signaling pathway)和先天性免疫反应(innate immune response),见表5。
Table 5. Enriched GO terms of hub genes in PPI network of BRCA
表5. BRCA数据集中的PPI网络枢纽基因的富集GO过程
GO terms |
Description |
p-value |
Genes |
GO:0002250 |
adaptive immune response |
1.58e−08 |
CD40LG, STAT5A, FLT3LG, CCL22, CD1C, LTB, FOXP3, TNFRSF17, SLAMF1, CLEC10A |
GO:0019221 |
cytokine-mediated signaling pathway |
2.23e−04 |
CCL22, STAT5A, TNFRSF17 |
GO:0045087 |
innate immune response |
2.11e−03 |
CCL22, SLAMF1, CLEC10A |
CD40LG是TNF超家族成员,能引起细胞死亡的细胞因子的一个超家族。据报道,CD40LG的表达水平与肿瘤微环境的免疫活性相关,CD40LG的低表达与BRCA肿瘤晚期和预后不良有关,高表达与BRCA患者的生存率呈正相关。因此,CD40LG可作为BRCA患者预后的潜在生物标志物[23]。
FOXP3是Forkhead-box-P蛋白家族成员,是调节性T细胞基因表达的关键调节因子。据报道,FOXP3在乳腺癌和前列腺癌患者样本中低表达,是一种潜在的肿瘤抑制基因。FOXP3的高表达与BRCA的良好预后有关,这可能与Toll受体通路、JAK/STAT通路、细胞周期和细胞凋亡有关[24]。
CCL22是趋化因子家族成员,主要由树突状细胞和巨噬细胞分泌产生,通过与细胞表面趋化因子受体CCR4相互作用,激发其对靶细胞产生影响。研究发现,CCL22在BRCA患者中高表达,CCL22肿瘤分泌与BRCA分子亚型的预后相关[25]。因此,CCL22是BRCA患者的预后预测因子,具有作为BRCA免疫治疗靶点的潜力。
STAT5A是信号转导和转录激活STAT家族成员,STAT相关信号传导在调节某些生物学过程起着至关重要的作用,包括增殖、转移、炎症和免疫反应。研究发现,与健康对照组相比,BRCA患者中STAT5A的表达显著下调,STAT5的高表达表征较好的BRCA预后[26]。
综合文献调查,这些枢纽基因参与了免疫反应和细胞因子介导相关的生物学过程,其表达水平与BRCA患者的生存预后有显著的相关性。进一步证实本文所提的GAN-1DCCox模型能够有效地提取具有生物学意义的预后特征基因。
4. 结论
本文提出了一种基于多组学融合和生成对抗网络GAN的生存分析GAN-1DCCox模型。通过生成对抗网络GAN实现多组学的特征提取,结合轻量级的1D-CNNCox生存分析架构,筛选出重要的预后特征基因。相比目前流行的生存分析模型,本文提出的GAN-1DCCox模型能够更有效地融合多组学数据,在8种不同癌症类型数据集上均取得了优越的生存预测性能。在未来的工作中,我们将探索该模型在蛋白质组和代谢组数据上的可行性,进一步改善模型的生存预测性能。
基金项目
陕西省自然科学基础研究计划面上项目(2022JM-026)。
NOTES
*通讯作者。