嘌呤代谢相关基因预测肺腺癌预后模型的构建及验证

李杰 李强

肺腺癌是一种常见的恶性肿瘤，发病率和死亡率都很高[1]，尽管治疗肺腺癌技术已取得了巨大进展，但肺腺癌患者的5 年总生存期（overall survival，OS）仍低于15%[2]。因此，探索影响肺腺癌患者预后的生物标志物非常有必要。嘌呤不仅能提供必要的能量和辅助因子以促进细胞存活和增殖，而且还参与免疫反应和宿主-肿瘤相互作用[3]。研究发现，嘌呤代谢与肿瘤的发生、发展密不可分[4]，其不仅在前列腺癌细胞的增殖和侵袭中发挥着重要的作用[5]，而且还影响着膀胱癌和乳腺癌的发生、发展[6]。但目前有关嘌呤代谢与肺腺癌关系的报道很少，因此，本研究基于癌症基因组图谱（the cancer genome atlas，TCGA）数据库中肺腺癌患者的数据，通过生物信息学技术建立预后模型，探讨嘌呤代谢相关基因与肺腺癌的预后关系。

1 材料和方法

1.1 数据收集 从TCGA 数据库官网（https://www.cancer.gov/tcga）中获取肺腺癌临床数据和mRNA 转录组数据，其中肺腺癌mRNA 转录组样本数据包含535个肺腺癌组织和59 个癌旁组织，从基因表达数据库（gene expression database，GEO）官网（http://www.ncbi.nlm.nih.gov/geo）中获取GSE26939 的临床数据和mRNA转录组数据，排除没有生存状态或随访时间＜1 d 的样本。此外，从人类基因（Genecards）数据库官网（https://www.genecards.org/）获取嘌呤代谢相关基因3 968 个。

1.2 差异基因筛选 使用R 4.1.5 软件中的“limma”包，以|log2差异倍数（fold change,FC）|＞2 且P＜0.05 为标准，筛选表达差异的嘌呤代谢相关基因，再使用“ggplot2”和“pheatmap”包，绘制成火山图和热图进行可视化。

1.3 预后差异基因富集分析 将上述差异基因与肺腺癌临床数据合并，筛选出与肺腺癌预后有关的嘌呤代谢相关基因，然后使用R 软件中的“clusterProfiler”和“enrichplot”包对预后相关的嘌呤代谢相关基因进行京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）和基因本体论（gene ontology，GO）富集分析，P＜0.05 为差异有统计学意义[7]。

1.4 预后模型的建立及分析 首先，使用Perl 软件将上述差异基因与肺腺癌临床数据进行整合，剔除无生存状态记录的患者；其次，使用R 软件中的“survival”和“forestplot”包进行单因素及多因素独立预后Cox 回归分析，筛选出P＜0.05 的嘌呤代谢差异基因；然后使用“survival”和“glmnet”包对上述筛出的差异基因进行套索（the least absolute shrinkage and selection operator，LASSO）回归分析，建立风险评分预后模型。风险评分公式为：风险评分与预后有关的嘌呤代谢相关基因的表达量（i）（Coef 为基因在多因素Cox 回归分析中的回归系数，n为与预后有关的嘌呤代谢相关基因的总数目），计算每例患者的风险评分，根据风险评分中位值，将患者分为高风险组和低风险组，再使用“survivalROC”包绘制ROC 曲线，并计算约登指数，从而预测模型的准确性[6]。

1.5 预后模型的评价 使用R 软件中的“survival”包绘制高、低风险组评分分布曲线、生存状态图及建模基因表达热图。使用“survminer”包绘制Kaplan-Meier生存曲线，评估高、低风险组患者的OS。使用“timeR‐OC”包绘制时间依赖性ROC 曲线，计算1、3、5 年肺腺癌患者OS 的AUC，以评价模型的准确性，AUC 越高代表模型越准确。

1.6 独立预后因素分析 将性别、年龄、肿瘤分期、肿瘤类型、吸烟和风险评分作为主要变量进行单因素及多因素Cox 回归分析，使用R 软件中的“forestplot”包绘制森林图来分析影响肺腺癌患者的预后因素。使用“rms”包建立列线图和决策曲线图来估计1、3、5 年肺腺癌患者的总复发率。

