基于TCGA数据库胆管癌自噬-临床预后模型的构建和评估

邹文强，符广华，杨艳梅，张新宇

（1.哈尔滨医科大学附属第二医院 普通外科，2.哈尔滨医科大学 肿瘤防治研究所，黑龙江 哈尔滨 150086）

胆管癌（cholangiocarcinoma，CCA）是起源于胆管上皮细胞的恶性肿瘤，根据病变的解剖部位可以分为肝内CCA、肝周CCA和远端CCA[1]。CCA的发病率约占消化道肿瘤的3%，我国及东南亚地区是最高发的地区，而且近几十年发病率呈上升趋势[2]。由于CCA患者早期发病症状隐匿且恶性程度较高，一旦发现，多数肿瘤已进入晚期。研究表明根治性手术是最有效的治疗方式，但由于淋巴结转移早，易侵犯血管及神经，术后远期预后较差[3]。CCA患者的预后评估方法主要是TNM分期和组织学诊断，但无法实现准确的个体化生存预测[4-6]；且临床常规的肿瘤标志物往往敏感度和特异性不高，因此，寻找可靠的预后标志物来优化预后评估系统是必不可少的。

近年随着对肿瘤的深入研究，发现自噬在肿瘤的发生发展中发挥着重要作用[7-9]。自噬是一种通过溶酶体自我分解细胞内的老化及受损的蛋白质和细胞器，实现细胞内稳态的途径[10]。细胞自噬在正常生理状态下可以促进机体抵御疾病，但当自噬异常时又会导致肿瘤、自身免疫性及神经退行性等疾病[11]，在肿瘤中的作用取决于组织类型、肿瘤分期和致癌突变类型等因素[12]。自噬对包括CCA在内的多种肿瘤，具有双重调节作用。自噬通过消除受损的线粒体和减少ROS介导的染色体损伤，在肿瘤初期发挥抑癌功能。自噬通过抑制p53 的诱导和维持线粒体的代谢功能，在肿瘤进展期促进肿瘤细胞存活[13-15]。Sasaki等[16]的研究表明，CCA组织样本中自噬相关蛋白（LC3、beclin-1）表达量显著增加。CCA进展及转移的重要驱动因素之一是上皮细胞向间质细胞的转变（EMT），Chen等[17]发现调节自噬的信号通路（PI3K/AKT、p53和beclin1）同时对EMT有显著调节的作用，促使CCA细胞存活。这些研究证实自噬和CCA存在密切关系，并表明ATGs作为CCA预后标志物的巨大潜力。既往研究表明依据自噬基因构建的预后风险模型，可以有效评估膀胱癌、乳腺癌及胃癌等患者的预后情况[18-20]。

本研究拟通过癌症基因组图谱数据库（The Cancer Genome Atlas，TCGA）分析与CCA预后密切相关的ATGs，同时构建自噬基因预后风险模型，观察其对CCA患者的预后价值，从而为临床诊断和指导治疗CCA提供个性化的模型。

1 材料和方法

1.1 数据下载及整理

在TCGA中获取CCA的转录组数据和临床数据，包括35 例CCA组织样本和9 例正常组织样本。使用Strawberry Perl软件将各个独立样本的基因表达量汇总在一个文件中，并将基因的“ID”转换为“Symbol”。自噬相关基因（ATGs）从人类自噬数据库（Human Autophagy Database，HADb）中下载，包括232个ATGs的信息。在TCGA数据库已经下载的CCA样本信息中，利用Perl脚本将232个ATGs的基因表达量提取出来。

1.2 差异表达的自噬相关基因（DEATGs）的识别

通过R语言（R 4.0.2）使用“Willcox.test”来筛选DEATGs，筛选条件为：FDR＜0.05，|logFC|＞1。通过R语言绘制DEATGs的火山图和箱线图。

1.3 DEATGs的富集分析

基因本体功能（gene ontology，GO）富集分析和京都基因与基因组百科全书通路（Kyoto Encyclopedia of Genes and Genomes pathway，KEGG pathway）分析，可以明确DEATGs的主要生物学属性。调用R语言中“org.Hs.eg.db”“enrichplot”和“gplot2”包，进行GO富集分析和KEGG通路分析；使用R语言中“GOplot”包，绘制GO气泡图和KEGG圈图。

1.4 预后相关ATGs的识别和基于ATGs的预后风险评分模型构建

