铁死亡相关的mRNAs特征预测胃癌预后的分析

谭琦焕,谢王亮,骆俞婷,蒋宁芳,马阿火

谭琦焕,谢王亮,骆俞婷,绍兴文理学院医学院 浙江省绍兴市 312000

蒋宁芳,绍兴市人民医院胃肠外科 浙江省绍兴市 312000

马阿火,绍兴文理学院附属第一医院消化内科 浙江省绍兴市 312000

0 引言

胃癌是最常见的胃肠道恶性肿瘤,占癌症发病率的第五位,但是其死亡率却居第三位.由于不良的生活习惯、幽门螺杆菌感染等,2020年全球新增病例1089103例,死亡768793例[1].虽然分子诊断技术,手术和辅助治疗技术等方面在不断进步,但胃癌患者的总生存率仍然很差,尤其是胃癌晚期患者,5年生存率低于20%[2,3].因此,更好地了解胃癌发生的分子机制,识别新的生物标志物,对提高胃癌的早期诊断和治疗水平至关重要.铁死亡是一种以脂质过氧化与活性氧的积累为特征的铁依赖程序性细胞死亡[4].近年来,铁死亡逐渐成为一种极具前景的诱导癌细胞死亡的治疗方法[5].有研究发现铁死亡与肿瘤的预后相关,如在肝细胞癌中已经建立一个铁死亡相关预后基因模型[6].在胃癌中,发现半胱氨酸双加氧酶1(cysteine dioxygenase 1,COD1)[7]和花生四烯酸15-脂氧合酶(arachidonic acid 15-lipoxygenase,ALOX15)[8]可促进人胃癌细胞铁死亡,而硬脂酰辅酶A去饱和酶1(stearoyl-coa desaturease1,SCD1)[9]、脂肪细胞分化相关蛋白2(adipocyte differentiation associated protein 2,PLIN2)[10]和生长分化因子15(growth and differentiation factor 15,GDF15)[11]则相反.但这些铁死亡相关基因是否与胃癌患者预后相关,仍值得我们探索.因此,在本研究中,我们基于TCGA数据库构建了铁死亡相关基因差异表达的预后模型,同时进行肿瘤免疫微环境分析,以探讨潜在的机制.

1 材料和方法

1.1 材料 获取基因表达数据及临床信息:从TCGA数据库下载胃癌相关的转录组数据及相关临床数据.我们获得32例癌旁组织及375例癌组织的mRNA数据.同时,我们获取了375例胃癌患者的临床信息,包括生存时间、生存状态、性别、年龄、肿瘤分级、病理分期.从FerrDb数据库获取铁死亡相关基因列表,最终得到245个铁死亡相关mRNA的表达数据[12].从GEO数据库下载GSE62254与GSE34942数据集,分别有300例胃癌组织的mRNA数据和56例胃癌组织的mRNA数据及相关临床生存数据.

1.2 方法

1.2.1 数据的预处理:首先对胃癌患者的数据进行筛选:(1)排除胃癌患者中缺乏生存时间或生存时间为0的数据;(2)排除胃癌患者中缺乏完整临床信息的病例;(3)排除病理诊断为非胃腺癌的病例;(4)排除同时有其他恶性肿瘤的病例.经过筛选我们在TCGA数据库中得到347例具有转录组数据及生存预后资料的胃癌样本,进行风险评分模型的构建.进一步选择TCGA数据库中306例具有完整临床相关资料的胃癌样本,用于相关预后临床风险因素的分析.

