自噬相关基因的表达量化及其对肾透明细胞癌的预后价值分析

伍艳婷,舒文莹,王艳平,周毅,简晓顺

在全球大多数国家,恶性肿瘤是导致死亡的主要原因[1],迄今为止,恶性肿瘤还没有绝对的治愈方法。肾细胞癌是泌尿系统最常见的肿瘤之一,具有隐匿性,对放化疗不敏感。其中肾透明细胞癌(kidney renal clear cell carcinoma,KIRC)是肾细胞癌最常见的亚型,占70%～75%,易发生转移[2]。针对多种激酶和雷帕霉素靶点的抑制剂可提高KIRC的总体存活率,但治疗反应因个体而异,大多数患者最终会产生耐药性[3]。因此,探索分析对个体化治疗或临床诊断有意义的生物标志物,对KIRC的早期诊断和提高生存率十分重要。

自噬在肿瘤的发生、耐药性、微环境相互作用等多个过程中发挥着至关重要的作用。自噬的发生是一个复杂的多步骤过程,受一系列自噬相关基因的严格控制,既可以抑制肿瘤的进展亦能促进肿瘤的进展[4],恶性肿瘤可以诱导微环境的非肿瘤细胞自噬,从而进一步促进肿瘤的生长和进展[5]。自噬相关的化疗耐药性是癌症治疗的一大挑战。目前,自噬相关基因在多药耐药中的作用越来越受重视,研究发现,自噬基因有望成为癌症的治疗靶点[6-8]。

但由于传统研究的局限性,大多集中在单基因对肿瘤影响的研究上,很少使用大规模基因表达队列研究自噬与KIRC预后的关系。近年来,随着自噬在肿瘤中的作用不断被研究,研究人员发现自噬相关基因可以作为预后生物标志物并成为潜在的治疗靶点[9]。为了更好地探索自噬相关基因对临床预后的影响,利用大型公共数据库筛选和验证KIRC自噬相关预后标志物已成为基础研究的重要补充和发展趋势。

本研究通过构建一个新的预后模型用于预测KIRC患者的总生存期(OS),并探索了由该模型基因定义的独特分子特征,包括免疫浸润和潜在的药物靶点。这为KIRC预后、免疫治疗和药物敏感性评价提供了一种新的思路和方法,报道如下。

1 资料与方法

