1. 引言

肺癌是全球范围内发病率和死亡率都较高的一种癌症[1]，其中以非小细胞肺癌(non-small cell lung cancer, NSCLC)最为常见，NSCLC按照组织学分类可以分为腺癌、鳞状细胞癌、腺鳞癌、大细胞癌、肉瘤样癌等病理亚型，其中又以肺腺癌(lung adenocarcinoma, LUAD)最为常见。

肺癌的传统治疗方法包括手术和化疗。近年来，随着医学研究的进展，免疫治疗逐渐成为一种有效的癌症治疗干预手段[2] [3]。然而，只有部分患者的预后通过免疫治疗后能得到极大改善，我们亟需预测癌症患者的预后及免疫治疗效果的手段。另一方面，肿瘤微环境(tumor microenvironment, TME)被证实在多种癌症的发生发展中起着关键作用，相关细胞如巨噬细胞、内皮细胞和成纤维细胞等，已被纳入各种研究[4]。已有的研究表明，TME确实对肺癌的发生、发展和治疗反应有重要影响[5]，相关基因可能成为未来预测肺癌预后和治疗反应的理想标志物。

在本研究中，我们从TCGA和GEO数据库中下载了相关数据，系统性分析了LUAD患者的TME相关基因特征。通过NMF方法将患者分为不同亚群，并建立了相应的预后模型。

2. 资料与方法

2.1. 资料

