基于自噬相关基因建立和验证口腔鳞状细胞癌的预后风险评分模型

齐延旭， 马 川， 朱 勇， 黄晟栋， 赵华强， 张风河

（山东大学齐鲁医学院口腔医学院，口腔医院口腔颌面外科，山东省口腔组织再生重点实验室，山东省口腔生物材料与组织再生工程实验室，山东 济南 250012）

自噬是一个高度保守的真核细胞循环过程，在细胞生存和维持中发挥着重要作用[1-2]。已有的研究表明，自噬在肿瘤发生、发展中是一把双刃剑，在肿瘤的形成和进展中起着相反的作用[3-6]。在口腔鳞状细胞癌（OSCC）中，自噬通过Beclin-1/Atg7/Atg12-Atg5通路抑制OSCC进展[7]；另一方面，也有研究显示，自噬促进了OSCC的干性和耐药性，线粒体自噬可能是导致OSCC化疗耐药性的原因[8]。目前已开发出许多促进OSCC自噬的药物，并取得了良好的治疗效果。然而，这些研究主要集中在自噬对OSCC进展和治疗的影响，很少有研究者研究自噬在OSCC预后中的作用[9]。

本研究通过分析OSCC中差异表达的自噬相关基因（ARGs），建立OSCC患者的风险评分模型，验证其独立预后价值，并构建列线图，以期对OSCC患者预后评估和临床治疗提供新思路。

1 材料和方法

1.1 数据获取

我们从人类自噬数据库 （human autophagy database，HADb）（http://www.autophagy.lu/）中 获 得232个自噬相关基因，然后从癌症基因组图谱（the cancer genome atlas，TCGA）数据库（https://portal.gdc.cancer.gov/）下载372个（340例OSCC患者和32例正常人）FPKM标准化的转录组测序（RNA-seq）数据和OSCC队列的临床病理和生存信息。从基因表达汇编（GEO）数据库（http://www.ncbi.nlm.nih.gov/geo/）下载GSE41613数据集（196例OSCC）作为验证集。

1.2 差异自噬相关基因筛选及功能富集

引用“limma”包对TCGA下载数据进行数据处理[10]，利用WilcoxTest统计分析OSCC样本与正常口腔组织样本中差异表达的自噬相关基因。之后利用错误发现率（false discovery rate，FDR）方法对每个基因通过WilcoxTest统计分析得到的P值进行校正，得到校正后的P值（adj.P）；并计算自噬基因在肿瘤组织中表达量的均值，将该均值除以正常组织中表达量的均值得到foldchange。设置筛选标准：adj.P＜0.05，|log2（foldchange）|＞1，即筛选出肿瘤中表达量是正常组织中表达量2倍以上的自噬基因。

使用“org.HS.eg.db”包对自噬相关基因进行Symbol-ID转换,使用“DOSE”“clusterProfiler”和“enrichplot”包进行基因本体论（gene ontology,GO）和京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）富集分析[11]，挑选出每个基因的ID所对应的GO功能和KEGG信号通路，将这些差异基因与上述挑选出的GO功能和KEGG信号通路在背景中所对应的基因进行比较，得到P值，检验这些差异自噬基因在这些GO功能和KEGG信号通路上是否有影响。最后使用“barplot”包进行可视化，设置以P＜0.05作为纳入标准，定义在差异表达基因中显著富集的GO功能和KEGG信号通路。

1.3 构建风险评分模型

在TCGA训练集中，利用log2（x+1）转换将OSCC队列中RNA-Seq数据进行标准化，将标准化后差异表达的ARGs表达量与生存信息合并，使用“survival”包进行单因素Cox比例风险回归分析，以确定与OS显著相关的ARGs。然后使用“glmnet”包进行套索（LASSO）回归分析[12]，LASSO根据每个样本的表达量，约束参数防止模型过拟合，选择最佳“λ”值，从而筛选重要预后基因，为每个预后基因分配相应的回归系数，构建风险评分模型。风险评分计算公式如下：

风险评分（risk score）=Gene1表达水平×β1+Gene2表达水平×β2+……+Gene n的表达水平×βn（β为回归系数，n为预后相关ARGs的总数目）。

