基于生物信息学和机器学习鉴定2型糖尿病肾病肾小管间质损伤相关基因

宿家铭，彭 景，陈海敏，周 盈，史 扬，董昭熙，温雅轩，林子萱，柳红芳

（1.北京中医药大学，北京 100029；2.北京中医药大学东直门医院，北京 100700）

糖尿病肾病（diabetic nephropathy，DN）作为糖尿病最主要的微血管病变之一，已成为全世界终末期肾病的最常见原因［1］，严重影响患者的预期寿命和生存质量。目前尚无特效的治疗方案，仍以治疗原发病延缓病情进展为主［2］。因此，对DN 发病与进展机制的探索是肾脏病学和内分泌学者亟需解决的难题。传统观点认为DN 主要以肾小球病变为主，但随着对DN 肾小管间质氧化应激损伤、细胞凋亡和自噬、血管生成异常以及过度炎症反应活化等肾脏损害机制研究的深入，研究结果提示肾小管病变和肾间质纤维化在DN 的发生和发展中起重要作用，甚至认为其损伤早于肾小球病变，并可作为预测DN 进展的相对独立因素［3］。

当前生物信息学方法被广泛应用于研发与诊断或预后相关的生物标志物，但考虑样本量有限造成的高假阳性率，可能难以在单芯片数据分析中获得可靠的结果［4］。遂本研究利用生物信息学多芯片联合分析方法全面搜索GEO 数据库中所有现存DN 患者肾小管间质组织基因测序芯片，寻找与DN肾小管间质损伤特异性相关的差异表达基因（differentially expressed genes，DEGs），并尝试联合机器学习算法对核心基因进行判断筛选，同时研究相关基因可能存在的具体富集途径和免疫浸润机制，结合临床特征相关性分析评估危险因素，从而揭示涉及DN 肾小管间质损伤相关病理、生理过程的分子机制，探索潜在的生物标志物，为DN 的早期诊断和靶向治疗提供理论参考和科学依据。

1 资料与方法

1.1 芯片数据集筛选

在GEO 数据库（http：//www.ncbi.nlm.nih.gov/geo）中以“diabetic nephropathy”为关键词检索基因芯片数据，纳入排除标准如下：（1）人类mRNA 表达数据集；（2）所有样本为肾小管间质组织；（3）数据集应包含对照组和实验组，且各组样本量大于2。筛选后获取人类DN 样本基因表达谱公共数据集GSE30122、GSE47185 和GSE99340 作 为 联 合 芯 片数据集，而GSE104954 作为独立验证数据集，下载各矩阵数据集文件及相应的平台文件，详细信息见表1。

表1 数据集详细信息Tab 1 Dataset details

1.2 数据质量与差异分析

基于平台注释信息，将各数据集中探针ID 转换为“Entrez ID”。 合 并GSE30122、GSE47185 和GSE99340 矩阵文件为联合芯片数据集，运用R 软件（版本4.2.0）Sva 包中函数Batch normalization 对数据集进行批次校正后，通过“limma”包对标准化的基因表达谱进行差异分析，并使用贝叶斯方程多重检验校正，以|log2FC|＞1 和校正后P＜0.05 作为标准筛选出DEGs，利用ggplot2 包绘制DEGs 的火山图，pheatmap 包绘制DEGs 的热图。

1.3 DEGs 功能注释和富集分析

为了探索DEGs 的主要功能和途径，通过R 软件中“Bioconductor”软件包对DEGs 进行基因本论（gene ontology，GO）富集和KEGG（kyotoencyclopedia of genes and genomes）通路富集分析，以P＜0.05为标准筛选出具统计学意义的生物过程（biological process，BP）、细胞成分（cellular component，CC）、分子功能（molecular function，MF）以及信号通路。为了得到更全面的结果，采用基因集富集分析（gene set enrichment analysis，GSEA）方 法，选 择c5.go.v7.4.symbols.gmt 和c2.cp.kegg.v7.4.symbols.gmt 作为参考基因集，使用GSEA 软件对联合芯片数据集中DN 患者的校正基因表达矩阵分别进行1 000 次模拟分析，获取GO 和KEGG 富集分析结果。

1.4 免疫浸润分析

利用R 软件对以上获得联合芯片数据集的校正基因表达矩阵进行免疫浸润分析，从CIBERSOFT 官网（https：//cibersort.stanford.edu/）下载22种免疫细胞基因表达数据，使用R 软件“e1071”包计算每个样品中22 种免疫细胞占比，绘制免疫浸润丰度图。分别使用“corrplot”包、“vioplot”包绘制相关性热图和小提琴图，用以分析免疫细胞浸润分布相关性及其差异，P＜0.05 代表两组间差异有显著性意义。

1.5 机器学习筛选核心基因

