原發(fā)性開角型青光眼鐵死亡相關(guān)基因表達譜的生物信息學分析

周日龍,廖 良

目的:通過生物信息學手段分析原發(fā)性開角型青光眼(POAG)進展過程中調(diào)控鐵死亡相關(guān)的關(guān)鍵基因,旨在進一步揭示鐵死亡在POAG中的生物學機制。

方法:從GEO數(shù)據(jù)庫中獲取小梁網(wǎng)來源的GSE27276數(shù)據(jù)集,其中包括19個POAG小梁網(wǎng)組織樣本和17個正常小梁網(wǎng)組織樣本;下載FerrDb數(shù)據(jù)庫整理的鐵死亡相關(guān)基因,將GSE27276數(shù)據(jù)集與鐵死亡基因集進行映射,篩選POAG中鐵死亡相關(guān)的預后差異表達基因(DE-FRGs)并進行相關(guān)性分析,進一步了解DE-FRGs的GO和KEGG通路富集。應用LASSO回歸模型與SVM-RFE模型兩種機器學習的算法篩選鐵死亡相關(guān)POAG的關(guān)鍵基因,將兩種模型的篩選結(jié)果取交集,獲得最佳特征基因,使用受試者特征曲線(ROC)評估臨床診斷能力;對最佳特征基因進行單基因基因組富集分析(GSEA)和變異分析(GSVA);借助視乳頭來源的GSE2378與GSE9944數(shù)據(jù)集驗證最佳特征基因的表達水平。

結(jié)果:與正常小梁網(wǎng)組織相比,POAG的小梁網(wǎng)組織有396個鐵死亡基因存在差異表達,其中39個為上調(diào)基因,64個為下調(diào)基因,Spearman相關(guān)性分析顯示上調(diào)基因和下調(diào)基因均有一定的相關(guān)性。GO功能和KEGG通路富集分析顯示,差異基因主要富集在氧化應激反應和鐵死亡通路;通過LASSO和SVM-RFE算法將18個DE-FRGs確定為關(guān)鍵基因,具有更高的診斷價值。GSEA和GSVA富集分析顯示GDF15、MFN2和OTUB1基因與谷胱甘肽代謝通路密切相關(guān),其中MFN2和OTUB1分別在高表達組和低表達組中激活谷胱甘肽代謝通路。GSE2378與GSE9944數(shù)據(jù)集交叉驗證明確視乳頭標本中CREB1的表達水平相對于正常視乳頭樣本顯著升高,這與GSE27276數(shù)據(jù)集小梁網(wǎng)樣本表達一致。

結(jié)論:基于生物信息學分析挖掘得到396個POAG的DE-FRGs,通過構(gòu)建機器篩選模型和外部數(shù)據(jù)集交叉驗證,篩選出CREB1有望成為潛在診斷生物標志物的最佳特征基因,為進一步深入闡明POAG鐵死亡相關(guān)的分子機制和診斷提供靶點。同時篩選的基因還需要進一步體內(nèi)、外實驗驗證,揭示鐵死亡在POAG中的生物學機制。

?KEYWORDS：primary open angle glaucoma; ferroptosis; bioinformatics

0 引言

原發(fā)性開角型青光眼(primary open angle glaucoma, POAG)是一種慢性進行性視神經(jīng)病變,其特征為視網(wǎng)膜神經(jīng)節(jié)細胞(retinal ganglion cell, RGC)進行性變性和視野缺失,好發(fā)于40～80歲人群[1-2]。POAG雖不危及生命,但具有早期隱匿性強、中后期致盲率高等特點[3],已成為導致中老年致盲最常見的原因之一。如果未得到及時有效的治療,POAG引起的視力喪失不僅影響生活質(zhì)量,而且造成重大的社會經(jīng)濟負擔[4]。目前,其發(fā)病機制尚不完全清楚,而遺傳學被證明在POAG的發(fā)病機制中起關(guān)鍵作用[3-5]。鐵死亡是2012年新發(fā)現(xiàn)的一種鐵依賴性程序性細胞死亡模式,其發(fā)生機制不同于細胞焦亡、凋亡、壞死和自噬[6]。鐵死亡的特征是線粒體萎縮和線粒體膜密度增加,鐵和脂質(zhì)活性氧的積累[7]。鐵死亡的發(fā)生機制是通過脂質(zhì)自由基的形成和谷胱甘肽(L-glutathione, GSH)的消耗或脂質(zhì)過氧化物酶4(glutathione peroxidase 4,GPX4)的失活來催化的[8]。鐵循環(huán)在鐵死亡的發(fā)展中起著關(guān)鍵作用。然而,POAG與鐵死亡的分子機制尚未完全明確,基于生物信息學的方法探索鐵死亡相關(guān)的最佳特征基因與通路,為闡明POAG的分子機制提供有力支持(詳細流程見圖1)。