1.1 数据收集 从UCSC xena数据库(https://xenabrowser.net/datapages/)下载24种常见癌种的表达数据log2(FPKM+1)(FPKM:每千个碱基的转录每百万映射读取的fragments)、表型数据和生存数据,通过R软件4.1.2将其表达数据转换为TPM(每千个碱基的转录每百万映射读取的Transcripts)数据;从人类自噬数据库(http://autophagy.lu/)中下载222个自噬相关基因,于GEPIA数据库(http://gepia.cancer-pku.cn/)中做与24种常见癌种的生存分析,最终得到159个与癌症预后相关的自噬相关基因(ARGs)。

1.2 自噬基因集在癌症中的表达情况 使用R包“GSVA”进行单样本基因集富集分析(ssGSEA)得到自噬评分(autophagy score)以量化ARGs在24种癌症中的表达水平,Kruskal-Wallis分析各癌症间的表达是否存在差异;Mann-WhitneyU检验评估23种癌症的肿瘤样本与正常样本间的差异[其中肾上腺皮质癌(ACC)无正常样本,故未纳入后面的分析];以自噬评分的中位数来定义高评分组与低评分组。

1.3 筛选KIRC差异表达基因 使用“edgeR”包比较KIRC患者的肿瘤样本和正常样本的基因表达差异(来自UCSC xena数据库的counts数据),筛选条件为|logFC|> 1,错误发现率(FDR)<0.05。

1.4 预后模型的构建 将531例KIRC患者数据随机分为训练集354个(占总KIRC数据的2/3)和测试集177个(占总KIRC数据的1/3);将159个ARGs与KIRC的差异表达基因取交集得到差异表达自噬相关基因(DEARGs);对训练集进行逐步回归分析、LASSO回归分析确定最终模型基因,根据多因素Cox回归分析的回归系数和基因的表达量,建立风险预后模型,建立预后模型公式为:风险评分(risk score)=∑(基因表达量×回归系数),测试集用于验证其预后能力。

1.5 预后模型的评估和验证 根据预后模型公式计算风险评分,以风险评分的中位数来定义高风险组和低风险组。绘制生存曲线和受试者工作特征曲线(ROC),并计算曲线下面积(AUC)以评估模型的特异度和敏感度。根据风险评分值进行排序,分布曲线、散点图和热图显示不同风险组的风险评分分布,患者生存状态及模型基因的表达情况。将风险评分与肿瘤大小(T)、淋巴结转移(N)、远处转移(M)、种族(Race)、年龄(Age)、性别(Sex)、临床分期(Stage)、病理分级(Grade)做单因素和多因素Cox回归分析,确定预后的独立影响因素。列线图、ROC曲线预测患者1、3、5年的生存率,并对临床特征进行分层做生存曲线分析,以进一步分析风险评分的临床潜力。

1.6 预后模型基因的表达验证 使用HPA数据库(https://www.proteinatlas.org/)验证模型基因在肿瘤组织和正常组织中的蛋白表达是否与差异表达分析的结果一致。

1.7 免疫细胞浸润分析 使用CIBERSORT计算免疫细胞在每个KIRC患者中的浸润情况。每种免疫细胞类型的总比例等于1。Mann-WhitneyU检验比较高低风险组患者免疫细胞的浸润差异情况。

1.8 药物敏感性分析 从CellMiner数据库(http://discover.nci.nih.gov/cellminer/)中下载NCI-60 细胞系的RNA表达数据和药物分析数据,选取经过临床试验和FDA批准的792个药物结果。提取模型基因的表达数据,计算药物敏感性(IC50),并进行Pearson相关性分析,最终得到每个模型基因的表达与药物敏感性之间的相关性。以P<0.01对分析结果进行筛选。

1.9 统计学分析 使用R软件4.1.2对数据进行统计分析和结果可视化。P<0.05为差异具有统计学意义。

2 结 果

2.1 ARGs在癌症中的表达差异情况 ssGSEA结果显示, ARGs在KIRC中表达最高,肾嫌色细胞癌(KICH)中表达最低(图1A)。Mann-WhitneyU检验结果发现,ARGs在14种癌症的肿瘤样本与正常样本之间表达存在差异(图1B),其中ARGs在多形性成胶质细胞瘤(GBM)、KIRC、肾乳头状细胞癌(KIRP)和甲状腺癌(THCA)的肿瘤样本中表达水平高于正常样本,在膀胱尿路上皮癌(BLCA)、乳腺浸润癌(BRCA)、结肠癌(COAD)、头颈鳞状细胞癌(HNSC)、肝细胞肝癌(LIHC)、肺腺癌(LUAD)、肺鳞癌(LUSC)、前列腺癌(PRAD)、直肠腺癌(READ)、胃癌(STAD)的肿瘤样本中表达水平低于正常样本(图1C),而在宫颈鳞癌和腺癌(CESC)、食管癌(ESCA)、KICH、胰腺癌(PAAD)、嗜铬细胞瘤和副神经节瘤(PCPG)、肉瘤(SARC)、皮肤黑色素瘤(SKCM)、胸腺癌(THYM)、子宫内膜癌(UCEC)中表达无差异性。对有差异表达的癌种进行单因素Cox回归分析,发现自噬评分与KIRC、LUSC的OS相关(图1D)。

注:A.ARGs在各癌症中的自噬评分;B.正常样本与肿瘤样本之间自噬评分存在差异的癌症;C.正常样本与肿瘤样本自噬评分的中位数;D.单因素Cox回归分析森林图。*P<0.05,**P<0.01,***P<0.001,****P<0.0001。BLCA.膀胱尿路上皮癌;BRCA.乳腺浸润癌;COAD.结肠癌;KIRP.肾乳头状细胞癌;HNSC.头颈鳞状细胞癌;KIRC.肾透明细胞癌;LIHC.肝细胞肝癌;READ.直肠腺癌;LUSC.肺鳞癌;GBM.多形性成胶质细胞瘤;STAD.胃癌;THCA.甲状腺癌;PRAD.前列腺癌;LUAD.肺腺癌。

2.2 自噬评分对KIRC患者预后的影响 生存曲线分析结果显示,KIRC低评分患者的预后比高评分患者差(P=0.039)(图2A),LUSC高评分患者的预后比低评分患者差(P=0.034)(图2D),这与已有文献结果一致[10-11]。单因素、多因素Cox回归分析结果发现,自噬评分可以作为KIRC患者预后的独立影响因素(图2B、2C),但不能作为LUSC患者的独立预后因素(图2E、2F)。

注:A.KIRC高低自噬评分组的生存曲线;B、C.影响KIRC预后的单、多因素Cox回归分析森林图;D.LUSC高低自噬评分组的生存曲线;E、F.影响LUSC预后的单、多因素Cox回归分析森林图。

2.3 确定KIRC差异表达的自噬相关基因并构建预后模型 KIRC的正常样本与肿瘤样本差异表达分析结果得到9 632个差异表达基因,将其与上述ARGs取交集,得到35个DEARGs(图3A)。其中有6个基因在肿瘤组织中下调,29个基因在肿瘤组织中上调(图3B)。对DEARGs进行逐步回归分析、LASSO回归分析,最终确定4个模型基因:BAG1、BIRC5、BNIP3、EIF4EBP1。基于多因素Cox回归分析结果(表1)构建预后模型,发现BIRC5和EIF4EBP1是高风险基因(HR>1),而BAG1和BNIP3是低风险基因(HR<1)。风险评分=(-0.037)×BAG1表达量+0.021×BIRC5表达量+(-0.003)×BNIP3表达量+0.003×EIF4EBP1表达量。免疫组化结果表明,与正常组织相比,KIRC组织中的BAG1(HPA018121)、BIRC5(HPA002830)、BNIP3(HPA003015)和EIF4EBP1(CAB005032)的蛋白水平上调(图3C),除BAG1外,其他3个基因的蛋白表达与DEARGs的表达结果一致。

注:A.差异表达基因与自噬相关基因维恩图;B.35个差异表达自噬相关基因热图;C.正常组织和肿瘤组织4个模型基因的蛋白水平比较。

表1 多因素Cox回归分析预后模型基因对预后的影响

2.4 模型对KIRC患者的临床预后价值 生存曲线结果显示,高风险组患者的预后比低风险组患者差(训练集P<0.001,AUC=0.755,见图4A、4B;测试集P<0.001,AUC=0.695,见图4C、4D)。风险因子关联图结果显示,风险评分值越高,患者死亡数越多,训练集(图4E)与测试集(图4F)结果基本一致。其中热图显示,BIRC5和EIF4EBP1在高风险组中上调,是危险因素,而BAG1和BNIP3在高风险组中下调,是有利因素,与多因素Cox回归结果一致。

注:A、B.训练集生存曲线和ROC曲线;C、D.测试集生存曲线和ROC曲线;E、F.训练集和测试集风险因子关联图,分布曲线(上)、散点图(中)和热图(下)。

单因素和多因素Cox回归分析结果发现,构建的模型风险评分、M分期和N分期均可作为KIRC患者预后的独立影响因素(P<0.05)(图5A、5B)。为了更好地预测KIRC患者1、3、5年的生存率,构建了列线图(图5C)并绘制ROC曲线,AUC分别为0.743、0.740、0.699(图5D)。为了进一步研究风险评分的临床预后潜力,基于临床特征(TNM分期、临床分期、肿瘤分级)进行生存曲线分析,发现风险评分似乎更适合用于预测临床分期、肿瘤分级、T期、M0期和N0期的总生存时间,除M1、N1期外,高风险评分患者的预后均比低风险评分患者差(P<0.05)(图5E)。

注:A、B.单因素、多因素Cox回归分析森林图;C、D.列线图和ROC曲线;E.生存曲线。

2.5 高低风险组免疫细胞浸润情况比较 基于自噬与免疫细胞浸润的相关性,计算了每个样本中肿瘤免疫微环境22种免疫细胞的浸润情况(图6A),评估了不同风险组中免疫微环境的组成。与低风险组比较,高风险组中CD8+T细胞、辅助性滤泡T细胞、调节性T细胞和M0巨噬细胞等比例显著升高,而静息杀伤细胞、单核细胞和静息肥大细胞等比例显著降低(图6B、C)。

注:A.KIRC患者22种免疫细胞浸润比例;B、C.训练集、测试集高低风险组免疫细胞浸润差异。*P<0.05,**P<0.01,***P<0.001,****P<0.0001。

2.6 模型基因的表达与多种药物敏感性的关系 表2显示每个基因与药物敏感性(IC50)最相关的前4个结果。在60个癌细胞系中BAG1与白屈菜红碱、AM-5992、7-羟基星形孢菌素呈正相关,与INK-128(mTOR抑制剂)呈负相关;BIRC5与SNS-314(极光激酶抑制剂)、5-氟脱氧尿苷10、博来霉素呈正相关,与司美替尼呈反比;BNIP3与P-529(mTOR抑制剂)、特拉替尼、IDH-C227(IDH1R132H抑制剂)呈正相关,与AM-5992呈负相关;EIF4EBP1与SNS-314、氟尿苷、奎扎替尼、维甲酰酚胺呈正相关。

表2 基因表达与药物敏感性的相关性

在8个肾癌细胞系中BAG1与丁硫氨酸—亚砜亚胺、PF-2771(CENP-E抑制剂)、三氧化二砷呈负相关,与帕博西尼呈正相关;BIRC5与恩曲替尼、Gandotinib(JAK2抑制剂)、AMD-070(CXCR4选择性抑制剂)、替沃扎尼呈负相关;BNIP3与Saridegib(Smo 特异性抑制剂)、阿扎胞苷、卡铂呈正相关,与CH-5132799(PI3K抑制剂)呈负相关;EIF4EBP1与AZD-5363、GSK-2141795、Afuresertib、Ipatasertib(AKT抑制剂)呈负相关,见表2。

3 讨 论

目前,尽管有大量研究表明自噬参与了KIRC的恶性进展,但以往的研究中仅探索了单个自噬相关基因与KIRC之间的相关性[12-16],缺乏对自噬相关基因在各癌症中的表达分析,且与多种癌症预后相关的自噬相关基因与KIRC患者预后之间的相关性尚未明确。虽然目前有文献报道了自噬相关基因对KIRC患者预后的影响[17-19],但涉及的基因多,临床经济实用性欠缺,且未曾发现有相关研究解释上述所有4个基因在KIRC中的预后价值。在本研究中,笔者发现自噬相关基因在不同癌症中的表达存在差异,其中在KIRC患者中表达最高。因此构建了由BAG1、BIRC5、BNIP3和EIF4EBP1组成的风险预后模型。发现在KIRC患者中BAG1和BNIP3是低风险基因,其表达越高,患者风险值越低;BIRC5和EIF4EBP1属于高风险基因,其表达越高,患者风险值越高。风险评分、M分期、N分期可以作为KIRC患者的独立预后因素。

肿瘤微环境中包含各种细胞,其中浸润性免疫细胞占很大比例且特定细胞对患者的生存有重大影响[20]。例如,肿瘤相关的巨噬细胞可通过多种方式帮助肿瘤细胞发生免疫逃逸、肿瘤血管生成和转移[21-23]。上述免疫浸润分析结果表明,与低风险组比较,高风险组中CD8+T细胞、辅助性滤泡T细胞、调节性T细胞和M0巨噬细胞浸润比例升高。高风险组的评分与患者预后呈反比,提示这些免疫细胞浸润增加可能与KIRC患者的不良预后有关。而静息杀伤细胞、单核细胞和静息肥大细胞在高风险组中浸润比例降低,说明这些免疫细胞的浸润增加可能与KIRC患者的良好预后有关。

近年来,自噬相关基因在多药耐药中的作用受到越来越多的关注。以往研究报道,在乳腺癌、结肠癌和卵巢癌细胞系中,星形孢菌素可以通过介导自噬蛋白p62的表达来逆转顺铂的耐药性[24],说明可以通过改变自噬水平或抑制自噬相关基因的表达,在一定程度上提高药效,从而为临床用药提供理论依据[25]。在本研究中构建的4个模型基因分别与多种癌细胞系的药物敏感性相关,其中在肾癌细胞系中,与BAG1表达呈正相关的药物2个,负相关的11个;与BIRC5表达呈正相关的3个,负相关的51个;与BNIP3表达呈正相关的3个,负相关的1个;与EIF4EBP1表达呈正相关的4个,负相关的21个。且相关系数可达到0.9以上。这可能为肾癌的预后提供潜在的药物治疗靶点,但仍需要更多的试验加以验证,以评估该模型基因对临床试验中新药的影响,为新药筛选提供选择。

综上所述,在本研究中,笔者从分析预后自噬相关基因在泛癌中的表达差异,聚焦到对KIRC患者的生存预后分析,确定了一种由4个差异自噬基因(BAG1、BIRC5、BNIP3和EIF4EBP1)组成的新预后标志物模型,同时证明了该模型与KIRC 患者的OS独立相关,并发现其在肿瘤免疫微环境、药物敏感性等方面也具有重要意义。但潜在的分子机制还需要进一步的实验研究。

利益冲突：所有作者声明无利益冲突

作者贡献声明

伍艳婷:提出研究方案并设计课题,实施研究过程,论文撰写;舒文莹:确定分析方法,论文修改;王艳平:文献调研,资料搜集和数据整理;周毅:课题设计审核;简晓顺:提出研究思路,论文审核