基于机器学习方法筛选IgA肾病铁死亡关键基因

杨柏新,付少杰,闫 冀,贾 冶,徐弘昭*

(1.武警吉林省总队医院 肾病科,吉林 长春130052;2.吉林大学第一医院 肾病科,吉林 长春130021;3.德惠市人民医院,吉林 德惠130300)

IgA肾病是最常见的原发性肾小球肾炎病理类型之一,病理生理过程复杂,预后多变,治疗也大多基于经验,临床上缺乏防止IgA肾病进展为终末期肾病的手段,给个人和社会带来沉重经济负担。铁死亡是新发现的由氧化还原反应驱动的细胞死亡程序,它构成了一种在遗传、形态和生化上都不同于凋亡和自噬的细胞死亡途径,包括脂质过氧化、氧化还原途径受损和铁过载。铁死亡处于新陈代谢、脂质过氧化和铁调节机制的交叉处,被定义为由脂质过氧化介导的铁依赖的细胞死亡形式,作为器官损伤的驱动因素存在于多种生物学系统中。铁死亡参与多种肾脏疾病,包括急性肾损伤、肾脏纤维化、多囊肾、糖尿病肾病等。本研究拟采用机器学习方法,探索IgA肾病铁死亡的关键基因,从铁死亡的角度为IgA肾病的防治提供新的靶点。

1 材料与方法

1.1 识别IgAN肾组织中差异表达的基因

