肺腺癌干性相关分子预后预测模型的构建及验证

蔡 杰， 蔡剑桥， 朱余明， 赵晓刚

(同济大学附属上海市肺科医院胸外科，上海 200433)

肺癌是世界第一大癌，每年全球约有200万的新发病例和170万的死亡病例，是恶性肿瘤致死的首要原因[1-3]。非小细胞肺癌(non-small-cell lung carcinoma, NSCLC)约占所有肺癌病例的85%，其中肺腺癌(lung adenocarcinoma, LAD)是最常见的病理类型，约占全部NSCLC病例的50%[2]。遗憾的是，尽管禁烟、空气治理、肺CT筛查和早期手术等一些早期干预措施的施行，超半数(57%)的NSCLC患者发现时已是中晚期，其5年生存率仅为6%[1]。过去20年，肺癌的治疗进入到传统治疗如手术、放化疗和新兴治疗如靶向治疗、免疫治疗相结合的个体化治疗阶段[4]，及早鉴别不良预后高风险肺癌患者人群并给予个体化综合治疗将极大改善其预后。

因此，本研究在获得TCGA肺腺癌数据集的干性指数mRNAsi基础上，应用加权基因共表达网络分析(WGCNA)筛选得干性相关基因，并进一步根据TCGA肺腺癌预后信息进行LASSO回归分析构建得一个由20个基因组成的风险模型预测肺腺癌的预后。最后，为更好地进行临床应用，建立了基于这一风险模型的列线图，该列线图能更准确有效地及早鉴别不良预后高风险肺腺癌患者，为患者带来临床获益。

1 资料与方法

1.1 TCGA肺腺癌数据集的准备及差异表达分析

在TCGA数据集(https:∥portal.gdc.cancer.gov/)中下载肺腺癌STAR-Counts原始数据，进行数据标准化和低质量数据清理，应用R语言包“DESeq2”进行肺腺癌和癌旁组织之间的差异表达分析，提取其中的差异表达mRNA(P<0.05且差异倍数>1.5)。R语言包“clusterProfiler”被用于进行GO和KEGG通路富集分析。

1.2 WGCNA分析筛选中枢干性相关差异表达基因

应用R语言包“WGCNA”对上诉差异表达mRNA进行WGCNA分析，该分析可根据基因间的加权共表达情况对基因进行富集分类。从原始发表论文获取肺腺癌干性指数mRNAsi数据，对各富集基因集与mRNAsi作相关分析，从3个与mRNAsi相关系数最高的基因集中进一步筛选中枢干性相关差异表达基因，要求符合模块身份(module membership, MM)>0.5且基因显著性(gene significance, GS)>0.5。

1.3 多基因LASSO回归模型的建立

将TCGA肺腺癌数据集随机分成试验组和内部验证组，在试验组中进行多基因LASSO回归模型的构建。首先，对每个候选中枢干性相关差异表达基因进行单因素Cox回归分析并计算P值，当P<0.2时，则纳入进一步分析以构建LASSO回归模型。R语言包“glmnet”返回一系列lambdas(λs)值和风险模型。进一步利用“glmnet”包中的“cv.glmnet”函数进行十折交叉验证，选择具有最小平均交叉验证误差的λ值，该值对应的LASSO模型即是所构建的理想LASSO风险模型。

1.4 模型预测有效性分析

分别在试验组、内部验证组和外部验证组中根据模型公式计算每个患者的风险评分，绘制ROC曲线，灵敏度及特异度之和最大时所对应的风险评分即是将肺腺癌患者分为预后低、高风险组的临界值。ROC曲线下面积，即AUC值，被用来评估不同变量的预测准确性。应用Kaplan-Meier方法绘制生存曲线，用对数秩检验比较不同组的生存差异。P<0.05被认为结果具有统计学意义。

1.5 列线图的构建