1.2.2 预后模型的建立:首先使用“limma”包进行差异分析,得到胃癌中铁死亡相关差异表达基因(differentially expressed genes,DEGs),筛选条件为:错误发现率(false discovery rate,FDR)＜0.05和|log2差异倍数(fold change,FC)|≥1.对347例有生存资料的胃癌患者进行单因素cox回归分析,评估铁死亡相关mRNA的预后价值(筛选条件为P＜0.05).我们把DEGs与具有预后价值的铁死亡相关mRNA取交集,得到具有预后价值的铁死亡相关DEGs.随后使用Lasso回归对有预后价值的铁死亡相关DEGs进行进一步筛选并构建预后模型.根据每个mRNA的表达水平及其相应的回归系数计算风险评分.风险评分=(mRNA1系数×mRNA1的表达量)+(mRNA2系数×mRNA2的表达量)+…+(mRNAn系数×mRNAn的表达量).根据风险评分的中位值,将胃癌患者分为高危组与低危组.应用Kaplan-Meier生存曲线来比较其之间的生存率(overall survival,OS).使用“survivalROC”包进行时间依赖的受试者工作特征曲线(receiver operating characteristic curve,ROC)分析来评估模型的可信度.最后我们使用单因素和多因素cox回归分析将风险评分与临床信息结合起来评估与患者预后的关系,识别独立预后因素.

1.2.3 免疫分析与基因表达:使用“gsva”包中的单样本基因集富集分析(single sample gene set enrichment analysis,ssGSEA)计算16种免疫细胞的浸润评分和13种免疫相关通路的活性[13].此外,还从先前的文献中检索了潜在的免疫检查点.

1.2.4 外部数据集的验证:使用“sva”包合并GSE62254与GSE34942数据集并进行去批次化,得到一个具有356例具有转录组数据及生存预后资料的胃癌样本.以同样的回归系数(coef值)进行风险评分模型的构建来进行验证.

统计学处理使用R语言(4.0.2版)、Bioconductor软件包、相关R语言软件包进行数据分析,模型构建等.正态和非正态分布变量分别采用不成对t检验和Wilcoxon检验进行分析.P＜0.05有统计学意义.

2 结果

2.1 胃癌中铁死亡相关基因的表达 首先我们把从TCGA数据库下载的转录组数据与FerrDb数据库下载的铁死亡相关基因取交集,得到245个铁死亡相关mRNA的表达数据.然后,进行胃癌与癌旁组织的差异分析,得到50个铁死亡相关DEGs(表1).接下来把347例TCGA数据集中245个铁死亡相关的mRNA进行单因素cox回归分析,得到25个与预后相关的铁死亡相关mRNA(表2).把得到的50个铁死亡相关DEGs与25个有预后价值的铁死亡相关mRNA取交集,得到7个具有预后价值的铁死亡相关DEGs(图1).

图1 TCGA队列中铁死亡相关基因的筛选. A:识别与生存相关的胃癌和癌旁组织之间差异表达基因的维恩图;B:森林图显示了具有预后价值的铁死亡相关DEGs单变量cox回归分析的结果;C:具有预后价值的铁死亡相关DEGs的相关网络,相关系数用不同的颜色表示.DEGs:差异表达基因.

表1 胃癌和癌旁组织差异表达的铁死亡相关基因

表2 与预后相关的铁死亡相关mRNA

2.2 构建胃癌铁死亡相关mRNA预后模型 我们把7个具有预后价值的铁死亡相关DEGs进行Lasso回归分析,通过交叉验证法进行迭代分析,构建出一个含有4个mRNA(DUSP1、MYB、CAV1、NOX4)的预后模型(图2).风险评分的公式:风险评分=(0.108×DUSP1的表达量)+(-0.089×MYB的表达量)+(0.036×CAV1的表达量)+(0.368×NOX4的表达量).使用上述方法我们计算了队列中每位患者的风险评分.根据风险评分(risk score,RS)的中位值(RS=0.8331),把患者分为高危组(RS≥0.8331)与低危组(RS＜0.8331).我们发现随着风险评分增加,高危组的患者死亡人数明显多于低危组(图3A和B).随后进行的主成分分析法(principal component analysis,PCA)和t-分布领域嵌入算法(t-distributed stochastic neighbor embedding,t-SNE)表明,不同风险组的患者分布在两个方向上(图3C和D).

图2 Lasso回归构建预后相关模型.A:Lasso筛选变量动态过程图;B:交叉验证参数λ的选择过程图.

