The DKD-related datasets GSE30528(11), GSE30529(11), GSE30122 (11,12) and GSE96804 (13,14) were downloaded from the GEO database via the R package GEOquery (15), in which the datasets GSE30528, GSE30529 and GSE30122 served as the test set and the dataset GSE96804 served as the validation set. All DKD samples and control samples were included in the present study. The characteristics of these datasets are listed in Table SI. The R package sva (16) was used to de-batch the datasets GSE30528, GSE30529 and GSE30122 to obtain the combined GEO datasets, which contained 19 DKD samples and 50 control samples. Finally, the combined GEO datasets and dataset GSE96804 were subjected to normalization using the R package limma (17), along with annotation of probes and other standardization and normalization procedures.

Ferroptosis-related genes (FRGs) were collected from the GeneCards database by searching for ‘ferroptosis’ and retaining only those categorized as ‘Protein Coding’, yielding 654 FRGs (18). Additionally, a systematic literature search on PubMed using the keyword ‘ferroptosis’ was performed, and FRGs with at least one in vitro or in vivo experimental validation reported in peer-reviewed publications were manually curated, yielding 60 FRGs (19). After merging and removing duplicates, a final set of 667 FRGs was obtained (Table SII). The data processing flowchart is depicted in Fig. 1.

The R package limma (17) was utilized to analyze the gene expression differences between the DKD and control samples. A threshold of |logFC|>0.0 and adj. P<0.05 was established to identify DEGs, with the intent of maximizing the capture of FRGs exhibiting statistically significant changes in DKD. This initial screening strategy was subsequently complemented by multi-step stringent filtering to ensure the robustness of the final diagnostic biomarkers. Specifically, genes with logFC >0.0 and adj. P<0.05 were considered upregulated, while genes with logFC <0.0 and adj. P<0.05 were deemed downregulated. Subsequently, the DEGs were intersected with FRGs to obtain FRDEGs. The R package ggplot2 demonstrated the results of the differential analysis for volcano mapping.

To determine the effect of expression levels of all genes in the combined GEO datasets on DKD, gene set enrichment analysis (GSEA) (20) and gene set variation analysis (GSVA) (21) were performed to explore the differences in pathways and biological processes of the samples in the combined GEO datasets. The c2.all.v7.5.1.symbols.gmt gene set for GSEA and GSVA was obtained from the Molecular Signatures Database (MSigDB) database (22). The screening criteria for GSEA were adj. P<0.05 and FDR value (q-value) <0.25, and the screening criteria for GSVA were |logFC|>0.50 and P<0.05.

GO analysis (23) is essential for exploring biological function, including biological process (BP), cellular component (CC) and molecular function (MF). The KEGG (24) is a widely used database storing information about genomes, biological pathways diseases and drugs and is often used to explore potential pathways. GO and KEGG analysis of FRDEGs were performed using the R package clusterProfiler (25), with screening criteria of adj. P<0.05 and FDR value (q-value) <0.25 considered statistically significant, and adj.p correction by Benjamini-Hochberg.

In order to screen the co-expression modules of the combined GEO datasets and analyze the core genes in the network, WGCNA was performed on the genes in the top 25% of the variance of all samples in the combined GEO datasets using the R package WGCNA package (26). During the process, the correlation of both DKD and control groups with different modules was measured, and the genes in each module were recorded, where the genes in each module were considered module signature genes. Finally, the modules with |r value| >0.30 were screened, and all the genes in the modules were intersected with FRDEGs respectively. All the intersected genes obtained from the different modules were the key genes.

In order to obtain the DKD diagnostic model of combined GEO datasets, the key genes were screened with P<0.05. Subsequently, a logistic regression model was constructed, and the molecular expression of the key genes included in this model was illustrated using a Forest Plot. Furthermore, Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis was performed on the selected key genes included in the logistic regression model using the R package glmnet (27), with the parameter set.seed ‘500’. The outcomes of the LASSO regression analysis were visually presented through diagnostic model plots and variable trajectory plots. In addition, the Support Vector Machine (SVM) model was constructed based on key genes by the SVM (28) algorithm and the key genes were screened based on the number of genes with the highest accuracy and lowest error rate. The key genes in the LASSO regression model and those in the SVM model were taken as intersections to get ferroptosis-related hub genes. R package rms drew a nomogram based on the results of LASSO regression analysis to show the interrelationships of ferroptosis-related hub genes.

To validate the accuracy and resolution of the DKD diagnostic model, the R package ggDCA was used to plot decision curve analysis (DCA) (29) based on the ferroptosis-related hub genes in the combined GEO datasets. Meanwhile, the R package pROC was used to plot the combined ROC curves in the combined GEO datasets and calculate the area under the curve (AUC) values to assess the diagnostic effect of LASSO risk score expression on the occurrence of DKD. The AUC of the ROC curves generally ranges between 0.5 and 1. The closer the AUC is to 1, the better the diagnostic effect. The AUC has a low accuracy at 0.5 to 0.7, a certain accuracy at 0.7 to 0.9, and a high accuracy at AUC above 0.9. The formula for calculating the LASSO risk score is as follows:

