0   引言

肺癌的发病率和死亡率在世界范围内居高不下，其中肺腺癌（Lung adenocarcinoma, LUAD）是最常见的病理类型，早期症状不明显，易被忽视，确诊时常处于晚期^[1]。近年来，尽管LUAD的诊断技术和放化疗、免疫治疗、靶向治疗等治疗手段不断进步，但LUAD患者的5年生存率约为20%，预后不佳^[2]。早期发现、及时干预、合理的预后评价是影响肺腺癌预后的关键因素。因此，探索影响LUAD患者预后的生物标志物，开发新的模型来更准确地预测患者的预后情况至关重要。

RNA存在100余种化学修饰，甲基化是主要的修饰形式，包括N1-甲基腺苷（N1-methyladenosine, m1A）、5-甲基胞嘧啶（5-methylcytosine, m5C）、N6-甲基腺苷（N6-methyladenosine, m6A）和7-甲基鸟苷（N7-methylguanosine, m7G）等^[3]。来自供体的一个活性甲基被添加到RNA中胞嘧啶碱基的碳5位置以形成m5C修饰，是一种广泛存在的RNA修饰，在mRNA和非编码RNA中都有检测到，包括tRNA、rRNA、snRNA和eRNA等。m5C修饰主要由甲基转移酶、去甲基化酶和识别蛋白介导^[4]。

tRNA中的m5C修饰主要参与维持稳态、优化密码子-反密码子配对、调节应激反应和控制翻译效率。mRNA中的m5C则与mRNA稳定、剪接和核质穿梭、DNA损伤修复、增殖和迁移、干细胞的发育、分化和重编程等多种生物学进程有关。mRNA m5C的异常修饰与多种疾病相关，包括动脉硬化、自身免疫性疾病和癌症等^[5]。近年来，有报道称mRNA m5C的异常修饰参与膀胱癌、胃癌和食管鳞状细胞癌的发生发展过程^[6-7]。目前m5C相关基因对LUAD预后影响的研究较少，本研究就此构建了LUAD的预后预测模型，评估了新模型的预后预测价值，同时进行了免疫浸润分析，以期为LUAD的临床治疗提供新的参考。

1   资料与方法

1.1   数据收集与整理

TCGA和GEO（GSE30219、GSE31210、GSE50081）数据库中LUAD患者的临床资料及转录谱数据获取同课题组前期研究^[8]。其中TCGA、GSE30219、GSE31210、GSE50081数据库中分别包含LUAD样本535、124、226、129例，正常样本59、14、20、0例。在进行预后分析时，均剔除随访资料不完整、随访时间≤30天的样本。从人类基因数据库官网（https://www.genecards.org/）中检索并下载获取136个m5C修饰相关基因。

1.2   差异表达基因分析

采用DESeq2包对TCGA数据库中的数据进行差异表达基因分析，筛选标准为|log₂（差异倍数）|>1且调整后P<0.05。

1.3   富集分析

通过基因本体功能（Gene ontology, GO）分析预后相关m5C修饰相关基因参与的生物过程（Biological process, BP）、细胞组分（Cellular component, CC）、分子功能（Molecular function, MF），通过京都基因和基因组百科全书（Kyoto Encyclopedia of Genes and Genomes, KEGG）分析m5C修饰相关基因参与的主要信号通路。利用提供的分子数值通过GOplot包计算每个富集条目对应的z得分值，使用ggplot 2包对富集分析结果进行可视化。

1.4   预后模型的建立

参照本研究课题组之前建立的流程与方法^[8]，使用caret包对TCGA数据集按照7∶3的比例进行随机分层抽样，获得训练集（346例）和验证集（148例）。使用survminer包对训练集行单因素Cox回归分析，初步筛选预后相关的m5C修饰相关基因，为避免遗漏，设定P<0.2做进一步分析。随后使用glmnet包进行最小绝对收缩和选择算子（Least absolute shrinkage and selection operator, Lasso）回归分析缩小候选基因范围。最后通过多因素Cox回归分析建立预后模型。在TCGA全集、GSE30219、GSE31210和GSE50081中根据预后模型分别计算患者的风险评分：风险评分= $ {\displaystyle\sum }_{i=1}^{n}exp_i\mathrm{*}Coef_i $（其中exp和Coef分别代表DRGs的表达量和风险系数），并以各数据集中风险评分的中位值为界，将患者分为高风险组和低风险组。