使用Strawberry Perl软件将DEATGs表达量和生存状态及时间，汇总在一个文件中。调用R软件的“survival”包，筛选标准为P＜0.05，对DEATGs进行单因素Cox回归分析，筛选出与CCA预后相关的自噬基因（ATGs）。再进行多因素Cox回归分析，删除相关性较高的基因来优化模型，最后得到用于构建预后风险评分模型的风险基因。使用R软件获取风险基因的评分系数（Coef）。模型的计算公式为：患者风险评分=（基因1的表达量×coef1+基因2的表达量×coef2+…）。按照患者风险评分，以风险值中位数为界限，将患者分为高风险和低风险组。

1.5 预后模型的验证

使用R 4.0.2软件进行生存分析、独立预后分析和临床相关性分析，绘制风险曲线和多指标ROC曲线。

2 结果

2.1 自噬相关基因的差异分析

从TCGA数据库中下载了35 个肿瘤组织样本和9 个正常组织样本的转录组数据和临床数据。在HADb数据库中提取了232 个自噬相关基因，并提取自噬相关基因表达量。按照FDR＜0.05和log2FC＞1 的筛选标准，过滤得到49 个自噬相关差异基因（DEATGs），其中上调的基因有48个，下调的基因有1个（图1）。

图1 49个差异自噬基因的箱线图

2.2 差异自噬基因的富集分析

对49个差异表达的自噬相关基因（DEATGs）进行GO富集分析和KEGG通路分析。经过GO富集分析，在生物学过程方面，基因主要富集在自噬、利用自噬机制的过程、宏观自噬、自噬调控、自噬信号通路的调控、细胞生长及内在凋亡信号通路等；在细胞成分方面，基因主要富集在晚期内涵体、细胞的焦点粘连、核被膜及液泡膜上；在分子功能方面上，基因主要富集在细胞黏附分子结合、泛素样蛋白连接酶结合、钙黏素结合及蛋白丝氨酸/苏氨酸激酶活性等（图2）。KEGG通路分析如图3，这些基因主要与自噬-动物、志贺氏症、人乳头瘤病毒感染、细胞凋亡、PI3K-AKT信号通路、铂类耐药、肝细胞癌、慢性粒细胞白血病、胰腺癌等通路显著相关。

图2 GO富集分析的气泡图

图3 KEGG通路分析的圈图和热图

2.3 预后相关自噬基因的筛选

将差异自噬基因的表达量与生存时间和生存状态合并，按照P＜0.05的筛选标准，通过单因素Cox回归分析，过滤得到5个与预后相关的差异自噬基因（RHEB、PPP1R15A、ATG101、BNIP3、NRG1），见图4。

图4 单因素Cox分析森林图

2.4 自噬基因的预后模型的构建

将这5个预后相关基因通过多因素Cox回归分析进行优化，删除相关性较高基因，结果这5个基因都与患者的预后显著相关，可以用来构建预后模型（见表1）。自噬基因的预后风险模型公式为（NRG1表达量×1.1207+BNIP3表达量×1.4002+ATG101表达量×1.3457-PPP1R15A表达量×1.1613-RHEB表达量×1.3279）。

表1 多因素Cox分析CCA预后影响因素

2.5 自噬预后模型的验证

根据模型计算每个患者风险值，18 例属于高风险组，18 例为低风险组。为了验证自噬预后模型在预测临床预后方面的可靠性，进行以下分析。Kaplan-Meier生存分析中P=8.52e-05（P＜0.05），因此高、低风险组的患者在生存率上有差异的可能性较大（图5），高风险组3 年生存率为18.7%，低风险组3年生存率为70.3%。图6A为风险曲线图，从左到右患者的风险值逐渐增大；图6B为生存状态图，描述了风险值与生存时间/生存状态的关系，随着风险值增大死亡患者增多，随着风险值减小死亡患者减少；图6C为风险热图，描述了风险值与基因表达量的关系，随着风险值的增加，4个基因的表达量增加（RHEB、ATG101、BNIP3、NRG1），表明这4个基因为高风险基因；随着风险值的增加，1个基因的表达量减少（PPP1R15A），表明这1个基因为低风险基因。在单因素独立预后分析中，风险值与生存时间及生存状态显著相关（HR1.382，95%CI1.160～1.645，P＜0.001）（图7A），在多因素独立预后分析中风险值与生存时间及生存状态也显著相关（HR1.427，95%CI1.161～1.756，P＜0.001）（图7B），表明风险值可以作为独立预后因子，进一步验证了Cox模型在预测临床预后方面的可靠性。在多指标ROC曲线中，风险值的AUC值最大（图8），表明风险得分有更好的预测临床预后的能力。

图5 生存曲线分析

