原发性开角型青光眼铁死亡相关基因表达谱的生物信息学分析

周日龙,廖 良

目的:通过生物信息学手段分析原发性开角型青光眼(POAG)进展过程中调控铁死亡相关的关键基因,旨在进一步揭示铁死亡在POAG中的生物学机制。

方法:从GEO数据库中获取小梁网来源的GSE27276数据集,其中包括19个POAG小梁网组织样本和17个正常小梁网组织样本;下载FerrDb数据库整理的铁死亡相关基因,将GSE27276数据集与铁死亡基因集进行映射,筛选POAG中铁死亡相关的预后差异表达基因(DE-FRGs)并进行相关性分析,进一步了解DE-FRGs的GO和KEGG通路富集。应用LASSO回归模型与SVM-RFE模型两种机器学习的算法筛选铁死亡相关POAG的关键基因,将两种模型的筛选结果取交集,获得最佳特征基因,使用受试者特征曲线(ROC)评估临床诊断能力;对最佳特征基因进行单基因基因组富集分析(GSEA)和变异分析(GSVA);借助视乳头来源的GSE2378与GSE9944数据集验证最佳特征基因的表达水平。

结果:与正常小梁网组织相比,POAG的小梁网组织有396个铁死亡基因存在差异表达,其中39个为上调基因,64个为下调基因,Spearman相关性分析显示上调基因和下调基因均有一定的相关性。GO功能和KEGG通路富集分析显示,差异基因主要富集在氧化应激反应和铁死亡通路;通过LASSO和SVM-RFE算法将18个DE-FRGs确定为关键基因,具有更高的诊断价值。GSEA和GSVA富集分析显示GDF15、MFN2和OTUB1基因与谷胱甘肽代谢通路密切相关,其中MFN2和OTUB1分别在高表达组和低表达组中激活谷胱甘肽代谢通路。GSE2378与GSE9944数据集交叉验证明确视乳头标本中CREB1的表达水平相对于正常视乳头样本显著升高,这与GSE27276数据集小梁网样本表达一致。

结论:基于生物信息学分析挖掘得到396个POAG的DE-FRGs,通过构建机器筛选模型和外部数据集交叉验证,筛选出CREB1有望成为潜在诊断生物标志物的最佳特征基因,为进一步深入阐明POAG铁死亡相关的分子机制和诊断提供靶点。同时筛选的基因还需要进一步体内、外实验验证,揭示铁死亡在POAG中的生物学机制。

•KEYWORDS：primary open angle glaucoma; ferroptosis; bioinformatics

0 引言

原发性开角型青光眼(primary open angle glaucoma, POAG)是一种慢性进行性视神经病变,其特征为视网膜神经节细胞(retinal ganglion cell, RGC)进行性变性和视野缺失,好发于40～80岁人群[1-2]。POAG虽不危及生命,但具有早期隐匿性强、中后期致盲率高等特点[3],已成为导致中老年致盲最常见的原因之一。如果未得到及时有效的治疗,POAG引起的视力丧失不仅影响生活质量,而且造成重大的社会经济负担[4]。目前,其发病机制尚不完全清楚,而遗传学被证明在POAG的发病机制中起关键作用[3-5]。铁死亡是2012年新发现的一种铁依赖性程序性细胞死亡模式,其发生机制不同于细胞焦亡、凋亡、坏死和自噬[6]。铁死亡的特征是线粒体萎缩和线粒体膜密度增加,铁和脂质活性氧的积累[7]。铁死亡的发生机制是通过脂质自由基的形成和谷胱甘肽(L-glutathione, GSH)的消耗或脂质过氧化物酶4(glutathione peroxidase 4,GPX4)的失活来催化的[8]。铁循环在铁死亡的发展中起着关键作用。然而,POAG与铁死亡的分子机制尚未完全明确,基于生物信息学的方法探索铁死亡相关的最佳特征基因与通路,为阐明POAG的分子机制提供有力支持(详细流程见图1)。

图1 原发性开角型青光眼铁死亡相关基因鉴定的详细流程。

1 材料和方法