机器学习算法可获得更精细的模型，已广泛应用于生物标志物的探索，本研究引入最小绝对值收敛和选择算子（LASSO）、支持向量机-递归特征消除（SVM-RFE）与随机森林（RF）3 种机器学习算法［5］。LASSO 由“glmnet”R 包构建，在拟合广义线性模型的同时进行变量筛选和复杂度调整；SVMRFE 是一种有监督的机器学习技术，可根据递归对特征进行排序［6］；RF 通过建立决策树分类器模型对分类变量进行反复迭代评分，产生高精确性分类特征，共同筛选联合芯片数据集中核心基因。

1.6 核心基因验证及诊断价值分析

将GSE104954 作为独立的验证数据集，采用非配对t检验，P＜0.05 为差异有统计学意义，验证筛选出的核心基因在两组间表达差异。随后，建立受试者工作曲线（ROC），计算ROC 曲线下面积（AUC）值，分别评估核心基因诊断DN 的效能。

1.7 预测模型的构建

整合核心基因在联合芯片数据集中表达矩阵，使用逻辑回归分析构建预测模型，利用R 软件可视化为列线图，预测DN 患者的肾小管间质损伤，采用ROC 曲线以辨别模型性能。

1.8 核心基因临床特征分析

Nephroseq 数 据 库（https：//www. nephroseq.org）是存储了肾脏疾病及其对照组的基因表达数据的临床数据库，被广泛用于肾脏病的研究［7］。本研究通过该数据库对筛选出的核心基因进行验证，采用Pearson 相关分析探索核心基因于DN 肾小管间质组织表达情况对DN 患者肾小球滤过率（GFR）、24 h 蛋白尿、血肌酐、尿素氮等临床指标的影响，P＜0.05 为差异有统计学意义。

2 结果

2.1 DN 肾小管间质组织相关DEGs 的识别

联合芯片数据集总计包含46 名DN 患者和37名健康对照者样本，批次校正后各数据集间的数据分布趋于一致，结果见图1。根据筛选标准与健康对照者相比，总计获得107 个DN 肾小管间质组织相关DEGs，其中26 个基因下调，81 个基因上调。为显示DEGs 的变化及聚类关系，绘制火山图和热图，见图2。

图1 芯片数据集的标准化Fig 1 Standardization of chip data set

图2 DEGs 的火山图（A）和热图（B）Fig 2 Volcano plots（A）and heatmap（B）of differentially expressed genes

2.2 富集分析结果

将筛选出的107 个DEGs 靶点导入R 软件进行GO 富集分析，其中生物学过程主要涉及到内肽酶活性的调节、肽酶活性的调节、水解酶活性的负调节、白细胞介导免疫和细胞因子产生的正调节等；细胞组分则包括了胶原蛋白-含有细胞外基质、分泌颗粒腔和细胞质泡腔等；分子功能主要富集在细胞外基质结构成分、酶抑制剂活性和肽酶调节活性等。KEGG 通路富集结果共34 条，主要相关通路为金黄色葡萄球菌感染、补体和凝血级联、细胞黏附分子、吞噬体、细胞外基原（extracellular matrix，ECM）-受体相互作用、Th1 和Th2 细胞分化等，见图3。进一步通过GSEA 富集分析发现，在DN 患者肾小管间质基因表达矩阵中，活跃的GO 功能主要富集于适应性免疫反应、淋巴细胞介导免疫、免疫效应器过程的调节等免疫相关过程；活跃的KEGG通路主要富集于细胞黏附分子、细胞因子-细胞因子受体相互作用、ECM-受体相互作用等相关通路，见图4。

图3 DEGs 的GO 分析（A）和KEGG 分析（B）Fig 3 GO analysis（A）and KEGG analysis（B）of differentially expressed genes

图4 DN 肾小管间质基因表达矩阵的GSEA 富集分析Fig 4 GSEA enrichment analysis of tubulointerstitial gene expression matrix in DN

2.3 免疫浸润分析

构建87 名健康对照者和71 名DN 患者肾小管间质组织的免疫细胞含量矩阵，显示了22 种免疫细胞在不同样本中浸润分布差异，见图5A。通过免疫细胞相关性分析发现，嗜酸性粒细胞与幼稚B 细胞之间相互作用最明显且呈正相关（r=0.66），见图5B。与健康对照者相比，DN 患者肾小管间质组织内22 种免疫细胞组间浸润显著（P＜0.05）的有5种，其中记忆性静息CD4 T 细胞、γδ T 细胞、静息肥大细胞和中性粒细胞上调，CD8 T 细胞下调，见图5。

图5 22 种免疫细胞浸润分析结果Fig 5 Infiltration of 22 kinds of immune cells

2.4 机器学习筛选核心基因

