The study process is shown in the flowchart in Figure 1. The datasets GSE30122, GSE30528, GSE30529, and GSE96804 datasets freely available through the GEO database (https://www.ncbi.nlm.nih.gov/geo/) were used in this study. The GSE30122 dataset belongs to the GPL571 platform, which includes a total of 50 normal samples and 19 samples with diabetic nephropathy. Of the 50 normal samples, 26 were glomerular-sourced, 12 were tubular-derived and the rest were microtubule-sourced, with 9 of the 19 diabetic nephropathy originating from the glomeruli and the rest from the microtubules. The GSE30122 dataset was used as a training cohort, with GSE30528, GSE30529, GSE96804, and other datasets used as external data to verify the diagnostic performance and expression of key genes in diabetic nephropathy. In addition, we retrieved the gene set related to disulphide stress from the GeneCards database and used data from the Nephroseq v5 database to verify the expression patterns of key genes.

Perl (Strawberry Edition) was used to annotate the data to obtain a gene expression matrix and the average was used if one gene corresponded to multiple probes. The expression profiles of disulfidptosis-related genes were obtained and screened to identify differentially expressed genes in diabetic nephropathy patients through limma packets. The top 20 differentially expressed genes were visualised through the PheatMap package and co-expression network analysis was performed using the WGCNA package.

The gene ontology (GO) enrichment analysis was performed to further understand the physiological mechanisms [biological process (BP), cellular component (CC) and molecular function (MF)] of the differentially expressed genes. KEGG enrichment analysis was also conducted to investigate the pathway enrichment using the enrichplot package in clusterProfiler (p-value = 0.05, q-value = 1). The enrichment results were visualised using the circlize package in ggplot2. Finally, GSEA enrichment analysis was performed for the disulfidptosis-related genes.

The minimum contraction operator method (Lasso) and the support vector machine recursive feature cancellation (SVM-RFE) were used to screen the key genes in diabetic nephropathy. LASSO is an algorithm that builds generalised linear models and filters variables. By introducing the penalty term λ, we can build a simpler model. SVM-RFE is a machine-learning algorithm that selects feature variables to gradually delete unimportant trait genes. Finally, the model genes obtained by the two algorithms were intersected as the key genes of diabetic nephropathy.

The basic principle of a nomogram is a multivariate regression model. Each gene was scored according to its contribution to diabetic nephropathy, then the scores were summed to obtain the probability of developing diabetic nephropathy. Nomograms are being increasingly used in clinical practice because they can transform complex regression equations into easy-to-understand graphs. The scale on the nomogram makes it very easy for clinicians to assess the patient’s risk of developing the disease.

The LASSO and SVM-RFE results were compared to identify the key genes in patients with diabetic nephropathy. Incorporating key genes into nomograms makes predictions easier to understand and increases the readability of the model. The diagnostic performance of the key genes was verified using three datasets (GSE30528, GSE30529, and GSE96804). The GSE30528 dataset has a total of 22 samples, including 9 diabetic nephropathy samples and 13 normal samples. All samples in the GSE30528 dataset were derived from glomerular cells. The GSE30529 dataset also included 22 samples, including 10 diabetic nephropathy samples and 12 normal samples. All samples in the GSE30529 dataset were derived from tubule cells. The GSE96804 dataset contains 41 glomerular samples from patients with diabetic nephropathy and 20 normal glomerular samples. A ROC curve was plotted for each key gene to infer its diagnostic performance and the differential expression of key genes in the different datasets was analysed. Finally, the expression patterns of key genes were verified in the Nephroseq v5 database.

Immunotherapy plays an important role in various diseases, so a comprehensive analysis of immune infiltration in patients with diabetic nephropathy was performed to understand their immune characteristics. The infiltration of 28 immune cells in all samples of the GSE30122 dataset was obtained by ssGSEA enrichment analysis. Next, we analysed the differential infiltration of immune cells between groups using R software and visualised the results in a violin diagram. The Cibersort algorithm was then used to estimate the abundance of immune cells. Next, the differences in immune cell abundance between different groups were analysed to determine the relative content of immune cells and the correlation between immune cells, as well as the correlation between genes characteristic of diabetic nephropathy and immune cells.

The identification of molecular subtypes is particularly important for the personalised treatment of patients with diabetic nephropathy. The ConsensusClusterPlus package was used to type patients with diabetic nephropathy for analysis. Principal component analysis was performed to distinguish patients with different subtypes before the expression patterns of key genes and the immune characteristics of different subtypes were analysed.

Serum creatinine and glomerular filtration rate are often used as clinical indicators to evaluate renal function, so the correlation between key gene expression and serum creatinine and glomerular filtration rate was analysed using the Nephroseq v5 platform.

All statistical analyses were performed in R software (version 4.2.2). A p-value<0.05 was considered significant (∗ = 0.05, ∗∗ = 0.01, ∗∗∗ = 0.001).

Ninety-one differentially expressed disulfidptosis-related genes in diabetic nephropathy patients were identified and the top twenty are presented in heat maps (Figure 2A). The WGCNA analysis identified eight modules relevant for diabetic nephropathy patients (Figure 2B), of which, the turquoise and the green modules were statistically significant (Figure 2C). There were 353 genes in the green module, of which, there are 39 key genes. The intersection of these key genes with the differentially expressed disulfidptosis-associated genes revealed 20 intersecting genes (Figure 2D). The location of the intersecting genes on the chromosomes is shown in Figure 2E. The correlation circle plot of these genes indicated that all genes were positively correlated (Figure 2F).

The circle plot of GO analysis provides a visual representation of the enrichment results (Figure 3A), with the outermost part representing the ID and different colours representing different ontologies. The second circle indicates the number of genes enriched, whereas the third circle indicates the number of differential genes enriched. The bubble chart presents the top ten functions of the biological processes (BP), cellular components (CC), and molecular functions (MF) (Figure 3B). The differentially expressed disulfidptosis-related genes are mainly involved in cellular immunity and cell death in various organelles as well as intracellular and intracellular tissue structures and bind to various functional molecules. The KEGG enrichment analysis revealed that these disulfidptosis-related genes were associated with immune activities such as primary immunodeficiency, neutrophil extracellular trap formation, and ECM-receptor interaction. There was also a correlation with COVID−19. The AGE-RAGE signalling pathway is also strongly associated with differential disulfidptosis-related genes (Figure 3C). GSEA enrichment analysis of the normal group showed that it was closely related to synaptic activity (Figure 3D). The GO analysis of the diabetic nephropathy group showed that the differentially expressed disulfidptosis-related genes were strongly associated with the activation of immune cells as well as the defence response (Figure 3E). The KEGG analysis showed that patients with diabetic nephropathy were associated with systemic lupus erythematosus as well as primary immunodeficiency (Figure 3F).

The top twenty intersecting genes were subjected to Lasso and SVM-RFE analysis. In the Lasso regression analysis, the least error was found when 5 genes were taken, and the model was the most concise when 4 genes were taken (Figures 4A, B), so five genes with the smallest error (CXCL6, LTF, CD48, C1QB, and COL6A3) were used for the subsequent analysis. In the SVM-RFE algorithm, 15 genes (CD48, C1QB, CCL5, CXCL6, IL10RA, TNC, VCAN, COL6A3, CD163, ITGB2, C1QA, COL1A2, CASP1, MMP7, and TLR7) gave the smallest error and the highest accuracy (Figures 4C, D). The intersection of the two algorithms obtained four key genes for diabetic nephropathy (CXCL6, CD48, C1QB, and COL6A3) (Figure 4E) which were used to construct the nomogram (Figure 4F). The DCA curve indicates that the nomogram has some application (Figure 4G) and the calibration curve shows the accuracy of the nomogram predictions (Figure 4H). In the volcano plot, CXCL6, CD48, C1QB, and COL6A3 are all highly expressed in diabetic nephropathy patients (Figure 4I), with areas under the ROC curve of 0.908, 0.928, 0.907, and 0.895, respectively (Figure 4J). The area under the ROC curve showed that all four genes had good diagnostic value, with CD48 having the greatest value, so the nomogram with these four key genes shows better diagnostic value (Figure 4K).

The AUC values of CXCL6, CD48, C1QB, and COL6A3 in the GSE30528 dataset were 0.838, 0.915, 0.974, and 0.889, respectively. The 95% confidence intervals were 0.615–1.00, 0.760–1.00, 0.897–1.00, and 0.709–1.00 (Figures 5A–D). The ROC curve areas under CXCL6, CD48, C1QB, and COL6A3 in the GSE30529 dataset were 1.00, 0.983, 0.942, and 0.942, respectively. The 95% confidence intervals were 1.00–1.00, 0.925–1.00, 0.825–1.00, 0.825–1.00 (Figures 5E–H). The ROC curve areas under CXCL6, CD48, C1QB, and COL6A3 in the GSE96804 dataset were 0.610, 0.692, 0.706, and 0.920, respectively. The 95% confidence intervals were 0.468–0.744, 0.541–0.830, 0.561–0.843, 0.844–0.977 (Figure 5I–L). The area under the ROC curve in the external verification data was greater than 0.6, indicating the satisfactory diagnostic performance of the key genes of diabetic nephropathy.

The volcano map in Figure 4I showed that CXCL6, CD48, C1QB, and COL6A3 were highly expressed in diabetic nephropathy patients, suggesting that these genes may be strongly related to the onset of diabetic nephropathy. This was verified in three external GEO datasets (Figure 6A–L) and the Nephroseq v5 database (Figure 6M–P).

The ssGSEA analysis revealed that almost all immune cells were present in higher amounts in samples from diabetic nephropathy patients than in controls, with the difference in 17 immune cells being statistically significant (Figure 7A). The Cibersort analysis showed that T cells gamma delta, resting NK cells, monocytes, M2 macrophages, resting dendritic cells and resting mast cells were significantly different between the two groups, with only monocytes being more abundant in the normal group (Figures 7B, C). The correlation heatmap shows the correlation between immune cells, where red represents a positive correlation and blue represents a negative correlation (Figure 7D).

The correlation of CXCL6, CD48, C1QB, and COL6A3 with immune cells is shown by lollipop plots, with the most strongly correlated immune cells with key genes shown separately by correlation scatter plots (Figures 8B, C, E, F, H, I, K, L). CXCL6 was positively correlated with resting mast cells, M2 macrophages, and T cells gamma delta and negatively correlated with monocytes and resting NK cells (Figure 8A). CD48 is positively correlated with T cells gamma delta, eosinophils, and M2 macrophages and negatively correlated with regulatory T cells and resting NK cells (Figure 8D). C1QB was positively correlated with M2 macrophages, T cells gamma delta, M1 macrophages and resting mast cells resting, and negatively correlated with neutrophils (Figure 8G). COL6A3 was positively correlated with T cells gamma delta, M2 macrophages, eosinophils, resting dendritic cells and resting mast cells but negatively correlated with monocytes, neutrophils, and resting NK cells (Figure 8J).

From the heat map of the classification, we can see that it is most suitable when classifying diabetic nephropathy patients into type 2 (Figures 9A–H). Principal component analysis also showed good differentiation between the two subtypes (Figure 9I). In addition, CXCL6 and COL6A3 expression was significantly different in the C1 and C2 subtypes (Figure 9J), with CXCL6 and COL6A3 highly expressed in the C1 subtype. In the immune infiltration of subtypes, resting NK cells and activated mast cells differed significantly between the two subtypes (Figure 9K), with extensive infiltration of both cells in the C2 subtype. Similarly, the proportion of immune cell infiltration in different subtypes of samples is evident in the scale plot (Figure 9L).

The glomerular filtration rate of patients gradually decreases as the expression of CXCL6, COL6A3, and C1QB increases (Figures 10A–C), that is, the patient’s kidney function gradually deteriorates with the expression of key genes. In addition, CXCL6, CD48, and C1QB expression were positively correlated with plasma creatinine (Figures 10D–F), indicating the higher CXCL6, CD48, and C1QB expression, the poorer the renal function. Thus, CXCL6, CD48, C1QB, and COL6A3 may be involved in kidney function deterioration in diabetic nephropathy patients. Since there are certain similarities between different kidney diseases, we explored the expression of key genes in different kidney diseases, finding that C1QB was highly expressed in lupus nephritis with the lowest expression in membranous glomerulonephropathy (Figure 10G). CXCL6 is highly expressed in diabetic nephropathy, second only to diabetic nephropathy in vasculitis, and approximately the same expression in other kidney diseases (Figure 10H). CD48 expression is similar to CXCL6 (Figure 10I). COL6A3 maintains high levels of expression in all kidney diseases (Figure 10J).