基于生物信息學和機器學習算法鑒定PDR鐵死亡關鍵基因

徐徑舟 黃晉 李智豪 張赟 陳鼎

作者單位：溫州醫科大學附屬眼視光醫院，溫州 325027

糖尿病視網膜病變（Diabetic retinopathy，DR）是糖尿病引起的常見眼部并發癥之一，近年來發病率逐年升高，已成為全球主要的致盲性眼病[1]。根據疾病發展的不同階段，DR可分為非增殖期和增殖期，其中增殖性糖尿病視網膜病變（Proliferative diabetic retinopathy，PDR）預后較差，常常造成不可逆的視力損傷。以往的研究認為，DR的發生機制主要涉及氧化應激[2]、慢性炎癥損傷[3]等。最新的研究表明，凋亡、焦亡、鐵死亡等調節性細胞死亡也與DR的發生發展密切相關[4]，這為我們探究DR背后的病理機制提供了新的線索。

鐵死亡是一種鐵離子依賴的、以細胞內脂質過氧化產物超量蓄積為主要特征的調節性細胞死亡方式[5]。在細胞形態學上，鐵死亡主要表現為線粒體膜皺縮、線粒體嵴減少或消失、細胞核無明顯變化等[6]。在分子水平，鐵死亡主要表現為細胞內Fe2+含量升高、谷胱甘肽含量降低、谷胱甘肽過氧化物酶4 （Glutathione peroxidase 4，GPX4）活性下降、脂質過氧化產物及活性氧大量蓄積等[7]。目前已知鐵死亡參與到多種疾病的病理機制中，包括神經退行性疾病、心血管疾病等，鐵死亡是當前疾病基礎機制研究的熱點問題[8-9]。

本研究的目的是綜合生物信息學和機器學習算法，從正常視網膜組織與PDR組織的轉錄組測序數據中篩選出與疾病關系最密切的鐵死亡差異基因（Differentially expressed ferroptosis-related genes，DEFRGs），進而為探究PDR潛在的分子機制與藥物治療靶點提供新的參考。

1 材料與方法

1.1 數據集的獲取

本次研究所用數據集來自美國國家生物信息中心的基因表達綜合數據庫（Gene expression omnibus，GEO），共獲得GSE102485和GSE60436 2個數據集。其中GSE102485數據集中包含3個正常視網膜組織樣本和22個PDR組織樣本，GSE60436數據集中包含3個正常視網膜組織樣本和6個PDR組織樣本。GSE102485數據集用于篩選疾病特征基因，GSE60436數據集用于后續驗證。

從FerrDb數據庫（www.zhounan.org/ferrdb/）中下載鐵死亡相關基因，包括驅動基因、抑制基因、標記基因，共484個鐵死亡相關基因。本研究遵守GEO數據庫和FerrDb數據庫使用規則。

1.2 鐵死亡相關差異基因篩選

使用R軟件limma軟件包，以|log2FC|>1 和P<0.05 為閾值標準對GSE102485 數據集中基因表達譜數據進行差異分析，篩選差異表達基因（Differentially expressed genes，DEGs），并使用pheatmap軟件包、ggplot 2軟件包繪制DEGs聚類熱圖和火山圖。將篩選出的DEGs與鐵死亡相關基因取交集，獲得DEFRGs。

1.3 功能富集分析

使用R軟件clusterProfilter軟件包對篩選出的DEFRGs進行基因本體（Gene ontology，GO）富集分析和京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）富集分析，以P<0.05為差異有統計學意義。GO富集分析主要明確DEFRGs參與的生物學過程、細胞組分、分子功能。KEGG富集分析主要明確DEFRGs參與的細胞信號通路。

1.4 加權基因共表達網絡分析

使用加權基因共表達網絡分析（Weighted gene co-expression network analysis，WGCNA）軟件包對標準化處理轉錄組測序數據進行樣本聚類后去除離群樣本，根據計算出的最佳軟閾值構建無尺度網絡后進行基因聚類分析，形成不同顏色表示的基因模塊。通過動態剪切樹法識別基因模塊，對相似模塊進行合并，隨后對不同模塊與臨床性狀行相關性分析，篩選出與疾病特征相關性較高的模塊基因。將篩選出的模塊基因與DEFRGs取交集，獲得與疾病特征相關性較高的鐵死亡相關差異基因。