图3 4个胃癌铁死亡相关mRNA预后模型.A:风险评分分布曲线图;B:胃癌患者生存状态的分布散点图;C:胃癌患者的PCA图;D:胃癌患者的t-SNE分析图.PCA:主成分分析法;t-SNE:t-分布领域嵌入算法.

2.3 生存结果分析 我们对预后模型进行Kaplan-Meier生存分析,发现高危组患者的OS明显低于低危组(P＜0.001)(图4A).接下来,利用单因素与多因素cox回归分析,确定风险评分是否是OS的独立预后预测因子.我们发现在单因素Cox回归分析中,风险评分与OS显著相关[风险比(hazard ratio,HR)=3.267,P＜0.001](表3).校正其他混杂因素后,在多变量cox回归分析中,发现风险评分仍然是OS的独立预测因子(HR=3.211,P＜0.001)(表4).接下来,为了验证预后模型的灵敏度与特异度,我们绘制了预后模型中临床特征的ROC曲线,胃癌患者的年龄、性别、肿瘤分级、病理分期的ROC曲线下的面积(area under curve,AUC)分别为0.589、0.520、0.563、0.596,而风险评分的AUC为0.644,说明在预测胃癌预后方面,预后模型的风险评分明显优于传统临床病理学特征(图4B).最后,我们分析了预后模型是否可用于判断不同临床特征下的生存情况,对年龄(≤65和＞65)、性别(男性和女性)、肿瘤分级(G1-2和G3)、病理分期(Stage Ⅰ-Ⅱ和Stage Ⅲ-Ⅳ)进行分层分析.我们发现,除了在女性组中,高危组与低危组之间的生存率无明显差别(P=0.329),其余各亚组中,低危组的OS明显优于高危组(P＜0.05),图5.

图4 4个胃癌铁死亡相关mRNA预后模型特征.A:不同风险评分胃癌患者的生存曲线;B:不同临床指标以及预后风险值的时间依赖性ROC曲线以及AUC值.ROC:受试者工作特征曲线;AUC:曲线下的面积.

图5 不同临床特征中胃癌患者的生存曲线. A:年轻患者组的Kaplan-Meier生存曲线;B:老年患者组的Kaplan-Meier生存曲线;C:女性组的Kaplan-Meier生存曲线;D:男性组的Kaplan-Meier生存曲线;E:低级别患者组的Kaplan-Meier生存曲线;F:高级别患者组的Kaplan-Meier生存曲线;G:早期患者组的Kaplan-Meier生存曲线;H:晚期患者组的Kaplan-Meier生存曲线.

表3 单因素cox回归结果

表4 多因素cox回归结果

2.4 免疫与基因表达 为了进一步探索风险评分与免疫状态之间的相关性,我们使用ssGSEA量化不同免疫细胞亚群、相关功能及途径的富集评分.我们发现在抗原呈递过程中像成熟树突状细胞、B细胞、CD8+T细胞、树突状细胞、未成熟树突状细胞、巨噬细胞、肥大细胞、中性粒细胞、NK细胞、浆细胞样树突状细胞、辅助性T细胞、滤泡辅助性T细胞、Th1细胞、肿瘤浸润淋巴细胞、调节性T细胞、APC共抑制、APC共刺激、细胞因子-细胞因子受体、检查点、人类白细胞抗原等的评分,在预后模型的高危组与低危组之间存在显着差异(P＜0.05)(图6A和B).随着检查点抑制剂相关的免疫疗法越来越受到大家的关注,我们进一步探索了两组之间免疫检查点表达的情况.我们发现两组之间在肿瘤坏死因子(tumor necrosis factor,TNF)超家族(如CD40LG、CD70、TNFSF4、CD27、CD40等)、B7-CD28家族(如CD274、PDCD1LG2、PDCD1、CD276、CTLA4、ICOS等)和其他免疫检查点(如LAG3、IDO2、HAVCR2等)存在差异,尤其是在高危组表达上调(图7).