1.7 预后模型验证 从GEO 数据库中选择GSE26939为外部检验集，根据风险评分公式，计算验证集中每个样本的风险评分。根据风险评分中位值，将样本分为高风险组和低风险组，使用Kaplan-Meier 生存曲线、ROC 曲线的AUC、生存状态图及建模基因表达热图来检验预后模型的准确性。

2 结果

2.1 差异基因分析 从3 968 个嘌呤代谢相关基因中，共筛选出355 个有表达差异的嘌呤代谢相关基因，其中155 个上调，200 个下调，见图1（插页）。

图1 嘌呤代谢相关基因差异分析（A：差异基因的热图；B：差异基因的火山图）

2.2 功能富集分析 KEGG 结果显示，嘌呤代谢相关基因主要富集在细胞周期、黄体酮介导的卵母细胞成熟、卵母细胞减数分裂、细胞衰老、p53 信号通路、神经活性配体-受体相互作用、人T 细胞白血病病毒1 感染、胆汁分泌、唾液分泌和造血细胞谱系等通路上（均P＜0.05），见图2A（插页）。分子功能富集分析显示，嘌呤代谢相关基因与蛋白质丝氨酸/苏氨酸激酶活性、作用于DNA 的催化酶活性、组蛋白激酶活性、蛋白酪氨酸激酶活性、蛋白质丝氨酸/苏氨酸/酪氨酸激酶活性、ATP 酶活性、单加氧酶活性、G 蛋白偶联肽受体活性、肽受体活性、DNA 依赖性ATP 酶活性等均有关（均P＜0.05），见图2B（插页）。生物学过程富集分析显示，嘌呤代谢相关基因与核分裂、细胞器分裂、有丝分裂核分裂、有丝分裂核分裂的调节、核分裂调节、有丝分裂中期/后期转变的调节、染色体分离的调节、细胞周期中期/后期转变的调节、有丝分裂细胞周期的中期/末期过渡、细胞周期的中期/后期过渡等均有关（均P＜0.05），见图2C（插页）。细胞定位富集分析显示，嘌呤代谢相关基因与染色体区域、染色体着丝粒区、浓缩染色体、浓缩染色体着丝粒区、浓缩染色体着丝粒、浓缩染色体外动着丝粒、凝聚核染色体、凝聚核染色体着丝粒区、着丝粒、纺锤体等均有关（均P＜0.05），见图2D（插页）。

图2 功能富集分析（A：通路富集分析；B：分子功能富集分析；C：生物学过程富集分析；D：细胞定位富集分析）

2.3 肺腺癌预后模型预测患者预后 对上述355 个有表达差异的嘌呤代谢相关基因进行单因素Cox 回归分析，得到24 个与肺腺癌患者预后相关的基因，见表1。使用LASSO 回归对这24 个基因进行交叉验证分析，得到7 个与预后密切相关的基因，见图3（插页）。将这7个基因进行多因素Cox 回归分析，最终得到5 个与肺腺癌患者预后显著相关的基因（CD19、CYP17A1、KH‐DRBS2、INHA、PLK1），其中INHA 和PLK1 基因HR值＞1，为高风险基因，表示INHA 和PLK1 高表达对肺腺癌的预后较差；CD19、CYP17A1 和KHDRBS2 基因HR值＜1，为低风险基因，表示CD19、CYP17A1 和KH‐DRBS2 低表达对肺腺癌的预后较差，见表2。根据5个嘌呤代谢相关基因的β值和基因的表达量计算每例患者的风险评分，风险评分=（-0.542×CD19 表达量）+（-0.542×CYP17A1 表达量）+（-0.463×KH‐DRBS2 表达量）+（0.134×INHA表达量）+（0.129×PLK1表达量）。

表1 单因素Cox 回归分析结果

表2 多因素Cox 回归分析结果

图3 预后风险评分模型评估（A、B：套索回归分析）

2.4 预后模型的预测价值 5 个基因的生存分析见图4（插页）。根据风险评分的中位值（1.77），将患者分为低风险组和高风险组。Kaplan-Meier 生存曲线（红色为低风险组，蓝色为高风险组）显示，低风险组患者OS高于高风险组（HR=3.85，95%CI：2.79～5.31，P＜0.01），见图4A（封三）。生存状态图及建模基因表达热图（红色为高风险组，蓝色为低风险组）显示，高风险组患者预后较差，见图4B（封三）。ROC 曲线分析显示，1、3、5年肺腺癌患者OS 的AUC 分别为0.76、0.74、0.77，见图4C（封三）。综上表明，此预后模型对肺腺癌患者的预后预测准确性很高。

