基于WGCNA分析鉴定糖尿病肾病新的诊断标志物

安 敏,杨晓玲,马 亮,付少杰,马福哲*

(1.唐山职业技术学院附属医院 肾内科,河北 唐山063000;2.唐山职业技术学院附属医院 内分泌科,河北 唐山063000;3.吉林大学第一医院 肾病内科,吉林 长春130012)

随着人们生活方式以及饮食习惯的改变,糖尿病的发病率逐年增高,而糖尿病肾病(diabetic kidney disease,DKD)作为糖尿病的主要微血管并发症,也已成为终末期肾脏病(end-stage renal disease,ESRD)的首要原因[1]。ESRD患者需要进行肾移植或长期的肾脏替代治疗,生活质量严重下降[2]。DKD早期症状隐匿,临床上易被忽视,疾病被发现时常已进展至中晚期,这是导致DKD患者预后不佳的原因之一。目前对DKD的治疗尚缺乏特异性的手段,主要是通过积极控制血压、血糖以及改善生活方式等对症治疗,且疾病一旦进展至中晚期则不可逆转[3]。因此早期诊断和及时干预对延缓DKD患者的肾功能衰竭以及ESRD的发生具有重要意义。

加权基因共表达网络分析(weighted correlation network analysis,WGCNA) 是一种将基因数据转换为共表达模块,从而探索基因与给定特征之间相关性的生物信息学方法[4]。近年来WGCNA分析已被广泛应用于对各种疾病生物标志物的鉴定,并展现出良好的潜力[4-5],但在DKD新型标志物鉴定方面的研究尚十分缺乏。因此本研究旨在通过对基因综合表达(gene expression omnibus,GEO)数据库进行挖掘,鉴定出DKD的差异表达基因(differentially expressed genes,DEGs),构建加权基因共表达网络,筛选出DKD的诊断标志物,并探究标志物与免疫细胞浸润之间的关系,以期为DKD的临床诊疗以及药物研发提供新的线索。

1 材料与方法

1.1 数据来源及DEGs的获取