图6 免疫细胞亚群、相关功能的关系.A:ssGSEA相关免疫细胞的表达;B:免疫相关功能的表达.aP＜0.05,bP＜0.01,cP＜0.001,与低危组比较.aDCs:抗原呈递过程中像成熟树突状细胞;B cells:B细胞;CD8+T cells:CD8+T细胞;DCs:树突状细胞;iDCs:未成熟树突状细胞;Macrophages:巨噬细胞;Mast cells:肥大细胞;Neutrophils:中性粒细胞;NK cells:NK细胞;pDCs:浆细胞样树突状细胞;T helper cells:辅助性T细胞;Tfh:滤泡辅助性T细胞;Th1 cells:Th1 细胞;Th2 cells:Th2 细胞;TIL:肿瘤浸润淋巴细胞;Treg:调节性T细胞;APC co inhibition:抗原提呈细胞共抑制;APC co stimulation:抗原提呈细胞共刺激;CCR:细胞因子-细胞因子受体;Check point:检查点;Cytolytic activity:细胞溶解活性;HLA:人类白细胞抗原;Inflamation promoting:炎症促进;MHC class Ⅰ:主要组织相容性复合体Ⅰ;Parainflamation:副炎症;T cell co inhibition:T细胞共抑制;T cell co stimulation:T细胞共刺激;Type I IFN Reponse:Ⅰ型干扰素反应;TypeⅡIFN Reponse:Ⅱ型干扰素反应.

图7 高危组与低危组免疫检查点的表达. aP＜0.05,bP＜0.01,cP＜0.001,与低危组比较.HAVCR2:甲型肝炎病毒细胞受体2;CD40LG:CD40配体基因;BTNL2:嗜乳脂蛋白样2;TNFRSF9:肿瘤坏死因子受体超家族成员9;TNFRSF9:肿瘤坏死因子受体超家族成员14;ADORA2A:腺苷A2a受体;TNFSF4:肿瘤坏死因子超家族第4成员;LAIR1:人白细胞相关免疫球蛋白样受体1;TNFSF14:肿瘤坏死因子超家族成员14;CD200R1:CD200受体1;PDCD1LG2:程序性死亡分子1配体2;KIR3DL1:杀伤细胞免疫球蛋白样受体基因单体分型模式3;TNFRSF4:肿瘤坏死因子受体超家族成员4;LGALS9:可溶性半乳糖苷结合凝集素9;TNFRSF8:肿瘤坏死因子受体超家族成员8;PDCD1:程序性细胞死亡蛋白1;IDO2:吲哚胺2,3-双加氧酶2;CTLA4:细胞毒T淋巴细胞相关抗原4;TNFRSF25:肿瘤坏死因子受体超家族成员25;TIGIT:含T细胞免疫球蛋白和ITIM 结构域蛋白;BTLA:B、T淋巴细胞衰减因子;NRP1:神经纤毛蛋白1;ICOS:可诱导共刺激分子;TMIGD2:受体穿膜和免疫球蛋白结构域2;HHLA2:人内源性逆转录病毒-H长末端重复关联蛋白2;TNFSF18:肿瘤坏死因子超家族成员18;LAG3:淋巴细胞活化基因3.

2.5 外部数据集的验证 我们使用相同的coef值在外部GEO数据集进行验证,同样根据风险评分将患者分为高危组与低危组.我们发现高危组患者的OS明显比低危组患者差(P=0.003),这与我们在TCGA数据集中得出的结果相同(图8A).我们进一步构建了预测患者1、2、3年生存的时间依赖性ROC曲线的AUC分别为0.620、0.590和0.578(图8B).

图8 GEO验证集的相关预后特征. A:验证集中不同风险评分胃癌患者的生存曲线;B:验证集中预后模型的1、2、3年生存的时间依赖性ROC曲线.ROC:受试者工作特征曲线;AUC:曲线下的面积.

3 讨论

在本次研究中,我们基于TCGA数据库构建了与铁死亡相关基因的预后模型.然后我们进一步探讨肿瘤微环境中免疫细胞浸润和免疫检查点抑制剂在胃癌的作用.

