基于TCGA数据库建立胆管癌自噬相关基因预后预测模型及其应用

史华帝，左瑜芳，钟富兰，易小琼，徐祖敏

广东医科大学附属医院，广东湛江524000

多数胆管癌患者确诊时已至晚期，失去手术机会，并且放化疗效果不佳。既往25 年的数据表明，胆管癌的发病率和病死率在增加，然而胆管癌的预后和治疗并没有明显改善，其5 年生存率仅为5%～10%，晚期胆管癌的中位生存期不超过12 个月［1］。因此，寻找影响胆管癌患者预后的相关因素具有重要的临床意义。美国癌症联合委员会（AJCC）的分期手册已经成为癌症患者分类、确定预后和确定最佳治疗方法的基准。然而随着对癌症研究的深入，人们发现TNM 分期系统并不足以评估预后，也不能反映癌症的生物学异质性。即使是同一分期的患者，其预后和治疗反应也有很强的异质性，其他因素如年龄、表现状况、肿瘤位置等也会影响患者的生存，因此它对临床预后提供的信息有限［2］。因此，我们迫切需要开发可靠的预后生物标志物，以提供更好的临床预后价值。自噬是通过吞噬细胞质蛋白或细胞器及其内含物进入囊泡，并与溶酶体融合，降解内含物的内容，实现细胞的代谢需要和某些细胞器的更新，从而形成自噬溶酶体的动态稳态过程［3］。这种分解过程涉及多种自噬相关基因（ARGs）［4］。在炎症、神经退行性变、衰老、肿瘤等病理条件下，自噬可以被抑制，以维持细胞内稳态［5］。基于ARGs构建的预后模型已应用于多种癌症类型［6-7］。2020年8月—2020 年9 月，本研究通过对TCGA 数据库中CHOL 队列的转录组和临床数据进行生物信息的综合分析，对胆管癌组织中差异表达的ARGs 进行筛选，并利用Lasso 和Cox 回归构建基于ARGs 的预后预测模型，为胆管癌患者的预后风险分层、预后评估及治疗策略的选择提供新的参考指标。

1 资料与方法

1.1 人类ARGs 集的获取 从人类自噬数据库（HADb，http：//autophagy. lu/clustering/index. html）中获取232 个ARGs，从分子特征数据库v6.2（MSig-DB，http//software. broadinstitute. org/gsea/msigdb）的GO_AUTOPHAGY 基因集中获取394 个ARGs。将两个数据库的ARGs 合并，删除重叠的ARGs，最后获得531个ARGs。

1.2 原始资料的下载与整理 在癌症基因组图谱数据库（TCGA，https：//portal. gdc. cancer. gov/）中选择CHOL 队列的转录组和临床数据进行下载。包含胆管癌组织36 例和正常胆管组织9 例。应用Perl5.28.1 软件将原始测序数据进行合并，转换基因的ID，并提取所有ARGs 的表达数据，若有多个探针检测同一个基因的表达量，则取该基因表达量的均值作为该基因的表达值。本研究纳入分析的临床数据完整，未删减患者的临床信息。

1.3 胆管癌组织中差异ARGs 的筛选及其GO 功能富集和KEGG 信号通路分析 应用R 3.6.1 软件的“limma”包对36 例胆管癌组织和9 例正常组织的ARGs 进行表达差异分析，筛选标准为：FDR＜0.05，∣logFC∣≥2，其中FDR 为错误发现率，FC 为差异倍数。利 用R 软 件 的“clusterprofiler”、“org. HS. eg.db”、“enrichplot”、“ggplot2”包进行基因本体（GO）功能富集，包括生物过程、细胞成分和分子功能。采用京都基因百科全书（KEGG）进行信号通路分析，并进行可视化。

1.4 胆管癌患者预后相关关键ARGs 的筛选 将胆管癌患者的生存信息与差异表达的ARGs 合并在一起，通过R 语言的“survival”包首先进行单因素Cox 回归分析，计算每个差异ARGs 与胆管癌患者生存的风险比（HR）和P 值，P＜0.05 的差异ARGs 被选出来进行下一步分析。应用Lasso 回归分析以减少基因之间共线性的影响，防止后续构建的预后风险模型变量过度拟合，通过构建惩罚函数对自变量的回归系数进行压缩，实现基因数据的降维，进而获得与胆管癌患者预后相关度更高的ARGs。最后，对Lasso 回归筛选得到的ARGs 进行多因素Cox 回归分析，计算每个ARGs 的多因素回归系数，得到与患者预后相关的关键ARGs。

