基于TGGA数据库构建肺腺癌自噬相关基因预后风险模型

王雪芹，刘亚锋，吴 静，周家伟，邢应如，张 鑫，李丹婷，谢 军，丁选胜，胡 东

自噬是维持哺乳动物体内稳态及细胞应激适应性反应的关键组成部分。在饥饿环境中，自噬捕获并降解溶酶体中的细胞器和蛋白质从而保护生物体免受疾病侵害，主要包括心脏病、感染和神经变性[1]。但有研究[2]表明，黑素瘤和乳腺癌在内的癌症会诱导自噬并依赖自噬生存。自噬介导的循环维持能量稳态和线粒体功能可以促进肿瘤细胞的生长和增殖。因此，抑制自噬功能可能对癌症治疗产生有益效果。

肺癌占世界新癌症发病率和病死率的第一位[3]。近年来肿瘤切除术、免疫治疗和放化疗综合应用于治疗肺腺癌，但是肺腺癌患者平均5年生存率仍低于20%[4]。因此，鉴定疾病特异性生物标志物，提高对肺腺癌分子机制的了解，开发新的治疗决策，以改善患者的预后是非常有必要的。

1 材料与方法

1.1 数据采集于2020年8月从癌症基因组图谱(the cancer genome atlas,TCGA)数据库(https://tcga-data.nci.nih.gov/tcga/)下载肺腺癌RNA-seq数据以及临床数据。从人类自噬数据库(human autophagy database,HADBb)(http：//www.autophagy.lu/)获得了232个自噬相关基因(autophagy related genes，ARGs)。TCGA数据库中肺腺癌患者的基本临床特征：女性213例(45.7%)，男性254例(54.3%)；65岁以上有227例(48%)，65岁及以下有240例(52%)；临床Ⅰ、Ⅱ、Ⅲ期分别为251例(54%)、108例(23%)、75例(16%)；T1、T2、T3、T4型分别为162例(34%)、244例(52%)、39例(8%)、19例(6%)；M0、M1型分别为314例(92%)、24例(8%)，N分型中N0、N1、N2～3型分别为300例(66%)、87例(19%)、68例(15%)。

1.2 鉴定差异表达自噬基因自TCGA数据库下载肺腺癌RNAs表达谱，并从中提取人自噬相关基因的表达数据，使用R软件“limma”软件包对正常样本和癌症样本鉴定差异表达自噬基因(DEARGs)，通过“pheatmap”软件包绘制热图。

1.3 自噬相关基因预后模型的构建用R语言“survival”程序包进行单因素Cox分析，并将初步筛选的基因与肺腺癌预后相关的自噬基因(P<0.05)纳入多因素Cox分析，以进一步验证基因对肺腺癌患者预后的价值。经过上述分析，获得风险评分模型。风险评分计算如下：

其中Coef是系数，x是每个选定自噬基因的表达水平。

1.4 高、低风险组之间免疫细胞浸润程度的比较为探讨风险评分与肺腺癌预后之间的关系，对风险评分与肿瘤微环境之间的关系进行分析。使用R软件ESTIMATE算法确定每个样品的免疫评分，通过Wilcoxon检验比较高危组和低危组之间免疫细胞浸润程度差异。此外使用CIBERSORT软件包评估22种免疫细胞亚型比例。

2 结果

2.1 筛选差异表达自噬基因使用R软件整理TCGA数据库中497个肿瘤样本及54个肿瘤患者的癌旁正常样本。转录组数据并从中提取232个自噬基因(HADb数据库获得)的表达数据，进行差异表达分析，筛选出30个差异自噬基因(FDR<0.05且|log2FC|>1)。其中有18个上调的ARGs[含有杆状病毒IAP重复序列5(baculoviral IAP repeat containing 5,BIRC5)、CDKN2A、TMEM74、内质网氧化还原酶1-α(endoplasmic reticulum oxidoreductase 1 alpha，ERO1A)、ATG9B、GAPDH、PTK6、EIF4EBP1、整合素4 (integrin subunit beta 4,ITGB4)、VMP1、IFNG、ATIC、BNIP3、IKBKE、P4HB、PARP1、ERBB2和HSPA5]和12个下调的ARGs[MAP1LC3C、GRID1、CCL2、PRKCQ、NRG1、FOS、PPP1R15A、DAPK2、HSPB8、NLR家族含有CARD结构域4 (NLR family CARD domain containing 4,NLRC4)、DLC1和NRG3]。见图1。

图1 肺腺癌和正常组织间差异表达的自噬基因

2.2 肺腺癌患者预后相关自噬基因鉴定及风险评分模型构建为初步筛选与肺腺癌预后相关的自噬基因，对差异表达的30个自噬mRNAs进行单因素Cox回归分析，共发现有9个自噬基因与患者预后相关(P<0.05)(图2A)。经过多变量Cox回归分析确定4个自噬基因用于构建预后风险模型(图2B)。风险评分计算方法如下：风险评分=(0.123 746×BIRC5表达量)+(0.324 086×ERO1A表达量)+(0.121 162×ITGB4表达量)+(-0.547 430×NLRC4表达量)。

通过风险预后模型来计算497例患者的风险值并以中位值将患者分为高、低风险两组。HR(风险比)>1的基因(BIRC5、ERO1A、ITGB4)被认为是危险基因，而HR<1的基因(NLRC4)被认为是保护性基因。随风险值升高，病例病死率呈升高趋势(图2C)。高风险组患者表达危险基因增多，表达保护基因减少(图2D)。

图2 预后相关自噬基因筛选及表达特征A：自噬相关基因的单变量Cox回归分析；B：自噬相关基因的多变量Cox回归分析；C：风险评分和患者生存时间以及生存/死亡状态的分布；黑色虚线：将患者分为低风险和高风险组的最佳临界值；D：低风险和高风险组自噬相关基因表达谱

2.3 风险评分模型中基因的蛋白表达水平为了进一步确定这4个自噬基因在肺腺癌中的表达，利用人类蛋白质图谱数据库(HPA)的免疫组织化学数据加以验证，结果显示，与正常肺组织相比，肺腺癌中的BIRC5、ERO1A、ITGB4表达增加，而NLRC4在肿瘤细胞中则未发现(图3)。肺腺癌中BIRC5、ERO1A、ITGB4、NLRC4蛋白水平的变化与其基因水平的变化一致。

图3 风险评分模型中4个基因的蛋白表达水平A:BIRC5；B:ERO1A；C:ITGB4；D:NLRC4

2.4 风险评分模型的性能检测Log-rank检验对低、高风险组进行Kaplan-Meier分析，结果表明高风险组的总生存期较差(P<0.000 1)，见图4A。ROC曲线下面积(AUC)用来评估风险评分模型预测效能，该风险评分模型预测患者1年生存率的AUC值为0.757(95%CI：0.727～0.787),3年生存率的AUC值为0.723(95%CI：0.667～0.779)，说明该模型具有良好的短期预测能力(图4B)。然而，5年生存率的AUC值仅为0.581(95%CI：0.516～0.646),提示该模型不能很好地预测患者长期(5年以上)生存状态。采用ROC曲线比较风险评分模型与四个单基因及临床参数的预测能力，风险评分模型的AUC为0.757(95%CI：0.727～0.787)；BIRC5的AUC为0.667(95%CI：0.609～0.725)；ERO1A的AUC为0.700(95%CI：0.649～0.751)；ITGB4的AUC为0.621(95%CI：0.555～0.687)；NLRC4的AUC为0.329(95%CI：0.258～0.400)；年龄的AUC为0.511(95%CI：0.442～0.580)；性别的AUC为0.578(95%CI：0.512～0.644)；疾病分期的AUC为0.721(95%CI：0.670～0.772)；T分型的AUC为0.671(95%CI：0.617～0.725)；M分型的AUC为0.504(95%CI：0.434～0.574)；N分型的AUC为0.674(95%CI：0.622～0.726)(图4C、4D)。Delong检验表明，风险评分模型的AUC高于年龄(P=0.033)、性别(P=0.039)、T分型(P=0.047)、M分型(P=0.031)及N分型(P=0.042)的AUC，但是风险评分模型与疾病分期的AUC差异无统计学意义(P=0.083)，其原因可能是样本量不足。由此可见，无论与单基因或者其他临床参数相比，该模型具有较好的预测能力和准确性。

图4 风险评分模型高、低风险组Kaplan-Meier曲线、ROC曲线及临床参数的相关性A:基于4个自噬基因的KM曲线；B:基于4个自噬基因ROC曲线；C:基于4个自噬基因与单独自噬基因的ROC曲线；D:基于4个自噬基因与临床参数的ROC曲线

2.5 高、低风险组之间的免疫浸润评分免疫评分显示高、低风险评分组之间存在差异(图5A)。Kaplan-Meier分析显示，肺腺癌中不同的免疫评分具有不同的生存期(图5B)。肺腺癌的22种免疫细胞比例，其中M0巨噬细胞、M2巨噬细胞和T细胞CD4、CD8为肺腺癌浸润主要免疫细胞亚群(图5C)。在高、低风险组中，有10群免疫细胞水平存在差异，其中CD4记忆静止细胞、CD4记忆激活细胞、单核细胞、M0巨噬细胞、静止树突状细胞、肥大静止细胞及肥大激活细胞差异最为明显(图5D)。

图5 风险评分模型与肿瘤免疫微环境的相关性A:免疫评分、基质评分和风险评分的关联；B:免疫评分对肺腺癌总体生存影响的KM分析；C:肺腺癌中22种免疫细胞的比例；D:高风险组和低风险组之间免疫细胞差异性表达；与高风险组比较：*P<0.05，**P<0.01，***P<0.001，****P<0.000 1

3 讨论

自噬是肿瘤起始和癌症进展的抑制因子，而在肿瘤晚期自噬作为促进生存、生长和侵略性的保护因子。越来越多的证据表明自噬与癌症密切相关。如NLRC5通过抑制自噬能够促进子宫内膜癌细胞的增殖[5]；孕酮和脂联素受体3(PAQR3)通过调节表皮生长因子受体(EGFR)介导的自噬抑制非小细胞肺癌细胞生长[6]；下调自噬通过抑制GLI系列锌指2(GLI2)的产生，从而促进Shh途径抑制的抗肺腺癌功效[7]。

本文基于TCGA数据库中的RNA-seq转录组数据及临床资料，结合HADBb进行生物信息学分析，表明BIRC5、ERO1A和ITGB4是肺腺癌危险因素，NLRC4是肺腺癌保护因素。BIRC5能够通过启动抗凋亡机制引起肺癌细胞增殖、肿瘤形成并可增强肿瘤细胞对化疗和放射治疗的耐药性[8-9]；ERO1A能够通过上调血管内皮生长因子，激活肿瘤血管生成信号从而促进肿瘤的生长[10]；ITGB4能够调控肿瘤、呼吸系统疾病和神经系统病变等多种疾病的进展[11]。NLRC4通过抑制细胞增殖和促进细胞死亡来预防结直肠癌[12]。本研究利用自噬相关基因结合风险评分构建肺腺癌预后预测模型，总结出高风险组的总生存期较差，且无论与单基因或其他临床参数相比，该模型能较好地预测肺腺癌患者预后。

目前，免疫疗法是肺腺癌的一种治疗方案，已有研究证明免疫细胞在肿瘤细胞上的作用。此外，肺腺癌组织中的免疫成分可用于评估患者的治疗效果和肺腺癌预后[13]。该研究结果表明，4个自噬基因组合的风险评分与免疫浸润评分具有相关性。使用CIBERSORT算法计算免疫细胞亚型，结果显示，高、低风险评分组表达不同的免疫细胞亚型，表明免疫细胞成分比例失衡与癌症患者的不良预后和低存活率相关。肺腺癌中的免疫细胞主要包括M0巨噬细胞、M2巨噬细胞、CD4+记忆静止T细胞和CD8+T细胞。其中，高风险组CD4+记忆静止T细胞和M2巨噬细胞水平较低，而M0巨噬细胞水平较高。这暗示着高危组中CD4+记忆静止T细胞与M0巨噬细胞比例的失衡可能降低患者的生存率。