从GEO数据库[6](http://www.ncbi.nlm.nih.gov/geo/)获取包含有DKD患者以及健康对照的肾小球微阵列数据集(GSE30528,GSE96804,GSE99339)作为训练集,其中共有DKD患者64例和健康对照44例。获取共包含有DKD患者33例和健康对照67例的数据集(GSE47185,GSE30122)作为验证集。使用R语言SVA包对各数据集进行批次校正后合并。使用R语言中的“limma”包对训练集中DKD患者以及健康对照肾小球中各基因mRNA的表达水平进行差异分析,以|log2差异倍数|>1以及矫正后的P值<0.05为筛选标准,从而鉴定出DEGs。

1.2 WGCNA网络的构建

使用R语言中的“WGCNA”包构建WGCNA网络。首先对训练集的基因表达数据进行聚类分析,根据基因间的无标度拓扑拟合指数(R2)以及平均连接度,确定软阈值β为3,使各基因调控关系符合无尺度分布。通过构建共表达网络绘制模块聚类树状图,并将所有基因分为10个模块,使用Pearson法分析基因间的相关性,获得关系矩阵,选择相关性最高的模块作为关键基因模块。

1.3 关键DEGs的筛选及富集分析

将WGCNA分析获得的关键模块基因与DEGs取交集,筛选出重叠基因作为关键DEGs,它们可能与DKD的发生发展高度相关。使用R语言 “clusterProfiler”包对关键DEGs分别进行基因本体论功能(Gene Ontology,GO)及基因组百科全书通路(Kyoto Encylopaedia of Genes and Genomes,KEGG)富集分析,矫正后的P值<0.05被认为具有统计意义。

1.4 构建关键DEGs的蛋白互作网络

将关键DEGs导入STRING数据库进行蛋白互作(protein-protein interaction,PPI)分析,并构建PPI网络。将PPI信息导入Cytoscape3.9.1内,计算网络的拓扑学参数,并利用cytoHubba插件进一步筛选出核心基因。

1.5 受试者操作特征曲线以及列线图的绘制

利用训练集中DKD和对照组样本的微阵列数据,根据每个核心基因的表达量分别绘制受试者操作特征(Receiver Operating Characteristic,ROC)曲线。为便于临床应用,将核心基因表达量纳入多因素Logistic回归分析,构建诊断模型,并绘制模型的校准曲线以及决策曲线。将训练集中ROC曲线下的面积(Area Under Curve,AUC)大于0.9的基因筛选出来作为DKD诊断潜在的生物标志物,并使用验证集数据再次绘制ROC曲线对其诊断效能进行再次验证。ROC曲线、校准曲线以及决策曲线分别是由R语言中的“pROC”,“rms”以及“rmda”包绘制。

1.6 生物标志物的基因集富集分析

为了进一步探究与生物标志物可能的相关机制,将DKD样本根据标志物基因表达量中位值划分为高低表达两组,进行基因集富集分析(Gene Set Enrichment Analysis,GSEA)。从分子特征数据库MSigDB中下载“c2.cp.kegg.symbols.gmt”的基因集,利用 GSEA算法来进行分析,矫正后的P值<0.05被认为富集出的结果具有统计学意义。

1.7 免疫浸润分析

采用CIBERSORT算法对DKD以及对照组肾组织中22种免疫细胞浸润的特征进行分析。使用R语言中的“ggplot2”包对免疫细胞浸润矩阵的主成分分析进行聚类,使用 R语言中的“vioplot”绘制小提琴图,将DKD患者和对照组样本肾组织间免疫细胞浸润的差异可视化。R语言中的Spearman等级相关分析被用来进一步分析生物标志物的表达与免疫细胞浸润间的关系,并使用R语言中的“ggstatsplot”和“ggplot2”包将生物标志物的表达与免疫细胞浸润间的关系可视化。

2 结果

2.1 DEGs的获取

从训练集的64例DKD患者肾小球组织以及44例健康对照肾小球组织中共筛选出87个DEGs,其中在DKD肾小球中上调表达的DEGs有39个,下调表达的有48个,见图1。

图1 DKD患者肾小球组织中的DEGs的鉴定

2.2 WGCNA网络的构建

使用R语言中的“WGCNA”包对训练集的基因表达数据进行聚类分析,选择3作为软阈值β,以符合无尺度网络分布,见图2A。利用构建共表达网络绘制模块聚类树状图,表达数据相似的基因被聚类在一起形成一个分支,见图2B。使用Pearson法分析基因模块与临床特性的相关性,绘制模块-临床特性关联图,结果发现粉橙色基因模块相关系数最大为0.67,且P值最小为2E-15,见图2C。提示该模块聚类的基因与DKD的发生发展最为密切,因此挑选粉橙色基因模块为关键基因模块。

图2 WGCNA聚类分析

2.3 关键DEGs的筛选及富集分析

将WGCNA分析获得的关键模块基因与DEGs取交集,筛选出14个重叠基因作为关键DEGs,可能与DKD的发生发展高度相关,见图3A。关键DEGs的GO富集分析结果提示,关键DEGs主要富集于白细胞趋化、趋化因子结合等条目,见图3B;KEGG富集分析结果提示,主要富集于白细胞介素17 信号通路、肿瘤坏死因子信号通路等通路,见图3C。

图3 关键DEGs的筛选及GO、KEGG富集分析

2.4 构建关键DEGs的PPI网络

将关键DEGs导入STRING数据库,进行PPI分析,构建出1个含有14个节点,48条边的PPI网络,见图4A。将网络信息导入Cytoscape3.9.1内,计算拓扑学参数,并借助cytoHubba插件使用“MCC”法对网络进一步筛选,获得10个打分最高基因,见图4B,将打分最高的5个基因筛选为核心基因,它们是PTGS2,FOS,FOSB,EGR1,DUSP1。

图4 构建关键差异表达基因的PPI网络

2.5 诊断标志物的确定

通过对训练集绘制ROC曲线对核心基因进行分析,结果发现5个核心基因的AUC均大于0.8,而DUSP1和FOSB的AUC大于0.9,提示在区分DKD以及正常样本方面均表现出良好的诊断能力,见图5A。随后,将5个核心基因表达量纳入多因素Logistic线性回归构建诊断模型,并绘制模型的决策曲线以及校准曲线,分别见图5B和图5C。对校准曲线进行Hosmer-Lemeshow拟合优度检验,结果表明理想曲线与校准曲线的差异不具有统计学意义(χ2=5.371,P=0.719),提示模型具有较高的校准度,其预测概率与实际接近。由列线图可知DUSP1以及FOSB的基因表达水平对是否发生DKD的影响最大,因此将DUSP1和FOSB筛选为DKD的潜在生物标志物,并在验证集中对它们的表达情况进行进一步验证,并分别绘制它们的ROC曲线,见图6。在验证集中DUSP1和FOSB的AUC也均大于0.75,进一步验证了作为DKD诊断标志物的诊断潜力。

图5 核心基因在训练集中的ROC曲线及诊断列线图

图6 诊断标志物在验证集中的差异表达情况及ROC曲线

2.6 生物标志物的GSEA富集分析

为了进一步探究与新鉴定出的生物标志物可能相关的机制,根据标志物基因表达量的中位值,将DKD样本划分为高表达和低表达组,并进行GSEA富集分析。结果发现随着DUSP1表达的增加,β丙氨酸代谢、脂肪酸代谢、色氨酸代谢等通路被抑制,见图7A。而随着FOSB表达的增加,色氨酸代谢、脂肪酸代谢、丁酸代谢等通路被抑制,见图7B。

图7 两个诊断标志物的GSEA富集分析

2.7 免疫浸润分析

免疫浸润结果发现,相较正常对照组,M2巨噬细胞、γδT细胞、记忆B细胞、静止的树突状细胞等在DKD肾小球组织中浸润显著增加,而辅助性T细胞、中性粒细胞以及活化的肥大细胞在DKD肾小球组织中浸润显著减少,见图8A。进一步对DUSP1以及FOSB两个标志物的表达与肾组织中免疫细胞浸润的关系进行相关性分析,结果发现DUSP1的表达与浆细胞、γδT细胞以及活化的CD4阳性的记忆T细胞的浸润呈负相关,见图8B;FOSB的表达与M1巨噬细胞浸润呈负相关,见图8C。

图8 免疫浸润分析及诊断标志物与免疫细胞浸润相关性分析

3 讨论

研究显示,我国糖尿病患者已高达1.14亿,且糖尿病导致的肾脏损害已超过慢性肾小球肾炎成为我国慢性肾脏病的首要病因[7]。目前DKD的早期诊断主要依赖于对尿蛋白的检测,但DKD的早期阶段的尿蛋白可以是阴性的[8]。因此发掘新的具有良好诊断效能的标志物,有助于DKD的早期诊断并改善患者的预后。本研究利用WGCNA分析确定了DUSP1和FOSB为DKD的潜在诊断标志物,并联合GSEA以及免疫浸润分析,进一步了解与它们表达相关的可能机制。结果发现DUSP1和FOSB对DKD具有良好的诊断效能,并且可能成为治疗DKD的潜在靶点,为DKD的诊断及治疗提供了新思路。

DUSP1是一种组蛋白磷酸酶,可以特异性地调节丝裂原活化蛋白激酶(Mitogen-Activated Protein Kinase,MAPK)的活性,从而发挥多种生物学功能[9]。本研究发现DUSP1在DKD肾小球组织中是低表达的。有研究报道lncRNA NR_038323可以通过靶向DUSP1抑制p38MAPK信号通路,在HK-2近端小管细胞系中发挥抗纤维化作用[10]。同时,DUSP1过表达可通过调节线粒体自噬减轻肾小管损伤,阻断线粒体裂变因子相关的线粒体过度裂变[11]。此外,在T细胞中抑制DUSP1的表达可能有利于多发性硬化自身免疫性疾病的治疗[12]。提示DUSP1有潜力成为包括DKD在内一些疾病的治疗靶点。

FOSB是FOS家族的一员,属于原癌基因的一种,其与Jun蛋白家族组成异源二聚体转录因子复合物AP-1,在肿瘤的进展中发挥着重要作用[13]。FOSB的表达增加可以抑制细胞的凋亡并促进细胞的增殖[14]。既往有研究表明,FOSB可作为IgA肾病的诊断标志物[15]。也有研究表明,MicroRNA-27a-3p可以靶向FOSB从而调节IgA肾病的炎症以及纤维化水平[16],但其在DKD中的相关研究尚少见。

近年来,越来越多的研究发现免疫紊乱在DKD的发病中发挥重要作用[17]。本研究也发现,在DKD肾小球中存在免疫细胞的异常浸润,且 DUSP1的表达与浆细胞、γδT细胞以及活化的CD4阳性的记忆T细胞的浸润呈负相关,FOSB的表达与M1巨噬细胞的浸润呈负相关。提示DUSP1和FOSB可能参与对肾脏免疫细胞浸润的调控,从而对DKD的发生及发展产生影响。但是,这一推断需进行更多的实验来进行验证。

综上所述,本研究发现DUSP1和FOSB对DKD具有良好的诊断效能,并且与DKD肾小球免疫细胞的浸润相关,有望成为DKD的新型诊断标志物以及治疗靶点。但本研究为生物信息学研究,具有一定局限性,需要更进一步的临床验证。