1.5 機器學習算法篩選疾病關鍵基因并驗證

對已得到的與疾病特征相關性較高的鐵死亡相關差異基因同時采用最小絕對收縮和選擇算子（Least absolute shrinkage and selection operator，LASSO）算法和支持向量機遞歸特征消除（Support vector machines-recursive feature elimination，SVMRFE）算法，將兩者結果取交集，最終篩選出3個鐵死亡關鍵基因。LASSO算法由R軟件中glmnet軟件包實現，SVM-RFE算法由e1071 軟件包實現。使用pROC軟件包繪制疾病關鍵基因受試者工作特征（Receiver operating characteristic，ROC）曲線，計算曲線下面積（Area under the curve，AUC），評估3個關鍵基因對疾病的診斷能力。進一步構建基于3個關鍵基因的邏輯回歸模型，通過ROC曲線評估邏輯回歸模型的診斷能力。使用GSE60436數據集對各關鍵基因的診斷能力進行驗證。

1.6 單基因基因集富集分析

以3 個關鍵基因在PDR組織樣本中基因表達量的中位數為標準，將疾病組樣本分為高、低表達組。使用clusterProfiler軟件包，選擇c2.cp.kegg.symbols.gmt作為參考基因集進行單基因基因集富集分析（Gene set enrichment analysis，GSEA），按照P<0.05、|NSE|>1、FPDR<0.25為閾值篩選并展示排名前6位的基因集富集結果。

1.7 構建基因-藥物調控網絡

藥物基因相互作用數據庫（D r u g-G e n e Interaction database，DGIdb）是一個提供基因與已知或潛在藥物相互作用關系的免費數據庫。在DGIdb中輸入篩選出的3個鐵死亡關鍵基因，獲得與關鍵基因相互作用的藥物信息，構建基因-藥物調控網絡。

1.8 免疫浸潤分析

使用preprocessCore、e1071、parallel軟件包，基于CIBERSORT算法分析常見的22 種免疫細胞在不同樣品中的浸潤比例。對正常組與疾病組免疫浸潤結果進行差異分析并使用vioplot軟件包繪制小提琴圖，以P<0.05 為篩選閾值。使用limma、reshape2、tidyverse、ggplot2軟件包進行鐵死亡關鍵基因與22 種免疫細胞的相關性分析，并繪制相關性熱圖。

1.9 統計學方法

所有統計分析結果由R v4.22軟件和Graphpad prism 8.0 軟件實現，計量資料的正態性檢驗采用Shapiro-Wilk檢驗，對符合正態分布的計量資料，組間比較采用獨立樣本t檢驗；對不符合正態分布的計量資料，組間比較采用Mann-WhitneyU檢驗。Pearson相關分析用于不同基因模塊與臨床性狀之間的相關性分析以及揭示3個鐵死亡關鍵基因與22種免疫細胞之間的關系。以P<0.05 為差異有統計學意義。

2 結果

2.1 DEFRGs篩選

從GSE102485數據集中共篩選出3 678個DEGs（|log2FC|>1且P<0.05），其中1 815個基因在PDR組中上調，1 863 個基因在PDR組中下調，差異基因的分布情況如圖1所示。從FerrDb數據庫中下載鐵死亡相關基因，包括369 個驅動基因、348 個抑制基因和11個標記基因，去除重復注釋基因，共484個鐵死亡相關基因。將484個鐵死亡相關基因與3 678個DEGs取交集，共獲得83個DEFRGs。

圖1.GSE102485數據集中PDR組織樣本與正常視網膜組織樣本行差異基因篩選結果Figure 1.Identification of DEGs between PDR samples and normal retina samples in GSE102485 dataset

2.2 功能富集分析

對篩選出的83個DEFRGs行GO和KEGG富集分析。GO富集分析結果顯示，83 個DEFRGs參與的生物學過程主要包括：細胞對化學應激的應答、細胞對氧化應激的應答、細胞對外部刺激的應答等。細胞組分方面，83個DEFRGs主要富集于囊泡、黏著斑等，見圖2A。KEGG富集分析顯示，83 個DEFRGs主要富集于鐵死亡、長壽調節通路、自噬、FOXO信號通路等，見圖2B。