首先,我们综合分析了铁死亡相关基因在胃癌中的表达并构建了可以预测患者生存状态的4-mRNA预后模型.通过预后模型,我们可以明显的发现高危组患者的生存率比低危组差.同时,基于单因素与多因素cox回归结果,4-mRNA预后模型(DUSP1、MYB、CAV1、NOX4)可以独立预测胃癌患者的预后.双特异性磷酸酶1(dual specificity phosphatase-1,DUSP1),也称丝裂原活化蛋白激酶磷酸酶-1 (mitogen-activated protein kinase phosphatase-1,MKP-1K),是双特异性磷酸酶的原型成员,在调节肿瘤细胞中丝裂原活化蛋白激酶(mitogenactivated protein kinase,MAPK)家族的功能方面发挥着作用,调节MAPK信号通路参与肿瘤细胞增殖、分化、转化、周期阻滞和细胞凋亡[14,15].在人胰腺癌细胞中,DUSP1的表达与活性通过erastin和RSL3以细胞外调节蛋白激酶(extracellular regulated protein kinases,ERK)依赖的方式增加,随后形成抑制自噬依赖性铁死亡反馈机制[16].DUSP1可能通过阻断unc-51样自噬激活蛋白激酶1(unc-51 like autophagy activating kinase 1,ULK1)或BECN1的磷酸化来限制自噬小体的形成.因此抑制DUSP1可能为铁死亡治疗诱导提供新的可能性.而且,DUSP1与肿瘤细胞的耐药性密切相关,在人卵巢癌细胞中,DUSP1过表达通过p38 MAPK介导的p-糖蛋白导致对紫杉醇的耐药[17].在人胃癌细胞中,DUSP1高表达会增加对阿帕替尼的耐药性[18].如果DUSP1表达降低,肿瘤细胞的化疗耐药性降低,使肿瘤细胞具有更高的敏感性.MYB(v-myb avian myeloblastosis viral oncogene homolog)是一种原癌基因蛋白,与活性氧(reactive oxygen species,ROS)的改变相关[19].在胃癌细胞中,抑制c-MYB可以抑制COD1的转录,从而增加谷胱甘肽(glutathione,GSH)的生成,阻止ROS的生成,最后抑制erastin诱导的铁死亡[7].且有研究表明c-Myb可以通调节髓性锌指1(一种调节铁转运蛋白表达的转录因子)来调节铁相关细胞活性和肿瘤细胞生长[20].窖蛋白-1(Caveolin-1,CAV1)是一种细胞表面的穴样内陷中的主要蛋白,在肿瘤发生和转移中起着重要作用,具有抗凋亡潜力,增强肿瘤细胞迁移和侵袭能力,增加化疗耐药性[21].在人胃癌细胞中,CAV1通过激活WNT信号通路和间充质上皮转化因子(mesenchymal to epithelial transition factor,MET)-人表皮生长因子受体-2(human epidermal growth factor receptor 2,HER-2)串扰来增加胃癌细胞对顺铂的耐药性[22].而且CAV1的表达与进行根治性胃切除术后患者的不良预后明显相关,可能靶向CAV1是治疗胃癌的潜在选择[23,24].还原型辅酶Ⅱ氧化酶4(nicotinamide adenine dinucleotide phosphate oxidase 4,NOX4)是NOX家族的成员,是一种致癌因子,参与免疫调节.NOX4的高表达会产生ROS并激活胶质瘤相关癌基因1(glioma-associated oncogene 1,GLI1)信号,促进胃癌细胞的增殖与转移[25].而且,高等人[26]发现NOX4在胃癌组织中表达升高,通过Janus激酶/信号转导器和转录信号通路激活剂参与胃癌细胞粘附和侵袭.与我们的结果一致,我们发现与正常组织相比,胃癌组织中的NOX4 mRNA显著上调.Lee等人[27]提出NOX4可以作为预测胃癌患者手术后复发的预后基因之一.同时,有研究表明抑制NOX4会增强CD8+T细胞在肿瘤中的浸润,提高免疫治疗的疗效[28].尽管这些基因与胃癌的发病机制密切相关,但将它们结合起来作为铁死亡相关基因在胃癌患者预后的相关研究较少,因此我们此次的发现为肿瘤的诊断与治疗提供一个建议.