在单因素Cox预后分析的基础上，将P<0.1的临床病理因素(包括干性LASSO模型)纳入以构建预测肺腺癌预后的列线图。校准曲线和一致性指数(C指数)被用于评估该列线图的有效性。决策曲线分析被用于评估列线图的临床获益度。R语言包“rms”被用于完成上述分析。

2 结 果

2.1 TCGA肺腺癌差异表达基因筛选

应用R语言包“DESeq2”筛选TCGA肺腺癌数据集中的差异表达mRNA，共获得上调表达基因4 111个和下调表达基因3 899个(图1A)。GO分析结果显示，上调基因主要富集在核分裂、DNA复制和ATP代谢等通路，下调基因则主要富集在细胞黏附、MAPK信号和免疫应答等通路(图1B、1C)。KEGG分析结果显示，上调表达基因主要富集在DNA复制、氨基酸合成和细胞周期等通路，下调表达基因主要富集在细胞黏附、Ca2+信号和cGMP-PKG信号等通路(图1D、1E)。

图1 TCGA肺腺癌差异表达基因筛选

2.2 中枢干性相关差异表达基因的筛选

图2 中枢干性相关差异表达基因的筛选

2.3 肺腺癌干性相关分子预后预测LASSO模型的构建

将TCGA肺腺癌数据集随机均分为试验组和内部验证组，在试验组中进行LASSO模型构建。对上述中枢干性相关差异表达基因进行单因素Cox回归分析，当P<0.2时，则纳入进一步分析以构建LASSO回归模型。R语言包“glmnet”返回一系列LASSO风险模型，每一条曲线即代表一个基因，每一个不同的λ值处所包含的具有非零系数的基因和相应的非零系数即构成该λ值处的一个LASSO风险模型。为进一步挑选最优的风险模型，进行十则交叉验证分析，结果显示当λ=0.038 79，该模型的交叉验证误差最小(图3)。