The dataset GSE96804 verified the accuracy and resolution of this DKD diagnostic model.

Group comparisons were plotted based on the expression of ferroptosis-related hub genes in the DKD and control groups in the training and validation sets. R package pROC was used to plot the ROC curves of ferroptosis-related hub genes in the training and validation sets and calculate the AUC to assess the diagnostic effect of the expression of ferroptosis-related hub genes on the occurrence of DKD.

CIBERSORT (30) is based on the principle of linear support vector regression for the back-convolution of transcriptome expression matrices to estimate the composition and abundance of immune cells in a mixture of cells. The CIBERSORT algorithm was used to combine the LM22 feature gene matrix and filter the output of data with immune cell enrichment scores greater than zero to finally obtain the specific results of the immune cell infiltration matrix. Finally, the R package ggplot2 was used to draw group comparison plots to show the expression differences of LM22 immune cells in the DKD and control groups in the combined GEO datasets. The R package pheatmap was used to draw correlation heatmaps to show the LM22 immune cells and the correlation analysis results between ferroptosis-related hub genes and LM22 immune cells.

In the present study, a total of 9 male db/db mice (BKS.Cg-Leprdb/Leprdb) and 3 male db/m littermates (BKS.Cg-Leprdb/+) were obtained from Jiangsu Ailingfei Biotechnology Co., Ltd. db/m mice were used in the control group and db/db mice in the DKD group as the DKD model. All mice were 8 weeks of age at the start of the study, with body weights of 40±2.5 g for db/db mice and 25±2 g for db/m mice. The animals were housed in a specific pathogen-free facility under controlled conditions: Temperature, 23±2˚C; humidity, 50-60%; and a 12/12 h light/dark cycle, with free access to standard chow and water. Mice were acclimated for 1 week prior to the experiment. Both rosiglitazone (RZ) and hyperoside (HPS) groups used db/db mice as DKD models, and the two groups were gavaged respectively with RZ (cat. no. 210310; Chengdu Hengrui Pharmaceutical Co., Ltd.) at a dose of 5 mg/kg/d and HPS (cat. no. B20631; Shanghai Yuanye Biotechnology Co., Ltd.) at a quantity of 20 mg/kg/d. The intervention period was 12 weeks. Total RNA was extracted from the kidney of mice using the Animal Total RNA Isolation Kit (cat. no. RE-03011; Chengdu Fuji Biotechnology Co., Ltd.) and quantified using a spectrophotometer (cat. no. N50Touch; Implen GmbH). The thermocycling conditions for qPCR were as follows: Initial denaturation at 95˚C for 2 min, followed by 40 cycles of denaturation at 95˚C for 10 sec and annealing/extension at 60˚C for 30 sec. A melting curve analysis was performed from 65 to 95˚C with a ramp rate of 0.5˚C/5 sec. Relative quantification was performed using the comparative quantification cycle (Cq) (2^-ΔΔCq) method (31). The iScript^™ cDNA Synthesis Kit (cat. no. 1708891; Bio-Rad Laboratories, Inc.) was used strictly according to the manufacturer's instructions. The standard reverse transcription temperature protocol includes primer annealing at 65˚C for 5 min followed by immediate ice incubation, cDNA synthesis at 42-55˚C for 20-60 min [42-45˚C for oligo(dT), 50-55˚C for thermostable reverse transcriptases and GC-rich RNA], and final enzyme inactivation at 70-85˚C for 5-15 min. Subsequently, qPCR was performed using the Sso Advanced^™ Universal SYBR^® Green Supermix kit (cat. no. 1725274; Bio-Rad Laboratories, Inc.) to detect the expression of lactotransferrin (LTF). The experiment was repeated at least three times with three sub-wells for each experiment. The primer sequences utilized in the investigation are presented in Table SIII. Anesthesia was performed when conducting animal sampling and surgical manipulation in the present study. Pentobarbital sodium as the anesthetic agent at a dosage of 50 mg/kg body weight, administered via intraperitoneal injection. All experimental mice were humanely euthanized by cervical dislocation after isoflurane inhalation anesthesia. Animal experiments were approved by the Animal Ethics Committee in Affiliated Hospital of Nanjing University of Traditional Chinese Medicine (approval no. 2023DW-039-01; Nanjing, China).