1.1材料本研究的POAG样本和正常样本的基因表达数据是从GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)中获得。其中将小梁网来源的GSE27276数据集[9]作为训练集,芯片信息:Sentrix Human-6 Expression BeadChip,平台是GPL2507;视乳头来源的GSE2378和GSE9944数据集[10-11]作为验证集,芯片信息:(1)Affymetrix Human Genome U95 Version 2 Array,GPL8300;(2)Affymetrix Human Genome U133A 2.0 Array,GPL571(表1)。此外,从FerrDb数据库(http://www.datjar.com:40013/bt2104/)获取铁死亡相关的调控基因FRG(n=728)。

表1 GEO数据芯片及其平台数据

1.2方法

1.2.1差异表达基因的筛选及相关性分析运用R语言(ver4.2.2)limma和pheatmap包对数据集基因与铁死亡基因进行映射,以P<0.05和|Log2FC|>1作为过滤条件,获取与铁死亡相关的预后差异表达基因(differentially expressed ferroptosis-related genes, DE-FRGs)以及表达量,并绘制DE-FRGs在正常(Control)组与原发性开角型青光眼(POAG)组的相关性热图。为了进一步了解DE-FRGs的相关程度,采用corrlot包中的Spearman相关对上调和下调的基因进行相关性分析。

1.2.2 GO功能富集分析和KEGG通路富集分析对DE-FRGs进行功能与通路富集化分析,以确定相关的信号通路,并揭示与生物过程相关的潜在分子机制。使用ClusterProfiler包进行基因本体(gene ontology, GO)富集分析和京都基因与基因组百科全书(kyoto encyclopedia of genes and genomes, KEGG)通路富集分析,两者均以P<0.05作为过滤条件。

1.2.3 POAG的最佳特征基因的筛选通过glmnet包中的LASSO算法对DE-FRGs进行筛选,鉴定POAG的特征基因。同时,用SVM包构建支持向量机-递归特征消除(SVM-RFE)模型,通过其10倍交叉验证的平均误判率进行比较,得到POAG的关键基因。将两种模型的筛选结果取交集,获得最佳特征基因。其诊断能力是通过ROC曲线进行判断,测量曲线下面积(area under curve, AUC)可评估特征基因的准确性、敏感性和特异性表达情况。在此基础上,利用GLM包构建了基于关键基因的Logistic回归模型,预测训练数据集中的样本类型,Logistic回归模型的诊断能力通过ROC曲线进行评估。最后,运用R包ggpubr将得到的最佳特征基因在POAG组和正常组中的表达进行可视化展示。

1.2.4单基因基因组富集分析为了进一步探索最佳特征基因的相关途径,运用R中的enrichplot包进行单基因基因组富集分析(gene set enrichment analysis, GSEA)。通过分析训练数据集中所有其他基因与最佳特征基因的相关度,并根据其相关度的大小进行排序,排序后的基因则被认为是待检测的基因集。最后,调用KEGG信号通路集作为预定义的基因集,检测其在基因集中的富集程度。

1.2.5单基因基因组变异分析以KEGG通路组作为背景基因组,使用gsva、ggpub和rlimma包对18个最佳特征基因进行单基因基因组变异分析(gene set variation analysis, GSVA),计算其样本的GSVA得分差异,以P<0.05和|t|>2作为差异过滤条件,获取最佳特征基因在高表达组和低表达组样本的差异情况。若t>0,则提示该通路在高表达组中被激活,反之,则提示该通路在低表达组中被激活。

1.2.6最佳特征基因外部数据集的交叉验证利用sav和limma包对3个数据集的基因进行合并以及重复基因的表达量取均值,并对整合的数据进行批次矫正,以相同的筛选标准分析整合的验证数据集中特征基因的表达水平。最后与GSE27276训练数据集中最佳特征基因的表达水平进行对比。

2 结果

2.1差异表达铁死亡相关基因的鉴定及相关性分析在728个铁死亡基因中确定了396个POAG相关的铁死亡基因,其中上调基因39个,下调基因64个(表2)。并绘制103个在POAG组和正常组之间差异表达DE-FRGs的热图(图2A)。为了进一步探索上调基因和下调基因之间的相关度,利用生物信息学的统计学方法进行相关性分析,结果表明上调基因和下调基因均有较高相关性(图2B)。

图2 差异表达基因的筛选 A:DE-FRGs的差异表达图;B:DE-FRGs的相关性热图。

表2 103个差异表达的铁死亡相关的上调基因和下调基因

2.2差异基因的GO功能富集分析和KEGG通路富集分析DE-FRGs主要参与对氧化应激的反应、细胞对化学应激的反应、对养分水平的反应、细胞对氧化应激的反应和对金属离子的反应等生物学过程(biological processes, BP);主要定位于外膜、细胞器外膜、线粒体外膜、自噬体、自噬体组装位点等细胞成分(cell component, CC);主要参与泛素蛋白连接酶结合、泛素样蛋白连接酶结合、血红素结合、双氧化活性和亚铁结合等分子功能(molecular function, MF)(图3A,表3)。KEGG富集分析主要集中在铁死亡、线粒体自噬、Kaposi肉瘤相关疱疹病毒感染、自噬和流体剪切力与动脉粥样硬化等相关信号通路(图3B)。

表3 GO功能富集分析BP、CC和MF排名前6通路

2.3 POAG最佳特征基因的筛选在GSE27276数据集中运用LASSO和SVM-RFE两种不同的机器学习算法,筛选出POAG铁死亡相关特征基因。一方面,LASSO逻辑回归算法筛选出18个与POAG相关的铁死亡特征基因(图4A、B)。另一方面,SVM-RFE算法确定103个DE-FRGs均为铁死亡特征基因(图4D、E)。LASSO和SVM-RFE模型得到的特征基因取交集处理,确定了18个最佳特征基因(CYBB、SLC1A5、HILPDA、AQP5、SLC25A28、MFN2、IFNA13、FADS2、SCD、OTUB1、ZFP36、PROM2、GCH1、GDF15、CREB1、FABP4、MPC1和LCN2)进行后续分析(图4F)。基于上述18个最佳特征基因的ROC曲线结果表明,Logistic回归模型可以区分正常和POAG样本,AUC=1,95%CI:1.000～1.000(图4C);为了区别单个基因的POAG样品和正常样本,生成了18个基因的ROC曲线,显示18个基因的AUC均大于0.7(图4G)。同时CYBB、FADS2、SCD、CREB1和MPC1在训练集中高表达;SLC1A5、HILPDA、AQP5、SLC25A28、MFN2、IFNA13、OTUB1、ZFP36、PROM2、GCH1、GDF15和FABP4在训练集中低表达(图5)。

图4 最佳特征基因的筛选 A、B:LASSO回归模型;C:用于识别疾病样本AUC的逻辑回归模型;D、E:SVM-RFE模型;F:LASSO和SVM-RFE模型交集基因的韦恩图;G:18个标记基因的ROC曲线。

2.4最佳特征基因与各种POAG的通路密切相关通过单基因GSEA-KEGG途径分析,旨在找到特征基因区分POAG样本和正常样本的潜在通路。提取每个标记基因富集的前6条通路(图6,表3)。发现这些基因富集在核糖体、氧化磷酸化、血管平滑肌收缩、抗原处理和提呈、溶酶体和各种疾病途径(肥厚型心肌病、系统性红斑狼疮、婴儿利什曼原虫、阿尔茨海默病、亨廷顿病、帕金森病、病毒性心肌炎)。特征基因主要富集在钙离子信号通道和趋化因子信号通路。此外,GDF15、MFN2和OTUB1基因与谷胱甘肽代谢通路密切相关。

图6 GSEA富集分析 A:AQP5;B:CREB1;C:FABP4;D:CYBB;E:FADS2;F:GCH1;G:GDF15;H:HILPDA;I:IFNA13;J:LCN2;K:MFN2;L:MPC1;M:OTUB1;N:PROM2;O:SCD;P:SLC1A5;Q:SLC25A28;R:ZFP36。

根据各特征基因的表达水平结合GSVA分析,观察高、低表达组之间的差异性激活途径。MFN2与POAG发病机制有关的各种途径在其高表达组中被富集(图7A)。如谷胱甘肽代谢、烟酸和烟酰胺代谢、哺乳动物的昼夜节律、核糖核酸聚合酶、类固醇生物合成、N-聚糖的生物合成、糖胺聚糖生物合成-硫酸软骨素、不饱和脂肪酸的生物合成、叶酸生物合成;而低表达组与刺猬信号通路、核黄素代谢、嗅觉转导、神经活性配体-受体相互作用、味觉转导密切相关。OTUB1在疾病中的低表达可能通过激活谷胱甘肽代谢、Notch信号通路、三羧酸循环、基础转录因子、不饱和脂肪酸的生物合成、哺乳动物的昼夜节律来诱导POAG;而OTUB1的过表达则激活了肾素-血管紧张素系统、α-亚油酸代谢、亚油酸代谢、嗅觉转导、抗坏血酸和醛酸代谢、神经活性配体-受体相互作用(图7B)。

图7 GSVA富集分析 A:MFN2;B:OTUB1。

2.5最佳特征基因的外部数据集的交叉验证为了进一步验证训练集的可靠性,通过视乳头来源的验证集进行上述18个特征基因的表达水平验证,结果显示,POAG视乳头标本中CREB1的表达水平相对于正常样本显著升高,与GSE27276数据集小梁网样本的表达一致;GCH1和FADS2的表达水平相对于正常样本显著降低,与GSE27276数据集的表达水平却相反;而9个基因包括CYBB、SLC1A5、AQP5、MFN2、SCD、ZFP36、GDF15、FABP4和LCN2则在验证集中表达差异无统计学意义(图8)。

图8 最佳特征基因外部数据集的交叉验证 A:AQP5;B:CREB1;C:CYBB;D:GCH1;E:FADS2;F:FABP4;G:LCN2;H:ZFP36;I:MFN2;J:SLC1A5;K:SCD;L:GDF15。

3 讨论

POAG作为一种进行性致盲性眼病,但其发病机制尚不完全清楚。虽然很多研究集中在眼压的作用上,但其他因素,如眼球血流异常、覆膜结构的异常易感性、低颅压、自身免疫和线粒体功能障碍也可能参与[12]。同时,已经确定多种危险因素可能与POAG的发病有关,例如年龄、眼压升高、家族史等[13-14]。GO富集分析和KEGG富集分析结果提示,103个DE-FRGs参与氧化应激的反应和铁死亡相关的通路等。氧化应激影响的主要组织包括小梁网(TM)[15]和RGC[16]。在POAG中,过量活性氧的积累可诱发TM损伤,从而导致房水传统流出途径受损[17],并加剧视乳头和RGC的损伤[18]。在过去的10a中,越来越多的证据表明铁死亡与青光眼密切相关[19-21]。蛋白质组学分析[22]表明,铁死亡参与了N-甲基-d-天冬氨酸(N-Methyl-D-aspartic acid ,NMDA)诱导的雄性Sprague-Dawley大鼠视网膜和视神经细胞死亡,并且发现ACSL3(Q63151)和Prnp(P13852)蛋白在NMDA损伤的视网膜中上调并与铁死亡有关。有实验研究证实[23],在高眼压的青光眼模型中,使用去铁酮(一种铁螯合剂)干预9wk后,视神经损失减少的RGC数量显著增加,这表明铁螯合对青光眼具有潜在的保护作用。以上研究均表明铁死亡参与了POAG发生的重要环节。

本研究共筛选出CYBB、SLC1A5、HILPDA、AQP5、SLC25A28、MFN2、IFNA13、FADS2、SCD、OTUB1、ZFP36、PROM2、GCH1、GDF15、CREB1、FABP4、MPC1和LCN2 18个与铁死亡有关的特征基因。18个基因的ROC曲线下面积所代表的AUC值均大于0.7,说明这18个基因对于区分POAG样本和正常样本有一定的准确性和特异性。在铁死亡中已经确定谷胱甘肽过氧化物酶-4-谷胱甘肽(glutathione peroxidase 4, GPX-4-GSH)作为主流途径,而谷胱甘肽代谢通路是GPX-4-GSH途径的重要机制之一[24]。GSEA分析的结果提示GDF15、MFN2和OTUB1基因均与谷胱甘肽代谢通路密切相关。GDF15是TGF-β超家族的成员[25],先前已被证实在视神经轴索损伤后视网膜中表达上调[26],而后被检测到房水(aqueous humor, AH)的GDF15表达增加[27],AH中的GDF15则来源于视网膜神经节细胞层(GCL),认为GCL分泌到AH中的GDF15可能是青光眼性神经变性的可量化指标。MFN2位于1号染色体短臂1p36.22位置上,编码一种线粒体膜蛋白,参与线粒体的聚集和融合,并有助于线粒体网络的维持和运行[28]。有研究基于多标记单倍型的关联检验证实了POAG与MFN2的关联[29]。MFN2的GSAV分析也证实在高表达组,谷胱甘肽代谢通路被激活。OTUB1以其去泛素化的特性,首先被确定为OTU家族的核心成员[30],位于人染色体11q13.1上,其基因产物是由271个氨基酸组成的蛋白质,OTUB1蛋白的大小为31kDa,是由一组保守的半胱氨酸、组氨酸、天冬氨酸残基组成的半胱氨酸蛋白酶[31-33]。目前对该基因的研究主要集中在免疫和癌症等疾病方面[34],其在POAG中的作用尚待进一步研究。OTUB1的GSVA分析发现,谷胱甘肽代谢通路在低表达组中高度表达,这可能是身体对疾病的一种保护机制,控制疾病的继续发展。

在POAG的发生发展过程中,早期多为小梁网组织变化,中后期多为视乳头改变及其神经节细胞凋亡,伴随神经纤维层丢失。故将小梁网组织作为训练集,借助视乳头样本进行外部数据集的交叉验证,一则为增加最佳特征基因的局部特异性,二则能有效识别出贯穿于POAG的疾病发生发展全过程的潜在生物标志物。外部数据集交叉验证的结果显示,POAG视乳头标本中CREB1的表达水平相对于正常样本显著升高,在训练集(GSE27276小梁网样本数据集)与验证集(GSE2378、GSE9944视乳头样本数据集)的表达一致,并且该基因先前均已在视神经相关疾病或铁死亡方面有过报道。CREB1也称为CREB,该转录调节因子响应有害应激(如缺氧和氧化应激)而被激活,并参与细胞对这些应激的防御[35],并促进神经元存活[36]。有研究表明CREB1是青光眼最重要的上游调节因子[37],在POAG的发病机制中存在潜在的关键保护作用[38]。Guo等[39]从在分子水平上,已经确定CREB是钙离子-钙调蛋白依赖性蛋白激酶Ⅱ(calcium-cam-dependent protein kinase Ⅱ, CaMKⅡ)的关键下游效应子,当RGC受到兴奋性毒性或视神经损伤时,CREB的转录活性急剧下降,CaMKⅡ再激活可以挽救受损的CREB活性,并且CREB活性对于CaMKⅡ介导的RGC起到充分保护,从而证实介导兴奋性毒性和在视神经损伤模型中起到RGC的保护作用。

综上所述,本研究基于生物信息学挖掘构建了POAG中铁死亡相关差异基因的模型,运用两种机器学习的模型和外部数据集交叉验证,筛选出CREB1表达一致的最佳疾病特征基因作为潜在的诊断生物标志物,为研究POAG在铁死亡方面的发病机制和治疗方法提供了新见解。但仍存在一定的局限性,主要体现在以下几个方面:(1)虽然本研究增加了外部数据集的交叉验证,但样本数量特别有限,尤其是训练集中同一样本来源的数据集;(2)缺乏进一步的体内外实验来验证POAG铁死亡相关特征基因的诊断价值及作用机制,后续将从分子生物学和临床试验等角度解决偏倚结果和结论方面可能存在的局限性。