1.4 风险评分模型评估

将所有TCGA训练集中的OSCC患者按照风险评分中位值分为高风险组和低风险组。利用“survival”包绘制高、低风险组的风险评分分布曲线、生存情况分布图及预后相关基因表达量热图。绘制Kaplan-Meier（K-M）生存曲线评价2组患者的生存差异，采用对数秩（log-rank）检验评估其显著性。为了评估预后模型的预测准确性，用“survival ROC”包构建了1、3、5年生存率的时间依赖受试者工作特征（receiver operating characteristic,ROC）曲线，并计算曲线下面积（area under the curve,AUC）。对GEO验证集（GSE41613）采用同样的方法检验风险评分模型的普遍性和可靠性。

1.5 独立预后因素鉴定

为确定风险评分是否为OSCC的独立预后因素，利用“survival”包进行单因素及多因素Cox回归分析，比较TCGA训练集中风险评分与其他临床病理因素（年龄、性别、病理分级、肿瘤分期、T分期、N分期）的预后相关性，若单因素及多因素Cox回归分析均显示与预后显著相关（P＜0.05），则表示该因素为OSCC独立预后因素。

1.6 列线图的建立及验证

为了预测OSCC患者1、3、5年生存的概率，我们利用“rms”包，结合年龄、性别、病理分级、肿瘤分期、T分期、N分期及风险评分，在训练集中构建列线图，然后绘制校准曲线，以图形方式评估实际生存率与预测生存率之间的一致性[13]。

1.7 肿瘤免疫浸润细胞特征分析

利用CIBERSORT算法绘制训练集中OSCC患者免疫细胞的比例[14]；随后通过“corrplot”包进行Pearson相关性分析，确定免疫细胞之间的共表达模式；最后通过非参数WilcoxTest进行检验比较高、低风险组免疫细胞浸润的差异，探究高、低风险组患者预后差异的潜在免疫机制。

1.8 统计学分析

所有统计分析均采用R4.0.5软件（https://www.r-project.org/）进行。采用Wilcoxon秩和检验评估风险评分与临床病理特征的相关性。在所有的统计检验中，以P＜0.05为差异有统计学意义。

2 结果

2.1 差异表达自噬相关基因的鉴定

我们首先从HADb数据库中获得232个ARGs，然后从TCGA数据库中下载OSCC队列中372例患者的RNA-seq数据和临床及预后数据。以FDR＜0.05，|log2（foldchange）|＞1为筛选标准，得到37个差异表达的ARGs，其中11个为下调ARGs,26个为上调ARGs（图1A、B）。肿瘤组织与正常组织间差异表达的ARGs的表达情况见图1C。

2.2 差异表达ARGs的功能注释

对37个差异表达的ARGs进行功能富集分析。GO富集分析表明，在生物学过程中，ARGs主要富集于神经元死亡、凋亡信号通路的调节，细胞对未折叠蛋白的反应等；在细胞成分中，ARGs主要富集于整合素复合物、参与细胞黏附的蛋白质复合物和自噬体膜；在分子功能方面，ARGs主要富集在细胞因子活性、受体配体活性、信号受体激活剂活性等方面（图1D）。在KEGG通路中，ARGs主要富集于细胞凋亡、铂类药物耐药、EGFR酪氨酸激酶抑制剂耐药、ErbB信号通路及癌症中的PD-L1表达和PD-1检查点通路等（图1E）。

图1 OSCC和正常对照样本中ARGs的表达差异和基因功能富集分析Figure 1 The expression difference and gene functional enrichment analysis of ARGs in OSCC and normal samples

2.3 基于预后相关的ARGs建立风险评分模型

经单因素Cox回归分析，在TCGA数据集中共筛 选 出6个ARGs（BID、NKX2-3、DDIT3、VEGFA、FADD、BIRC5）作为预后相关基因（P＜0.05），森林图显示了每个与预后相关的ARG的危险比（harzard ratio,HR）和相应的置信区间（图2A），这些预后基因被纳入进一步的LASSO回归分析（图2B）以筛选最佳建模基因，图2C为6个基因的系数分布图，以此构建风险评分模型。在这些ARGs中，BID、DDIT3、VEGFA、FADD、BIRC5的HR＞1，提示高表达与患者预后不良有关，而NKX2-3的HR＜1，其低表达预示着患者的不良预后。风险评分计算如下：

