在全球范围内癌症相关死亡中，肺癌仍是其主要原因[1]。肺癌主要包括两种的组织学类型：非小细胞肺癌（non-small cell lung cancer，NSCLC）和小细胞肺癌，NSCLC约占所有肺癌的80%[2]。其中，NSCLC的主要组织学亚型是肺腺癌（lung adenocarcinoma，LUAD），约占肺癌发病率的40%[3]。根据2011年国际肺癌研究协会（International Association for the Study of Lung Cancer，IASLC）/美国胸科学会（American Thoracic Society，ATS）/欧洲呼吸学会（European Respiratory Society，ERS）和2015年世界卫生组织（World Health Organization，WHO）肺肿瘤分类的组织学分型，浸润性肺腺癌（invasive lung adenocarcinoma，IAC）主要被分为5种病理亚型：贴壁型、腺泡型、乳头型、微乳头型以及实体型，并基于其最主要的生长亚型进行了级别的划分：低级别为贴壁为主；中级别为腺泡或乳头为主；高级别为实体或微乳头为主，级别越高，侵袭性越高，预后越差[4-5]。目前，肺癌的治疗主要依赖于组织类型和临床分期，但由于IAC的高度异质性，即使是同样级别组织类型和临床分期的IAC患者预后也不同。近年来，高通量技术的进步使医学研究人员能够从基因层面发现癌症中的遗传改变，并被广泛用于筛选潜在的肿瘤生物标志物[6-7]。同时，机器学习算法已被引入且应用于遗传和基因组数据，并从生物医学数据集中，通过构建不同的预测模型来预测临床结果[8]。因此基于高通量测序数据来检测IAC不同级别病理亚型的基因表达情况，以及结合生物信息学方法系统分析相关差异表达基因及其调控机制，然后运用机器学习算法构建临床风险预后预测模型，对LUAD患者进行个体化、精准化治疗具有重要意义。
在本研究中，我们筛选出与LUAD患者预后相关的特征基因，然后构建基于基因特征的风险预测模型，预测患者的生存状况，为LUAD患者的个体化诊疗提供参考。
1 资料与方法
1.1 获取数据
于2023年1月3日从UCSC Xena网站（https://xenabrowser.net/）下载了全部的LUAD RNA测序（RNA-Seq）数据和临床病理数据，包括性别、年龄、肿瘤TNM分期、病理亚型、生存时间及生存状态等。经过筛选后共纳入具有明确IAC病理亚型的180例样本（中级别包括98例样本，高级别包括82例样本），用以筛选DEGs和分析基因功能。按照2∶1随机分层抽样的方法分为训练组（120例样本）和验证组（60例样本）。通过训练组样本构建风险预测模型，验证组样本评价模型的准确性（图1）。
图1
研究流程
1.2 筛选差异表达基因
使用R-4.2.1软件中的“DESeq2”包筛选中级别与高级别病理亚型间的DEGs，根据P≤0.05和| log2(FoldChange)| >1的阈值对DEGs进行筛选。使用“pheatmap”和“ggplot2”包绘制火山图和热图。
1.3 基因功能和基因集富集分析
使用R-4.2.1软件的“enrichplot”、 “org.Hs.eg.db”和“clusterProfiler”等包对DEGs进行基因本体功能富集分析（Gene Ontology，GO）、京都基因与基因组百科全书通路分析（Kyoto Encyclopedia of Genes and Genomes pathway，KEGG）和基因集富集分析（Gene Set Enrichment Analysis，GSEA）。P≤0.05具有统计学意义，使用“ggplot2”包对结果进行可视化。
1.4 构建风险预测模型
将训练组中120例样本的DEGs表达量和IAC样本生存信息合并，剔除没有生存时间（overall survival，OS）及生存状态（Censor）的患者后，剩余115例IAC样本纳入Cox分析。把通过P≤0.05和| log2(FoldChange)| >1筛选出来的DEGs纳入单因素Cox分析，通过R-4.2.1软件的“survival”包得到具有独立预后相关的基因。用“glmnet”包把独立预后相关的DEGs行LASSO回归分析，得到的结果纳入多因素Cox分析，最终构建临床风险预测评分模型。风险评分（RiskScore）预测模型计算公式：
RiskScore=Exp1×C1+ Exp2×C2+…+ Expn×Cn（Exp：基因的表达量，C：LASSO Cox回归分析得到的回归系数，n：DEGs在多因素Cox分析中筛选出的基因总数目）。
根据上面的公式计算每例样本的RiskScore，以RiskScore中位数为截断值（cut-off）将IAC患者分为低风险组和高风险组，P≤0.05差异具有统计学意义。
1.5 评价预测风险模型和绘制列线图
利用R-4.2.1软件的“pheatmap”包绘制低风险组和高风险组样本的临床病理因素和候选基因表达量，同时应用“survival”、“survminer”等包绘制Kaplan-Meier生存曲线来比较两组的OS，最后通过“timeROC”包绘制时间依赖性的ROC曲线，算出训练组样本OS在1年、3年、5年的AUC值，来评价预测风险模型的准确性。将两组的临床病理因素和RiskScore纳入Cox回归分析，加载“rms”等包把具有独立危险因素的变量绘制成列线图，上部为评分系统，下部为预测系统。患者的1年、2年、3年、5年、10年生存率均可通过评分系统预测，并且利用校准曲线和c指数来显示生存预测的准确性。
1.6 验证预测风险模型
用验证组的样本数据来验证预测模型，通过上述的公式计算出验证组各个样本的RiskScore，以RiskScore中位数为cut-off值分为低风险组和高风险组，应用Kaplan-Meier生存曲线比较两组的OS。同样的，绘制1年、3年、5年的ROC曲线，计算AUC值评估模型预测预后的能力。
2 结果
2.1 筛选差异基因和功能富集分析
从UCSC Xena网站上收集了全部的LUAD mRNA表达数据和临床病理数据，提取了IAC中级别病理亚型98例样本和高级别病理亚型82例样本，以P≤0.05和| log2(FoldChange)| >1为阈值，筛选出782个DEGs，包括327个上调基因和455个下调基因，可视化火山图、热图结果所示（图2）。
图2
不同级别病理亚型相关基因的差异分析
a：IAC中级别和高级别病理亚型差异表达基因火山图（327个上调基因和455个下调基因）；b：IAC中级别和高级别病理亚型差异表达基因热图（|log2(FoldChange)| >3的差异表达基因）
为了探索不同级别病理亚型相关的生物学功能和途径，进一步缩小筛选标准，以P≤0.05和| log2(FoldChange)| >1.5得到280个DEGs，对上述280个DEGs进行GO及KEGG分析（图3a～d）。GO分析显示，DEGs主要富集在免疫反应、物质代谢及酶活性调节等相关的生物学功能；KEGG分析主要富集在细胞色素P450调节物质代谢、神经活性配体-受体相互作用通路等相关生物学功能。把P≤0.05的DEGs纳入GSEA分析（图3e～f），其分析结果主要富集于细胞色素P450相关的物质代谢、抗原的呈递与自然杀伤细胞介导的免疫反应等相关生物学过程，进一步验证了中级别与高级别病理亚型DEGs与免疫反应和物质代谢之间的密切关系。以上功能富集分析的结果解释了IAC不同级别病理亚型之间基因层面的关系，也为以后寻找与预后相关生物学标志物的潜在靶点提供了不同的研究思路。
图3
功能富集分析
a和c：GO富集分析 主要富集在免疫反应、物质代谢及酶活性调节等相关的生物学功能；b和d：KEGG通路富集分析 主要富集在细胞色素P450调节物质代谢、神经活性配体-受体相互作用通路等相关生物学功能；e和f：GSEA富集分析 主要富集于细胞色素P450相关的物质代谢、抗原的呈递与自然杀伤细胞介导的免疫反应等相关生物学过程
2.2 Cox分析、LASSO回归及风险模型构建
把782个DEGs纳入单因素Cox分析，得到66个与独立预后显著相关的基因（P＜0.05），随后进行LASSO回归的降维分析（图4），虚线标注处为logλ最小值，筛选出15个基因，经进一步经多因素Cox分析，最终得到最佳构建风险模型的5个基因（MELTF，C=0.1851、MAGEA1，C=0.0475、FGF19，C=0.0269、DKK4，C=0.2300和C14ORF105，C=0.0758）。根据5个基因的基因表达量和风险系数计算每个样本的RiskScore：
图4
LASSO回归分析
通过 LASSO 回归分析进一步确定与IAC的 OS 相关的15个关键基因
RiskScore=（MELTF表达量×0.185 1）+（MAGEA1表达量×0.047 5）+（FGF19表达量×0.026 9）+（DKK4表达量×0.230 0）+（C14ORF105表达量×0.075 8）。
2.3 检验预测风险模型性能
以训练组IAC患者的RiskScore中位数为cut-off值分为低风险组（L=57）和高风险组（H=58）。不同危险组患者表现出不同的相关临床病理特征（图5a）。绘制Kaplan-Meier生存曲线（蓝色代表低风险组，红色代表高风险组），结果显示两组OS差异具有统计学意义（P<0.01），表明高风险组比低风险组更易具有不良预后（图6a）。以训练组RiskScore绘制ROC曲线和时间依赖性ROC曲线结果所示（图6b～c），RiskScore的ROC曲线AUC值为0.675，时间依赖性ROC曲线1年AUC值为0.893、3年AUC值为0.713、5年AUC值为0.632。表明该预测风险模型具有良好的敏感性和特异性。单因素和多因素Cox分析显示，RiskScore是一个独立的预后因素，可作为IAC患者的独立预后预测因子，将低风险组和高风险组的临床病理因素及RiskScore纳入Cox回归分析，识别出具有独立危险因素的变量，绘制得到的列线图结果所示（图7a），该个体化预测模型可估计IAC患者（1年、2年、3年、5年和10年）生存率。校准曲线中的列线图和实际观测结果在训练组和验证组中显示出准确的重叠，表明具有良好的预测精度（图7b）。并且该列线图模型的c指数为0.78，高于其他单一因素的预测模型（图7c）。
图5
风险评分与临床特征的关系
a～b：分别是训练组和验证组每个患者临床特征（Grade、M、N、T、Gender、Age、RiskScore）和5个代表性基因（MELTF、MAGEA1、FGF19、DKK4、C14ORF105），按风险评分升序排列
图6
预测模型的评估
a：训练组K-M生存曲线（P<0.01）；b：训练组ROC曲线，AUC值为0.675；c：训练组时间依赖ROC曲线，1年AUC值为0.893，3年AUC值为0.713，5年AUC值为0.632；d：验证组K-M生存曲线（P<0.01）；e：验证组ROC曲线，AUC值为0.626；f：验证组时间依赖ROC曲线，1年AUC值为0.705，3年AUC值为0.700，5年AUC值为0.684
图7
个体化预测模型
a：Nomogram图可以准确预测患者1年、2年、3年、5年和10年生存率；b：校准图显示了训练组和验证组中1年、2年、3年和5年预测OS与实际OS之间生存概率的比较；c：采用C-Index评价个体化预测模型的预测效果
2.4 验证组数据集验证模型
将验证组中60例样本的DEGs表达量和临床信息合并，经剔除无OS和Censor的患者后，计算剩余56例患者的RiskScore，以中位数为cut-off值分为低风险组（L=28）和高风险组（H=28），验证组的相关临床病理特征可视化结果所示（图5b）。然后绘制Kaplan-Meier生存曲线、ROC曲线及时间依赖性ROC曲线，二者OS仍具有显著统计学差异（P<0.01），时间依赖性ROC曲线1、3、5年AUC值为0.705、0.700、0.684（图6d～f）。上述结果说明该模型在验证组中仍具有较好的预测性能。
3 讨论
尽管免疫治疗和分子靶向治疗已经取得了很大的进展，但由于缺乏准确的早期诊断和预后生物学标志物，肺癌患者的OS仍然较差，5年生存率仅约15%[9]。不同肺癌亚型具有不同的临床特征和预后，确定特异的预后生物学标志物和新的治疗靶点仍然是至关重要的。
在本研究中，我们首先筛选出IAC患者中级别和高级别病理亚型间DEGs，包括327个上调基因和455个下调基因，并对其进行了功能富集分析，经GO、KEGG和GSEA富集分析表明这些基因主要富集在细胞色素P450相关物质代谢、自然杀伤细胞介导的免疫反应、抗原的呈递和酶活性调节等生物学过程。细胞色素P450（CYP）是存在于动物、植物和微生物中的一个参与物质代谢的酶超家族，其密切参与了肺癌、肝癌、乳腺癌等癌症的相关生物学功能[10]。自然杀伤（NK）细胞是免疫系统中不可或缺的肿瘤杀伤效应细胞，通过依赖的细胞毒性和产生细胞因子进行免疫调节，相关证据表明NK细胞可以有效地防止肺癌的发生和发展[11]。NK细胞介导的免疫反应在肿瘤形成的免疫监测中具有至关重要的作用[12]。通过对中级别和高级别两组病理亚型DEGs的分析，解释了两组之间生物侵袭性和预后等因素基因层面的差异，可能为更深层次的研究提供了一定的科研思路。其次，通过对DEGs进行单因素Cox分析，得到66个与独立预后显著相关的基因，再行LASSO回归分析和多因素Cox分析后，我们最终确定了由5个特征基因（MELTF、MAGEA1、FGF19、DKK4、C14ORF105）组成的风险评估模型，并根据中位RiskScore将样本分为低风险组和高风险组。随后，我们对RiskScore与临床病理特征的相关性进行了分析，结果显示随着复发评分的增加临床病理呈现不对称分布特征。进一步对两组患者的预后进行研究，结果显示，高风险组患者的OS明显短于低风险组（P<0.01）。使用ROC曲线和时间依赖性ROC曲线评价预后预测性能，结果显示，ROC曲线AUC为0.675，时间依赖性ROC曲线1年、3年和5年生存率的AUC值分别为0.893、0.713、0.632，表明模型具有良好的准确性，在验证集（P<0.01，1年AUC=0.705、3年AUC=0.700、5年AUC=0.684）中也得到验证。最后，用RiskScore和两组具有独立预后的临床病理因素绘制出列线图，该模型能精确预测患者1年、2年、3年、5年和10年的生存率，校准曲线显示出准确的重叠，c指数也高于其他单一因素（RiskScore、Grade）的预测模型。
上述5个特征基因已被报道与多种癌症的发生发展密切相关。MELTF又称MTf（Melanotransferrin）或MTF1（Metal regulatory transcription factor 1），是铁结合转铁蛋白的同源物，可作为LUAD和胃癌的预后标志物[13-14]。研究[14-15]表明MTf能促进癌细胞侵袭、迁移、增殖和EMT进展，是一种具有很大潜力的靶点。MAGEA1是属于X染色体簇的癌症/睾丸抗原，在多种肿瘤中表达，如肺癌，它可能是免疫治疗的理想抗原靶点[16]。研究[17]表明，随着肿瘤的进展，MAGEA1的表达水平增加，说明了MAGEA1的检测有可能作为肺癌早期诊断的辅助检查方法。成纤维细胞生长因子（FGFs）家族由22个成员组成，分为7个亚家族，调节着一系列生物学功能和多种发育过程[18]。作为FGF家族内分泌亚群的一员，FGF19在维持蛋白质和葡萄糖代谢方面都起着关键作用[19]。FGF19具有与4种类型的酪氨酸激酶FGF受体（FGFR 1-4）结合的潜力，但亲和力相当低，需要β-klotho （KLB）作为辅受体，一旦FGF19结合到FGFRs和KLB，它就会启动二聚作用和自身磷酸化，最终激活细胞质信号转导通路[19-21]。相关研究[19]表明FGF19在NSCLC中表达明显上调，且FGF19过表达与不良预后显著相关，可能作为预测NSCLC患者不良预后的潜在生物学标志物。DKK4作为DKK家族中的一员，是Wnt/β-catenin信号通路的抑制剂[22]。研究[23]发现DKK4在肾细胞癌、结直肠癌、胃癌、NSCLC和上皮性卵巢癌中表达上调，在肝细胞癌中表达下调，DKK4也参与肿瘤生长、侵袭、迁移和化疗耐药等生物功能。目前已经有很多关于DKK4的研究，然而，DKK4在癌症或非癌症疾病中的作用机制和分子机制仍不明确，许多前沿领域有待探索，所以本文为进行更深层次的研究提供了可能的方向。研究发现C14ORF105在耐紫杉醇的卵巢癌细胞中过表达，提示其可能参与了紫杉醇耐药[24]。到目前为止，C14ORF105在肺癌中的作用尚未被广泛研究。综上所述，MELTF、MAGEA1、FGF19、DKK4和C14ORF105均参与了肿瘤的发生或进展。
总的来说，我们的研究通过Cox分析和LASSO回归分析构建了一个准确有效的5基因组预测模型。基于该5基因组的RiskScore可用于确定LUAD患者的OS。将RiskScore和临床参数相结合的列线图可用于预测LUAD患者1年、2年、3年、5年和10年的生存率，所以，它将有助于LUAD患者的预后和随访监测，为LUAD患者的个体化诊疗提供参考。但需要指出的是，我们的研究还存在一定的局限性。首先，我们的研究数据主要来自TCGA数据集，有必要在大型独立临床队列中评估其预测效能。其次，我们纳入的分组仅包括了中级别和高级别病理亚型，排除了样本量较少的低级别病理亚型，并且需要进一步外部数据集的验证。最后，这5个基因在LUAD发病中的生物学机制有待进一步通过功能研究来阐明。
利益冲突：无。
作者贡献：崔严奇负责病例筛选、数据整理与论文设计、初稿撰写等；曾志勇、廖毅、倪琳、连铎煌、杨鲸蓉、叶仕新和张锦灿负责论文审阅与修改。
在全球范围内癌症相关死亡中，肺癌仍是其主要原因[1]。肺癌主要包括两种的组织学类型：非小细胞肺癌（non-small cell lung cancer，NSCLC）和小细胞肺癌，NSCLC约占所有肺癌的80%[2]。其中，NSCLC的主要组织学亚型是肺腺癌（lung adenocarcinoma，LUAD），约占肺癌发病率的40%[3]。根据2011年国际肺癌研究协会（International Association for the Study of Lung Cancer，IASLC）/美国胸科学会（American Thoracic Society，ATS）/欧洲呼吸学会（European Respiratory Society，ERS）和2015年世界卫生组织（World Health Organization，WHO）肺肿瘤分类的组织学分型，浸润性肺腺癌（invasive lung adenocarcinoma，IAC）主要被分为5种病理亚型：贴壁型、腺泡型、乳头型、微乳头型以及实体型，并基于其最主要的生长亚型进行了级别的划分：低级别为贴壁为主；中级别为腺泡或乳头为主；高级别为实体或微乳头为主，级别越高，侵袭性越高，预后越差[4-5]。目前，肺癌的治疗主要依赖于组织类型和临床分期，但由于IAC的高度异质性，即使是同样级别组织类型和临床分期的IAC患者预后也不同。近年来，高通量技术的进步使医学研究人员能够从基因层面发现癌症中的遗传改变，并被广泛用于筛选潜在的肿瘤生物标志物[6-7]。同时，机器学习算法已被引入且应用于遗传和基因组数据，并从生物医学数据集中，通过构建不同的预测模型来预测临床结果[8]。因此基于高通量测序数据来检测IAC不同级别病理亚型的基因表达情况，以及结合生物信息学方法系统分析相关差异表达基因及其调控机制，然后运用机器学习算法构建临床风险预后预测模型，对LUAD患者进行个体化、精准化治疗具有重要意义。
在本研究中，我们筛选出与LUAD患者预后相关的特征基因，然后构建基于基因特征的风险预测模型，预测患者的生存状况，为LUAD患者的个体化诊疗提供参考。
1 资料与方法
1.1 获取数据
于2023年1月3日从UCSC Xena网站（https://xenabrowser.net/）下载了全部的LUAD RNA测序（RNA-Seq）数据和临床病理数据，包括性别、年龄、肿瘤TNM分期、病理亚型、生存时间及生存状态等。经过筛选后共纳入具有明确IAC病理亚型的180例样本（中级别包括98例样本，高级别包括82例样本），用以筛选DEGs和分析基因功能。按照2∶1随机分层抽样的方法分为训练组（120例样本）和验证组（60例样本）。通过训练组样本构建风险预测模型，验证组样本评价模型的准确性（图1）。
图1
研究流程
1.2 筛选差异表达基因
使用R-4.2.1软件中的“DESeq2”包筛选中级别与高级别病理亚型间的DEGs，根据P≤0.05和| log2(FoldChange)| >1的阈值对DEGs进行筛选。使用“pheatmap”和“ggplot2”包绘制火山图和热图。
1.3 基因功能和基因集富集分析
使用R-4.2.1软件的“enrichplot”、 “org.Hs.eg.db”和“clusterProfiler”等包对DEGs进行基因本体功能富集分析（Gene Ontology，GO）、京都基因与基因组百科全书通路分析（Kyoto Encyclopedia of Genes and Genomes pathway，KEGG）和基因集富集分析（Gene Set Enrichment Analysis，GSEA）。P≤0.05具有统计学意义，使用“ggplot2”包对结果进行可视化。
1.4 构建风险预测模型
将训练组中120例样本的DEGs表达量和IAC样本生存信息合并，剔除没有生存时间（overall survival，OS）及生存状态（Censor）的患者后，剩余115例IAC样本纳入Cox分析。把通过P≤0.05和| log2(FoldChange)| >1筛选出来的DEGs纳入单因素Cox分析，通过R-4.2.1软件的“survival”包得到具有独立预后相关的基因。用“glmnet”包把独立预后相关的DEGs行LASSO回归分析，得到的结果纳入多因素Cox分析，最终构建临床风险预测评分模型。风险评分（RiskScore）预测模型计算公式：
RiskScore=Exp1×C1+ Exp2×C2+…+ Expn×Cn（Exp：基因的表达量，C：LASSO Cox回归分析得到的回归系数，n：DEGs在多因素Cox分析中筛选出的基因总数目）。
根据上面的公式计算每例样本的RiskScore，以RiskScore中位数为截断值（cut-off）将IAC患者分为低风险组和高风险组，P≤0.05差异具有统计学意义。
1.5 评价预测风险模型和绘制列线图
利用R-4.2.1软件的“pheatmap”包绘制低风险组和高风险组样本的临床病理因素和候选基因表达量，同时应用“survival”、“survminer”等包绘制Kaplan-Meier生存曲线来比较两组的OS，最后通过“timeROC”包绘制时间依赖性的ROC曲线，算出训练组样本OS在1年、3年、5年的AUC值，来评价预测风险模型的准确性。将两组的临床病理因素和RiskScore纳入Cox回归分析，加载“rms”等包把具有独立危险因素的变量绘制成列线图，上部为评分系统，下部为预测系统。患者的1年、2年、3年、5年、10年生存率均可通过评分系统预测，并且利用校准曲线和c指数来显示生存预测的准确性。
1.6 验证预测风险模型
用验证组的样本数据来验证预测模型，通过上述的公式计算出验证组各个样本的RiskScore，以RiskScore中位数为cut-off值分为低风险组和高风险组，应用Kaplan-Meier生存曲线比较两组的OS。同样的，绘制1年、3年、5年的ROC曲线，计算AUC值评估模型预测预后的能力。
2 结果
2.1 筛选差异基因和功能富集分析
从UCSC Xena网站上收集了全部的LUAD mRNA表达数据和临床病理数据，提取了IAC中级别病理亚型98例样本和高级别病理亚型82例样本，以P≤0.05和| log2(FoldChange)| >1为阈值，筛选出782个DEGs，包括327个上调基因和455个下调基因，可视化火山图、热图结果所示（图2）。
图2
不同级别病理亚型相关基因的差异分析
a：IAC中级别和高级别病理亚型差异表达基因火山图（327个上调基因和455个下调基因）；b：IAC中级别和高级别病理亚型差异表达基因热图（|log2(FoldChange)| >3的差异表达基因）
为了探索不同级别病理亚型相关的生物学功能和途径，进一步缩小筛选标准，以P≤0.05和| log2(FoldChange)| >1.5得到280个DEGs，对上述280个DEGs进行GO及KEGG分析（图3a～d）。GO分析显示，DEGs主要富集在免疫反应、物质代谢及酶活性调节等相关的生物学功能；KEGG分析主要富集在细胞色素P450调节物质代谢、神经活性配体-受体相互作用通路等相关生物学功能。把P≤0.05的DEGs纳入GSEA分析（图3e～f），其分析结果主要富集于细胞色素P450相关的物质代谢、抗原的呈递与自然杀伤细胞介导的免疫反应等相关生物学过程，进一步验证了中级别与高级别病理亚型DEGs与免疫反应和物质代谢之间的密切关系。以上功能富集分析的结果解释了IAC不同级别病理亚型之间基因层面的关系，也为以后寻找与预后相关生物学标志物的潜在靶点提供了不同的研究思路。
图3
功能富集分析
a和c：GO富集分析 主要富集在免疫反应、物质代谢及酶活性调节等相关的生物学功能；b和d：KEGG通路富集分析 主要富集在细胞色素P450调节物质代谢、神经活性配体-受体相互作用通路等相关生物学功能；e和f：GSEA富集分析 主要富集于细胞色素P450相关的物质代谢、抗原的呈递与自然杀伤细胞介导的免疫反应等相关生物学过程
2.2 Cox分析、LASSO回归及风险模型构建
把782个DEGs纳入单因素Cox分析，得到66个与独立预后显著相关的基因（P＜0.05），随后进行LASSO回归的降维分析（图4），虚线标注处为logλ最小值，筛选出15个基因，经进一步经多因素Cox分析，最终得到最佳构建风险模型的5个基因（MELTF，C=0.1851、MAGEA1，C=0.0475、FGF19，C=0.0269、DKK4，C=0.2300和C14ORF105，C=0.0758）。根据5个基因的基因表达量和风险系数计算每个样本的RiskScore：
图4
LASSO回归分析
通过 LASSO 回归分析进一步确定与IAC的 OS 相关的15个关键基因
RiskScore=（MELTF表达量×0.185 1）+（MAGEA1表达量×0.047 5）+（FGF19表达量×0.026 9）+（DKK4表达量×0.230 0）+（C14ORF105表达量×0.075 8）。
2.3 检验预测风险模型性能
以训练组IAC患者的RiskScore中位数为cut-off值分为低风险组（L=57）和高风险组（H=58）。不同危险组患者表现出不同的相关临床病理特征（图5a）。绘制Kaplan-Meier生存曲线（蓝色代表低风险组，红色代表高风险组），结果显示两组OS差异具有统计学意义（P<0.01），表明高风险组比低风险组更易具有不良预后（图6a）。以训练组RiskScore绘制ROC曲线和时间依赖性ROC曲线结果所示（图6b～c），RiskScore的ROC曲线AUC值为0.675，时间依赖性ROC曲线1年AUC值为0.893、3年AUC值为0.713、5年AUC值为0.632。表明该预测风险模型具有良好的敏感性和特异性。单因素和多因素Cox分析显示，RiskScore是一个独立的预后因素，可作为IAC患者的独立预后预测因子，将低风险组和高风险组的临床病理因素及RiskScore纳入Cox回归分析，识别出具有独立危险因素的变量，绘制得到的列线图结果所示（图7a），该个体化预测模型可估计IAC患者（1年、2年、3年、5年和10年）生存率。校准曲线中的列线图和实际观测结果在训练组和验证组中显示出准确的重叠，表明具有良好的预测精度（图7b）。并且该列线图模型的c指数为0.78，高于其他单一因素的预测模型（图7c）。
图5
风险评分与临床特征的关系
a～b：分别是训练组和验证组每个患者临床特征（Grade、M、N、T、Gender、Age、RiskScore）和5个代表性基因（MELTF、MAGEA1、FGF19、DKK4、C14ORF105），按风险评分升序排列
图6
预测模型的评估
a：训练组K-M生存曲线（P<0.01）；b：训练组ROC曲线，AUC值为0.675；c：训练组时间依赖ROC曲线，1年AUC值为0.893，3年AUC值为0.713，5年AUC值为0.632；d：验证组K-M生存曲线（P<0.01）；e：验证组ROC曲线，AUC值为0.626；f：验证组时间依赖ROC曲线，1年AUC值为0.705，3年AUC值为0.700，5年AUC值为0.684
图7
个体化预测模型
a：Nomogram图可以准确预测患者1年、2年、3年、5年和10年生存率；b：校准图显示了训练组和验证组中1年、2年、3年和5年预测OS与实际OS之间生存概率的比较；c：采用C-Index评价个体化预测模型的预测效果
2.4 验证组数据集验证模型
将验证组中60例样本的DEGs表达量和临床信息合并，经剔除无OS和Censor的患者后，计算剩余56例患者的RiskScore，以中位数为cut-off值分为低风险组（L=28）和高风险组（H=28），验证组的相关临床病理特征可视化结果所示（图5b）。然后绘制Kaplan-Meier生存曲线、ROC曲线及时间依赖性ROC曲线，二者OS仍具有显著统计学差异（P<0.01），时间依赖性ROC曲线1、3、5年AUC值为0.705、0.700、0.684（图6d～f）。上述结果说明该模型在验证组中仍具有较好的预测性能。
3 讨论
尽管免疫治疗和分子靶向治疗已经取得了很大的进展，但由于缺乏准确的早期诊断和预后生物学标志物，肺癌患者的OS仍然较差，5年生存率仅约15%[9]。不同肺癌亚型具有不同的临床特征和预后，确定特异的预后生物学标志物和新的治疗靶点仍然是至关重要的。
在本研究中，我们首先筛选出IAC患者中级别和高级别病理亚型间DEGs，包括327个上调基因和455个下调基因，并对其进行了功能富集分析，经GO、KEGG和GSEA富集分析表明这些基因主要富集在细胞色素P450相关物质代谢、自然杀伤细胞介导的免疫反应、抗原的呈递和酶活性调节等生物学过程。细胞色素P450（CYP）是存在于动物、植物和微生物中的一个参与物质代谢的酶超家族，其密切参与了肺癌、肝癌、乳腺癌等癌症的相关生物学功能[10]。自然杀伤（NK）细胞是免疫系统中不可或缺的肿瘤杀伤效应细胞，通过依赖的细胞毒性和产生细胞因子进行免疫调节，相关证据表明NK细胞可以有效地防止肺癌的发生和发展[11]。NK细胞介导的免疫反应在肿瘤形成的免疫监测中具有至关重要的作用[12]。通过对中级别和高级别两组病理亚型DEGs的分析，解释了两组之间生物侵袭性和预后等因素基因层面的差异，可能为更深层次的研究提供了一定的科研思路。其次，通过对DEGs进行单因素Cox分析，得到66个与独立预后显著相关的基因，再行LASSO回归分析和多因素Cox分析后，我们最终确定了由5个特征基因（MELTF、MAGEA1、FGF19、DKK4、C14ORF105）组成的风险评估模型，并根据中位RiskScore将样本分为低风险组和高风险组。随后，我们对RiskScore与临床病理特征的相关性进行了分析，结果显示随着复发评分的增加临床病理呈现不对称分布特征。进一步对两组患者的预后进行研究，结果显示，高风险组患者的OS明显短于低风险组（P<0.01）。使用ROC曲线和时间依赖性ROC曲线评价预后预测性能，结果显示，ROC曲线AUC为0.675，时间依赖性ROC曲线1年、3年和5年生存率的AUC值分别为0.893、0.713、0.632，表明模型具有良好的准确性，在验证集（P<0.01，1年AUC=0.705、3年AUC=0.700、5年AUC=0.684）中也得到验证。最后，用RiskScore和两组具有独立预后的临床病理因素绘制出列线图，该模型能精确预测患者1年、2年、3年、5年和10年的生存率，校准曲线显示出准确的重叠，c指数也高于其他单一因素（RiskScore、Grade）的预测模型。
上述5个特征基因已被报道与多种癌症的发生发展密切相关。MELTF又称MTf（Melanotransferrin）或MTF1（Metal regulatory transcription factor 1），是铁结合转铁蛋白的同源物，可作为LUAD和胃癌的预后标志物[13-14]。研究[14-15]表明MTf能促进癌细胞侵袭、迁移、增殖和EMT进展，是一种具有很大潜力的靶点。MAGEA1是属于X染色体簇的癌症/睾丸抗原，在多种肿瘤中表达，如肺癌，它可能是免疫治疗的理想抗原靶点[16]。研究[17]表明，随着肿瘤的进展，MAGEA1的表达水平增加，说明了MAGEA1的检测有可能作为肺癌早期诊断的辅助检查方法。成纤维细胞生长因子（FGFs）家族由22个成员组成，分为7个亚家族，调节着一系列生物学功能和多种发育过程[18]。作为FGF家族内分泌亚群的一员，FGF19在维持蛋白质和葡萄糖代谢方面都起着关键作用[19]。FGF19具有与4种类型的酪氨酸激酶FGF受体（FGFR 1-4）结合的潜力，但亲和力相当低，需要β-klotho （KLB）作为辅受体，一旦FGF19结合到FGFRs和KLB，它就会启动二聚作用和自身磷酸化，最终激活细胞质信号转导通路[19-21]。相关研究[19]表明FGF19在NSCLC中表达明显上调，且FGF19过表达与不良预后显著相关，可能作为预测NSCLC患者不良预后的潜在生物学标志物。DKK4作为DKK家族中的一员，是Wnt/β-catenin信号通路的抑制剂[22]。研究[23]发现DKK4在肾细胞癌、结直肠癌、胃癌、NSCLC和上皮性卵巢癌中表达上调，在肝细胞癌中表达下调，DKK4也参与肿瘤生长、侵袭、迁移和化疗耐药等生物功能。目前已经有很多关于DKK4的研究，然而，DKK4在癌症或非癌症疾病中的作用机制和分子机制仍不明确，许多前沿领域有待探索，所以本文为进行更深层次的研究提供了可能的方向。研究发现C14ORF105在耐紫杉醇的卵巢癌细胞中过表达，提示其可能参与了紫杉醇耐药[24]。到目前为止，C14ORF105在肺癌中的作用尚未被广泛研究。综上所述，MELTF、MAGEA1、FGF19、DKK4和C14ORF105均参与了肿瘤的发生或进展。
总的来说，我们的研究通过Cox分析和LASSO回归分析构建了一个准确有效的5基因组预测模型。基于该5基因组的RiskScore可用于确定LUAD患者的OS。将RiskScore和临床参数相结合的列线图可用于预测LUAD患者1年、2年、3年、5年和10年的生存率，所以，它将有助于LUAD患者的预后和随访监测，为LUAD患者的个体化诊疗提供参考。但需要指出的是，我们的研究还存在一定的局限性。首先，我们的研究数据主要来自TCGA数据集，有必要在大型独立临床队列中评估其预测效能。其次，我们纳入的分组仅包括了中级别和高级别病理亚型，排除了样本量较少的低级别病理亚型，并且需要进一步外部数据集的验证。最后，这5个基因在LUAD发病中的生物学机制有待进一步通过功能研究来阐明。
利益冲突：无。
作者贡献：崔严奇负责病例筛选、数据整理与论文设计、初稿撰写等；曾志勇、廖毅、倪琳、连铎煌、杨鲸蓉、叶仕新和张锦灿负责论文审阅与修改。