分别利用LASSO 回归、SVM-RFE 算法和RF算法对DN 肾小管间质组织DEGs 进一步筛选，其中构建LASSO 回归模型并进行交叉验证，误差最小值对应16 种特征基因，SVM-RFE 算法通过5 折交叉验证后挑选出8 种特征基因，RF 算法鉴定了10个特征基因。取交集得到MARCKSL1、CX3CR1、FSTL1、AGR2、GADD45B5 个核心基因，见图6。

图6 机器学习算法筛选DN 肾小管间质损伤核心基因Fig 6 The core genes of DN with tubulointerstitial injury were screened by machine learning algorithm

2.5 核心基因表达和诊断能力的验证

使用GSE104954 验证数据集进行外部交叉验证，结果显示：与健康对照者相比，DN 患者肾小管间 质 内MARCKSL1、CX3CR1和FSTL1基 因 表 达显著上调，而GADD45B基因表达显著下调，同联合芯片数据集内表达趋势一致，但AGR2 表达量差异无统计学意义，见图7A～E。同时，ROC 曲线显示筛选出的5 个核心基因在验证数据集（图7F～J）内对区分DN 患者与健康对照者具有较高的诊断效能（AUC＞0.7）。

图7 核心基因验证的结果Fig 7 The verification results of core genes

2.6 列线图预测模型

基于联合芯片数据集的核心基因表达矩阵，通过逻辑回归构建预测模型并可视化为列线图，预测模型的C 指数为0.994，具有较高的关联度。此外，ROC 曲线表明，与其他单一核心基因模型相比，组合列线图模型在预测DN 患者肾小管间质损伤方面的性能最高，见图8。

图8 列线图预测模型Fig 8 Model of prediction nomogram

2.7 核心基因临床特征分析结果

核心基因表达与DN 临床特征的相关性通过已有的肾病临床数据库（Nephroseq）进行验证，检索得到ERCB Nephrotic Syndrome TubInt（10 名DN 患者）、Ju CKD TubInt（17 名DN 患者）、Schmid Diabetes TubInt（9 名DN 患 者）、Woroniecka Diabetes TubInt（10 名DN 患者）等4 个DN 肾小管间质组织基因表达和临床特征数据集，Pearson 相关分析发现，基因MARCKSL1、CX3CR1和AGR2的表达与DN 患者GFR（mL/min）呈负相关，见图9A～C；基因MARCKSL1、FSTL1和GADD45B的 表 达 与DN 患者尿蛋白（g/d）水平呈正相关，见图9D～F。余结果无统计学意义，故未在图中列出。

图9 关键基因临床特征分析结果Fig 9 The clinical characteristics results of core genes

3 讨论

DN 肾小管间质损伤并非继发于肾小球病变，而在DN 早期存在且于疾病进展中发挥重要作用，所以相关生物标志物的研究成为DN 早期诊断的突破口之一［8，9］。本研究基于GEO 数据库中2022 年1月前的DN 患者肾小管间质组织基因测序芯片，利用生物信息学方法联合机器学习技术探索潜在分子机制。

在DN 患者与健康对照者样本之间共获得107个DEGs 与DN 肾小管间质病变特异性相关，综合DEGs 富集分析结果发现DN 肾小管间质损伤中免疫失调和炎性反应、细胞因子作用、ECM 沉积等机制尤为突出，深入探究发现各病理损伤之间存在紧密联系，并最终指向肾间质纤维化。肾间质纤维化是DN 进行性肾功能衰竭的重要原因，也是DN 进展至终末期肾病的最终结局［10］。ECM 过度沉积决定了肾间质纤维化的程度和进展［11］，DN 早期局部浸润的免疫炎性细胞释放出多种细胞因子、生长因子和血管活性物质造成肾脏固有细胞损伤，激活纤维蛋白原、纤连蛋白等多种ECM 蛋白转录并渗入受损部位，同时在细胞黏附分子的介导下将成纤维细胞和免疫细胞黏附于受损处以促进修复或清除病原体［12，13］。但随着免疫微炎症状态的长期浸润和高糖刺激持续存在，损伤后修复延伸为病理改变［14］，肾小管间质中ECM 过度沉积，且部分水解成具有生物活性的片段，刺激周围细胞转化为难以降解的胶原蛋白和纤连蛋白等纤维化ECM，导致正常肾脏组织结构和功能丧失，形成肾间质纤维化损伤［15，16］。遂基于本研究基因富集途径推测免疫失调和炎性反应为DN 肾小管间质损伤使动因素，而ECM 沉积造成的肾间质纤维化为最终结局。