在该λ值处，共获得20个带有非零风险系数的基因，组成的风险模型公式如下: 风险评分=(-0.821 95×PREX1表达状态)+(-0.475 02×DPH2表达状态)+(-0.465 33×SLC15A2表达状态) +(-0.301 49×ADRB2表达状态)+(-0.310 36×TFB2M表达状态)+(-0.242 76×THSD1表达状态)+(-0.163 58×PRKCE表达状态)+(-0.158 38×STXBP6表达状态)+(-0.145 52×ZNF25表达状态+(-0.143 29×IL11RA表达状态)+(-0.101 92×MEST表达状态)+(-0.075 26×GJC2表达状态)+(-0.048 07×IL34表达状态)+(-0.0128 0×CA5B表达状态)+ (-0.011 42×PREB表达状态)+(0.004 20×DHFR表达状态)+(0.075 43×FAM111B表达状态)+(0.241 86×CDC7表达状态)+(0.392 36×RBMS2表达状态)+(0.752 92×IQSEC3表达状态)(图3)。在该公式中，基因的低表达状态等于0，高表达状态等于1。

图3 肺腺癌干性相关分子预后预测LASSO模型的构建

2.4 该LASSO模型的有效性验证

图4 该LASSO模型的有效性验证

2.5 基于该LASSO模型的列线图构建及其临床获益评估

图5 基于该LASSO模型的列线图构建及其临床获益评估

3 讨 论

肺癌是目前世界上第一高发和致死肿瘤，其中肺腺癌是肺癌最常见病理类型[12]。TNM分期系统是一种确定肿瘤病变范围的分类方法，被广泛用于评估患者的肿瘤分期，指导治疗并进行预后评判[13]。然而，由于该分期系统仅关注临床病理特征而忽视肿瘤分子遗传背景的差异，其并不能保证准确可靠地预测肺腺癌患者的预后[14]。

肿瘤干性是引起肿瘤恶性进展的重要因素，既往研究已表明部分基因可以影响肺腺癌的干性维持，促进肺腺癌恶性进展，影响肺腺癌术后预后。Liu等[15]发现，EHD1结合CD44形成EHD1/CD44/Hippo/SP1正反馈环路，促进肺腺癌细胞干性维持和迁移侵袭，导致肺腺癌患者的不良预后。此外，Xu等[16]发现，SIRPγ连接MST1和PP2A以促进MST1去磷酸化，激活Hippo/YAP信号和CD47表达，导致肺腺癌细胞干性维持和免疫逃脱，高表达SIRPγ预示肺腺癌的不良预后。由于肿瘤分子遗传背景的多样性和系统性，单一分子标志物因为敏感性和特异性不足，往往无法准确预测肿瘤的预后。因此，应用特定的分析方法建立多分子标志物的预测模型，以完善目前的肿瘤分期系统，可以实现及早鉴定肺腺癌不良预后患者，并进行更为密切的随访和早期干预，以提高肺腺癌患者的总体生存。

干性指数mRNAsi由Malta等[8]在2018年提出，其基于TCGA中各肿瘤基因表达数据集，应用机器学习方法计算获得，该指数能有效评估肿瘤的干性程度。本研究首先筛选得TCGA肺腺癌-癌旁间的差异表达基因，分析差异表达基因与mRNAsi的相关性以进一步筛选得中枢干性相关差异表达基因。最后利用R语言包“glmnet”构建一系列LASSO风险模型，十折交叉验证后选择了平均交叉验证误差最小的最优风险模型，该风险模型公式由20个基因和对应的非零系数组成。该风险模型中部分基因已经被报道与细胞干性或分化有关。CDC7被发现可以促进TGF-β介导的平滑肌细胞分化[17]。Fawal等[18]发现，DHFR在脑肿瘤特别是脑肿瘤始动细胞中显著高表达，靶向DHFR可以抑制脑肿瘤始动细胞的干性维持，促进脑肿瘤始动细胞的良性分化，减少脑肿瘤形成。Huang等[19]发现，PRKCE通过抑制ABCB1的转运及PI3K/Akt、MAPK/ERK信号通路激活，进而抑制肾癌细胞的干性维持。牙周韧带干细胞是牙周组织再生的主要细胞来源，Hasegawa等[20]发现，MEST在低-高分化牙周韧带干细胞间显著差异表达，敲低MEST可显著抑制牙周韧带干细胞的增殖及其干细胞标志物如CD105、CD146、p75NTR、N-cadherin和NANOG的表达。

进一步的分析显示，该LASSO风险模型被证明在试验组、内部验证组和外部验证组这3个队列中均能显著区分肺腺癌生存，并且其对肺腺癌1、3和5年预后均有良好预测价值。为更好地实现临床应用，本研究基于风险模型进一步构建列线图，校准曲线表明整合风险模型后的列线图的预测准确度更高，决策曲线分析证实整合风险模型后的列线图能使肺腺癌患者的总体临床收益度显著提升。

不可否认，本研究还存在一定局限性。首先，除了mRNA之外，lncRNA、microRNA、CpG和单核苷酸多态性等其他标志物已被证实有一定的肿瘤预后预测价值[21-24]，更多维数据的整合分析将进一步提高风险模型的预测效率。其次，一些重要的临床病理指标如病理亚型、胸膜浸润、脉管内癌栓以及影像特征等数据无法在TCGA中获得，这可能会降低整合风险模型后列线图的预后预测价值。最后，本部分研究基于公共数据集的分析，未在实际肺腺癌队列中得到验证，下一步应进行前瞻性临床试验以评估该风险模型对肺腺癌预后的实际预测价值。

总之，本研究构建了一个实用的风险模型能够预测肺腺癌的预后，可能是对目前肺腺癌TNM临床分期的有效补充。医生或者患者可以通过这一风险模型公式对肺腺癌预后进行个体化预测，并制定下一步的诊疗计划，使社会和国家的医疗资源得到更为恰当的利用，符合当前个性化医疗的大趋势。