图2 差异表达的ARGs的单因素Cox回归分析及LASSO回归分析Figure 2 Univariate Cox regression and LASSO regression analysis of ARGs with differential expression

风 险 评 分=（0.197×BID表 达 量）+（-0.409×NKX2-3表达量）+（0.091×DDIT3表达量）+（0.226×VEGFA表达量）+（0.052×FADD表达量）+（0.072 7×BIRC5表达量）。

2.4 风险评分模型评估

我们在TCGA数据集内部和GEO验证集中验证了ARGs风险评分模型的性能。将TCGA数据集中的OSCC患者按照风险评分中位数分为高风险组和低风险组。低风险组NKX2-3的表达水平高于高风险组，另一方面，高风险组中其他自噬相关预后基因表达水平较高，将2组的风险评分与生存数据进行可视化（图3A）。时间依赖性ROC曲线显示，ARGs预后模型对OS的预测效果较好，1、3、5年的AUC值分别为0.61、0.66、0.67（图3B）。生存分析结果显示，高风险组的总体生存率明显低于低风险组（图3C）。在GEO验证集中用同样的方法进行验证，得到的结果与TCGA测试集相似。在GEO数据集中，风险评分分布、生存状态和基因表达模式如图4A所示，除保护性基因NKX2-3外，其余5个ARGs在高风险组的表达水平均高于低风险组。GEO数据集中风险评分1、3、5年的AUC值分别为0.65、0.67、0.68（图4B）。高风险组的总体生存率显著低于低风险组（图4C）。

图3 基于ARGs的风险评分模型在TCGA数据集显示出良好的预测能力Figure 3 ARGs-based risk score model displays favorable prediction ability in the TCGA datasets

图4 基于ARGs的风险评分模型在GEO验证集显示出良好的预测能力Figure 4 ARGs-based risk score model displays favorable prediction ability in the GEO verification datasets

2.5 鉴定独立预后因素

风险评分在单因素（HR=2.522,95%CI=1.496～4.249,P＜0.001）及多因素（HR=3.131,95%CI=1.933～5.073,P＜0.001）Cox回归分析中均显示，其与OSCC患者的预后显著相关。另外，在其他临床病理因素中，年龄在单因素（HR=1.028,95%CI=1.010～1.046,P＜0.01）及多因素（HR=1.023,95%CI=1.006～1.039,P＜0.01）Cox回归分析中与患者OS显著相关（图5A、B）。以上结果说明风险评分和年龄这两个因素具有独立预后价值，可作为OSCC的独立预后因素。ROC曲线验证了ARGs风险评分模型的预测准确性，风险评分AUC值为0.627，高于其他临床特征（图5C）。

图5 独立预后因素评估Figure 5 Evaluation of independent prognostic factors

2.6 风险评分与临床病理变量之间的关系

为了确定自噬相关风险评分模型是否与OSCC进展有关，我们分析了TCGA数据集中风险评分与临床病理变量之间的相关性。≤65岁与＞65岁2组间（P＞0.05，图6A），男性与女性性别间（P＞0.05，图6B）以及病理分级G1～2与病理分级G3～4 2组间（P＞0.05，图6C）风险评分差异无统计学意义。Ⅲ、Ⅳ期的风险评分高于Ⅰ、Ⅱ期，且差异具有统计学意义（P＜0.001，图6D），T3～4期的风险评分显著高于T1～2期（P＜0.01，图6E），N1～3期的风险评分显著高于N0期（P＜0.05，图6F）。上述结果表明，风险评分纳入的ARGs与OSCC患者肿瘤分期、T分期及N分期显著相关。

图6 OSCC患者组间风险评分分布Figure 6 Risk score distribution among OSCC patients

2.7 构建并验证列线图