1.5 风险预后模型的建立与分析 根据Cox 和Lasso回归筛选得到的关键ARGs，构建基于ARGs表达的风险评分方程［8-9］：风险评分Coefi×Xi，其中Coef是指基因在多因素Cox 回归分析中的回归系数，X 为基因的表达量，n 为预后相关ARGs 的数量。根据公式计算出每个患者的风险评分，取其中位数作为截断值，将胆管癌患者分为低风险评分组和高风险评分组。采用R 软件“survival”包进行Kaplan-Meier 生存曲线分析，应用“timeROC”包绘制模型的ROC 曲线，计算曲线下面积（AUC），评价模型的敏感度和特异度。利用“rms”包绘制列线图，为评估列线图实际生存和预测生存的一致性，进一步绘制校准曲线以评价模型的可靠性。

2 结果

2.1 胆管癌组织差异表达ARGs 的筛选 与正常胆管组织比较，胆管癌组织中有324 个存在差异表达的ARGs，其中表达上调的基因311 个、表达下调的基因13个。

2.2 胆管癌组织差异表达ARGs 功能富集和信号通路分析结果 GO功能富集分析显示，差异表达的ARGs 在生物学过程中包括自噬、利用自噬机制的过程、自噬的调节、线粒体的拆卸等，在细胞组分中包括内体膜、泡膜、晚期内体、自噬小体、溶酶体膜等，在分子功能中包括蛋白丝氨酸/苏氨酸激酶活性、泛素样蛋白连接酶结合、细胞黏附分子结合、钙黏着蛋白的结合、蛋白激酶调节活性、蛋白磷酸酶的结合等，见表1。KEGG 信号通路分析显示，差异表达的ARGs 主要涉及自噬—动物、线粒体自噬—动物、自噬—其他、志贺菌病、凋亡等信号通路，见表2。

表1 GO功能富集分析（前10位）

表2 KEGG信号通路分析（前10位）

2.3 胆管癌患者预后相关关键ARGs 的筛选结果将324个ARGs纳入单因素Cox回归分析，以P＜0.05为筛选条件，共发现12个基因可能与胆管癌患者的生存相关，分别是VPS25、VPS11、EVA1A、BNIP3、FAM13B、PPP1R15A、GABARAP、GMIP、VPS4B、UBC、FXR2、ATG101。Lasso 回归将上述单因素Cox筛选出来的ARGs 进行再次筛选，结果见图1A，并使用交叉验证建立模型，见图1B。最后纳入下一步 分 析 模 型 的ARGs 有VPS11、EVA1A、BNIP3、GABARAP、VPS4B。

2.4 胆管癌患者预后预测模型的构建结果 将上述筛选出来的5 个ARGs 纳入多因素Cox 回归分析，最终得到这5 个关键ARGs 组成的胆管癌患者预后预测模型，根据5 个ARGs 的mRNA 表达水平和风险系数计算每个患者的风险评分。风险评分=（-3.739×VPS11）+（1.691×EVA1A）+（1.734×BNIP3）+（5.776×GABARAP）+（-1.310×VPS4B）。取风险评分的中位数0.96为截断值，将36例患者分为低风险组（风险评分＜0.96，18 例）和高风险组（风险评分＞0.96，18例）。

2.5 胆管癌患者预后预测模型的预测价值 Kaplan-Meier生存分析显示，高风险组总的生存时间低于低风险组（P＜0.05），见图2A。ROC 曲线分析显示，预测1、2、3年生存率的ROC曲线下面积（AUC）分别为0.973、0.989、0.941，均大于0.9。见图2B。

图1 Lasso筛选ARGs和交叉验证结果

图2 生存曲线与ROC曲线

2.6 胆管癌患者预后预测模型的可靠性 利用R软件构建基于5 个ARGs 的列线图，见图3A。该列线图模型经评价得到C指数为0.822，95%可信区间为0.721～0.924，表明该模型具有较好的区分能力。为评估列线图实际生存和预测生存的一致性，我们绘制预测1年、2年、3年生存率的校准曲线几乎都落在了45°的对角线上（见图3B、3C、3D），提示该模型准确性较高。

3 讨论

