鉴定铁死亡相关基因作为狼疮性肾炎诊断生物标志物的研究

吴昱升，林栩▲，王蓉，黄慕源，徐璐瑶，张洁，杨岚茵，蓝梦麟

(1.右江民族医学院附属医院肾内科，广西百色533000；2.广西免疫相关性疾病医学科研基础保障重点实验室，广西百色 533000)

系统性红斑狼疮(systemic lupus erythematosus，SLE)是一种涉及对内源性核颗粒的不适当免疫反应[1]，狼疮性肾炎(lupus nephritis，LN)是SLE最常见和最严重的靶器官表现之一[2]。10%～30%的LN患者在确诊后15年内发展为终末期肾病(end-stage renal disease，ESRD)，这是SLE死亡的主要原因[3]。最近的研究表明，LN易感基因(破坏免疫耐受)可以增强先天免疫信号通路，促进淋巴细胞活化，从而导致肾损伤[4]。铁死亡是一种新发现的程序性细胞死亡，其特点是产生脂质活性氧(reactive oxygen species，ROS)和铁超载，导致胱天蛋白酶和坏死体非依赖性细胞死亡。中性粒细胞铁死亡在狼疮发病机制中起到关键作用，铁死亡抑制剂治疗可以显著改善狼疮小鼠的疾病严重程度[5]。本研究以铁死亡相关基因(ferroptosis-related genes，FRGs)为目的基因筛选LN的生物标志物，旨在为LN的治疗提供新的策略。

1 材料与方法