1.5   预后模型的验证

预后模型的验证方法同李梦涵等^[9]的研究。使用survival包对TCGA训练集、验证集及全集中LUAD患者进行Kaplan-Meier生存分析。使用timeROC包绘制ROC曲线，评估该模型对LUAD患者1、3、5 年生存期的预测效能。在GSE30219、GSE31210和GSE50081这3个测试集中重复上述操作进行外部验证。

1.6   列线图、校准曲线和决策曲线绘制

在TCGA数据集中，纳入年龄、性别、吸烟、T分期、N分期、肿瘤分期和风险评分进行多因素Cox回归分析。纳入LUAD患者独立的预后因素，使用rms包绘制列线图，同时绘制校准曲线用于评估实际生存率和预测生存率的一致性。使用dcurves包绘制决策曲线（Decision curve analysis, DCA），观察列线图模型预测患者预后的优越性。

1.7   免疫浸润分析

在TCGA数据集中，通过ESTIMATE算法评估高风险组和低风险组肿瘤纯度、基质评分、免疫评分和综合评分的差异，通过TIMER算法评估高风险组和低风险组6种免疫细胞丰度的差异。从文献^[10]中收集常见免疫检查点，比较高风险组和低风险组两组间36种免疫检查点表达的差异。

1.8   统计学方法

采用R 4.3.1软件进行数据分析及绘图。两组间计量资料比较采用Wilcoxon秩和检验，多组间比较采用Kruskal-Wallis法。P<0.05为差异有统计学意义。

2   结果

2.1   肺腺癌特异性m5C修饰相关基因筛选

在TCGA数据库中筛选出9 878个基因在LUAD中表达上调，2 843个基因在LUAD中表达下调（P<0.05）。 在GSE31210数据集中筛选出329个基因在LUAD中表达上调，559个基因在LUAD中表达下调。将上述基因与136个m5C修饰相关基因取交集，在TCGA数据集中，共获得75个表达上调的m5C修饰相关基因和21个表达下调的m5C修饰相关基因，在GSE31210数据集中共获得34个表达上调的m5C修饰相关基因和18个表达下调的m5C修饰相关基因。两个数据集中共同表达上调的m5C修饰相关基因有30个，共同表达下调的有13个。

为进一步探索筛选出的差异化表达基因在LUAD中潜在的生物学功能和信号通路，我们进行了GO/KEGG分析。在生物过程方面，主要与高分子甲基化、甲基化、RNA甲基化、ncRNA加工、RNA修改等有关；在细胞组分方面，主要与内溶酶体膜、丝氨酸苏氨酸蛋白激酶复合物、内溶酶体、蛋白激酶复合物等有关；在分子功能方面，主要与s-腺苷蛋氨酸依赖的甲基转移酶活性、甲基转移酶活性、转移酶活性、转移单碳基团、RNA甲基转移酶活性、作用于RNA的催化活性等有关。同时也与p53信号通路、肿瘤中的microRNA、卡波西肉瘤相关疱疹病毒感染、半胱氨酸和蛋氨酸代谢、乙型肝炎等通路相关，见图1。

图  1  肺腺癌特异性表达m5C修饰相关基因的筛选

Figure  1  Screening of specifically-expressed m5C modification-related genes in LUAD

A: volcano plots revealed differentially expressed genes between LUAD tissues and normal tissues in TCGA and GSE31210 cohorts; B: heatmap revealed the expression level of m5C modification-related genes in TCGA and GSE31210 cohorts; C, D: GO and KEGG enrichment analyses for m5C modification-related genes in TCGA cohort.