圖1 原發(fā)性開角型青光眼鐵死亡相關(guān)基因鑒定的詳細流程。

1 材料和方法

1.1材料本研究的POAG樣本和正常樣本的基因表達數(shù)據(jù)是從GEO數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/geo/)中獲得。其中將小梁網(wǎng)來源的GSE27276數(shù)據(jù)集[9]作為訓練集,芯片信息:Sentrix Human-6 Expression BeadChip,平臺是GPL2507;視乳頭來源的GSE2378和GSE9944數(shù)據(jù)集[10-11]作為驗證集,芯片信息:(1)Affymetrix Human Genome U95 Version 2 Array,GPL8300;(2)Affymetrix Human Genome U133A 2.0 Array,GPL571(表1)。此外,從FerrDb數(shù)據(jù)庫(http://www.datjar.com:40013/bt2104/)獲取鐵死亡相關(guān)的調(diào)控基因FRG(n=728)。

表1 GEO數(shù)據(jù)芯片及其平臺數(shù)據(jù)

1.2方法

1.2.1差異表達基因的篩選及相關(guān)性分析運用R語言(ver4.2.2)limma和pheatmap包對數(shù)據(jù)集基因與鐵死亡基因進行映射,以P<0.05和|Log2FC|>1作為過濾條件,獲取與鐵死亡相關(guān)的預后差異表達基因(differentially expressed ferroptosis-related genes, DE-FRGs)以及表達量,并繪制DE-FRGs在正常(Control)組與原發(fā)性開角型青光眼(POAG)組的相關(guān)性熱圖。為了進一步了解DE-FRGs的相關(guān)程度,采用corrlot包中的Spearman相關(guān)對上調(diào)和下調(diào)的基因進行相關(guān)性分析。

1.2.2 GO功能富集分析和KEGG通路富集分析對DE-FRGs進行功能與通路富集化分析,以確定相關(guān)的信號通路,并揭示與生物過程相關(guān)的潛在分子機制。使用ClusterProfiler包進行基因本體(gene ontology, GO)富集分析和京都基因與基因組百科全書(kyoto encyclopedia of genes and genomes, KEGG)通路富集分析,兩者均以P<0.05作為過濾條件。

1.2.3 POAG的最佳特征基因的篩選通過glmnet包中的LASSO算法對DE-FRGs進行篩選,鑒定POAG的特征基因。同時,用SVM包構(gòu)建支持向量機-遞歸特征消除(SVM-RFE)模型,通過其10倍交叉驗證的平均誤判率進行比較,得到POAG的關(guān)鍵基因。將兩種模型的篩選結(jié)果取交集,獲得最佳特征基因。其診斷能力是通過ROC曲線進行判斷,測量曲線下面積(area under curve, AUC)可評估特征基因的準確性、敏感性和特異性表達情況。在此基礎上,利用GLM包構(gòu)建了基于關(guān)鍵基因的Logistic回歸模型,預測訓練數(shù)據(jù)集中的樣本類型,Logistic回歸模型的診斷能力通過ROC曲線進行評估。最后,運用R包ggpubr將得到的最佳特征基因在POAG組和正常組中的表達進行可視化展示。

1.2.4單基因基因組富集分析為了進一步探索最佳特征基因的相關(guān)途徑,運用R中的enrichplot包進行單基因基因組富集分析(gene set enrichment analysis, GSEA)。通過分析訓練數(shù)據(jù)集中所有其他基因與最佳特征基因的相關(guān)度,并根據(jù)其相關(guān)度的大小進行排序,排序后的基因則被認為是待檢測的基因集。最后,調(diào)用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最佳特征基因外部數(shù)據(jù)集的交叉驗證利用sav和limma包對3個數(shù)據(jù)集的基因進行合并以及重復基因的表達量取均值,并對整合的數(shù)據(jù)進行批次矯正,以相同的篩選標準分析整合的驗證數(shù)據(jù)集中特征基因的表達水平。最后與GSE27276訓練數(shù)據(jù)集中最佳特征基因的表達水平進行對比。

2 結(jié)果

2.1差異表達鐵死亡相關(guān)基因的鑒定及相關(guān)性分析在728個鐵死亡基因中確定了396個POAG相關(guān)的鐵死亡基因,其中上調(diào)基因39個,下調(diào)基因64個(表2)。并繪制103個在POAG組和正常組之間差異表達DE-FRGs的熱圖(圖2A)。為了進一步探索上調(diào)基因和下調(diào)基因之間的相關(guān)度,利用生物信息學的統(tǒng)計學方法進行相關(guān)性分析,結(jié)果表明上調(diào)基因和下調(diào)基因均有較高相關(guān)性(圖2B)。

圖2 差異表達基因的篩選 A:DE-FRGs的差異表達圖;B:DE-FRGs的相關(guān)性熱圖。

表2 103個差異表達的鐵死亡相關(guān)的上調(diào)基因和下調(diào)基因

2.2差異基因的GO功能富集分析和KEGG通路富集分析DE-FRGs主要參與對氧化應激的反應、細胞對化學應激的反應、對養(yǎng)分水平的反應、細胞對氧化應激的反應和對金屬離子的反應等生物學過程(biological processes, BP);主要定位于外膜、細胞器外膜、線粒體外膜、自噬體、自噬體組裝位點等細胞成分(cell component, CC);主要參與泛素蛋白連接酶結(jié)合、泛素樣蛋白連接酶結(jié)合、血紅素結(jié)合、雙氧化活性和亞鐵結(jié)合等分子功能(molecular function, MF)(圖3A,表3)。KEGG富集分析主要集中在鐵死亡、線粒體自噬、Kaposi肉瘤相關(guān)皰疹病毒感染、自噬和流體剪切力與動脈粥樣硬化等相關(guān)信號通路(圖3B)。

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

2.3 POAG最佳特征基因的篩選在GSE27276數(shù)據(jù)集中運用LASSO和SVM-RFE兩種不同的機器學習算法,篩選出POAG鐵死亡相關(guān)特征基因。一方面,LASSO邏輯回歸算法篩選出18個與POAG相關(guān)的鐵死亡特征基因(圖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)進行后續(xù)分析(圖4F)。基于上述18個最佳特征基因的ROC曲線結(jié)果表明,Logistic回歸模型可以區(qū)分正常和POAG樣本,AUC=1,95%CI:1.000～1.000(圖4C);為了區(qū)別單個基因的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的通路密切相關(guān)通過單基因GSEA-KEGG途徑分析,旨在找到特征基因區(qū)分POAG樣本和正常樣本的潛在通路。提取每個標記基因富集的前6條通路(圖6,表3)。發(fā)現(xiàn)這些基因富集在核糖體、氧化磷酸化、血管平滑肌收縮、抗原處理和提呈、溶酶體和各種疾病途徑(肥厚型心肌病、系統(tǒng)性紅斑狼瘡、嬰兒利什曼原蟲、阿爾茨海默病、亨廷頓病、帕金森病、病毒性心肌炎)。特征基因主要富集在鈣離子信號通道和趨化因子信號通路。此外,GDF15、MFN2和OTUB1基因與谷胱甘肽代謝通路密切相關(guān)。

圖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。

根據(jù)各特征基因的表達水平結(jié)合GSVA分析,觀察高、低表達組之間的差異性激活途徑。MFN2與POAG發(fā)病機制有關(guān)的各種途徑在其高表達組中被富集(圖7A)。如谷胱甘肽代謝、煙酸和煙酰胺代謝、哺乳動物的晝夜節(jié)律、核糖核酸聚合酶、類固醇生物合成、N-聚糖的生物合成、糖胺聚糖生物合成-硫酸軟骨素、不飽和脂肪酸的生物合成、葉酸生物合成;而低表達組與刺猬信號通路、核黃素代謝、嗅覺轉(zhuǎn)導、神經(jīng)活性配體-受體相互作用、味覺轉(zhuǎn)導密切相關(guān)。OTUB1在疾病中的低表達可能通過激活谷胱甘肽代謝、Notch信號通路、三羧酸循環(huán)、基礎轉(zhuǎn)錄因子、不飽和脂肪酸的生物合成、哺乳動物的晝夜節(jié)律來誘導POAG;而OTUB1的過表達則激活了腎素-血管緊張素系統(tǒng)、α-亞油酸代謝、亞油酸代謝、嗅覺轉(zhuǎn)導、抗壞血酸和醛酸代謝、神經(jīng)活性配體-受體相互作用(圖7B)。

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

2.5最佳特征基因的外部數(shù)據(jù)集的交叉驗證為了進一步驗證訓練集的可靠性,通過視乳頭來源的驗證集進行上述18個特征基因的表達水平驗證,結(jié)果顯示,POAG視乳頭標本中CREB1的表達水平相對于正常樣本顯著升高,與GSE27276數(shù)據(jù)集小梁網(wǎng)樣本的表達一致;GCH1和FADS2的表達水平相對于正常樣本顯著降低,與GSE27276數(shù)據(jù)集的表達水平卻相反;而9個基因包括CYBB、SLC1A5、AQP5、MFN2、SCD、ZFP36、GDF15、FABP4和LCN2則在驗證集中表達差異無統(tǒng)計學意義(圖8)。

圖8 最佳特征基因外部數(shù)據(jù)集的交叉驗證 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作為一種進行性致盲性眼病,但其發(fā)病機制尚不完全清楚。雖然很多研究集中在眼壓的作用上,但其他因素,如眼球血流異常、覆膜結(jié)構(gòu)的異常易感性、低顱壓、自身免疫和線粒體功能障礙也可能參與[12]。同時,已經(jīng)確定多種危險因素可能與POAG的發(fā)病有關(guān),例如年齡、眼壓升高、家族史等[13-14]。GO富集分析和KEGG富集分析結(jié)果提示,103個DE-FRGs參與氧化應激的反應和鐵死亡相關(guān)的通路等。氧化應激影響的主要組織包括小梁網(wǎng)(TM)[15]和RGC[16]。在POAG中,過量活性氧的積累可誘發(fā)TM損傷,從而導致房水傳統(tǒng)流出途徑受損[17],并加劇視乳頭和RGC的損傷[18]。在過去的10a中,越來越多的證據(jù)表明鐵死亡與青光眼密切相關(guān)[19-21]。蛋白質(zhì)組學分析[22]表明,鐵死亡參與了N-甲基-d-天冬氨酸(N-Methyl-D-aspartic acid ,NMDA)誘導的雄性Sprague-Dawley大鼠視網(wǎng)膜和視神經(jīng)細胞死亡,并且發(fā)現(xiàn)ACSL3(Q63151)和Prnp(P13852)蛋白在NMDA損傷的視網(wǎng)膜中上調(diào)并與鐵死亡有關(guān)。有實驗研究證實[23],在高眼壓的青光眼模型中,使用去鐵酮(一種鐵螯合劑)干預9wk后,視神經(jīng)損失減少的RGC數(shù)量顯著增加,這表明鐵螯合對青光眼具有潛在的保護作用。以上研究均表明鐵死亡參與了POAG發(fā)生的重要環(huán)節(jié)。

本研究共篩選出CYBB、SLC1A5、HILPDA、AQP5、SLC25A28、MFN2、IFNA13、FADS2、SCD、OTUB1、ZFP36、PROM2、GCH1、GDF15、CREB1、FABP4、MPC1和LCN2 18個與鐵死亡有關(guān)的特征基因。18個基因的ROC曲線下面積所代表的AUC值均大于0.7,說明這18個基因?qū)τ趨^(qū)分POAG樣本和正常樣本有一定的準確性和特異性。在鐵死亡中已經(jīng)確定谷胱甘肽過氧化物酶-4-谷胱甘肽(glutathione peroxidase 4, GPX-4-GSH)作為主流途徑,而谷胱甘肽代謝通路是GPX-4-GSH途徑的重要機制之一[24]。GSEA分析的結(jié)果提示GDF15、MFN2和OTUB1基因均與谷胱甘肽代謝通路密切相關(guān)。GDF15是TGF-β超家族的成員[25],先前已被證實在視神經(jīng)軸索損傷后視網(wǎng)膜中表達上調(diào)[26],而后被檢測到房水(aqueous humor, AH)的GDF15表達增加[27],AH中的GDF15則來源于視網(wǎng)膜神經(jīng)節(jié)細胞層(GCL),認為GCL分泌到AH中的GDF15可能是青光眼性神經(jīng)變性的可量化指標。MFN2位于1號染色體短臂1p36.22位置上,編碼一種線粒體膜蛋白,參與線粒體的聚集和融合,并有助于線粒體網(wǎng)絡的維持和運行[28]。有研究基于多標記單倍型的關(guān)聯(lián)檢驗證實了POAG與MFN2的關(guān)聯(lián)[29]。MFN2的GSAV分析也證實在高表達組,谷胱甘肽代謝通路被激活。OTUB1以其去泛素化的特性,首先被確定為OTU家族的核心成員[30],位于人染色體11q13.1上,其基因產(chǎn)物是由271個氨基酸組成的蛋白質(zhì),OTUB1蛋白的大小為31kDa,是由一組保守的半胱氨酸、組氨酸、天冬氨酸殘基組成的半胱氨酸蛋白酶[31-33]。目前對該基因的研究主要集中在免疫和癌癥等疾病方面[34],其在POAG中的作用尚待進一步研究。OTUB1的GSVA分析發(fā)現(xiàn),谷胱甘肽代謝通路在低表達組中高度表達,這可能是身體對疾病的一種保護機制,控制疾病的繼續(xù)發(fā)展。

