We screened three diabetic nephropathy datasets: GSE30528 (GPL571) contained nine cases of diabetic nephropathy and thirteen controls; GSE104948 (GPL22945) served as a validation set and contained seven cases of diabetic nephropathy and eighteen controls; and GSE131882 (GPL24676) contained three early diabetic nephropathy and three control samples for single nucleus RNA sequencing. Additionally, using a relevance score of greater than 7 as a screening criterion, we were able to extract 855 genes associated with oxidative stress from the Genecard database. Table 1 displays the pertinent details.

With |log2 fold change (FC)| > 0.5 and p < 0.05 as screening criteria, differentially expressed genes (DEGs) from GSE30528 were identified utilizing “Limma” R package, where log FC > 0.5, p < 0.05 was Up, log FC < -0.5, p < 0.05 was Down. The heat map and volcano map of DEG were plotted using the “Pheatmap” R package and “ggplot2” R package, respectively.

Subsequently, the obtained DEGs were intersected with 855 oxidative stress-related genes to obtain differentially expressed genes related to oxidative stress (DEOSGs).

CIBERSORT employs a deconvolution algorithm to estimate the composition and abundance of immune cells in a mixture of cells based on transcriptome data. In the present study, we first assessed the proportion of 22 immune cell species in normal and diabetic nephropathy samples in GSE30528 using the CIBERSORT algorithm (7).

Weighted Gene Go-expression Network Analysis (WGCNA) is performed to identify modules of highly correlated genes, summarize the interconnections between modules and associations with external sample traits, and identify candidate biomarkers or therapeutic targets. In our research, WGCNA was constructed by the R package “WGCNA” to identify the modules with the highest relevance to immune cells in diabetic nephropathy patients (8). Specifically, we preprocessed the sample data and removed the outliers. Subsequently, the correlation matrix was constructed by the “WGCNA” software package. The optimal soft threshold was chosen to convert the correlation matrix into an adjacency matrix, and a topological overlap matrix (TOM) was created from the adjacency matrix. The TOM-based phase dissimilarity metric was utilized to categorize genes with similar expression patterns into gene modules using average linkage hierarchical clustering. The two modules with the strongest relevance to immune cells were selected as key modules for subsequent analysis.

Finally, the genes in DEOSGs and key modules were intersected, and the intersected genes were described as differentially expressed immune-related oxidative stress genes (DEIOSGs) for further study.

In this research, the “clusterProfiler” R package was implemented to conduct GO and KEGG functional enrichment analysis in R to assess gene-related biological processes (BP), molecular functions (MF), cellular components (CC), and gene-related signaling pathways.

Least Absolute Shrinkage and Selection Operator (LASSO) logistic regression analysis is a data mining method that sets the coefficients of less important variables to zero by applying the L1-penalty (lambda) in order to filter out the significant variables and construct the best classification model (9). Support Vector Machine-Recursive Feature Elimination (SVM-RFE) analysis is a supervised machine learning technique for identifying the optimal core genes by dropping the feature vectors generated by SVM (10). Random Forest (RF) analysis is a decision tree-based machine learning method that focuses on evaluating the significance of variables by scoring the importance of each variable (11). In combination with machine learning algorithms, the cytoHubba plugin is frequently applied for the identification of key genes. On the one hand, diagnostic genes from DEIOSG were assessed using the three machine learning algorithms separately (12). After that, the intersection of the three machine learning algorithms was established.

On the other hand, the STRING database was exploited to establish protein-protein interaction (PPI) networks, which Cytoscape then visualized. The differential genes were then evaluated using 12 algorithms in the cytoHubba plugin, and finally the top 10 genes for each algorithm were taken as intersection and visualized through the ImageGP platform (13).

Ultimately, the genes obtained by both methods in total were identified as hub genes.