进一步通过免疫浸润分析发现记忆性静息CD4 T 细 胞、γδ T 细 胞、静 息 肥 大 细 胞、中 性 粒 细胞、CD8 T 细胞等许多免疫系统成分参与DN 肾小管间质损伤，尽管DN 并非“免疫介导”为主的肾脏疾病，但大量研究证实先天免疫和适应性免疫异常造成的炎症反应失调导致DN 患者进行性肾功能损害［17］，包括多种免疫和炎性细胞（T 淋巴细胞、单核/巨噬细胞、中性粒细胞等）、炎性因子［血管内皮生长因子、肿瘤坏死因子-α（TNF-α）、转化生长因子-β 1、白介素（IL）、C 反应蛋白、NF-κB、结缔组织生长因子CTGF、单核细胞趋化蛋白MCP-1 等］通过多种信号通路和相互交叉作用，加剧机体微炎症状态与氧化应激反应，共同促进了ECM 沉积、细胞凋亡、肾小管硬化以及肾脏血流动力学改变，造成肾功能持续恶化，故近年来认为DN 是一种免疫炎症性疾病［18，19］。

为了实现DN 的早期诊断与评估预后，利用机器学习方法发现了DN 肾小管间质损伤新的分子特征，包括MARCKSL1、CX3CR1、FSTL1、AGR2和GADD45B，并通过ROC 曲线和列线图预测模型证实了出色的诊断效能。其中MARCKSL1在多种组织中表达［20］，参与调控细胞迁移、分泌、增殖和分化等多种生理活动，现有研究表明MARCKSL1在免疫系统中发挥重要的调节作用，包括促进炎症细胞的迁移以及细胞因子的分泌，可通过抑制其磷酸化进而抑制p38、JNK MAPKs 和NF-κB，降低TNF-α和IL-6 等炎性细胞因子水平，同时影响巨噬细胞、中性粒细胞的迁移和黏附［21］，但在DN 肾小管间质免疫浸润与微炎症状态中的功能和作用有待深入挖掘；多项临床研究已证实CX3CR1在DN 患者的肾脏中表达上调［22］，且Song 等［23］研究认为CX3CR1通过上调ECM 合成而在DN 中发挥重要作用，抑制CX3CR1可减少DN 小鼠模型的ECM 沉积，改善肾脏巨噬细胞浸润与纤维化，因此可能成为防治DN的有效靶点；FSTL1是一种成纤维细胞衍生的细胞因子，与多种组织器官（肾、肝、肺等）纤维化密切相关［24，25］，经临床和基础研究证实FSTL1在慢性肾脏病患者体内表达上调，且促进肾脏纤维化、炎症和细胞凋亡过程，可能为慢性肾脏病新的治疗靶点；AGR2可通过多种途径参与细胞增殖和生长，目前在肿瘤发生、发展及靶向治疗中的作用日益得到认可［26］，而在验证数据集中并未表现出显著差异，但Zhou 等［27］研 究 结 果 同 样 认 为AGR2为DN 肾 小 管间质病变关键枢纽基因之一；GADD45B参与细胞周期阻滞、细胞存活或凋亡及DNA 损伤修复等过程，少数研究报道高糖可刺激db/db 小鼠肾脏组织和人近端肾小管上皮细胞中GADD45B的高表达，促进肾小管上皮-间充质转化和细胞凋亡，但也指出GADD45B可能在其他细胞系和不同疾病模型中具有抗凋亡作用［28］，而本研究发现其在DN 患者肾小管间质组织中表达下调显著。

另一方面，DN 临床特征表现为GFR 恶化、进行性蛋白尿、血清肌酐和尿素氮水平升高等进行性肾功能损害，研究发现相对于肾小球病变，肾小管间质损伤对肾功能恶化的判断更加灵敏［29］。本研究基于核心基因在Nephroseq 数据库中进行临床特征相关性分析，发现MARCKSL1、CX3CR1和AGR2的表达与DN 患者GFR、尿蛋白水平具有相关性，提示对判断DN 患者病情进展具有一定意义。由此本研究结果充分表明MARCKSL1、CX3CR1、FSTL1在DN 肾小管间质损伤中具有良好的诊断与预测效能，甚至有足够潜力作为DN 治疗靶点，而AGR2和GADD45B还有待进一步深入研究。

综上所述，基于整合GEO 数据库中DN 肾小管间质组织差异表达的基因图谱，借助生物信息学联合机器学习方法，阐释了相关生物标志物在DN 肾小管间质损伤中的生物学意义，也为诊断DN 的生物标志物和治疗靶点提供新的思路，但相关检测方法以及具体应用的范围和准确性仍需进一步开展实验验证。

作者贡献说明：

宿家铭、彭景、董昭熙：统计分析数据并进行结果的分析与解释，负责撰写论文；宿家铭、陈海敏、史扬、董昭熙：共同完成论文相关数据资料的获取和处理；周盈、温雅轩、林子轩、彭景：协助完成数据分析和讨论；柳红芳：提供本文思路并指导写作。

所有作者声明不存在利益冲突关系。