Data processing and analysis were based on R software (version 4.2.0) and GraphPad Prism (version 9.5.0; Dotmatics). The continuous variables were presented as the mean ± standard deviation. Comparisons between the two groups were made using the Wilcoxon rank sum test. If not specified, results were calculated by Spearman's correlation analysis of the correlation coefficients between the different molecules. P<0.05 was used as a criterion for statistically significant difference.

Firstly, the distribution box plots indicated the successful elimination of the batch effect in the DKD dataset following the batch removal procedure (Fig. 2A and B). Among the combined GEO datasets, 1,237 DEGs were identified, with 731 genes exhibiting upregulation and 506 genes displaying down-regulation (Fig. 2C). By intersecting the 1,237 DEGs with the 667 FRGs, a subset of 47 FRDEGs was obtained, which were LTF, CD44, CYBB, chemokine ligand 5 (CCL5), IGKC, RBBP4, ENO1, ALOX5, PRKCB, LCN2, TPM1, TUBA1A, PDIA4, CASP8, XRCC5, LASP1, ACTB, RPA2, HELLS, PTEN, SIAH2, ACSL4, GFRA1, KRT19, ANXA4, HNRNPL, EEF1A1, MUC1, CTNNB1, FHL2, PGK1, LPIN2, TP53, PDIA6, CP, CAD, VDAC1, SNAP29, SOX15, NFE2L1, G0S2, ALOX12, DLX2, CA9, CLTB, Forkhead box C1 (FOXC1) and TYRO3 (Fig. 2D). Then, based on the intersection results, the results of the expression differences of FRDEGs among different sample subgroups in the combined GEO datasets were shown in Fig. 2E.

To avoid the one-sidedness caused by only using the intersection gene enrichment, we also utilized GSEA on all the genes of the combined GEO datasets. The findings of the GSEA indicated that most genes exerted a significant impact on the inflammatory response pathway IL18 signaling pathway, IL-12 pathway and apoptosis, among others (Fig. 3A-E). The specific outcomes are presented in Table SIV.

The GSVA results showed that a total of 24 pathways exhibited statistical significance (P<0.05) in both the DKD and control groups (Fig. 4A and B). The DKD group exhibited a notable enrichment in the regulation of pathways associated with immunity, phosphorylation, apoptosis, methylation and deacetylation, such as the T cell receptor signaling pathway, monocyte signaling pathway, neutrophil signaling pathway, as well as the regulation of runt-related transcription factors and histone deacetylases. Conversely, pathways linked to lipid synthesis and accumulation were found to be suppressed. The details can be found in Table SV.

A total of 47 FRDEGs were performed for enrichment analysis (Fig. 5A). In GO-BP analysis (Fig. 5B), the FRDEGs were mainly related to apoptosis regulation and cell proliferation. In GO-CC analysis (Fig. 5C), cell cortex, cell membrane and macromolecular compounds were significantly enriched. The results of enrichment analysis in GO-MF (Fig. 5D) indicated that these FRDEGs were mainly related to protein, enzyme and transcription factors. According to the results of the KEGG analysis, these FRDEGs were mainly associated with ferroptosis, HIF-1 signaling pathway, leukocyte trans-endothelial migration and Influenza A (Fig. 5E).

Firstly, the constructed network exhibited higher conformity with the scale-free topology, as evidenced by the computation and presentation of scale-free fitting indices (Fig. 6A) across various soft threshold values. The findings indicated that a scale-free fitting index of 0.85 corresponds to the minimum soft threshold value of 11, which fulfilled the criteria for constructing an optimal scale-free network. The co-expression network was constructed using the optimal soft threshold, and genes with the top 25% variance are clustered using a clustering tree (Fig. 6B) and labeled with the grouping information. The results showed that genes in the top 25% of variance were clustered into 10 modules when the screening criterion was 0.2, and then the relationship between the genes and the merged modules was visualized (Fig. 6C). The correlation between all genes in the ten modules and the DKD group and control group in the combined GEO datasets was determined based on the expression patterns of the module genes (Fig. 6D). Genes within the MEyellow (|r value|=0.45) and MEblue (|r value|=0.32) modules were screened using a criterion of |r value|>0.30, and their intersections with the 47 FRDEGs were plotted using Venn diagrams (Fig. 6E and F). Ultimately, 10 key genes were identified: FOXC1, G0S2, TYRO3, CCL5, FHL2, IGKC, LCN2, LTF, MUC1 and PRKCB.