同时,在我们的研究中,高危组的免疫细胞浸润水平较高,如巨噬细胞、肥大细胞和中性粒细胞.有研究表明[29],肿瘤相关巨噬细胞与胃癌患者不良预后相关,会促进肿瘤细胞增殖[30]、侵袭[31]、转移、免疫逃逸[32]等.在肿瘤相关巨噬细胞通过TGFβ2/NF-κB/Kindlin-2 axis轴参与胃癌的进展[33].巨噬细胞来源的miR-21可以通过下调第10号染色体同源丢失性磷酸酶张力蛋白基因(gene of phosphate and tension homology deleted on chromosome ten,PTEN)抑制细胞凋亡及胃癌细胞顺铂耐药,增强PI3K/AKT信号通路的激活[34].肥大细胞水平升高会抑制T细胞免疫及促进肿瘤的生长,并且可以独立预测胃癌患者的生存率[35].肿瘤浸润性中性粒细胞也可以抑制正常T细胞免疫,与胃癌患者较差的生存相关[36].而且,高危组人群有着高水平的Ⅱ型干扰素反应,会导致抗肿瘤免疫的减弱.我们推测,胃癌患者的免疫细胞浸润水平较高导致的肿瘤免疫增强和抗肿瘤免疫的减弱导致预后不良.铁死亡的诱导与免疫检查点抑制剂相结合显示出协同增强的抗肿瘤活性,即使在免疫检查点抑制剂抗性肿瘤中[37].如B和T淋巴细胞衰减蛋白(B and T lymphocyte attenuator,BTLA)是TNF超家族的成员,其表达与癌症侵袭性有关[38].甲型肝炎病毒细胞受体2(hepatitis A virus cellular receptor 2,HAVCR2)主要位于自然杀伤细胞和巨噬细胞,抑制抗肿瘤免疫的激活[39].本次的研究为未来进一步研究提供了方向.

4 结论

铁死亡是一种新的细胞死亡形式,有望为肿瘤治疗提供一种新的途径.这项研究有一些局限性.我们的模型是通过公共数据库的回顾性数据构建和验证的,需要进一步的前瞻性研究来验证其临床价值.其次,应该进行铁死亡实验性研究来进一步揭示潜在机制.我们筛选了用于胃癌预后的铁死亡相关生物标志物,为胃癌的治疗提供了一定的理论基础.

文章亮点

实验背景

胃癌是最常见的消化道恶性肿瘤,我国胃癌发病率呈逐年升高趋势.虽然随着手术、化疗、靶向及免疫治疗等方式取得了进展,胃癌患者的生存期有所延长,但大多数患者就诊时已处于晚期,预后差.因此,寻找新的胃癌的标志物和潜在治疗靶点,研究其在胃癌中的分子机制具有重要意义.

实验动机

铁死亡是一种新型的调节细胞死亡形式,探讨铁死亡相关基因在胃癌中的表达及其与患者预后的关系.

实验目标

采用生物信息学的方法来构建铁死亡相关基因的预后模型,评估其作为胃癌患者预后分子标志物的可行性.

实验方法

通过分析TCGA数据库中的铁死亡相关基因表达数据与胃癌生存预后之间的关系,构建预后模型,应用Kaplan-Meier生存曲线分析不同风险组之间对胃癌患者的影响,同时进行全面科学的验证和评价.

实验结果

构建了一个含有4个mRNA(DUSP1、MYB、CAV1、NOX4)的预后模型,生存分析显示高危组与胃癌预后不良相关,且风险评分是生存的独立预后指标,同时找到了潜在的治疗靶点.基于铁死亡相关基因构建的预后模型具有良好的预测价值.

实验结论

铁死亡相关mRNAs预后模型与胃癌的预后密切相关,对胃癌的预后具有一定的预测意义并提供了潜在的治疗靶点.

展望前景

铁死亡相关mRNAs预后模型可以作为一种潜在的生物标记物来识别高危胃癌患者,有希望成为胃癌预后监测指标和新的治疗靶点.