图6 风险值与生存时间/生存状态、基因表达量的关系

图7 单因素（A）和多因素（B）分析患者预后影响因素

3 讨论

CCA作为临床常见的恶性肿瘤，预后极差，可能与肝血管的侵犯、远处转移、较早淋巴结的转移及术后早期复发有关，因此，寻找可靠的早期诊断标志物来完善预后评估系统，对于CCA患者预后有极大帮助。本研究在TCGA数据库中筛选出49个差异表达的自噬相关基因DEATGs，且通过富集分析证实自噬参与CCA发生发展及耐药的过程，本研究应用单因素及多因素Cox回归分析最终筛选出5个与CCA预后显著相关的自噬基因（RHEB、PPP1R15A、ATG101、BNIP3、NRG1）。

RHEB基因属于RAS家族的一员，编码的蛋白属于GTP结合蛋白，在与GTP结合的情况下Rheb被激活，当GTP被水解时Rheb受到抑制。在TSC1-TSC2-TBC1D7/RHEB/mTORC1经典的信号通路调节中，扮演重要角色。在参与调节自噬的发生发展、调控细胞分化及抑制肿瘤发面发挥重要作用[21]。在结肠癌组织中，p27KIP1受到RHEB的调控，RHEB/p27 激活mTORC1 信号通路来增强自噬，促进肿瘤细胞的存活，加速结肠癌的进展[22]。He等[23]的研究表明，FADD基因在乳腺癌中显著上调，FADD基因敲减降低了乳腺癌细胞系中Rheb的表达，抑制了mTORC1的活性，抑制肿瘤细胞的存活。Chen等[24]在宫颈癌组织样本中的研究表明，circMYLK与miR-1301-3p结合可以促宫颈癌细胞生长，其机制是RHEB/mTORC1信号通路的激活。

PPP1R15A基因又称为GADD34，在应激生长停止条件和DNA损伤剂治疗后其转录水平增加。PPP1R15A编码的蛋白质可以调节电离辐射后的凋亡，在促进细胞死亡和未折叠蛋白反应中起着重要的作用。Holczer等[25]在肝癌细胞中的研究发现，GADD34的表达被抑制后，导致mTOR的快速激活和自噬水平下调，随后凋亡细胞死亡。GADD34表达的减少显著抑制肿瘤，并导致MDSCs和T细胞的积累减少，抑制GADD34减少了MDSCs分泌的血管上皮生长因子α和转化生长因子β[26]。

ATG101基因位于染色体12q13.13，编码的蛋白是ULK1 自噬复合物的重要组成部分，与ATG13结合维持ULK1复合物的稳定性，进而调控自噬[27]。研究表明，PTCH1是已知的肿瘤抑制因子，抑制肿瘤细胞增殖和诱导细胞凋亡，减少PTCH1的表达量会增加ATG101介导的自噬通量来促进癌细胞增殖[28]。

BNIP3基因的全称为BCL2相互作用蛋白3基因，该基因编码一种线粒体蛋白，含有BH3 结构域，可以充当促凋亡因子。细胞内存在很多蛋白可以诱导线粒体自噬和凋亡，但BNIP3能将线粒体自噬和凋亡串联在一起，功能更加特殊[29]。m6A去甲基化酶（FTO）在乳腺癌中显著高表达，BNIP3是其下游靶点，促进肿瘤细胞增殖和转移[30]。Li等[31]在胰腺癌的研究中发现，BNIP3甲基化抑制了线粒体介导的胰腺癌细胞凋亡的诱导。

NRG1基因编码一种膜糖蛋白，介导多种细胞信号通路，该基因有多种异构体，该基因的失调与癌症等疾病密切相关。关于前列腺癌中一项研究发现，肿瘤微环境下产生的NRG1，可以促进肿瘤细胞的耐药性。研究表明，阻断NRG1会减少其介导的HER3的激活，抑制NRG1信号通路的诱导，起到抑制胰腺肿瘤细胞（PC）和肿瘤相关成纤维细胞（CAFS）相关的促肿瘤细胞增殖和转移的作用[32]。Trombetta 等[33]研究发现，肺腺癌中NRG1的高表达会通过PI3K-AKT和MAPK信号通路异常激活ErbB2/ErbB3，使肺腺癌细胞侵袭性增强。

本研究基于上述5 个自噬基因构建预后风险模型，发现风险模型的AUC值（0.906）要明显优于其他指标，因此可以作为CCA的独立预后指标。当然本研究还存在一定的局限性，譬如本研究属于纯数据库分析，缺少相关的功能实验进行验证。