圖2.差異鐵死亡基因功能富集分析柱狀圖Figure 2.Histograms of functional enrichment analysis of DEFRGs

2.3 WGCNA分析

使用W G C N A 軟件包對標準化處理后的GSE102485 轉錄組測序數據進行分析處理，使用pickSoftThreshold函數尋找合適的軟閾值，當軟閾值＝6 時，滿足無標度拓撲擬合指數R2=0.9，以構建共表達網絡。使用動態剪切樹法識別并合并相似的基因模塊，最終獲得16 個基因模塊，其中灰色模塊為無法聚類的基因模塊，見圖3A。計算各個基因模塊與臨床性狀之間的Pearson相關系數和P值，見圖3B，篩選出相關性較高的基因模塊：yellow（r=-0.89）、turquoise（r=0.85）、midnightblue（r=0.61）、black（r=0.60），共獲得疾病特征基因1 105個。將1 105 個疾病特征基因與83 個DEFRGs取交集，進而篩選出17個鐵死亡相關疾病特征基因。

2.4 機器學習算法篩選疾病鐵死亡關鍵基因并評估其診斷價值

基于篩選出的17 個鐵死亡相關疾病特征基因，使用LASSO算法，篩選出3 個鐵死亡關鍵基因：過氧化物酶體增殖物激活受體α（Peroxisome proliferator activated receptor alpha，PPARA）、ATP結合盒轉運蛋白C5（ATP binding cassette subfamily C member 5，ABCC5）、結節性硬化癥復合體亞單位1（TSC complex subunit 1，TSC1），見圖4A。同時使用SVM-RFE算法，從17 個鐵死亡相關疾病特征基因中篩選出13 個鐵死亡關鍵基因：TSC1、RBMS1、FTL、PPARA、CD44、ANO6、ABCC5、TMSB4X、TGFBR1、KDM4A、CTSB、SLC38A11、CIRBP，見圖4B。將兩者取交集，最終得到3個鐵死亡關鍵基因：PPARA、ABCC5、TSC1。分別繪制3 個基因的ROC曲線并計算AUC，結果顯示，3 個基因的AUC均超過0.95，見圖4C?；谏鲜?個鐵死亡基因構建邏輯回歸模型，繪制邏輯回歸模型ROC曲線并計算AUC，結果顯示，AUC=1.00，見圖4C。

圖4.綜合LASSO和SVM-RFE算法篩選鐵死亡關鍵基因Figure 4.Identification of ferroptosis-related hub genes by integrating LASSO algorithm and SVM-RFE algorithm

在獨立的驗證數據集GSE60436 中，基因PPARA、ABCC5、TSC1的表達量在正常組織樣本與疾病組織樣本中差異均具有統計學意義（Z=-2.32，P=0.024；Z=-2.12，P=0.036；Z=-2.24，P=0.025）。同時ROC曲線顯示，3個基因的AUC分別為0.917、0.972、0.944，見圖5。

圖5.外部數據集GSE60436驗證鐵死亡關鍵基因診斷效能Figure 5.Diagnostic efficiency of each ferroptosis-related hub gene in validation set GSE60436

2.5 單基因GSEA分析

為進一步探究篩選出的3個鐵死亡基因在PDR發生發展中可能參與的分子機制，我們基于KEGG數據集對3 個基因行GSEA分析。結果顯示，糖酵解和糖異生途徑、趨化因子信號通路、果糖和甘露糖代謝、肌萎縮側索硬化、近端小管對碳酸氫鹽重吸收途徑在PPARA高表達樣本中富集程度較高；嗅覺傳導通路在PPARA低表達樣本中富集程度較高。果糖和甘露糖代謝、淀粉與蔗糖的代謝、近端小管對碳酸氫鹽重吸收途徑、溶酶體、ABC轉運蛋白途徑在ABCC5高表達樣本中富集程度較高；嗅覺傳導通路在ABCC5低表達樣本中富集程度較高。溶酶體、黏著斑、致心律失常性右室心肌病、癌癥相關信號通路在TSC1高表達樣本中富集程度較高；神經配體-受體相互作用信號通路、嗅覺傳導通路在TSC1低表達樣本中富集程度較高。

