基于TCGA 数据库建立胆管癌自噬相关基因预后预测模型及其应用 史华帝，左瑜芳，钟富兰，易小琼，徐祖敏广东医科大学附属医院，广东湛江524000
摘要：目的利用癌症基因组图谱（TCGA ）数据库构建胆管癌自噬相关基因（ARGs ）预后预测模型。方法通
过人类自噬数据库和分子特征数据库获得531个胆管癌ARGs 。从TCGA 数据库中选择CHOL 队列的转录组和临床数据进行下载，包含胆管癌组织36例、正常胆管组织9例。采用Perl 软件将原始测序数据进行合并，提取所有ARGs 的表达数据；利用R 软件对胆管癌组织和正常胆管组织中的ARGs 进行差异表达分析，筛选出胆管癌组织中表达失调的ARGs ，并进行GO 功能富集和KEGG 信号通路分析。利用单因素Cox 及Lasso 回归模型筛选关键ARGs ，多因素Cox 回归模型建立预后预测模型，根据关键ARGs 的mRNA 表达水平和风险系数计算每个患者的风险评分，按其中位数将胆管癌患者分为高风险组和低风险组，绘制生存曲线，比较两组的生存期。绘制预后预测模型的ROC 曲线，评价其预测预后的敏感度和特异性。最后利用R 软件构建基于关键ARGs 的列线图，绘制校准曲线评估实际生存和预测生存的一致性。结果 与正常胆管组织比较，胆管癌组织中有324个表达失调的ARGs 。这些ARGs 主要涉
及自噬、利用自噬机制的过程、大自噬、自噬的调节、凋亡等生物学过程和信号通路。经过单因素Cox 和Lasso 回归分析，筛选出5个关键ARGs ，即VPS11、EVA1A 、BNIP3、GABARAP 、VPS4B 。以这5个关键ARGs 建立预测胆管癌患者预后的风险模型，风险评分=（-3.739×VPS11）+（1.691×EVA1A ）+（1.734×BNIP3）+（5.776×GABARAP ）+（-1.310×VPS4B ）。生存分析显示，高风险组的总生存时间低于低风险组，预测1年、2年、3年生存率的ROC 曲线下面积均大于0.9。构建了基于5个ARGs 的列线图（C 指数为0.822，95%CI 为0.721~0.924），另外，绘制预测1年、2年、3年生存率的校准曲线几乎都落在了45°的对角线上，提示该模型的准确性和区分能力较好。结论成功构建
了基于VPS11、EVA1A 、BNIP3、GABARAP 、VPS4B 共5个关键ARGs 表达的胆管癌预后风险预测模型，该模型可有效
预测胆管癌患者的预后。
关键词：胆管癌；预后风险模型；列线图；癌症基因组图谱数据库；自噬相关基因doi ：10.3969/j.issn.1002-266X.2021.02.002中图分类号：R735.8
文献标志码：A
文章编号：1002-266X （2021）02-0006-06
Construction of prognosis prediction model for cholangiocarcinoma based on autophagy -related genes of TCGA database and its application
SHI Huadi ，ZUO Yufang ，ZHONG Fulan ，YI Xiaoqiong ，XU Zumin
The Affiliated Hospital of Guangdong Medical University ，Zhanjiang 524000，China
Abstract ：Objective
To construct a survival model for predicting the prognosis of patients with cholangiocarcinoma
（CCA ）based on autophagy -related genes （ARGs ）of The Cancer Genome Atlas （TCGA ）database.Methods
A total of
531ARGs were obtained from the Human Autophagy Database and Molecular Signatures Database.The original expres‑
sion profiles and corresponding clinical data of CCA patients were downloaded from the CHOL cohort of the TCGA data‑base.There were 36cases of cholangiocarcinoma and 9cases of normal bile duct.Perl software was used to merge the orig‑inal sequencing data and extract the expression data of all ARGs.Differential expression analysis of CCA tissue and normal tissue was performed using R software to screen out ARGs with aberrant expression.The GO functional enrichment and
KEGG signaling pathway analysis were carried out using R software.ARGs were submitted to Lasso and univariate Cox re‑gression analyses to remove the genes which might not be related to the prognosis of CCA patients.Multivariate Cox regres‑sion model was used to establish the prognostic model and we calculated the risk score of each patient according to the
基金项目：广东省自然科学基金（2020A1515010048）；广东省中医药管理局科研项目（20201179）；湛江市非资助科技攻关项目（2019B01021）。
第一作者简介：史华帝（1986-），男，主治医师，主要研究方向为恶性肿瘤的基础和临床。E -mail ：shi_ 通信作者简介：徐祖敏（1982-），男，博士后，副主任医师，主要研究方向为恶性肿瘤的基础和临床。E -mail ：zuminxu@163.
com
开放科学（资源服务）
标识码（OSID ）
6
mRNA expression and risk coefficient of key ARGs.CCA patients were divided into the high-risk group and low-risk group according to the median of risk scores.The Kaplan Meier survival curve was plotted to analyze the median survival time. Time-dependent receiver operating characteristic（ROC）curve was drawn to investigate the sensitivity and specificity of the model.Finally，R software was used to construct the nomogram based on the key ARGs.The calibration curve was plotted to evaluate the consistency between the actual survival and the predicted survival.Results Compared with the normal tissues，there were324ARGs with aberrant expression in CCA tissues，
and these ARGs were mainly involved in autophagy，process utilizing autophagic mechanism，macroautophagy，regulation of autophagy，apoptosis，etc.After uni‑variate Cox and Lasso regression analysis，five key ARGs were selected，namely VPS11，EVA1A，BNIP3，GABARAP，and VPS4B.A risk model based on5ARGs was conducted for predicting the prognosis of CCA patients.The risk score was quantified by the following formula：risk score=（-3.739×VPS11）+（1.691×EVA1A）+（1.734×BNIP3）+（5.776×GABARAP）+（-1.310×VPS4B）.The K-M survival curve showed that the overall survival of CCA patients was significant‑ly lower in the high-risk group than in the low-risk group.The area under the ROC curve（AUC）for predicting1-year，2-year and3-year survival rates was all greater than0.9.A nomograph based on these5ARGs was constructed（C-index：0.822，95%CI：0.721-0.924）.In addition，the calibration curve for predicting1-year，2-year and3-year survival rates almost all fell along the45-degree line.It indicated that the model had good accuracy and distinguishing ability.Conclu⁃sion A risk prediction model based on VPS11，EVA1A，BNIP3，GABARAP，and VPS4B was successfully constructed，which could effectively predict the prognosis of CCA.
Key words：cholangiocarcinoma；prognostic risk model；nomogram；The Cancer Genome Atlas；autophagy-relat‑ed genes
多数胆管癌患者确诊时已至晚期，失去手术机会，并且放化疗效果不佳。既往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///gsea/msigdb）的GO_AUTOPHAGY基因集中获取394个ARGs。将两个数据库的ARGs合并，删除重叠的ARGs，最后获得531个ARGs。
1.2原始资料的下载与整理在癌症基因组图谱数据库（TCGA，https：//portal.v/）中选择CHOL队列的转录组和临床数据进行下载。包含胆管癌组织36例和正常胆管组织9例。应用Perl5.28.1软件将原始测序数据进行合并，转换基因的ID，并提取所有ARGs的表达数据，若有多个探针检测同一个基因的表达量，则取该基因表达量的均值作为该基因的表达值。本研究纳入分析的临床数据完整，未删减患者的临床信息。
1.3胆管癌组织中差异ARGs的筛选及其GO功能富集和KEGG信号通路分析应用R3.6.1软件的“limma”包对36例胆管癌组织和9例正常组织的ARGs进行表达差异分析，筛选标准为：FDR<0.05，∣logFC∣≥2，其中FDR为错误发现率，FC为差异倍数。利用R软件的“clusterprofiler”、“
7
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］：风险评分=∑i=1n Coef i×X i，其中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。
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）。
表1GO功能富集分析（前10位）
基因集名称
自噬
利用自噬机制的过程
大自噬
自噬的调节
液泡组织
大自噬的调节
自噬体组织
线粒体分解
自噬体装配
细胞器分解
基因数目
158
158
119
75
54
53
40
37
37
39
校正后的
P值
2.74E-167
2.74E-167
3.83E-136
1.84E-61
8.20E-53
4.46E-50
2.45E-43
3.02E-43
3.02E-43
1.73E-42
错误发现率
1.85E-167
1.85E-167
2.58E-136
1.24E-61
5.52E-53
3.00E-50
1.65E-43
2.03E-43
2.03E-43
1.16E-42
表2KEGG信号通路分析（前10位）
基因集名称
自噬—动物
线粒体自噬—动物
自噬—其他
志贺菌病
卡波西肉瘤相关疱疹病毒感染
寿命调节通路
沙门氏菌感染
凋亡
胰腺癌
肌萎缩侧索硬化
基因
数目
58
27
20
43
31
21
30
24
18
36
校正后的
P值
9.22E-53
1.07E-22
6.64E-22
3.00E-21
4.78E-14
1.22E-12
6.89E-12
1.21E-11
5.62E-11
7.68E-10
错误发现率
4.60E-53
5.36E-23
3.32E-22
1.50E-21
2.39E-14
6.11E-13
3.44E-12
6.06E-12
2.81E-11
3.84E-10
8
取风险评分的中位数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。
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，均被报道与肿瘤的发生发展相关，涉及多种肿瘤细胞的自噬、凋亡、铂类药物耐药等生物学过
注：实线代表高风险组，虚线代表低风险组。
图2生存曲线与ROC曲线
图1Lasso筛选ARGs和交叉验证结果
9
程。PETERSON等［18］报道，VPS11基因与黑素瘤细胞的顺铂耐药有关。SHEN等［19］报道，EVA1A诱导的自噬和细胞凋亡在抑制GBM的发展中起作用，在3种脑胶质瘤细胞系（U251、U87和SHG44）中，
EVA1A过表达通过激活自噬和诱导细胞凋亡，EVA1A过表达也成功抑制了NOD/SCID小鼠的肿瘤生长。EVA1A也被证实与肝癌、多发性骨髓瘤、乳腺癌有关［20-22］。LI等［23］报道，lncRNA DGC
R5/miR-27a-3p/BNIP3轴可通过调节p38/MAPK通路促进胰腺癌细胞凋亡。BORTNIK等［24］报道，GABARAP高表达与所有乳腺癌亚型的侵袭性疾病表型的不良预后和临床病理特征相关。SZYMANSKA等［25］发现，结直肠癌组织中VPS4B mRNA和蛋白表达下调，敲除VPS4B基因可诱导细胞死亡。以上研究表明，这些关键ARGs在肿瘤发生发展中起重要作用。然而，我们对这些ARGs在胆管癌中的了解仍然很有限，对于自噬如何影响胆管癌患者预后的分子机制及其对胆管癌患者临床的意义有待进一步研究。
准确判断癌症患者的预后有助于指导临床医师进行决策。近年来，列线图作为肿瘤预后评估的工具在医学研究和实践中得到了广泛的应用［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表达的风险预测模型，该模型具有较高敏感度和特异
图3列线图和校准曲线10
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议