目前对于胆管癌发生发展的分子机制尚不明确。近年研究发现，自噬参与了胆管癌的发生发展［10-11］。然而自噬在肿瘤发生中的作用目前仍存在争议，对于自噬是否影响胆管癌患者的预后，目前尚无定论。近年来，基因签名经常用于预测各种肿瘤的预后，在一定程度上甚至优于TNM 分期和组织病理学诊断［12-13］。基于ARGs 表达的预后模型已被报道用于多种癌症，如肺癌、胃癌、肠癌、乳腺癌、胶质瘤等［6，14-17］。本研究中，我们通过Cox和Lasso回归模型分析了531 个ARGs，最终筛选出5 个预后相关的关键ARGs，并成功构建了胆管癌预后预测模型，结果显示高风险组总的生存时间明显低于低风险组，预测胆管癌患者1 年、2 年、3 年生存率的ROC 曲线下面积AUC 均大于0.9，提示该模型是具有较好敏感性和特异性的理想模型。

本研究中，我们在TCGA 数据库中对胆管癌的ARGs 进行差异表达分析，通过GO 功能富集和KEGG 信号通路分析发现，差异表达的ARGs主要涉及的生物学过程和信号通路有自噬、凋亡、线粒体、志贺菌病等。我们筛选出的5 个预后相关的关键ARGs，即 VPS11、EVA1A、BNIP3、GABARAP、VPS4B，均被报道与肿瘤的发生发展相关，涉及多种肿瘤细胞的自噬、凋亡、铂类药物耐药等生物学过程。PETERSON 等［18］报道，VPS11 基因与黑色素瘤细胞的顺铂耐药有关。SHEN 等［19］报道，EVA1A 诱导的自噬和细胞凋亡在抑制GBM 的发展中起作用，在3 种脑胶质瘤细胞系（U251、U87 和SHG44）中，EVA1A 过表达通过激活自噬和诱导细胞凋亡，EVA1A 过表达也成功抑制了NOD/SCID 小鼠的肿瘤生长。EVA1A 也被证实与肝癌、多发性骨髓瘤、乳腺癌有关［20-22］。LI 等［23］报道，lncRNA DGCR5/miR-27a-3p/BNIP3 轴可通过调节p38/MAPK 通路促进胰腺癌细胞凋亡。BORTNIK 等［24］报道，GABARAP 高表达与所有乳腺癌亚型的侵袭性疾病表型的不良预后和临床病理特征相关。SZYMANSKA 等［25］发现，结直肠癌组织中VPS4B mRNA 和蛋白表达下调，敲除VPS4B 基因可诱导细胞死亡。以上研究表明，这些关键ARGs 在肿瘤发生发展中起重要作用。然而，我们对这些ARGs 在胆管癌中的了解仍然很有限，对于自噬如何影响胆管癌患者预后的分子机制及其对胆管癌患者临床治疗的意义有待进一步研究。

图3 列线图和校准曲线

准确判断癌症患者的预后有助于指导临床医师进行治疗决策。近年来，列线图作为肿瘤预后评估的工具在医学研究和实践中得到了广泛的应用［26-28］。列线图的原理是通过多因素回归模型（例如Cox、Logistic 回归等），把复杂的回归方程转变为可视化的图形后，临床医师可以根据模型中每个风险变量的评分，然后相加后得到总评分，很方便的计算出疾病发生的概率及判断患者预后的好坏。模型的区分能力常常用C 指数表示，其小于0.65 表示区分能力较差，0.65～0.75 表示模型的区分能力一般，大于0.75 则表示模型的区分能力较好［29］。我们模型的C指数为0.822，提示该模型对预后有较好的区分能力。除了区分能力外，校准曲线也是评价模型准确性的重要指标，它能够反映模型的预测风险与实际发病风险的一致程度，理想模型的校准曲线刚好落在45°度的对角线上，校准曲线越匹配对角线，预测准确性越高［30］。从图3B可以看出，我们构建模型的校准曲线几乎都落在了对角线上，提示该模型的准确性较高。然而，本研究也存在一定的局限性，我们的模型缺少独立的外部数据集进行验证。另外，模型真实的评估性能需要在实践中进行验证。

综上所述，胆管癌中差异表达的ARGs 主要涉及自噬、凋亡、mTOR 信号通路、铂耐药、PI3K-Akt 信号通路、HIF-1 信号通路等癌症相关通路，通过筛选与预后相关的ARGs，我们成功构建了基于VPS11、EVA1A、BNIP3、GABARAP、VPS4B 共5 个ARGs 表达的风险预测模型，该模型具有较高敏感度和特异度，可有效预测胆管癌患者的预后，对指导临床治疗、制定精准治疗方案具有重要价值。