为了定量估计临床环境下OSCC患者的生存概率，我们在TCGA训练集中构建了一个列线图，整合ARGs风险评分模型和临床病理特征（年龄、性别、病理分级、肿瘤分期、T分期、N分期）来预测患者1、3、5年的生存概率。通过图7A，我们能够定量得出列线图的预测效果，样本分数越高，预后越差。此外，我们构建了校准曲线，以图形化评估列线图预测与实际生存率之间的一致性。在数据集的1年（图7B）、3年（图7C）、5年（图7D）的生存概率中，实际生存和预测生存是吻合的。

图7 构建列线图预测OSCC患者1、3、5年生存率Figure 7 Nomogram is constructed to predict 1-,3-,and 5-year survival rate of OSCC patients

2.8 鉴定肿瘤免疫浸润细胞特征

如图8A所示，利用CIBERSORT算法检测到OSCC患者中22种免疫细胞的比例。通过Wilcoxon秩和检验比较高、低风险组免疫细胞浸润的差异，发现高风险组M0期巨噬细胞显著增加（P＜0.001）,CD8+T细胞、M1期巨噬细胞显著低于低风险组（P＜0.001）。其中CD8+T细胞在高、低风险组中浸润差异如图8B所示,高风险组CD8+T细胞浸润比例中位值为4.26%，低风险组中位值为7.94%（P＜0.001）。

图8 免疫细胞浸润模式Figure 8 Immune cell infiltration pattern

3 讨论

生物标志物一直被认为是癌症诊断和靶向治疗的重要工具。在过去的几十年里，大量的生物标志物被用于OSCC的早期检测或治疗过程中[15]。目前已经应用的OSCC患者靶向治疗方案包括表皮生长 因 子受 体 （epidermal growth factor receptor，EGFR）抑制剂和西妥昔单抗（cetuximab,Cet），但对EGFR抑制剂治疗的耐药性是临床使用该药物的一大阻碍。此外，低成本和高产量的标准生物标志物诊断模型仍然缺乏[16]。越来越多的研究表明，自噬在OSCC的肿瘤发生、进展和耐药中起着重要作用[17-21]。一些ARGs被认为是OSCC的潜在生物标志物，如Liang等[22]指出，高表达的细胞质自噬蛋白p62可能是OSCC患者预后不良的标志物。Hou等[23]构建了一个基于13个ARGs的OSCC预后模型，包括4个风险基因和9个保护性基因，在不同数据集的验证中均展现了良好的预测性能。本研究在保证预测准确性的基础上，纳入了更少的ARGs构建预后模型，可能在未来的前瞻性研究和临床实际应用中有更好的可操作性；此外，本研究还引入了免疫细胞浸润分析，深层次地挖掘该模型影响预后的潜在机制。

