基于生物信息学和机器学习算法鉴定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。基于上述3个铁死亡基因构建逻辑回归模型，绘制逻辑回归模型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、酰基辅酶A氧化酶、脂肪酰基辅酶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的病理分子机制，寻找检测指标及药物治疗靶点提供了新的方向。然而，本次研究缺乏进一步的实验验证，需要在后续的研究中改进完善。

利益冲突申明 本研究无任何利益冲突

作者贡献声明 陈鼎、徐径舟：参与课题设计；数据分析；撰写论文；根据编辑部的修改意见进行修改。黄晋、李智豪、张赟：参与数据收集；数据分析；参与修改论文中关键性结果、结论