2.6 疾病關鍵基因-藥物調控網絡構建

基于DGIdb 數據庫，挖掘與3 個鐵死亡基因潛在相關的靶向藥物，并分析藥物與基因之間的相互作用關系。通過檢索數據庫，篩選出10 個PPARA相關藥物，包括苯扎貝特、肉豆蔻酸、CHEMBL1089501、培馬貝特、氯馬扎利、CHEMBL107518、非諾貝特酸、LY-518674、GW590735、阿格列扎。其中培馬貝特、苯扎貝特、LY-518674、阿格列扎、非諾貝特酸、GW590735對PPARA有激動作用。對于ABCC5基因，共篩選出7個相關藥物：伊立替康、奧沙利鉑、脫氧氟尿苷、齊多夫定、氟尿嘧啶、甲酰四氫葉酸、格列本脲。對于TSC1基因，共篩選出3個相關性藥物：西羅莫司、阿司匹林、依維莫司。

2.7 免疫浸潤分析

免疫細胞浸潤差異分析結果顯示，M1 巨噬細胞、中性粒細胞、活化的記憶CD4+T細胞、γδ T細胞在疾病組織樣本中浸潤程度較高，差異具有統計學意義（t=2.62，P=0.016；t=3.10，P=0.005；t=2.45，P=0.023；t=2.28，P=0.034），見圖6A。22種免疫細胞與3個鐵死亡基因相關性分析結果顯示，ABCC5與調節性T細胞浸潤呈正相關（r=0.57，P=0.006），與單核細胞浸潤呈正相關（r=0.51，P=0.016），與γδ T細胞浸潤呈負相關（r=-0.45，P=0.035）；PPARA與單核細胞浸潤呈正相關（r=0.50，P=0.017）；TSC1與M2巨噬細胞浸潤呈正相關（r=0.51，P=0.014），見圖6B。

圖6.GSE102485數據集中PDR組織樣本與正常視網膜組織樣本免疫細胞浸潤分析Figure 6.Immune cell infiltration analysis between PDR samples and normal retina samples in GSE102485 dataset

3 討論

隨著人們生活方式、飲食結構的改變，糖尿病在全世界范圍內的發病率逐年增高，DR作為糖尿病常見的眼部并發癥，至2030 年，全球患病人數預計將達1.9億人[10]。PDR以新生血管為主要特征，是DR發展的中晚期階段，常伴有視網膜前出血、玻璃體積血、牽拉性視網膜脫離等嚴重并發癥[11]。因此，研究PDR發生的分子機制、鑒定新的特異性生物學標記物、篩選新的治療靶點是當前研究的重點。

鐵死亡是區別于凋亡、自噬、焦亡等的一種全新的調節性細胞死亡形式，以依賴鐵的細胞內脂質過氧化產物大量蓄積為主要特征，最早由Dixon于2012年提出[12]。在眼部疾病方面，已有研究表明，鐵死亡與視網膜色素變性、年齡相關性黃斑變性等密切相關[13-14]。有研究發現，高葡萄糖會刺激視網膜色素上皮細胞上調miR-338-3p表達，進而通過降解SLC1A5誘發氧化應激介導的鐵死亡通路[15]。此外，鐵死亡會強化高葡萄糖環境對視網膜毛細血管內皮細胞生長的抑制作用，這可能與GPX4的泛素化有關[16]?？傊?，鐵死亡與DR的病理機制關系密切，更多的潛在機制尚不明確。

本研究通過LASSO算法和SVM-RFE算法雙重篩選，共從17 個鐵死亡相關疾病特征基因中篩選出3個鐵死亡相關PDR關鍵基因：PPARA、ABCC5、TSC1。這3 個關鍵基因ROC曲線的AUC分別為：0.955、0.970、0.985，表明這3 個基因在區分疾病組織樣本與正常組織樣本方面具有很高的準確性。同時，基于3個關鍵基因構建的邏輯回歸模型同樣展現出很高的診斷效能（AUC=1.00）。上述結果表明，3個PDR鐵死亡關鍵基因可作為潛在的有預測價值的生物學標記物。