本研究分析了OSCC和正常口腔组织中差异表达的ARGs，GO和KEGG富集分析表明，差异表达的ARGs主要富集在凋亡中，提示这些ARGs与凋亡显著相关，与之前的研究一致，自噬介导的细胞凋亡通过AKT/mTOR通路促进OSCC进展[24-25]。这些ARGs主要参与铂耐药途径，其可能在OSCC治疗策略中发挥重要作用[26]。另外，这些ARGs富集于ErbB信号通路和EGFR酪氨酸激酶抑制剂耐药，这些都与OSCC进展相关[27]，提示这些ARGs可能为OSCC提供潜在的治疗靶点。通过单因素Cox回归分析及LASSO回归分析，将筛选出的6个ARGs（BID、NKX2-3、DDIT3、VEGFA、FADD、BIRC5）纳 入风险评分模型，并在TCGA数据集内部和GEO验证集中进行了评价。这6个ARGs中大多与OSCC或其他恶性肿瘤的发生和预后密切相关，有研究表明，在二乙基亚硝胺（diethylnitrosamine，DEN）诱导的原发性肝癌模型中，通过BID基因缺失抑制肝细胞死亡途径可以保护动物免受肿瘤发生[28]；BID蛋白是接受辅助治疗的结肠癌患者的独立预后变量[29]。恶性转化风险高的口腔白斑与NKX2-3基因的启动子异常甲基化有关[30]；在人黑色素瘤细胞系中能够观察到NKX2-3基因的高甲基化[31]。DDIT3通过上调胃癌中的CEBPβ促进肿瘤干细胞（cancer stem cells，CSCs）的干性，并与预后不良有关[32];DDIT3基因在内质网应激诱导的细胞凋亡和自噬中起重要作用[33-34]。VEGFA基因在OSCC组织中过表达,并且与患者较差的总生存率相关[35]；同时，血管生成增加往往会导致远处转移和预后不良[36]。FADD基因拷贝数和蛋白表达与中国台湾地区OSCC患者的淋巴结转移密切相关[37]。BIRC5/Survivin蛋白的细胞质表达与OSCC患者的总体低生存率有关，而核表达与较高的增殖率相关[38]。综上所述，风险评分模型中纳入的预后ARGs主要通过调节自噬和凋亡、内质网应激、血管生成等生物学过程参与肿瘤的发生、发展。同时，单因素、多因素Cox回归分析提示风险评分可作为独立的预后因素，ROC分析也表明，与其他临床危险因素相比，ARGs风险评分模型具有更好的预测性能。本研究还分析了风险评分与临床病理变量之间的相关性，风险评分的高低与肿瘤分期、T分期、N分期密切相关，但它们之间的因果关系仍需进一步探索。为了准确预测OSCC患者的生存率，本研究将风险评分与多种临床风险因素（年龄、性别、病理分级、肿瘤分期、T分期、N分期）结合，建立列线图用于生存概率的定量评估，校准曲线也证明了预测生存率与实际生存率之间的一致性。

最后，本研究分析了高、低风险组间免疫浸润细胞的差异，揭示高风险组预后差的潜在免疫机制。如，高风险组CD8+T细胞显著低于低风险组（P＜0.001），CD8+T细胞对于细胞内病原体和肿瘤的保护性免疫很重要，并与患者免疫治疗效果相关，CD8+T细胞减少往往预示着预后不佳[39]。结合图1E中KEGG富集的信息，OSCC中差异表达的ARGs与PD-L1表达和PD-1检查点通路相关，T细胞表面表达的PD-1与肿瘤细胞表面表达的PD-L1结合会向T细胞传递负向调控信号，抑制CD8+T细胞的产生和其破坏肿瘤细胞的能力[40]。因此，在本研究中，我们推测高风险组差异表达的ARGs可能参与调控PD-1和PD-L1的免疫检查点通路进而影响了CD8+T细胞免疫浸润，最终导致了不良的预后。另外，Carleton等[41]发现自噬失活Atg5-/-小鼠的肿瘤浸润淋巴细胞 （tumor infiltrating lymphocytes,TIL）中超过80%的CD8+TIL获得效应记忆表型，产生更多的干扰素γ（interferon-gamma,IFN-γ）和肿瘤坏死因子（tumor necrosis factor,TNF），而CD4+TIL效应记忆无变化，提示自噬缺失增强了CD8+T细胞效应记忆及功能；并发现自噬失活可阻止代谢下调H3K4me3水平，从而使免疫应答基因启动子的H3K4me3表达量得以维持。据此我们也可大胆假设本研究纳入的预后相关基因可能通过同样或类似的表观遗传修饰的生物学过程影响OSCC患者的免疫应答及CD8+T的产生。这些猜测需要进一步的实验加以论证，同时也为接下来研究自噬调控肿瘤发生、发展的潜在机制提供指引。

综上所述，本研究综合分析了ARGs，并且成功构建了基于6个预后相关ARGs的风险评分模型和列线图，其可以为OSCC患者的预后提供准确可靠的预测方法，帮助临床医生优化个性化治疗策略，同时也为OSCC肿瘤进展的潜在分子机制提供了一些思路。然而，本研究也存在一些局限性，首先，此次研究是回顾性的，因此风险评分模型和列线图需要在前瞻性研究和多中心临床试验中进一步验证；另外，本研究所鉴定的ARGs特异性生物学特性及其在OSCC肿瘤发生、发展中的潜在机制有待进一步的实验研究。