The Nephroseq v5 database (http://v5.nephroseq.org) (14) is a comprehensive information platform for evaluating the correlation between gene expression levels and clinical characteristics of kidney diseases. To explore the correlation between the expression of hub genes and clinical features, we mined the Nephroseq v5 database.

We performed a single-gene GSEA analysis to investigate the possible roles of hub genes.

The JASPAR database (15) and the TarBase database (16) were accessed by the NetworkAnalyst (https://www.networkanalyst.ca/) (17) to predict transcription factors (TFs) and miRNAs, respectively. Subsequently, the results were visualized using Cytoscape software.

We used the Enrichr platform (https://amp.pharm.mssm.edu/Enrichr/) (18) to access the DSigDB database (19) for potential drug prediction.

A single-cell sequencing database for kidney disease called the Kidney Integrative Transcriptomics (K.I.T.) database was developed by Ben Humphrey’s lab at Washington University (http://humphreyslab.com/SingleCell/) (20). To explore the distribution of hub genes in cell groups, we applied the database for analysis and visualization of the results. In one of them, we used single nucleus RNA sequencing data from diabetic nephropathy that was initially taken from the GSE131882 dataset.

GraphPad Prism 8.0 (GraphPad Software, CA, USA) was implemented to conduct the statistical analysis. The diagnostic value of hub genes was evaluated with ROC curve analysis. Hub genes were analyzed for correlation with clinical features via Pearson analysis. An unpaired t-test was performed for the assessment of hub gene differential expression. P < 0.05 was defined as statistically significant.

A total of 1696 DEGs were acquired from GSE30528, and another 855 oxidative stress-related genes were mined from the Genecard database, and 111 DEOSGs were generated by taking the intersection of the two ( Figures 2A–C ).

Five immune cell types, including T cells CD4 naive, T cells gamma delta, NK cells resting, Dendritic cells resting, and mast cells resting, were demonstrated to be comparable between DN and control samples using the CIBERSORT algorithm ( Figure 3A ).

The soft-threshold power in this research was calibrated to 14 (scale-free R² = 0.85) ( Figure 3B ). Last but not least, a sum of 11 modules was revealed by the WGCNA analysis ( Figure 3C ). In particular, the green module and the magenta module had strong positive correlations with T cell CD4 naive and gamma delta subsets, respectively. Due to their significance in association with immunological infiltrating cells, the green and magenta modules were considered for additional investigation.

DEIOSGs are the genes that overlap DEOSGs with the magenta and green modules generated by WGCNA, and a total of 24 DEIOSGs were identified ( Figure 4A ).

Furthermore, we performed the functional enrichment of 24 DEIOSGs via GO and KEGG. In the BP assessment, DEIOSGs were mostly engaged in superoxide metabolic processes, neutrophil activation, and other functions. DEIOSGs have been localized to the external side of the plasma membrane, endocytic vesicle, and other structures in CC. DEIOSG changes associated with MF include amide binding, integrin binding, and superoxide-generating NAD(P)H oxidase activity ( Figure 4B ). According to KEGG analysis, DEIOSGs are particularly abundant in leukocyte transendothelial migration, neutrophil extracellular trap formation, lipid and atherosclerosis, diabetic cardiomyopathy, natural killer cell mediated cytotoxicity and other pathways ( Figures 4C, D ).

Firstly, 6 genes were extracted from DEIOSGs using the LASSO regression algorithm ( Figure 5A ). Secondly, the SVM-RFE algorithm identified 6 genes ( Figure 5B ). Then, 7 genes were selected by the RF algorithm ( Figure 5C ). Subsequently, the three were overlapped by the Venn diagram and finally two genes were obtained, namely CD36 and SLC1A3 ( Figure 5D ). Meanwhile, from the PPI network, we obtained a gene, namely ITGB2, through the cytoHubba plugin ( Figures 6A, B ). Finally, a total of 3 hub genes were identified by both methods, all of which were up-regulated.

When compared to the normal control sample, we discovered in the GSE30528 dataset that these genes were expressed more highly in DN ( Figures 7A–C ). We next confirmed the expression of these genes using another dataset, and the results revealed that these genes were likewise more strongly expressed in DN than control in GSE104948, and they were all statistically significant ( Figures 7D–F ).

To explore the diagnostic efficacy of the 3 hub genes, we implemented a ROC curve analysis in which hub genes with an AUC value > 0.7 were used as diagnostic markers. In the GSE30528 dataset, the AUC values were 0.8215 for CD36, 0.9402 for SLC1A3, and 0.9060 for ITGB2 ( Figures 8A–C ).

In the GSE104948 dataset, the AUC values of CD36 were 1.000 (95% CI: 1.000-1.000), AUC values of SLC1A3 were 0.7937 (95% CI: 0.5244-1.000), AUC values of ITGB2 were 0.9921 (95% CI: 0.9669-1.000) ( Figures 8D–F ).

According to GSEA findings, the CD36 high expression group was highly enriched for primary immunodeficiency and viral protein interaction with cytokines and cytokine receptors ( Figure 9A ). The ITGB2 high expression group was mostly concentrated in the citrate cycle (TCA cycle) and proteasome ( Figure 9B ). Allograft rejection, primary immunodeficiency, and systemic lupus erythematosuswere all associated with increased SLC1A3 expression ( Figure 9C ).

In DN patients, correlation analysis revealed a negative correlation between CD36 expression and glomerular filtration rate (GFR) (r = -0.860, p < 0.001) and a positive correlation between CD36 expression and serum creatinine (r = 0.887, p < 0.001) ( Figures 10A, B ). ITGB2 expression was negatively correlated with glomerular filtration rate (GFR) (r = -0.2031, p = 0.6002) but not statistically different, whereas ITGB2 expression was positively correlated with serum creatinine (r = 0.5590, p = 0.020) ( Figures 10C, D ).

Using the JASPAR database, 31 TFs were finally obtained, among which, there were 9 TFs with degree≥2, and they were FOXC1, FOXL1, YY1, PPARG, STAT3, HINFP, MAX, USF1, USF2 ( Figure 11A ). Possible miRNAs were predicted by the TarBase database with 10 miRNAs of degree≥2 ( Figure 11B ).

Eighty-seven potential therapeutic agents were screened in the DSigDB database with a cut-off value of Adjusted p-value < 0.05 ( Supplementary Table 1 ).

By single nucleus RNA sequencing, we determined the distribution of CD36, ITGB2 and SLC1A3 in 12 cell groups ( Figure 12A ), among which CD36 was mainly distributed in endothelium and ITGB2 and SLC1A3 were highly expressed in leukocyte ( Figures 12B-D ).