在POAG的發(fā)生發(fā)展過程中,早期多為小梁網(wǎng)組織變化,中后期多為視乳頭改變及其神經(jīng)節(jié)細胞凋亡,伴隨神經(jīng)纖維層丟失。故將小梁網(wǎng)組織作為訓練集,借助視乳頭樣本進行外部數(shù)據(jù)集的交叉驗證,一則為增加最佳特征基因的局部特異性,二則能有效識別出貫穿于POAG的疾病發(fā)生發(fā)展全過程的潛在生物標志物。外部數(shù)據(jù)集交叉驗證的結(jié)果顯示,POAG視乳頭標本中CREB1的表達水平相對于正常樣本顯著升高,在訓練集(GSE27276小梁網(wǎng)樣本數(shù)據(jù)集)與驗證集(GSE2378、GSE9944視乳頭樣本數(shù)據(jù)集)的表達一致,并且該基因先前均已在視神經(jīng)相關(guān)疾病或鐵死亡方面有過報道。CREB1也稱為CREB,該轉(zhuǎn)錄調(diào)節(jié)因子響應有害應激(如缺氧和氧化應激)而被激活,并參與細胞對這些應激的防御[35],并促進神經(jīng)元存活[36]。有研究表明CREB1是青光眼最重要的上游調(diào)節(jié)因子[37],在POAG的發(fā)病機制中存在潛在的關(guān)鍵保護作用[38]。Guo等[39]從在分子水平上,已經(jīng)確定CREB是鈣離子-鈣調(diào)蛋白依賴性蛋白激酶Ⅱ(calcium-cam-dependent protein kinase Ⅱ, CaMKⅡ)的關(guān)鍵下游效應子,當RGC受到興奮性毒性或視神經(jīng)損傷時,CREB的轉(zhuǎn)錄活性急劇下降,CaMKⅡ再激活可以挽救受損的CREB活性,并且CREB活性對于CaMKⅡ介導的RGC起到充分保護,從而證實介導興奮性毒性和在視神經(jīng)損傷模型中起到RGC的保護作用。

綜上所述,本研究基于生物信息學挖掘構(gòu)建了POAG中鐵死亡相關(guān)差異基因的模型,運用兩種機器學習的模型和外部數(shù)據(jù)集交叉驗證,篩選出CREB1表達一致的最佳疾病特征基因作為潛在的診斷生物標志物,為研究POAG在鐵死亡方面的發(fā)病機制和治療方法提供了新見解。但仍存在一定的局限性,主要體現(xiàn)在以下幾個方面:(1)雖然本研究增加了外部數(shù)據(jù)集的交叉驗證,但樣本數(shù)量特別有限,尤其是訓練集中同一樣本來源的數(shù)據(jù)集;(2)缺乏進一步的體內(nèi)外實驗來驗證POAG鐵死亡相關(guān)特征基因的診斷價值及作用機制,后續(xù)將從分子生物學和臨床試驗等角度解決偏倚結(jié)果和結(jié)論方面可能存在的局限性。