Logistic regression was performed based on the ten key genes, and a logistic regression model was constructed and visualized by Forest Plot (Fig. 7A). The results showed that all 10 key genes included in the logistic regression model were statistically significant (P<0.05). Then, the LASSO regression model (Fig. 7B and C) was constructed based on 10 key genes, and the results contained four key genes, which were CCL5, FOXC1, G0S2, and LTF. Meanwhile, the SVM model was constructed based on 10 key genes and the SVM algorithm to get the number of genes with the lowest error rate (Fig. 7D) and the highest accuracy (Fig. 7E). The results showed that the SVM model has the highest accuracy when the number of genes is 5. The five key genes were LTF, PRKCB, CCL5, MUC1 and FOXC1. The intersection of the key genes of the two models was received and three ferroptosis-related hub genes were obtained, CCL5, FOXC1 and LTF (Fig. 7F). Finally, to further validate the value of the DKD diagnostic model, a nomogram was developed to illustrate the interrelationships of ferroptosis-related hub genes, as depicted in Fig. 7G. The findings indicated that the expression of LTF exhibited notably greater efficacy in the DKD diagnostic model compared with other variables, while the expression of CCL5 demonstrated significantly lower efficacy in the DKD diagnostic model than the other variables.

In order to evaluate and validate the accuracy and discriminatory capacity of the DKD diagnostic model, the role of the DKD diagnostic model was assessed in clinical utility through DCA based on ferroptosis-related hub genes in the combined GEO datasets. The results showed that the line of the DKD model was stably higher than all positive and all negative in a specific range, resulting in a substantial net gain. Consequently, the model demonstrated effectiveness in accurately diagnosing DKD (Fig. 8A). Furthermore, ROC curves were generated using the LASSO risk score derived from the combined GEO datasets (Fig. 8B). These ROC curves demonstrated the robustness of the LASSO risk score in accurately predicting the expression level within various subgroups, as evidenced by an AUC exceeding 0.9. The calculation of the LASSO risk score was performed as follows:

The DKD diagnostic model's accuracy and discriminatory capacity were validated using the GSE96804 dataset. The results of the DCA demonstrated that the model consistently outperformed both the ‘All positive’ and ‘All negative’ scenarios within a specific range, indicating a high net gain and overall effectiveness (Fig. 8C). Additionally, the ROC curves revealed that the LASSO risk score expression level in the GSE96804 dataset exhibited varying accuracy across different subgroups, with AUC values ranging from 0.7 to 0.9 (Fig. 8D). Finally, a nomogram was constructed to illustrate the interplay among ferroptosis-related hub genes in dataset GSE96804 (Fig. 8E). The results showed that the expression of FOXC1, a ferroptosis-related hub gene, exhibited notably greater efficacy in the DKD diagnostic model compared with other variables. Conversely, the expression of LTF demonstrated significantly diminished utility in the DKD diagnostic model compared with the other variables.

The expression of CCL5, FOXC1, and LTF were all statistically significant (P<0.05) in the DKD and control groups in the combined GEO datasets (Fig. 9A). Whereas in dataset GSE96804 (Fig. 9B), only the expression of FOXC1 and LTF were statistically significant (P<0.05). Subsequently, the accuracy of FOXC1 expression levels in the combined GEO datasets was confirmed a certain level of accuracy (0.7< AUC <0.9), as well as higher accuracy (AUC >0.9) in the GSE96804 dataset. Meanwhile, LTF in the combined GEO datasets exhibited high accuracy (AUC >0.9) in the expression levels across different groups while also demonstrating a certain level of accuracy (0.7< AUC <0.9) in the GSE96804 dataset among different subgroups (Fig. 9C-F).

The CIBERSORT algorithm was employed to calculate the correlation between the 22 immune cells in the DKD group and the control group, and a histogram illustrating the percentage of immune cells in the combined GEO datasets was plotted (Fig. 10A). The box line plot (Fig. 10B) provided evidence of the disparity in immune cell infiltration abundance between DKD and control groups in the combined GEO datasets. The statistical analysis revealed that the expression levels of eight immune cell types, namely B cells naive, plasma cells, T cells gamma delta, NK cells resting, macrophages M1, macrophages M2, mast cells resting and mast cells activated, were significantly different (P<0.05) between the two groups. Subsequently, the correlation between the infiltration abundance of the eight immune cells was illustrated through correlation dot plots (Fig. 10C). The findings revealed that T cells gamma delta and macrophages M2, as well as T cells gamma delta and macrophages M1, exhibited the highest level of positive correlation (r=0.38). Conversely, mast cells resting and mast cells activated displayed the most substantial negative correlation (r=-0.45). Finally, the correlation heatmap analysis revealed that LTF and CCL5, two of the ferroptosis-related hub genes, significantly correlated with immune cells, and both showed a significant positive correlation with immune cell T cells gamma delta (Fig. 10D).

In the DKD group, LTF and FOXC1 exhibited higher expression levels compared with the control group, while CCL5 expression was lower. Upon intervention with RZ and HPS on the DKD model, LTF and CCL5 expression levels were significantly increased compared with the DKD group. Similarly, FOXC1 expression was significantly elevated following interventions with HPS, although FOXC1 exhibited a trend to increase in response to RZ intervention, the difference was not statistically significant (Fig. 11).