下载: 全尺寸图片 幻灯片

2.2   基于m5C修饰相关基因的预后模型构建

通过单因素Cox回归分析，初步筛选出13个m5C修饰预后相关基因，见图2A。然后通过Lasso回归分析，设制十折交叉验证获取最优模型，最终获得8个与肺腺癌预后独立相关的m5C修饰相关基因：细胞周期蛋白依赖性激酶 1（Cyclin-dependent kinase 1, CDK1）、细胞周期蛋白依赖性激酶抑制因子 1A（Cyclin-dependent kinase inhibitor 1A, CDKN1A）、核仁蛋白2（Nucleolar protein 2, NOP2）、核糖核苷酸还原酶调节亚基M2（Ribonucleotide reductase regulatory subunit M2, RRM2）、T细胞白血病/淋巴瘤基因6（T cell leukemia/lymphoma 6, TCL6）、Toll样受体8（Toll-like receptor 8, TLR8）、tRNA天冬氨酸甲基转移酶（tRNA aspartic acid methyltransferase 1, TRDMT1）、YTH N6-甲基腺苷RNA结合蛋白2（YTH N6-methyladenosine RNA binding protein 2, YTHDF2），它们共同用于构建预后模型，其中CDK1、CDKN1A、NOP2、RRM2是预后的危险因素，TCL6、TLR8、TRDMT1、YTHDF2是预后的保护因素，见图2。风险评分=0.117*expCDK1+0.034*expCDKN1A+0.045*expNOP2+0.051*expRRM2−0.069*expTCL6−0.045*expTLR8−0.147*expTRDMT1−0.034*expYTHDF2。

图  2  预后模型构建及风险评分分布

Figure  2  Construction of the prognostic model and distribution of risk scores

A, B: univariate Cox and Lasso regression analyses of prognostic m5C modification-related genes in TCGA cohort; C: distribution of risk scores, survival status, and expression of m5C modification-related genes for each patient in TCGA cohort.

下载: 全尺寸图片 幻灯片

2.3   预后模型效果评估

使用Kaplan-Meier生存曲线分析高、低风险组的预后情况。结果显示，高风险组的总体生存率在TCGA训练集、验证集、全集及GSE30219、GSE31210、GSE50081三个测试集中，均比低风险组短，预后较差（均P<0.05），见图3。

图  3  TCGA训练集、验证集、全集及GEO 3个测试集预后模型的评估（Kaplan-Meier曲线）

Figure  3  Evaluation of prognostic model using TCGA training, validation, full, and three GEO test cohorts shown by Kaplan-Meier curve

A–F: survival analysis of high- and low-risk patients in TCGA training, validation, all, GSE30219, GSE31210, and GSE50081 cohorts.

下载: 全尺寸图片 幻灯片

使用ROC曲线评估模型的预测效能，结果提示该模型预测LUAD患者的1年、3年、5年存活率的敏感度和准确性较高，具有较好的预测效能，见图4。

图  4  TCGA训练集、验证集、全集及GEO 3个测试集对预后模型的评估（ROC曲线）

Figure  4  Evaluation of prognostic model using TCGA training, validation, full, and three GEO test cohorts shown by ROC curve

A–F: 1, 3, and 5-year ROC curves in the TCGA training, validation, all, GSE30219, GSE31210, and GSE50081 cohorts.

下载: 全尺寸图片 幻灯片

2.4   列线图、校准曲线和决策曲线分析

在TCGA数据集中，多因素Cox回归分析结果显示T分期、N分期、风险评分是LUAD患者独立的预后危险因素（均P<0.05），见图5A。纳入上述因素构建列线图模型，见图5B。校准曲线显示，该列线图模型能准确预测患者的3年生存预后情况，见图5C。决策曲线显示，该列线图模型的预测效能，远远高于T分期、N分期及风险评分等因素的预测效能，见图5D。

图  5  列线图的构建和验证

Figure  5  Construction and validation of nomogram