从GEO数据库(http://www.ncbi.nlm.nih.gov/geo/)中下载IgAN患者和对照组肾小球组织的微阵列数据集(GSE93798,GSE104948),两组共包含47例IgA肾病及48例正常对照组肾小球组织的芯片数据,通过FerrDb数据库获取全部铁死亡相关基因数据,将研究中的铁死亡相关基因筛选出来,从而获得IgA肾病和对照组中铁死亡相关基因的表达情况,对铁死亡相关基因在IgA肾病和对照组进行差异分析,从而获得IgAN组和对照组中有差异表达的铁死亡相关基因(P<0.05)。

1.2 差异表达基因的富集分析

为了解获得的IgAN肾小球组织的差异表达基因分析(DEGs)可能与哪些生物过程(Biological Process,BP)、细胞成分(Cellular Component,CC)以及分子功能(Molecular Function,MF)相关,对这些DEGs进行了基因本体论(Gene Ontology,GO)功能富集分析。GO富集分析是通过R语言中的“clusterProfiler”包来实现,矫正后的P值<0.05被认为富集出的结果具有统计学意义。

1.3 筛选潜在的诊断性生物标志物

为了确定IgA肾病铁死亡诊断生物标志物,使用最小绝对收缩和选择算子(Least Absolute Shrinkage and Selection Operator,LASSO),支持向量机递归特征消除(Support Vector Machine-Recursive Feature Elimination,SVM-RFE),以及随机森林(Random Forests,RF)三种机器学习算法来预测样本的疾病状态。LASSO是一种回归分析算法,又称之为套索回归,它通过执行变量选择和正则化来提高预测的精度和结果统计模型的可解释性,本研究中LASSO回归是通过使用R语言中的“glmnet”包来实现的。SVM是一种监督下的机器学习算法,广泛被应用于回归和分类研究,RFE算法经常被应用于避免过拟合的出现,在本研究中SVM-RFE是通过R语言中的“e1071”包实现的,其中交叉验证的参数设置为5倍。RF是一种有监督的,以决策树为基学习器的集成学习算法,本研究使用R语言中的“randomForest”包来实现,将重要性得分大于2.0为标准对基因进行筛选。最后,将三种分类模型均鉴定出的基因挑选出来作为潜在的基因生物标志物进行进一步的分析。

1.4 评估鉴定出的生物标志物对IgAN的诊断效能

为了评估所鉴定的基因生物标志物对IgAN的诊断效能,利用训练集中的47个IgAN患者样本和48个对照组样本的mRNA微阵列数据,根据每个基因标志物的表达量分别绘制了它们的接受者操作特征(Receiver Operating Characteristic,ROC)曲线。然后将GSE99339数据集作为验证集,再次绘制每个鉴定的基因生物标志物的ROC曲线,以更进一步的评估它们的诊断功效。计算每个ROC曲线下的面积(Area Under Curve,AUC)以对诊断功效进行定量评估。ROC曲线是用R语言中的“pROC”包绘制的,双侧P值<0.05表示差异具有统计学意义。

1.5 鉴定出的生物标志物的基因集变异分析(Gene Set Variation Analysis,GSVA)和基因集富集分析(Gene Set Enrichment Analysis,GSEA)

为了进一步探究与鉴定出的生物标志物相关的生物学过程,根据鉴定出的生物标志物基因表达量的中位值,将IgAN样本分别划分为高表达组和低表达组,并进行GSVA和 GSEA。从分子特征数据库MSigDB中下载“c5.go.bp.v7.4.symbols”的基因集,GSVA和GSEA算法分别被用来分析数据集,以探究与标志物可能相关的生物学过程,以找到样品间有显著差异的生物过程通路。矫正后的P值<0.05被认为富集出的结果具有统计学意义。

1.6 免疫细胞浸润分析

CIBERSORT算法被用来分析IgAN和对照组肾组织中免疫细胞浸润的特征。CIBERSORT是目前最广泛使用的免疫细胞浸润算法,它使用线性支持向量回归的方法对组织表达矩阵进行解构,可以准确地量化每个样本的22种免疫细胞的丰度分数。使用R语言中的“ggplot2”包进行免疫细胞浸润矩阵的主成分分析(Principal Components Analysis,PCA)聚类研究,使用“pheatmap”包绘制了不同分析方法获得的免疫细胞浸润矩阵的热图,以显示不同免疫细胞在各个样本肾组织中的浸润情况。还使用“pheatmap”包以热图的形式绘制了识别出的浸润的免疫细胞之间的相关性。R语言中的 “vioplot”包被用来绘制小提琴图,以帮助将DKD患者和对照组样本肾组织之间免疫细胞浸润的差异可视化。

2 结果

2.1 识别IgAN肾组织中差异表达的基因

本研究的流程图见图1。对从两个GEO数据集(GSE93798,GSE104948)中获得的47个IgAN患者和48个对照组的肾小球组织的微阵列数据进行了回顾性分析,并通过FerrDb数据库获取全部铁死亡相关基因数据。在去除了批次效应后,最终共鉴定出了157个IgAN组中的DEGs;其中64个基因上调,93个基因下调(P<0.01,见图2)。IgAN组和对照组中有统计学意义的差异表达的铁死亡相关基因52个(P<0.0001,见图3、图4)。

图1 本研究流程示意图

图2 获得IgA肾病组和对照组中有差异表达的铁死亡相关基因157个(P<0.05)

图3 IgAN组和对照组中有特别显著差异表达的铁死亡相关基因52个(P<0.0001)

图4 IgAN组和对照组中有统计学意义的差异表达的铁死亡相关基因之间表达的相关性分析(P<0.0001)

2.2 差异表达基因的富集分析

通过对DEGs进行GO功能富集分析来探究IgAN可能的发病机制。GO功能富集分析确定出了DEGs最相关的10个生物学过程、细胞成分和分子功能(见图5)。

图5 IgAN与对照组间有差异表达的铁死亡基因的GO功能富集分析

2.3 IgAN诊断标志物的鉴定

三种不同的机器学习算法被用来鉴定IgAN的潜在诊断标志物。DEGs中共有12个基因被LASSO回归算法构建的诊断模型纳入进来(见图6A);共有3个基因被SVM-RFE算法构建的诊断模型纳入进来(见图6B);共有1个基因被随机森林构建的诊断模型纳入进来(图6C)。而只有锌指蛋白36(ZFP36)是被3种算法都纳入诊断模型之中(见图6D),因此ZFP36被鉴定为IgAN的潜在诊断标志物。

图6 通过机器学习方法筛选关键铁死亡基因(A:LASSO回归;B:SVM-RFE;C:随机森林模型;D:3种方法取交集)

2.4 IgAN诊断标志物诊断效能的验证

如图7A所示,在训练集中这两个基因的表达能够有效地将IgAN和对照组样本区分开来,ZFP36的AUC为0.874(95%CI 0.795～0.945)。为了进一步评估这两个生物标志物的诊断效能,GSE99339数据集为验证集,其中包括11个IgAN患者样本和26个对照组样本。如图7B所示,验证集中也绘制了ZFP36的ROC曲线,ZFP36的AUC为0.874(95%CI 0.793～0.940)。这些结果表明,ZFP36的表达在诊断IgAN方面显示出良好的潜力。

图7 关键铁死亡基因的诊断效能评价及验证(A:以ZFP36的表达量为参数在训练集中构建诊断模型;B:以ZFP36的表达量为参数在验证集中构建诊断模型)

2.5 生物标志物的GSVA和GSEA分析

利用GSVA和GSEA探索与ZFP36相关的生物学过程的结果见图8;GSEA分析结果如图8A所示,ZFP36参与细胞死亡、DNA结合转录激活因子活性等生物过程。ZFP36参与调节趋化因子信号通路,细胞因子-细胞因子受体相互作用,利什曼原虫感染,MAPK信号通路和NOD样受体信号通路(图8B)。GSVA分析结果显示,ZFP36负调控成纤维细胞生长因子,Toll样受体4结合,细胞肝素反应以及转录因子AP1等生物过程(图8C)。ZFP36下调哺乳动物昼夜节律,NOD样受体信号通路等,且上调牛磺酸和次牛磺酸代谢(图8D)。

图8 GSEA和GSVA分析探究关键铁死亡基因相关的生物过程和信号通路

2.6 免疫细胞浸润的分析

基于CIBERSORT算法探究IgAN和对照组肾组织间差异浸润的免疫细胞,见图8。在CIBERSORT算法中,不同浸润的免疫细胞的小提琴图显示,激活的NK细胞、静止的树突细胞在IgAN样本中的浸润数量显著大于对照组样本,而静止的CD4记忆T细胞,静止的自然杀伤(Natural Killer,NK)细胞和中性粒细胞的浸润数量则较少(见图9A)。22种免疫细胞亚型浸润的相关性热图显示,静止树突细胞浸润正相关,而中性粒细胞、单核细胞和静止的NK细胞这三种免疫细胞的浸润负相关(见图9B)。ZFP36表达与中性粒细胞、静止的NK细胞浸润呈正相关,与静止的树突细胞、激活的NK细胞、巨噬细胞M1浸润呈负相关(见图9C)。

图9 IgA肾病和对照组肾小球组织中差异浸润的免疫细胞情况

3 讨论

目前IgAN的病理机制还未完全明确,但是被广泛接受的理论认为,IgA肾病是多种机制共同导致的,被称为“多重打击假说”,包括4个环节:①异常糖基化IgA1,进而生成半乳糖缺乏IgA1(galactose-deficient IgA1,Gd-IgA1)[1];②被循环系统中抗聚糖的自身抗体识别[2];③导致肾炎的免疫复合物沉积在系膜区[3];④激活系膜细胞,导致补体激活、进一步的免疫反应和肾损伤。IgA肾病患者的循环系统中,异常糖基化IgA1水平较高。异常糖基化IgA1与基因缺陷或遗传特征有关,但不足以导致疾病的发生;IgA肾病患者有血缘关系的人,也可能出现异常糖基化IgA1水平升高,但是他们不一定都会罹患IgA肾病[4]。因此,异常糖基化IgA1如何导致IgA肾病仍然值得研究。异常糖基化IgA1可能会产生Gd-IgA1。高水平的Gd-IgA1一般认为是由扁桃体淋巴细胞、骨髓和肠黏膜的浆细胞产生[5-6]。Gd-IgA1的大量合成是IgA肾病中IgA1免疫复合物形成的基础,而IgA1免疫复合物是IgA肾病发病的关键因素。在IgA肾病患者中,可检测到这些免疫复合物。虽然单独的IgA通常不能激活补体,但这些复合物却可通过自聚化形成聚合体,并被补体特异性识别、结合,形成循环免疫复合物,在肾小球系膜区异常沉积,最终导致IgA肾病[7]。而最近的研究表明,黏膜和补体免疫系统在IgA肾病中发挥重要作用。血尿的出现和黏膜免疫相关,而循环免疫复合物中常常可以发现补体,例如C3补体和补体因子H[8]。IgAN是一种慢性炎症性疾病,其发病和进展与许多铁代谢指标相关[9-10],此外,IgAN以血尿和蛋白尿为主要临床特征,肾小管上皮细胞可以吞噬和降解红细胞产生二价铁,与H2O2相互作用产生活性氧(ROS)导致肾脏细胞损伤[11-12]。

肾脏影响铁稳态的调节,而铁失衡会加重肾损伤[13-15]。铁死亡是2012年发现的一种独特的内源性修复和保护系统控制的古老细胞死亡形式[16]。铁死亡作为器官损伤的驱动因素存在,与多种生物学过程密切相关,包括发育、衰老、免疫与癌症等。细胞代谢机制、ROS的调控和铁的调节在铁死亡中交叉错集,构成了铁死亡的机制和生物学意义。铁死亡的并不是一般类型的氧化应激,而是发生于细胞膜的致命脂质过氧化物积累。PUFA活化并将其掺入膜脂(如PLs)中,以便产生致命性的过氧化物,最终驱动铁死亡。鉴于ACSL4参与活化和整合 PUFA,如花生四烯酸(AA)进入膜定位的脂质,说明PUFA需要存在于其膜结合环境中以表现过氧化后的杀伤力。ACSL4 失活是一种在不同情况下抑制铁死亡的关键机制,ACSL4的过表达对铁死亡敏感[17]。此外,ACSL4还积极参与正反馈循环以执行铁死亡[18]。因此,ACSL4可能更类似于细胞凋亡的caspase-3,它是细胞死亡的执行者,而不是管家蛋白。GPX4和谷胱甘肽抑制脂质ROS的积累,是一种铁死亡的中枢抑制因子[19]。GPX4的降解,而非抑制可以促进铁死亡。铁死亡抑制蛋白1(FSP1)/(A:与ZFP36相关的生物过程的GSEA分析;B:与ZFP36相关的信号通路的GSEA分析;C:与ZFP36相关的生物过程的GSVA分析;D:与ZFP36相关的信号通路的GSVA分析)CoQ10、二氢乳清酸脱氢酶(DHODH) 和环GTP水解酶1(GCH1)/四氢生物蝶呤(BH4),通过独立于GPX4的方式来抑制铁死亡[20-22]。顾名思义,铁死亡依赖于铁。非酶催化的、铁依赖的Fenton链式反应可能对于铁死亡是必要的:PLOOHs在GPX4被抑制时会持续存在更长时间,从而起始Fenton反应,使铁死亡标志物PLOOHs的含量迅速增加,PLOOHs能够与亚铁离子和铁离子反应,分别产生自由基PLO·和PLOO·引发破坏性的过氧化链式反应[23]。此外,LOXs和POR发挥催化活性需要铁的参与,铁对于许多与细胞内ROS产生相关的基于氧化还原的代谢过程也是必不可少的。由于铁在控制细胞生存和死亡中占据中心地位,细胞内的铁稳态受到精细调控,许多胞内通路通过改变细胞内不稳定铁的含量改变细胞对于铁死亡的敏感性[24]。IgAN时铁代谢和转运、ROS的产生和抗氧化防御系统也受到影响。转铁蛋白受体(transferrin receptor,TfR)介导的铁输入和自噬介导的铁蛋白(ferritin)降解等过程都促进铁死亡[25]。一项队列研究显示与健康对照组相比,IgAN患者中Gd-IgA1更倾向于与TfR结合,表明免疫复合物的形成有利于系膜 TfR-IgA1相互作用,同时可以通过正反馈增加TfR的表达[10,26]。此外,TfR与IgAN疾病进展显著相关并可能在IgAN中发挥生物学作用[27]。在一项纳入514例IgAN患者的回顾性研究发现尿转铁蛋白(Transferrin,Tf)与IgAN系膜细胞增生、内皮细胞增生以及间质纤维化相关。与对照组相比,IgAN患者有更低的SOD活性和更高的MDA水平[28]。GPX4表达在IgAN大鼠肾脏组织中表达明显降低[29]。鉴于大量的研究表明IgAN与铁死亡有关,可以通过调控铁死亡来治疗IgAN。通过机器学习,发现ZFP36可作为IgAN铁死亡新的诊断标记物,对诊断IgAN表现出良好潜力。

ZFP36转录后调控参与了许多炎症基因的调控,而IgAN时循环免疫复合物在系膜区的沉积可以刺激细胞外基质聚积、促炎性和促增殖因子的释放,从而导致系膜细胞增殖和促炎性细胞进入肾小球。此外,促炎性介质同时改变足细胞、肾小管上皮细胞和内皮细胞的基因表达,从而导致细胞损伤,说明IgAN是一种慢性炎症性疾病[30-32]。ZFP36家族蛋白是参与信使RNA(mRNA)代谢途径的RNA结合蛋白。ZFP36家族由ZFP36(也称为tristetraprolin,TTP)、ZFP36L1、ZFP36L2和ZFP36L3(仅在啮齿动物中)组成。ZFP36家族蛋白包含两个串联重复的CCCH型锌指基序,结合到特定mRNA的3′-非翻译区(3′UTR)中富含腺嘌呤尿苷的元件,并导致目标mRNA的衰变。虽然ZFP36家族成员在结构上相似,但已知它们发挥不同的功能并调节不同的靶mRNA,这可能是由于它们的细胞类型特异性表达模式所致。例如,ZFP36在全身炎症性疾病的小鼠模型中,通过下调各种促炎细胞因子(包括TNF-α)的产生,发挥抗炎调节剂的作用。

目前,免疫浸润细胞被认为在IgAN的发生和发展中起着重要作用。因此,揭示IgAN肾脏组织中免疫细胞的浸润情况,并探讨所发现的新的生物标志物与免疫浸润细胞之间的关系,具有重要意义。ZFP36表达与中性粒细胞、静止的NK细胞浸润呈正相关,与静止的树突细胞,激活的NK细胞,巨噬细胞M1浸润呈负相关。研究表明IgAN患者的树突细胞诱导IgA产生的能力受损[33]。NK淋巴细胞亚群在IgAN中的重新分布可能在免疫调节网络的损伤中起重要作用,可引起血尿[34-35]。中性粒细胞可引起IgAN肾小球损伤[36]。尽管本研究发现IgAN部分免疫细胞的失调,但相关免疫细胞的详细功能与IgAN的确切联系仍不清楚。

综上,铁死亡与IgAN发生、发展密切相关。ZFP36作为IgAN铁死亡的诊断标志物,将为IgAN的诊断、治疗及预防提供新的靶点。