殷清燕 车露美 刘星宇
(1.西安建筑科技大学 西安 710055)(2.西安邮电大学 西安 710121)
生存分析(Survival Analysis)是研究生存现象和响应时间的一门统计学分支,已被广泛运用于医学领域,如临床试验、疾病诊断和预后分析等[1~4]。生存数据中经常出现删失数据,导致生存时间信息不完全。传统的统计学习方法无法很好地解决生存删失数据的分析。设计适用于生存分析的机器学习算法,具有重要的研究意义和应用价值。
TCGA 是由美国国家癌症研究所和人类基因组研究所合作建立的可公开获取的癌症组学数据库[5]。TCGA 收录了33 种癌症数据,包含四种类型的 生 存 事 件:OS(Overall Survial),PFI(Progression-Free Interval),DFI(Disease-Free Interval),DSS(Disease-Specific Survival)。OS 指从观察开始到死亡的时间;PFI 指从观察开始到肿瘤进展的时间;DSS是指从诊断或治疗开始到死亡的时间;DFI指治疗后无病状态到肿瘤进展的时间,进展事件可能是局部复发和远处转移。本研究首先选择OS和PFI 事件进行初步Kaplan-Meier 估计,OS 事件进行深入的生存分析。
尽管Cox 比例风险模型在生存分析领域非常流行,但是比例风险假设在高维生存数据上难以成立。为了缓解这个问题,研究者们提出了具有较少的约束模型假设的非参数方法。随机生存森林(Random Survival Forset,RSF)是Ishwaran提出的一种适用于生存数据的决策树集成方法[9]。RSF 算法提出了CART 决策树框架下的生存树分裂规则,即将结点分裂成具有不相似生存时间的左右子结点,选取最大化结点间生存差异的特征作为最佳分裂特征。
RSF 算法通过Bootstrap 随机采样数据生成生存树;分割节点时使用随机选择的特征子集;生存树完全生长不进行剪枝;将多个生存树的预测结果进行多数投票得到最终的预测。
表1 梯度提升生存树的算法流程
成分梯度提升(CW-GBM)是2016 年提出的基于boosting和稳定性选择的梯度提升方法[12]。可以提高计算可行性,在稳定性选择中引入随机排列来达到控制假阳性率。
原始的TCGA 数据集包含33 种癌症类型的10915 个患者信息。本研究对原始数据进行数据清洗,得到预处理后的5384 个患者生存数据。生存信息包括OS,OS.time,PFI,PFI.time,临床特征包括年龄、性别、肿瘤状态和病理分期等。
Kaplan-Meier 估计通常是进行生存分析的第一步。根据离散的生存时间估计生存概率,绘制Kaplan-Meier 生存曲线,观察和比较不同组患者的生存状况。根据肿瘤状态特征将患者分为有肿瘤和无肿瘤,分别绘制生存曲线(图1)。上图表示OS事件的生存曲线,无肿瘤的最低生存概率约为0.55,有肿瘤的生存曲线下降非常快;下图表示PFI事件的生存曲线,无肿瘤的最低生存率约为0.72,有肿瘤的生存曲线下降得最快。有肿瘤的生存概率均在2000 天后低于0.1,有肿瘤和无肿瘤的生存曲线之间的差别非常大。有肿瘤的OS生存曲线和PFI 生存曲线具有相同的趋势,PFI 生存曲线下降速度更快。
图1 不同Tumor status下的K-M曲线
根据性别特征将患者分为男性(Male)和女性(Female),绘制生存曲线(图2)。男性患者的OS和PFI 生存曲线均低于女性患者,男性患者的生存概率低于女性患者。男性和女性患者在PFI 生存曲线上的差别更大。
图2 不同性别的K-M生存曲线
根据病理分期特征将患者分为六组。根据癌症分期将患者分为六期,分别绘制曲线(图3)。I/II/III 的最低生存概率分别约为0.51、0.28、0.21,癌症四期(Stage IV)在t=2000 时的生存概率小于0.3。在t=3000 时不同病理分期的生存曲线之间有明显的差别,病理分期等级越高的患者生存概率越低。
图3 不同Stage下的K-M生存曲线
根据不同癌症类型将患者进行分组,绘制生存曲线比较不同癌症类型患者的生存状况(图4)。所有癌症类型的患者在t=2000天(约5.5年)时达到最低生存概率。癌症患者的生存概率与癌症类型,性别,肿瘤状况(Tumor status)和病理分期(Stage)相关。
图4 不同癌症类型的K-M生存曲线
本文使用OrdinalEncoder 对“cancer type”,“tumor_stage”,“tumor_status”进 行 编 码;使 用One-HotEncoder 对“gender”特征进行编码,使用StandardScaler 对“age”进行标准化。采用贝叶斯超参数优化方法对每个算法的超参数进行调整[13]。横坐标表示参数的调参范围,纵坐标表示一致性指数(C-index)值,绘制各个算法模型的超参数优化过程图。
图5 是CW-GBM 算法的调参结果图。随着“n_estimators”的增大,模型的C-index 保持增大。n_estimators 达到300 之后,C-index 指数趋于稳定,选择最优n_estimators=493。图6 和图7 是GBM 和RSF算法的调参结果图。“min_samples_leaf”的调参范围是(1,25),min_samples_split 的调参范围是(2,25),max_features 的调参范围是(1,5),C-index 取值(0.80,0.82)。GBM 选取的最优超参数为min_samples_leaf=9,max_features=2,min_samples_split=11,n_estimators=188。RSF 算法选取的最优超 参 数 为min_samples_leaf=3,max_features=1,min_samples_split=16,n_estimators=21。
图5 CW-GBM模型的超参数优化
图6 GBM模型的超参数优化
图7 RSF模型的超参数优化
C-index 在0.5~1 之间,0.5 为完全不一致,说明该模型没有预测作用,1 为完全一致,说明该模型预测结果与实际完全一致。一般情况下C-index在0.50~0.70 为准确度较低:在0.71~0.90 之间为准确度中等;而高于0.90 则为高准确度。各算法在TCGA 测试集上的C-index 见图8。RSF、GBM、COX、CW-GBM 在测试集上的一致性指数依次为0.87、0.85、0.79 和0.79。RSF 模型在预测测试集上的表现较好,具有较高的准确性。
图8 各算法在TCGA测试集上的C-index
图9 是各算法的time-dependent AUC 曲线图,虚线表示AUC平均值[14]。Cox PH模型的AUC平均值为0.84,最低值为0.65。观察初期的AUC值增长较块,约1000 天后达到AUC 平均值。RSF 和GBM算法的AUC 曲线位于CoxPH 曲线上方,RSF 和GBM算法的性能优于CoxPH。因此,基于贝叶斯超参数优化的RSF和GBM具有较好的性能。
图9 各算法在TCGA数据集上的AUC
选择性能表现优异的随机生存森林和梯度提升树进行癌症患者的生存风险分层。根据预测的风险评分值的四分位数,将患者分为四个风险组:低危组、中低危组、中高危组和高危组[15]。在随机生存森林RSF中选取的风险评分分界点为0~1.84,1.84~8.45,8.45~22.53,22.53~100。在梯度提升树GBM 算法中选取的风险评分分界点为0~15.96,15.96~29.33,29.33~43.54,43.54~100。
RSF 风险评分的四分位数分别为1.84,8.45,22.53,据此划分的低危组1459 例、中低危组1458例、中高危组1167 例,高危组1750 例。GBM 风险评分的四分位数分别为15.96,29.33,43.54,据此划分的低危组1461例、中低危组1458例、中高危组1165 例,高危组1750 例(表2)。RSF 和GBM 四个风险组的生存曲线如图10 所示。可见两个算法中相邻风险组的生存曲线均有显著性的差异。风险评分较高的风险组的生存曲线也较低,表明生存风险分层方法的有效性。
图10 四个风险组的KM生存曲线
表2 TCGA患者的风险分层
结合表3 各风险组在性别、年龄、肿瘤状态等方面的差异。在性别方面,我们可以看到各风险组中女性的数量普遍较低,尤其在RSF 低风险组中,女性人数仅为9 人,而男性人数为484 人。随着风险程度的提高,男性人数逐渐增多,而在GBM 高风险组中,男性人数高达1078 人,女性人数也有667人。这表明,随着年龄的增长,男性患脑肿瘤的风险较高。在年龄方面,我们可以看到各风险组的年龄分布有所不同。在RSF低风险组中,平均年龄为47.98岁,随着年龄的增加,风险组平均年龄也逐渐上升。在GBM 高风险组中,平均年龄为64.7 岁,这说明随着年龄的增长,患脑肿瘤的风险也在增加。
表3 各风险组的信息
针对不同风险组的特征信息进行分析,结果表明肿瘤状态对患者生存风险的作用最显著。随着年龄的增长,男性患脑肿瘤的风险较高;患者的肿瘤状态为“With tumor”时,患者的生存风险较高;其次年龄特征对患者生存风险的影响显著。相比高风险组,低风险组患者的平均年龄较低。因此,针对不同风险组的人群,应采取有针对性的预防措施,以降低脑肿瘤的发病风险。
本研究针对TCGA 癌症患者生存数据,比较CoxPH 模型、随机生存森林RSF、梯度提升树GBM和CW-GBM 算法的生存分析性能。对各个算法进行贝叶斯超参数优化,获得C-index 和time-dependent AUC性能评价值。结果表明,随机生存森林和梯度提升树优于其他模型,在测试集上C-index 为0.87,time-dependent AUC 为0.90。其次,基于RSF和GBM 的生存风险评分,进行癌症患者的生存风险分层,划分为低危组、中低危组、中高危组和高危组。KM 生存曲线的实验结果表明,随机生存森林和梯度提升树算法对TCGA 癌症患者的生存风险具有很强的预测能力,在识别高危和低危患者群体方面有着显著的判别作用。