图4 预后模型的预测价值（A：Kaplan-Meier 生存风险曲线；B：生存状态图和建模基因表达热图；C：时间依赖性ROC 曲线）

2.5 独立预后因素分析 单因素Cox 回归分析显示肿瘤分期和风险评分均是肺腺癌患者预后的独立危险因素（均P＜0.05），见图5A。 将上述两个危险因素纳入多因素Cox 回归分析，结果显示肿瘤分期和风险评分均是肺腺癌患者预后的独立危险因素（均P＜0.05），见图5B。列线图得分可推测患者未来1、3、5 年的生存率，见图5C。决策曲线显示其预测效能较好，见图5D。

图5 独立预后生存分析（A：单因素Cox 回归分析；B：多因素Cox 回归分析；C：列线图；D：决策曲线）

2.6 外部数据集验证 从GEO数据库中选择GSE26939为外部验证集，GSE26939 数据集的风险评分中位值为0.99，通过验证集的Kaplan-Meier 生存曲线（红色为低风险组，蓝色为高风险组）显示，低风险组患者OS 高于高风险组（HR=8.94，95%CI：2.96～27.01，P＜0.01），见图6A（封三）。生存状态图及建模基因表达热图（红色为高风险组，蓝色为低风险组）显示，高风险组患者的预后较差，见图6B（封三）。验证集的ROC 曲线1、3、5 年肺腺癌患者OS 的AUC 分别为0.96、0.82、0.84，见图6C（封三）。综上表明该模型在外部验证集中有很好的预测性能。

图6 预后模型的性能评估（A：Kaplan-Meier 生存分析曲线；B：生存状态图和建模基因表达热图；C：时间依赖性ROC 曲线）

3 讨论

本研究以嘌呤代谢相关基因作为背景，最终建立了由5 个嘌呤代谢相关基因（CD19、CYP17A1、KH‐DRBS2、INHA、PLK1）组成的预后模型。结果显示INHA 和PLK1 高表达对肺腺癌的预后较差，CD19、CYP17A1 和KHDRBS2 低表达对肺腺癌的预后较差。GO 和KEGG 富集分析表明这些基因主要富集在细胞周期和细胞分裂信号通路上，由此可大胆推测，嘌呤代谢相关基因可能通过促进细胞增殖来影响肺腺癌患者的预后。

INHA 是重要的代谢相关基因之一，据报道，INHA的功能不仅与前列腺癌雄激素非依赖性转移和卵巢肿瘤的血管生成有关[8]，还与肺腺癌的免疫浸润有关[9]。PLK1 参与了多种细胞周期调节途径，其可作为G2/M 检查点并负责中心体、纺锤体组装和染色体分离的调节[10]。既往研究发现，PLK1 的过表达与多种癌症的发生、发展密切相关[10-11]。较高的PLK1 转录和蛋白水平对胃癌患者的预后有不良影响[12]。PLK1 高表达也与卵巢透明细胞癌密切相关[13]。还有研究证实，PLK1 的表达可以预测转移性非小细胞肺癌患者的生存[14]，此外，研究发现PLK1 还与肺腺癌的免疫微环境有关[15]。CD19 是重要的B 细胞表面标志物，目前被广泛应用于血液系统的肿瘤治疗中[16]。据报道CD19 可能在肺腺癌的不同阶段发挥双重作用，并且还与免疫预后相关[17]。CYP17A1 可以将睾酮转化为雌二醇，是非小细胞肺癌易感基因之一，但其多态性与亚洲人群中的非小细胞肺癌发展无关[18]。有研究表明CYP17A1较高的基因突变和拷贝数变异可能通过影响B 细胞功能来影响肺腺癌患者的易感性[19]。KHDRBS2 可以作为抑癌基因，既往研究报道，KHDRBS2 表达水平可以预测肺癌患者OS[20]，这与本研究结果一致。

综上所述，嘌呤代谢相关基因可影响肺腺癌患者的预后，有潜在的临床应用价值。但本研究仍有不足之处，如本研究的全部数据均来源于公共数据库，未进行实验研究及临床数据验证，因此，未来的研究需要多中心、大样本来进一步验证嘌呤代谢相关基因预测肺腺癌患者预后的准确性，为今后的治疗提供潜在的治疗靶点。