PPARA屬于核激素受體超家族的成員，是配體調節的轉錄因子。PPARA已被證實是脂代謝調控中的關鍵因子，在維持脂肪酸氧化分解、胰島素敏感性、血糖穩定性、細胞分化等方面發揮重要的調控作用[17]。當PPARA被激活后，其可通過上調肉堿棕櫚酰轉移酶1、?；o酶A氧化酶、脂肪?；o酶A合酶等基因的表達，增強脂肪酸β氧化的水平，發揮調節脂代謝異常的作用[18]。PPARA被發現可通過抑制小膠質細胞誘導的神經功能紊亂和神經血管復合體損傷在DR中發揮保護性作用[19]。在鐵死亡方面，有研究發現，在小鼠肝臟模型中，PPARA可以通過增加GPX4的表達水平發揮抑制鐵死亡的作用[20]。從DGIdb數據庫，我們篩選出10個PPARA相關性藥物，包含7個PPARA激動劑，其中培馬貝特，是一種新型的PPARA激動劑，在III期臨床試驗中被證明可以顯著降低甘油三酯，提高高密度脂蛋白水平，降低低密度脂蛋白水平，培馬貝特對胰島素抵抗也有改善作用[21]。有研究發現，在小鼠的DR模型中，培馬貝特在DR的早期階段對視網膜功能具有保護作用[22]。非諾貝特酸作為PPARA的經典激動劑，可顯著改善高甘油三酯血脂，降低心血管疾病的發生率[23]。一項非諾貝特酸的體外研究發現，非諾貝特酸可降低高糖環境中人視網膜色素上皮細胞中纖維蛋白連接水平和IV型膠原蛋白的過表達[24]。同時，在一項大樣本隨機對照研究中發現，與安慰劑組相比，使用非諾貝特酸治療的DR患者眼底激光治療需求下降37%[25]。有研究發現，ABCC5可通過穩定SLC7A11蛋白的方式減少細胞內脂質過氧化產物的蓄積，發揮對鐵死亡的抑制作用[26]。TSC1對鐵死亡的調控可能與Rheb1 的激活和線粒體穩態相關[27]。有研究表明，抑制TSC1介導的哺乳動物雷帕霉素靶蛋白復合物1信號通路可減緩糖尿病腎病中纖維化的發生，TSC1在上皮細胞的間質轉化中發揮重要作用[28]。DR與糖尿病腎病同為糖尿病常見的微血管并發癥，兩者在發病機制與病理變化中存在很多相近之處，TSC1在DR中發揮的作用有待進一步研究探索。

在比較PDR組織樣本與正常視網膜組織樣本的免疫細胞浸潤程度時，我們發現多種免疫細胞存在差異，且鐵死亡關鍵基因與某些免疫細胞相關。已有研究表明，在DR的發展過程中，伴隨著高糖及高脂質環境刺激，巨噬細胞、小膠質細胞等免疫細胞功能與代謝發生紊亂[29]。巨噬細胞在組織增生、修復、纖維化等生理過程中發揮重要的調節作用，其功能紊亂會導致大量炎癥因子、血管內皮生長因子釋放，參與纖維血管膜的形成[30]。目前，不同免疫細胞在PDR中發揮的作用尚未完全明確，更多的潛在細胞機制與免疫相關療法有待進一步研究。

總之，本次研究通過生物信息學分析方法結合機器學習算法挖掘PDR轉錄組測序數據中鐵死亡關鍵基因、潛在的治療藥物及免疫細胞浸潤分布情況，為進一步研究PDR的病理分子機制，尋找檢測指標及藥物治療靶點提供了新的方向。然而，本次研究缺乏進一步的實驗驗證，需要在后續的研究中改進完善。

利益沖突申明 本研究無任何利益沖突

作者貢獻聲明 陳鼎、徐徑舟：參與課題設計；數據分析；撰寫論文；根據編輯部的修改意見進行修改。黃晉、李智豪、張赟：參與數據收集；數據分析；參與修改論文中關鍵性結果、結論