A: multivariate Cox analyses of the prognostic risk score and clinical characteristics in the TCGA cohort; B: nomograms for predicting the 1-, 3-, and 5-year survival rates of patients; C: calibration curve to verify the accuracy of the nomogram; D: DCA curves of the OS nomograms of patients.

下载: 全尺寸图片 幻灯片

2.5   免疫浸润分析

在TCGA数据集中，评估高风险组和低风险组间肿瘤微环境、6种免疫细胞丰度和36种免疫基因表达的差异。高风险组患者的基质评分、免疫评分和ESTIMATE评分明显低于低风险组，说明高风险组患者肿瘤微环境中基质细胞和免疫细胞含量更低，肿瘤细胞含量更高，见图6A。高风险组中巨噬细胞、B细胞、中性粒细胞、CD4+T细胞、CD8+T细胞、树突细胞浸润显著减少，见图6B。两组免疫检查点表达差异结果表明，CD160、TNFSF14、CD200R1、CD80、ICOS、TNFSF15、CTLA4、CD28、CD200、ADORA2A、ICOSLG、CD48、CD86、HAVCR2、LAIR1、TNFRSF14、LGALS9、CD44在低风险组患者中高表达，而TNFSF4、TNFSF9、TNFRSF18、CD276在高风险组患者中高表达，见图6C。

图  6  高、低风险组免疫浸润分析

Figure  6  Comparison of immune infiltration between high- and low-risk groups

*: P<0.05, **: P<0.01, ***: P<0.001, ****: P<0.000 1, ns: P>0.05. A, B: ESTIMATE and TIMER analysis showed the differences of immune cell infiltration between high- and low-risk groups; C: difference analysis of immune checkpoints between high- and low-risk groups.

下载: 全尺寸图片 幻灯片

2.6   风险评分与临床病理特征的关系

在TCGA和GSE30219数据集中的分析结果表明，LUAD组织中的风险评分均显著高于正常组织，LUAD患者的T分期、N分期越晚风险评分越高，见图7A~C。TCGA和GSE31210数据集的结果表明，吸烟患者风险评分高于不吸烟患者，EGFR野生型患者的风险评分显著高于EGFR突变型患者，见图7D~E。

图  7  风险评分与患者临床病理特征、吸烟、EGFR基因突变的联系

Figure  7  Relationship between risk score and clinicopathological features, smoking, and EGFR mutation

A-C: relationship between risk score and clinicopathological features in TCGA and GSE30219 cohorts; D: relationship between risk score and smoking in TCGA and GSE31210 cohorts; E: relationship between risk score and EGFR mutation in TCGA and GSE31210 cohorts.

下载: 全尺寸图片 幻灯片

3   讨论

LUAD是肺癌的主要亚型之一，具有高发病率和高死亡率的特点。随着检验、影像等诊断技术的不断进步和治疗手段的改善，LUAD患者的整体预后尽管较前有所提升，但其5年生存率仍然较低^[11]。LUAD作为一种异质性较强的肿瘤，其发病机制复杂，因此，寻找新的分子标志物和潜在的治疗靶点及建立评估其预后的模型成为当前LUAD研究的重要方向之一。

本研究以m5C修饰相关基因为基础，最终构建了基于CDK1、CDKN1A、NOP2、RRM2、TCL6、TLR8、TRDMT1、YTHDF2等8个基因的预后风险模型，并在3个独立的外部数据集上进行了验证。Cox回归分析显示，风险评分是LUAD患者独立的预后危险因素；纳入肿瘤大小、淋巴结转移及风险评分的列线图模型，可以更准确和更稳定地预测患者的预后状态。