1.1 数据来源从GEO数据库(Gene Expression Omnibus，https://www.ncbi.nlm.nih.gov/geo/query)下载LN患者的转录组数据和临床数据。GSE32591数据集包含93例肾小管间质样本和肾小球样本，提取肾小球样本为本研究的研究对象，其中包含14例正常样本和32例LN患者。GSE157293数据集包含3例正常样本和3例LN患者的长链非编码RNA(long noncoding RNA，lncRNA)和miRNA表达数据。此外，从FerrDB数据库(http://www.zhounan.org/ferrdb)中获得了349个FRGs，包括激活基因、抑制基因和标记基因。

1.2 方法

1.2.1 差异分析使用“limma”R软件筛选GSE32591数据集中的14例正常样本和32例LN样本之间的差异表达基因(differentially expressed genes，DEGs)，筛选标准为P<0.05, |Log2FC|>0.5，利用“ggplot”和“pheatmap”R包绘制DEGs的火山图和热图。

1.2.2 差异表达的LN相关基因(LN-related DEGs，LN-DEGs)的筛选本研究使用“WGCNA”R包进行加权基因共表达网络分析(weighted gene co-expression network analysis，WGCNA)，用good Samples Genes函数来过滤离群样本。为了确保基因之间的相互作用最大限度地符合无尺度分布，首先确定了软阈值。设置每个模块的最小基因数为100，通过动态剪切树算法合并模块，设置MEDissThres为0.2来合并相似的模块。最后分析LN与模块之间的相关性来寻找LN相关基因(LN-related genes，LNGs)，并将DEGs和LNGs相交得到LN-DEGs。

1.2.3 功能富集分析为了寻找LN-DEGs的功能和相关通路，使用“clusterProfiler”包进行GO和KEGG富集分析，P<0.05为显著富集。

1.2.4 蛋白质互作网络(protein-protein interaction,PPI)为了研究LN-DEGs之间的相互作用，使用STRING(https://string-db.org)网站构建蛋白质相互作用网络，并选择置信水平为0.4。

1.2.5 筛选诊断基因将上述LN-DEGs与349个FRGs取交集，得到狼疮性肾炎-铁死亡相关差异表达基因(LN-FR DEGs)。使用“glmnet”R包对LN-FR DEGs进行LASSO回归分析。同时，利用“e1071”R包对LN-FR DEGs进行支持向量机递归特征消除(support vector machine recursive feature elimination feature，SVM-RFE)分析，获得特征基因。将LASSO和SVM得到的特征基因取交集，并绘制受试者工作特征(ROC)曲线，得到诊断基因。

1.2.6 免疫细胞与诊断基因的相关性分析用“GSVA”R包计算LN和正常样品中23种免疫细胞的含量，并通过Wilcoxon秩和检验比较两组样品中的免疫细胞含量。并通过Pearson相关分析探索23种免疫细胞与诊断基因的关系。

1.2.7 CeRNA网络分析使用“DESeq2”R包筛选GSE157293数据集中的差异miRNA和差异lncRNA，筛选条件为P<0.05和|Log2FC|>0.5。然后我们使用miRWalk网站，设置置信水平为1，预测与上述交叉miRNA结合的lncRNA，与差异lncRNA取交集。用上述交叉miRNA，lncRNA和诊断基因构建ceRNA网络。

1.2.8 诊断基因表达量的验证为了验证诊断基因的表达量，分别在GSE32591和GSE157293数据集中比较了LN样本和正常样本的同源性磷酸酶-张力蛋白(phosphatase and tensin homolog,PTEN)和孤束核受体4A1(NR4A1)的表达。

2 结 果

2.1 差异分析本研究共筛选出了1221个DEGs，包括694个上调基因，527个下调基因。热图展示了DEGs的表达(图1)。

2.2 差异表达的LN相关基因的筛选样本和性状聚类树图显示不需要对基因进行过滤(图2)。当软阈值为5时，基因间的相互作用符合无标度分布(图3)。通过动态混合树剪切算法得到13个模块(图4)，合并后得到11个模块(图5)。LN与模块之间的相关性显示Blue模块与LN显著相关(图6)。Blue模块与LN相关性散点图如图7。1221个DEGs与Blue模块中的2278个LN相关基因相交，共得到628个LN-DEGs。

注：每个小方格代表每个样本，每行表示每个基因，表达量越高颜色越红，越少越蓝图1 差异基因表达热图 图2 样本和性状树形图

图3 无尺度软阈值的筛选

图4 聚类模块树形图

注：基因通过层次聚类被分为各种模块，不同的颜色代表不同的模块，其中灰色默认是无法归类于任何模块的基因 注：纵坐标为不同模块，横坐标为临床性状，每一个方块表示某模块和某性状的相关性系数图5 模块的识别与合并 图6 模块与临床性状相关性热图

注：横坐标表示Blue模块内的连通度，纵坐标表示临床性状图7 Blue模块与LN性状相关性散点图

2.3 功能富集分析628个LN-DEGs共富集到1619个GO条目(图8)，主要富集到免疫相关通路，包括白细胞黏附和信息、激活免疫反应、白细胞增生、免疫反应调节信号通路和淋巴细胞增殖等。KEGG富集结果显示富集到了98条KEGG通路，图9显示了15个KEGG条目，包括病毒性心肌炎、病毒感染、肺结核、金黄色葡萄球菌感染、甲型流感、破骨细胞分化、抗原加工和呈递等通路。

注：横轴表示GO词条包含的目标基因数，纵轴表示GO词条的名称，颜色表征-log10(P-value)图8 LN-DEGs的GO富集图

注：横轴表示KEGG词条包含的目标基因数，纵轴表示KEGG词条的名称，颜色表征-log10(P-value)图9 LN-DEGs的KEGG 富集图

2.4 PPI网络蛋白质通过相互作用构成网络来参与生物信号传递、基因表达调节及细胞周期调控等生命过程的各个环节。628个LN-DEGs构建的PPI网络如图10所示。有557个基因有相互作用，共有5028个相互作用关系。

图10 蛋白质相互作用网络

2.5 诊断基因的筛选628个LN-DEGs与349个FRGs相交共得到21个LN-FR DEGs。通过LASSO回归分析得到基因系数图和交叉验证误差图(图11)，lambdamin为0.0049时筛出7个特征基因。SVM-RFE模型的特征基因排名如表1所示，SVM准确率和泛化误差与特征数的关系图表明，当基因从1～21变化时，预测LN样本和正常样本最佳点的错误率为0.0167，精确率为0.983，共纳入了5个特征基因(图12)。然后将LASSO和SVM得到的特征基因取交集，得到关键基因PTEN和NR4A1，绘制ROC曲线验证诊断模型的预测能力，发现诊断模型具有良好的诊断能力(AUC=0.98)，PTEN(AUC=0.90)和NR4A1为诊断基因 (AUC=0.79)(图13)。

图11 LASSO回归分析筛选诊断基因

注：横坐标代表的是特征基因的个数，左图纵坐标表示5折交叉验证下的准确性，右纵坐标表示5折交叉验证下的泛化误差。折线图线的趋势代表特征基因个数与准确度和泛化误差的关系图12 SVM准确率和泛化误差与特征数的关系图

图13 LASSO、SVM筛选预后诊断基因的ROC曲线

表1 SVM-RFE模型特征基因排名

2.6 免疫细胞与诊断基因的相关性分析LN和正常样本中23个免疫细胞的小提琴图如图14所示。LN样本中eoslinophils、NKCD56bright、NK细胞、TFH和Th17细胞含量低，aDC、毒性细胞、iDC、巨噬细胞、Tcm和Tgd含量高。相关性结果表明PTEN与aDC、iDC、巨噬细胞、中性粒细胞、Tcm、Tgd、Th1细胞、Th2细胞呈正相关，与NK CD56dim细胞、NK细胞、TFH、Th17细胞呈负相关。NR4A1与嗜酸性粒细胞、NKCD56bright、NK CD56dim细胞、Th1细胞、Th17细胞呈正相关，与aDC、iDC、巨噬细胞、T辅助细胞、Th2细胞呈负相关(图15)。

注：“*”表示P-value<0.05；“**”表示P-value<0.01；“***”表示P-value<0.001；“****”表示P-value<0.0001图14 ssGSEA算法免疫细胞类型小提琴图展示

图15 ssGSEA算法免疫细胞类型在各亚型中的免疫评分柱状堆叠图展示

2.7 CeRNA网络分析本研究共筛选出36个差异miRNAs，包括12个上调的差异miRNAs和24个下调的差异miRNAs。共筛选出992个差异lncRNAs，包括412个上调的差异lncRNAs和580个下调的差异lncRNAs。将差异miRNA和lncRNA与miRWalk网站预测到的miRNA和lncRNA取交集得到4个miRNA(hsa-miR-183-5p，hsa-miR-129-5p，hsa-miR-1269b，hsa-miR-642b-5p)，3个lncRNA(MIR497HG，LINC01963，TBX2-AS1)，构建的ceRNA网络见图16。结果表明LINC01963可能通过hsa-miR-129-5p调节PTEN，TBX2-AS1和MIR497HG可能通过hsa-miR-642b-5p调节NR4A1。

注：绿色菱形表示miRNA，红色圆形表示mRNA，黄色矩形表示lncRNA图16 核心基因ceRNA网络构建

2.8 诊断基因的表达验证在GSE32591数据集中，PTEN在LN样本中高表达，NR4A1在LN样本中低表达(图17)。在GSE157293数据集中，PTEN和NR4A1的表达趋势与GSE32591数据集相同(图18)，说明获得的诊断基因具有可靠性。

注：此处统计方法为Wilcoxon.test，“****”表示表示P-value<0.0001，“**”表示P-value<0.01 注：此处统计方法为Wilcoxon.test，“****”P-value<0.0001，“**”表示P-value<0.01图17 诊断基因表达量展示 图18 诊断基因外部数据集的验证

3 讨 论

LN是SLE病人的常见并发症，目前对于LN的诊断和治疗标准还不完善[6]。铁死亡是新近发现的一种伴随着大量铁积累和脂质过氧化的细胞死亡形式。铁死亡与许多疾病的病理生理过程密切相关，包括肿瘤、神经系统疾病、缺血再灌注损伤、肾损伤和铁代谢疾病等[7]。研究LN中铁死亡所涉及的机制有助于为LN的治疗提供新的策略。

本实验通过机器学习得到了PTEN、NR4A1两个诊断基因。PTEN基因是一种磷酸酶和张力蛋白的同系物,能够编码一个同时具有双特异性蛋白和磷脂酸化酶功能的肿瘤抑制蛋白(PTEN蛋白)。PTEN蛋白广泛表达并介导黏附、迁移、细胞存活和凋亡等细胞过程[8]。有研究报告指出，PTEN缺乏在体内外以多种方式加重高血糖状态下的肾脏足细胞损伤，包括足细胞骨架重排、焦亡、细胞自噬和上皮细胞-间充质转化[9]。足细胞中PTEN的表达升高可能通过代偿性改善自噬以及抑制凋亡保护肾脏免受高血糖的影响，特异性敲除足细胞中PTEN基因可以引起尿白蛋白排泄增加，中度肾小球硬化等症状[10-12]。这些研究结果增强了PTEN成为LN新的治疗干预手段的可能性。

孤束核受体NR4A1(nuclear receptor subfamily 4，group A，member 1)，也称为TR3、Nur77m或NGF-IB，属于类固醇/甲状腺激素受体家族[13]，参与很多细胞活动，如葡萄糖和脂质代谢、凋亡和血管内稳态等[14]。NR4A1的表达可由多种炎症刺激快速诱导并通过单核细胞和巨噬细胞中的核NF-κB途径激活，活化的NR4A1反过来可以通过阻断p65与DNA的结合以及直接诱导其他NF-κB抑制剂表达来抑制NF-κB的活化。NR4A1在生理条件下是炎症反应的天然对应物，但在炎症和慢性疾病(如神经炎症和纤维化等)中，NR4A1常常表达为下调或失活[15]。NR4A1作为细胞凋亡传感器和组织稳态协调器，维持自身免疫耐受，有望成为自身免疫性疾病治疗中新的靶点。

本文GO条目主要涉及多种细胞生物学活动，如白细胞黏附和信息、激活免疫反应、白细胞增生、免疫反应调节信号通路、白细胞游走等。致病性自身抗体的产生、免疫复合物(IC)的沉积和补体级联的激活都是狼疮患者引发肾炎的原因[16]；巨噬细胞和淋巴细胞(主要是T细胞)的间质浸润，常常引起肾脏细胞损伤、间质纤维化和肾小管萎缩；细胞免疫在肾小球肾炎发病机制中也起到重要作用，白细胞从循环中迁移到周围组织在发炎的血管内皮上初始滚动、趋化因子对白细胞的激活、整合素及其配体相互作用介导的血管壁附着以及跨内皮迁移等一系列活动贯穿LN的发展[17]。在LN开始阶段，免疫沉淀物和自身抗体上调，从而引起炎性细胞因子和趋化因子表达及白细胞浸润和激活，活化的白细胞继而产生增强炎症反应的细胞因子，多种触发因素持续产生的细胞因子与LN的进展相关[18]。因此，从LN的起始阶段到进展阶段，细胞活动以及细胞因子是必不可少的。

KEGG富集结果表明，LN-DEGs主要富集于病毒性心肌炎、病毒感染、肺结核、金黄色葡萄球菌感染、甲型流感、破骨细胞分化、抗原加工和呈递、EB病毒感染、补体和凝血级联反应、自然杀伤细胞介导的细胞毒性等条目。获得性免疫缺陷的易感因素包括补体缺陷(尤其是C1q和C4)、细胞因子调节受损、T细胞增殖和B细胞功能改变，与外部因素(主要是感染)共同引发了SLE的发展[19-20]。有研究表明，狼疮性肾炎与EBV抗原潜伏膜蛋白(LMP)1之间存在关联，与健康对照组相比，SLE患者中EBV定向抗体的频率更高，滴度更高，表明EBV频繁被再激活[21]。EB病毒可能参与LN患者体内各种自身抗体的形成，并有望成为LN治疗的新策略。

ssGSEA结果显示LN样本中eoslinophils、NKCD56bright细胞、NK细胞、TFH和Th17细胞含量低，aDC、毒性细胞、iDC、巨噬细胞、Tcm和Tgd含量显著升高。肾小球中的免疫复合物沉积，以及T细胞、B细胞和髓样细胞(主要位于肾小球外)参与的炎症过程是LN的主要致病特征。有研究显示，与健康人相比，SLE患者外周血中NK细胞的比例和总数显著降低，细胞毒性降低，并最终引起免疫失调[22]，这与本文的免疫浸润分析结果一致。此外，也有实验结果证实了处于活动期SLE病人体内的IFN-α浓度显著超过非活动期SLE病人[23]，因此IFN-α是SLE发病机制中的一个关键细胞因子，重要的是，研究表明血清IFN-α水平与NK细胞中IFN-γ的产生直接相关。NK细胞产生IFN-γ的能力和IFN-α血清水平之间的直接相关性以及在先天免疫和适应性免疫之间的联系可能成为LN新型免疫治疗策略的发展。

ceRNA全称内源竞争RNA(competing endogenous RNA)，一般常用的分析结构有lncRNA-miRNA-mRNA分析或circRNA-miRNA-mRNA分析，其中miRNA处于调控的核心地位，当miRNA被lncRNA或circRNA这类ceRNA竞争结合时，受miRNA家族调控的mRNA转录水平会上升。本文ceRNA网络结果表明LINC01963可能通过hsa-miR-129-5p调节PTEN，TBX2-AS1和MIR497HG可能通过hsa-miR-642b-5p调节NR4A1，这在以往的文献中还没有报道过，可能成为LN治疗的新的靶向策略。

综上所述，本研究通过GEO数据库中LN患者转录组数据和临床数据，以及349个铁死亡(激活、抑制、标记基因)相关基因进行生物信息学分析，得到PTEN、NR4A1两个诊断基因，其表达量在外部数据集中的验证结果均有显著差异，这为LN的发病机制研究及治疗提供了新的思路。

利益冲突：所有作者都声明不存在利益冲突。