通过从TCGA数据库(https://portal.gdc.cancer.gov/)获得了539例肺腺癌组织样本和59例正常肺组织样本的基因表达数据。剔除临床信息不完整的样本，获得507个数据完整的样本进行研究。从GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)下载验证队列GSE13213作为验证数据集。本研究中使用的TME相关基因来自于已发表的高质量相关研究[6]-[12]。

2.2. 方法

2.2.1. DEGs的筛选及NMF分型

采用R软件中的“limma”R包，筛选出LUAD组织中与TME相关的差异表达基因(differentially expressed genes, DEGs) [13]，将筛选出的基因及其表达水平融合为基因表达矩阵。采用非负矩阵分解(nonnegative matrix factorization, NMF)算法对所有LUAD样本进行50次迭代，以提取生物相关系数来预测基因表达矩阵的内部特征结构[14]。

2.2.2. 预后模型的构建

将507个LUAD样本按照最优比例7:3分为训练集和测试集。使用“survival”R包对所有表达差异的TME相关基因进行单因素Cox回归的赋值分析，以得到显著的预后基因来构建模型。根据Lasso回归分析法，使用“glmnet”和“caret”R包评估LUAD患者的相关预后风险特征。基于筛选出的基因，成功构建了一个包含5个TME相关基因的LUAD患者预后模型。根据差异基因表达量值和回归分析系数计算风险评分，模型的风险计算公式如下：

$riskscore={\displaystyle {\sum }_{i=0}^{n}coefi\ast Xi}$

其中，coef为风险系数，X为基因表达水平。

所有样本按中位风险评分分为高危组和低危组。采用受试者工作特征曲线(receiver operating characteristic, ROC)评价风险评分的预后价值。绘制相应的K-M (Kaplan-Meier)曲线，并使用Log-rank方法评估总生存率(OS)。

2.2.3. 预后模型的验证与评价

使用受试者工作特征曲线(ROC)、校准曲线(calibration curve, CC)和决策曲线分析(decision curve analysis, DCA)来评估模型的性能。将模型与其他已建立的同类型模型进行比较。通过受限平均生存(restricted mean survival, RMS)时间评估模型的预后准确性。此外，还计算并比较了所有模型的一致性指数(C指数)。使用R中的“limma”包来研究风险评分与各种临床因素之间的相关性，使用“ggpubr”包进行可视化。

2.2.4. GSEA和免疫相关性分析

使用GSEA软件进行基因集富集分析(gene set enrichment analysis, GSEA)，错误发现率(false discovery rate, FDR) < 0.05的基因集被认为有统计学意义。并进一步研究风险评分与免疫细胞浸润水平的相关性，风险评分与免疫检查点、DNA复制、错配修复、上皮–间质转化(epithelial-mesenchymal transition, EMT)等特征性基因表达水平的相关性。

2.3. 统计学处理

本研究采用R软件(4.1.3)进行相关分析和绘图。使用了“limma”包、“ggplot2”包、“ggpubr”包、“NMF”包、“survival”包、“glmnet”包、“caret”包、“Survminer”包、“timeRoc”包、“RMS”包、“MCPcounter”包和“corrplot”等软件包。采用Wicoxon秩和检验进行两组之间的比较，两组及以上的比较采用Kruskal-Wallis秩和检验。采用Kapaln-Meier法进行生存分析。采用Pearson相关性分析法，以P < 0.05为有统计学意义。

3. 结果

3.1. 通过NMF法构建LUAD亚型

从TCGA数据库下载相关数据，计算肿瘤和正常组织之间的差异表达基因(DEGs)，并筛选出其中与TME相关的基因，共鉴定出2854个基因，其中903个在LUAD组织中扩增，见图1(A)。然后，采用NMF聚类分析法，将所有样本划分为两个不同的亚型(C1和C2)，发现当k = 2时，聚类结果最佳，一致性矩阵热图的边界最为清晰，见图1(B)。C2组患者比C1组患者具有更好的总生存期(OS)和无进展生存期(PFS)，见图1(C)、图1(D)。

Figure 1. (A) Volcano map of TME-related differentially expressed genes; (B) The results of NMF clustering; ((C) (D)) Survival analysis of patients in Subtype 1 (C1) and Subtype 2 (C2). LUAD: Lung adenocarcinoma; TME: Tumor microenvironment; NMF: Nonnegative matrix factorization

图1. (A) TME相关差异表达基因的火山图；(B) NMF法聚类结果；((C) (D)) 亚型1 (C1)和亚型2 (C2)中患者的生存分析。LUAD：肺腺癌；TME：肿瘤微环境；NMF：非负矩阵分解

已知免疫细胞在TME中起着重要作用，我们评估了两种亚型中特异性免疫细胞的丰度，其中中性粒细胞、NK细胞、T细胞、细胞毒性淋巴细胞、内皮细胞、成纤维细胞和髓系树突状细胞在两个亚型间的免疫评分差异有统计学差异(P < 0.05)。见图2(A)~(G)。之前的研究已经确定了6种与TME相关的肿瘤免疫亚型[15]，与我们得到的LUAD亚型进行比较，见图2(H)。

Figure 2. (A)~(G) Immune scores of cells of the tumor microenvironment (TME) in different subtypes; (H) Plot of the relationship between established molecular subtypes and C1 and C2

图2. (A)~(G) 不同亚型中肿瘤微环境(TME)细胞免疫评分；(H) 已建立的分子亚型分类与C1和C2的关系图

3.2. 肺癌预后的基因特征

在7:3的最优比例下，将507个样本随机分为训练集(n = 356)和测试集(n = 151)。采用单因素Cox分析(P < 0.05)，在保持高精度的同时，减少研究基因的数量，为了使过拟合的风险最小化，采用Lasso回归算法，见图3(A)、图3(B)，最后通过Cox比例风险模型对差异显著的TME相关基因进行分析，构建出以下风险评分公式：

$\begin{array}{c}Riskscore=\left(\text{C1QTNF}6*0.365178453387755\right)+\left(\text{PLEK2}*0.284495496551624\right)\\ +\left(\text{FURIN}*0.154361274455861\right)-\left(\text{TM6SF1}*0.396920071856817\right)\\ +\left(\text{IGF2BP}1*0.208860117937174\right)\end{array}$

通过得到的风险评分公式发现，高危基因(C1QTNF6、PLEK2、FURIN和IGF2BP1)提示患者预后较差。相比之下，TM6SF1基因的高表达则与LUAD患者更长的生存期相关。

我们分别研究了模型在训练集和测试集中的预测效果。在训练集中，高危组的LUAD患者的生存期更短，见图3(C)。为了进一步评价模型的有效性和准确性，我们还绘制了时间相关的ROC曲线图，见图3(E)，1年、3年和5年的曲线下面积(AUC值)分别为0.698、0.713和0.710，AUC值可以在一定程度上反映模型预测的准确性，且与模型准确性呈正相关。与训练集一致，测试集中高危组和低危组的患者也表现出明显不同的预后，见图3(D)，且1年、3年和5年的AUC值均高于0.65，较为理想，见图3(F)。

最后，为了进一步验证模型的有效性，我们从GEO数据库下载了GSE13213 (n = 117)进行进一步验证。在GSE13213数据集中，高危患者的生存时间仍低于低危患者，见图3(G)，1年、3年和5年的AUC值分别为0.845、0.602和0.622，见图3(H)。这些结果与训练集和测试集获得的结果高度一致，进一步证明了我们的预测模型的有效性和稳定性。

Figure 3. (A) Lasso coefficient path plot; (B) Cross-validation plot; (C) Survival analysis of training set; (D) Survival analysis of test set; (E) ROC curve plot of training set; (F) ROC curve plot of test set; (G) Survival analysis of GSE13213 dataset; (H) ROC curve plot of GSE13213 dataset

图3. (A) Lasso系数路径图；(B) 交叉验证图；(C) 训练集的生存分析；(D) 测试集的生存分析；(E) 训练集的ROC曲线图；(F) 测试集的ROC曲线图；(G) GSE13213数据集的生存分析；(H) GSE13213数据集的ROC曲线图

3.3. 列线图的建立

Figure 4. (A) Nomogram for predicting overall survival; (B) Calibration plot of survival rate (relationship between risk score and clinicopathological factors); (C) Multi-indicator ROC curve of risk score and other parameters; (D) Relationship between risk score and tumor staging; (E) Relationship between risk score and T staging (primary tumor extent); (F) Relationship between risk score and N staging (regional lymph node metastasis status)

图4. (A) 用于预测总生存期的列线图；(B) 生存率校准图(风险评分与临床病理因素的关系)；(C) 风险评分与其他指标的多指标ROC曲线；(D) 风险评分与分期的关系；(E) 风险评分与T分期(原发肿瘤范围)的关系；(F) 风险评分与N分期(区域淋巴结转移状态)的关系

为整合各种临床危险因素，我们绘制了列线图来量化患者的风险，结合性别、年龄、分期(TNM分期)、T (原发肿瘤大小及范围)、N (淋巴结受累情况)和风险评分来预测LUAD患者的1年、3年和5年的总生存(OS)率，见图4(A)，校准曲线显示，预测的生存率与实际的生存率密切相关，见图4(B)。

我们使用多指标ROC曲线进一步评估模型预测的准确性，列线图和风险评分相较于其他单一评估指标显示出更好的预测准确性，见图4(C)。此外，我们还研究了风险评分与性别及年龄之间的关系。发现风险评分会随着LUAD的进展而增加，见图4(D)~(F)，且不受性别(女性和男性)和年龄(≤65和>65)的影响，见图5(A)~(D)。以上结果证明预后模型的准确性和稳定性。

Figure 5. ((A) (B)) Relationship between risk score and gender (female and male); ((C) (D)) Relationship between risk score and age (≤65 years and >65 years)

图5. ((A) (B)) 风险评分与性别的关系(女性和男性)；((C) (D)) 风险评分与年龄的关系(≤65和>65)

3.4. 模型的比较分析

将我们的预后模型与其他相同类型的模型进行比较，以进一步验证其预测性能，本研究中共涉及其他三个研究团队建立的预后模型(Feng等[16]、Chen等[17]和Wu等[18])。所有模型中，低危组的预后均明显优于高危组(P < 0.05)，见图6(A)~(D)。相较于其他模型，我们的模型在第五年的AUC值更大，说明我们的预后模型具有更好的预测性能，见图6(E)~(H)。进一步计算了一致性指数(C指数)，我们的模型的值最高，为0.676，见图6(I)。最后，通过受限平均生存(RMS)时间评估模型在不同时间点的预测效果，在>80个月时，我们的模型表现较好，见图6(J)。以上结果均能够提示我们的模型预测的准确性和稳定性。

Figure 6. (A)~(D) Kaplan-Meier curves of the four prognostic models; (E)~(H) ROC curves of the four prognostic models; (I) C-indices of the four prognostic models; (J) Restricted mean survival (RMS) time curves of the four prognostic models

图6. (A)~(D) 四个预后模型的K-M曲线；(E)~(H) 四个预后模型的ROC曲线；(I) 四个预后模型的C指数；(J) 四个预后模型的受限平均生存时间曲线

3.5. GSEA和免疫相关性分析

最后，我们对来自两个风险组的样本进行了基因集富集分析(GSEA)。结果显示，高危组的LUAD患者的基因主要在以下生化通路和活性方面，见图7(A)、图7(B)：

KEGG_CELL_CYCLE、KEGG_DNA_REPLICATION、KEGG_ECM_RECEPTOR_INTERACTION、KEGG_FOCAL_ADHESION、KEGG_PROTEASOME、GOBP_CHROMOSOME_SEGREGATION、GOBP_CORNIFICATION、GOBP_DNA_DEPENDENT_DNA_REPLICATION、GOBP_DNA_REPLICATION、GOBP_EPIDERMAL_CELL_DIFFERENTIATION。

由于肿瘤免疫微环境是影响免疫治疗应答的重要因素之一[19] [20]。因此，我们采用微环境细胞种群计数器算法(MCP-counter algorithm)对计算患者的免疫浸润细胞评分，以研究其与风险评分的相关性分析。结果显示，风险评分与成纤维细胞呈正相关，与T细胞、单核细胞谱系、髓系树突状细胞、中性粒细胞和内皮细胞呈负相关，见图7(D)；与DNA复制相关基因POLE2、FEN1和MCM6，上皮–间质转化(EMT)相关基因LOXL2的表达水平呈显著正相关，见图7(C)。

Figure 7. (A) KEGG analysis of the high-risk group; (B) GO analysis of the high-risk group; (C) Relationship between risk score and expression of representative genes in oncology; (D) Relationship between the risk score and immune infiltrating cells

图7. (A) 高危组KEGG分析；(B) 高危组GO分析；(C) 风险评分与肿瘤学中具有代表性的基因的表达之间的相关性；(D) 风险评分与免疫浸润细胞之间的相关性

4. 讨论

组学技术的快速发展极大地提高了我们对TME的认识，我们逐渐开始发现它在肿瘤免疫治疗中的重要作用[21] [22]。进一步加深相关理解将有助于提高免疫治疗的有效性和稳定性，因此，近年来关于TME的相关研究日益增多。通常肿瘤组织内部和周围会有大量免疫细胞浸润，免疫细胞和肿瘤细胞之间的关系值得进一步研究探讨。分析肿瘤组织免疫细胞的组成和比例已成为TME相关研究的重要组成部分之一[23]，既往的研究已经证明，通过生物信息学方法评估肿瘤的纯度具有重要意义，值得进一步研究[24]-[26]。

近年来，免疫治疗已成为一种重要的肿瘤治疗方法，且相关研究表明其可以极大地改善部分肺癌特别是非小细胞肺癌(NSCLC)患者的预后[27] [28]。但免疫治疗也有较大的不足，其中最突出的是其高昂的使用成本及明显的临床副作用，因此确定哪些患者可以从免疫治疗中获得更大的益处尤为重要，识别肿瘤中的免疫细胞并进一步解析TME在肿瘤发生发展中的作用在临床医学和基础研究中越来越重要。

本研究旨在建立一个基于TME相关基因的LUAD预后模型，以研究TME的异质性在预测肺癌患者的预后和治疗反应中的重要作用[29]-[31]。非负矩阵分解(NMF)方法，可以提取基因表达矩阵中数据的生物相关系数，捕获内部结构特征，并对样本进行分组，目前在疾病分类中得到了广泛应用。本研究通过NMF法将研究样本分为两组，结果显示，不同亚型间TME免疫浸润细胞的评分存在显著差异。此外，预后指标总生存期和无进展生存期的不同也证实了两组间的差异，以上结果清楚地提示了LUAD中TME的异质性。

通过验证，证实了我们的预后模型能够较为准确稳定地预测LUAD患者的预后，并得到了5个TME相关基因以构建风险评分公式。其中，C1QTNF6基因在许多类型的癌症中都存在过表达现象，研究表明它可以独立作为许多肿瘤预后不良的指标；此外，已经证实C1QTNF6的表达水平与肿瘤标志物通路评分、肿瘤微环境相关通路评分和肿瘤细胞免疫特征呈正相关；在药敏分析方面，C1QTNF6的高表达水平也与较高的耐药性相关[32]。相关研究表明PLEK2基因的过表达可能是LUAD患者无进展生存期(PFS)预后不良的特异性生物标志物，其表达水平与肿瘤细胞侵袭、细胞周期、DNA损伤和DNA修复呈正相关[33]。另外，GF2BP1基因的高表达提示鼻咽癌患者的预后不良[34]。而在Zhong等[35]建立的肺腺癌预后模型中，TM6SF1的高表达提示预后良好。

通过进一步GSEA富集分析发现，高危组LUAD患者的富集分析结果主要集中在免疫相关的生物过程。免疫相关性分析结果也表明，风险值与成纤维细胞的免疫评分呈正相关，与T细胞、单核细胞谱系、髓系树突状细胞、中性粒细胞和内皮细胞的免疫评分呈负相关。在以前的研究中也获得了类似的结果，可能表明高风险组LUAD患者对免疫治疗的反应较差[36]。

目前免疫检查点相关基因受到了越来越多研究者的关注[37] [38]。本研究中，风险评分与DNA复制相关基因(POLE2、FEN1、MCM6)和上皮–间质转化(EMT)相关基因(LOXL2)的表达水平呈显著正相关。DNA复制相关基因的表达水平对DNA复制和细胞周期调控有很大的影响[39]，其错误表达将极大地促进细胞增殖和肿瘤发生[40] [41]。之前的研究也表明，POLE2和MCM6在一些肿瘤组织中过表达，而其高表达与不良预后相关[42] [43]。EMT基因阳性可能导致肿瘤患者有更高的转移倾向和较差的预后[44]。我们得到的研究结果提示，一些特殊基因参与了肺癌的发生发展，在未来可能为其诊断和治疗提供强有力的依据。

我们的研究仍有一些不足之处。首先，我们的原始研究数据都来自TCGA和GEO数据库，样本量不够充足，可能会导致结果的偏倚。其次，我们的结论缺乏实验验证，因此，未来有必要开展多中心、大样本、前瞻性双盲试验进行进一步验证。

5. 结论

本研究成功构建了TME相关的预后模型，可以较为准确稳定地预测肺癌患者的预后。可为未来进一步研究肺癌患者的个体化治疗提供方向和依据。

NOTES

^*通讯作者。

参考文献