此外，我们还进行了免疫浸润和临床特征相关性分析，以研究该模型对LUAD增殖、转移和免疫微环境的影响。CDK1主要负责G₂/M期的转换，与细胞增殖密切相关^[12]，m5C修饰通过调节CDK1 mRNA的稳定性，促进其在癌细胞中的表达，导致细胞周期加速，增加肿瘤细胞的增殖能力。既往研究表明，CDK1的异常表达与肺癌、乳腺癌、胃癌等多种癌症的恶性进展相关，其高表达往往与不良预后相关^[13]。抑制CDK1可能有助于减缓肿瘤的生长，因此CDK1被认为是一个潜在的抗癌靶点。而CDKN1A是细胞周期抑制因子，通常通过抑制CDK1和其他CDKs来阻止细胞从G₁期进入S期，在细胞应对DNA损伤时发挥作用^[14]。m5C修饰可以调节CDKN1A基因的表达，导致其功能失效或表达降低，从而促使肿瘤细胞逃避细胞周期检查点，加速癌症进展。NOP2属于m5C甲基转移酶家族，能够在rRNA中催化m5C修饰。研究表明NOP2在肝癌中高表达，其可以通过催化c-Myc mRNA的m5C甲基化上调c-Myc的表达，激活肝癌细胞的糖酵解过程，为肿瘤细胞提供生长优势，促进肝癌的发生发展^[12]。NOP2敲除可抑制肝癌细胞的增殖和转移，为肝癌的靶向治疗提供潜在新策略。核苷酸还原酶（RR）参与核苷酸代谢，是DNA生物合成、修复和复制所必需的限速酶，哺乳动物的RR包含两种不同的同源二聚体亚基，即大亚基RRM1和小亚基RRM2。RRM2的表达受到m5C修饰的调控，在多种癌症中表现为上调，其过表达不仅促进了肿瘤细胞的DNA合成和增殖能力，还与肿瘤免疫、化疗耐药等密切相关^[15]。RRM2在前列腺癌、肺腺癌等癌症中被认为是铁摄取相关代谢基因^[16-17]，抑制RRM2的活性有望作为一种抗癌治疗策略。TCL6是一种相对较新的癌症相关基因，最早在T细胞白血病患者中被发现，尽管其具体功能尚不完全明确，但有研究发现其可能在细胞周期和细胞分化过程中发挥作用^[18]。TLR8是一种重要的模式识别受体，在多种免疫细胞中均有表达，包括单核细胞、巨噬细胞和调节性T细胞等。TLR8能够逆转调节性T细胞免疫抑制功能，具有较强的抗肿瘤作用。TLR8激动剂可导致多种促炎细胞因子的释放，如IL-6、IL-12、TNF-α等，激活先天性免疫^[19]，m5C修饰可能通过调控TLR8的表达或活性，影响肿瘤细胞的免疫逃逸和微环境变化。根据TLR8作用机制，可与免疫检查点抑制剂和化疗药物等联合使用来提高疗效。既往研究中，TRDMT1被认为是直肠癌患者的保护因子^[20]。TRDMT1是m5C修饰主要的甲基转移酶，负责tRNA中的m5C修饰，通过稳定tRNA分子，维持蛋白质合成的正常进行，从而促进肿瘤细胞的快速增殖。YTHDF2是m6A（N6-甲基腺嘌呤）的主要阅读蛋白之一，能够识别并结合m6A修饰的蛋白，与胃癌、肺癌、肝癌等多种肿瘤密切相关，其在不同肿瘤细胞中调控机制不同，亦可能发挥双重作用，既能促进肿瘤细胞增殖和转移，也能抑制肿瘤细胞生长^[21]。

本研究加深了我们对LUAD发病机制中RNA甲基化相关调节机制的认识，构建出基于m5C修饰相关基因的LUAD预后预测模型，并通过内部和GEO外部验证集验证了该风险模型的良好预后评估效能。该模型中关键的m5C修饰相关基因在LUAD的发生和发展中发挥了重要的作用，有望成为潜在的治疗靶点。然而本研究仍存在一定的局限性：首先，数据来源于公共数据库，可能存在数据质量和样本代表性的问题；其次，仅针对LUAD进行了分析，未涉及其他的肺癌类型，模型的适用范围有限。未来，我们将收集临床样本进行验证，并通过分子生物学、细胞学等方法进一步探讨m5C修饰相关基因影响肿瘤发生、发展的作用与机制。