The GEO database (http://www.ncbi.nlm.nih.gov/GEO) was accessed in order to retrieve DN microarray datasets. There were a total of five glomerular DN datasets (GSE30122, GSE47185, GSE96804, GSE99340, and GSE104948), one tubulointerstitial DN dataset (GSE104954), and one kidney biopsy DN dataset (GSE142025) downloaded, with the specifics of each dataset shown in Table 1 and Supplementary Table 1 . Each cohort’s probe ID was initially annotated and converted to “Entrez ID” according to platform annotation documents. Within five glomerular DN datasets, R (version 4.2.0) was used to construct a multi-chip dataset containing 100 healthy controls and 90 DN patients. Then, the R package “SVA” (version 3.46.0) containing the “Combat” function was applied to the multi-chip dataset to correct the batch effects. Using principal component analysis (PCA) to evaluate if the batch effect has been abolished, with GSE104954 and GSE142025 serving as the validation data set at last.

Utilizing the R package “limma” (version 3.54.2) (15), a differential expression analysis was conducted between the DN and control groups using the multi-chip dataset. Then, genes with a |log-fold change (FC)| > 0.5 and a Benjamini and Hochberg adjusted with p-value < 0.05 were selected as differentially expressed genes (DEGs).

To further illuminate the biofunction of the selected DEGs, Gene Ontology (GO) enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, and Disease Ontology (DO) enrichment analyses were done on DEGs using “ClusterProfiler” (version 4.6.2) within R and then visualized using “ggplot2” (version 3.4.2). And the biological process (BP), cellular component (CC), molecular function (MF), signal pathway, and disease type were screened with a threshold of p-value < 0.05. In addition, gene set enrichment analysis (GSEA) was utilized to determine the functional words that were most significant between the DN and control groups (16). The hallmark gene sets “c5.go.v7.4.symbols.gmt” and “c2.cp.kegg.v7.4.symbols.gmt” were obtained from the Molecular Signatures Database (MSigDB). These gene sets summarize and represent well-defined biological states or processes and have consistent expression. To obtain a normalized enrichment score for each analysis, 1000 times gene set permutations were performed. After 1000 permutations, a false discovery rate (FDR) < 0.25 and p-value < 0.05 was considered highly enriched.

Weighted gene co-expression network analysis (WGCNA) is a system biology technique used to examine gene connection patterns between samples. It can be used to identify highly synergistic gene expression matrices, and to identify candidate core genes based on the interconnectedness of expression matrices and the degree of association between genes and phenotypes (17). In the current investigation, the “WGCNA” (version 1.72-1) function within R was used to discover related gene modules and to screen gene sets that may be especially associated with glomerular damage in DN. The module analysis soft threshold β is obtained by assessing the scale independence and average connectivity of modules under various weighted coefficients. After determining the soft threshold β, the scale-free topological distribution network was constructed, and the correlation matrix was transformed according to the pearson correlation coefficient between genes into the adjacency matrix, then into the topological overlap matrix (TOM) to obtain the differences between genes (1-TOM). In addition, the modules incorporated hierarchical clustering and dynamic tree cutting function detection. To classify genes with similar expression profiles into gene modules, the average linkage hierarchical clustering of the gene tree was performed using a “TOMbased” difference measurement method with a minimum genome size of 50, and the module membership (MM, correlation between specific genes and module characteristic genes) and gene significance (GS, correlation between specific genes and clinical variables) were calculated. Finally, the network of feature genes were shown.

First, we extracted 1406 ERSRGs from the MSigDB (GOBP response to endoplasmic reticulum stress and GOBP regulation of response to endoplasmic reticulum stress) and GeneCards database (with relevance scores ≥ 10) (18), which are exhaustive datasets that curate ERSRGs from research articles ( Supplementary Table 2 ). Second, we intersected these ERSRGs with the DEGs derived from the multi-chip dataset and genes in the major modules of WGCNA to obtain the gene expression profile of hub ERSRGs. Using the “pheatmap” (version 1.0.12) R package, the heat map depicting the expression of hub ERSRGs was created. As previously noted, GO enrichment analysis and KEGG pathway analysis were performed using the “ClusterProfiler” tool to show the biofunction of hub ERSRGs.

The discovery of biomarkers has made extensive use of machine learning algorithms, which can produce more detailed models. Three machine learning classifiers were used to filter feature genes in the current study: least absolute shrinkage and selection operator (LASSO), support vector machine recursive feature elimination (SVM-RFE), and random forest classifier (RF) (19). Using the “cv.glmnet” function in the R package “glmnet” (version 4.1-7), a LASSO regression prediction model was developed that could fit the generalized linear model, while variables were filtered and complexity was simultaneously changed. The R package “e1071” (version 1.7-13) was used to do SVM-RFE, the 10-fold cross-validation algorithm was used as the resampling method for SVM-RFE, and the features were sorted by recursion. Moreover, using the “RandomForest” (version 4.7-1.1) function, a RF analysis was conducted, and decision tree classifier models were configured to score the classification variables iteratively, thereby finding features with high classification accuracy. The genes within the intersection of three subsets were then selected for further study as characteristic genes.

At last, to verify the selected intersecting genes reliability, gradient boosting decision tree (GBDT) model was used as well. GBDT is one of boosting methods, which uses the method of gradient boosting to carry on each iteration and finally builds a strong model (20). And the algorithm often needs to generate a certain number of decision trees to achieve the accuracy of satisfaction. So we used GBDT based on recursive feature elimination with 10-fold cross-validation to test our result in above three models. This part was carried out in Python (version3.9).

A co-expression pattern network diagram was developed using the R package “corrplot” (version 0.92) based on the correlation between gene expression levels in order to elucidate the interaction between the characteristic genes. Besides, the tubulointerstitium (GSE104954) and kidney biopsy (GSE142025) datasets were used for external validation of the ability of the characteristic genes to differentiate DN from healthy control. Using the unpaired t-test, it was determined if there was a difference in the expression of specific genes between the two groups at p-value < 0.05.

The R package “rms” (version 6.6-0) was used to combine a total of the characteristic genes into a logistic regression (LR) model, which was displayed as a nomogram. The area under the receiver operating characteristic (ROC) curve (AUC) was recognized as the quantitative evaluation criterion for determining the discrimination capacity of each characteristic gene and the nomogram. The ROC analysis was performed using the R package “pROC” (version 1.18.0) (21).

We used the CIBERSORT method to identify the ratios of various immune cells in multi-chip data sets of DN and control groups and to depict the abundance of immune cells based on the gene expression matrix of multi-chip datasets. The “corrplot” program was used to construct the heat map illustrating the quantitative relationship between distinct immune cells, and p-value < 0.05 indicates a statistically significant difference between the two groups. In addition, “ggplot2” was employed to analyze the correlation between the expression of characteristic genes and the proportions of immune cells.

Using the k-means algorithm of unsupervised clustering analysis (“ConsensusClusterPlus” R package version 1.62.0) (22), we classified 90 DN samples into distinct clusters based on the expression level of the characteristic genes with a total of 1,000 iterations. The best number of categories was found through cumulative distribution function (CDF) curves, a consensus matrix, and a cluster score that was greater than 0.90. Then, PCA analysis was done to illustrate the distributional distinction between immune subtypes, which was shown graphically.

GSVA, which is a non-parametric and unsupervised method for assessing the variation of gene set enrichment across a sample population (23), is utilized to elucidate the molecular mechanisms features between various ER stress subtypes. As reference sets, “c5.go.v7.4.symbols.gmt” and “c2.cp.kegg.v7.4.symbols.gmt” were selected from the MSigDB. The R package “GSVA” (version 1.46.0) and the “ssGSEA” function were utilized to determine the GSVA score for each gene set. The GSVA score indicates the enrichment of each genome in absolute terms. And the “limma” package was utilized to compare differences in GSVA score between subtypes for each genome with a threshold of p-value < 0.05.

The CMap (https://portals.broadinstitute.org/cmap) is a public database that links diseases, genes, and medications based on similar or opposite gene expression profiles (24). It was utilized to identify prospective DN-targeting therapies. Upregulation and downregulation ERSRGs were translated into chip-specific probe sets for querying the CMap. Using norm cs and FDR, we created a list of CMap instances predicted to invert ERSRGs.

Computer is running the Windows 10 64-bit system, with Intel Core i5-13500H processor and RTX 4050 graphics card, and the highest Rui frequency is 4.7 GHz. All data calculations and statistical analysis were performed with the R (version 4.2.0) or Origin (version 9.1) software. ROC curves and AUC values were utilized to assess the predictive performance of the diagnostic model. Pearson’s analysis was used to conduct correlation analysis. In addition, an unpaired t-test was utilized to assess the differential expression levels of the DN-specific genes. All p-values were bilateral, and p-value < 0.05 was deemed statistically significant.

The multi-chip dataset included 100 healthy controls and 90 DN patients. After batch correction, PCA analysis revealed that the data distribution of each data set tended to be uniform ( Figures 1A, B ), implying that normalization was likely done properly. According to the screening criterion (|logFC| > 0.5, p-value < 0.05), 497 DN-specific related DEGs were found across five databases, 242 of which were downregulated and 255 of which were upregulated compared to healthy controls ( Supplementary Figure 1A ).

Further functional enrichment analysis uncovered the following GO terms to be the most significantly enriched: BP: urogenital system development; CC: collagen-containing extracellular matrix; MF: glycosaminoglycan binding ( Supplementary Figure 1B ). Further significantly enriched KEGG pathways were screened: PI3K-Akt signaling pathway, complement and coagulation cascades, focal adhesion, phagosome, AGE-RAGE signaling pathway in diabetic complications, protein digestion and absorption, ECM-receptor interaction, and renin-angiotensin system et al. ( Supplementary Figure 1C ). DO analysis revealed that majority of the genes were involved in urinary system and kidney diseases ( Supplementary Figure 1D ). In addition, the GSEA results demonstrated that in the renal tissue gene expression matrix of DN patients, active GO functions are mainly enriched in biological processes including activation of immune response, adaptive immune response, and external encapsulating structure organization and extracellular matrix-related functions ( Figure 1C ). The active KEGG pathway were mainly enriched in cell adhesion molecules cams, chemokine signaling pathway, ECM receptor interaction, and focal adhesion ( Figure 1D ).

This work did a WGCNA analysis on the multi-chip dataset, and 100 healthy controls and 90 DN samples were chosen to cluster the samples and remove the obviously aberrant samples by establishing a threshold, as depicted in Figure 2A . Then, as shown in Figure 2B , when R2 is greater than 0.90 and the average connectivity is high, we set the soft threshold to 6. After combining the strongly linked modules using a 0.25 clustering height restriction, 13 modules were chosen for further investigation and displayed beneath the clustering tree ( Figure 2C ). Using the frontal correlations between ME values and clinical characteristics, the relationship between modules and clinical symptoms was investigated. The brown module was positively correlated with control (r = 0.54, p = 6e-16) and negatively linked with DN (r = -0.54, p = 6e-16), and the salmon module was positively correlated with control (r = 0.51, p = 5e-14) and negatively linked with DN (r = -0.51, p = 5e-14), while the black module was negatively connected with control (r = -0.6, p = 9e-20) and positively correlated with DN (r = 0.6, p = 9e-20), and the green module was negatively connected with control (r = -0.51, p = 9e-14) and positively correlated with DN (r = 0.51, p = 9e-14) ( Figure 2D ). Transcriptional correlation study within modules validated the reliability of module delineation by revealing no meaningful relationship between modules ( Figure 2E ). The findings of an examination of the correlation between modules revealed that there was no significant relationship between them ( Figure 2F ). Clinically relevant modules have been discovered. The findings revealed that the brown, black, green, and salmon modules were strongly associated with DN ( Figure 2G ). All genes within the four modules were investigated in greater detail.

After overlapping major module genes from WGCNA, DEGs, and ERSRGs using a Venn diagram, we discovered 49 overlapping genes ( Figure 3A ; Supplementary Table 3 ). A heat map was then created to illustrate the change and cluster relationship of hub ERSRGs ( Figure 3B ). In addition, we performed functional analysis in order to acquire a deeper comprehension of the biological functions of the hub ERSRGs. These genes were discovered to be linked with oxidative stress, collagen-containing extracellular matrix, endoplasmic reticulum lumen, extracellular matrix structural constituent, Toll-like receptor binding, and additional GO terms ( Figure 3C ). Further KEGG analysis showed that hub ERSRGs mainly participated in AGE-RAGE signaling pathway in diabetic complications, PI3K-Akt signaling pathway, focal adhesion, fluid shear stress and atherosclerosis, lipid and atherosclerosis, ECM-receptor interaction, endocrine resistance, phagosome, IL-17 signaling pathway, Fc gamma R-mediated phagocytosis, and HIF-1 signaling pathway et al. ( Figure 3D ).

Separately, the LASSO, SVM-RFE, and RF algorithms were used to identify feature genes. When the construction of LASSO based on 10-fold cross-validation was applied, the minimal error value corresponded to twenty-one characteristic genes, including CCL2, CDKN1B, COL1A1, COL4A1, COMP, EGR1, ELOVL4, EXT1, FKBP5, FOS, GDF15, IGF1, KPNA2, LPL, MARCKS, NQO1, NUPR1, PLA2G4A, RSAD2, and S100A9 ( Figure 4A ; Supplementary Table 4 ). The SVM-RFE algorithm was validated by 10-fold cross-validation as well, and when the algorithm was most accurate and the estimation error was smallest, nine genes, including CCL2, CDKN1B, EGR1, FKBP5, GDF15, MARCKS, NQO1, PLA2G4A, and PLN, were determined ( Figure 4B ; Supplementary Table 5 ). RF in combination with feature selection was used to determine the association between error rate, number of classification trees, and the top 20 genes by weight were selected ( Figure 4C ). Consequently, five overlapping genes (CDKN1B, EGR1, FKBP5, GDF15, and MARCKS) were recognized based on the findings of above machine learning models ( Figure 4D ). Moreover, when testing the five overlapping genes reliability, these genes were selected by GBDT in the same way ( Figure 4E ).

Based on the microarray expression matrix, the expression patterns of characteristic genes were investigated and validated. A heat map of co-expression correlations among the five genes were produced and displayed a significant interaction relationship ( Figure 5A ). Then we verified the expression of these five genes in the multi-chip dataset and found that CDKN1B, EGR1, FKBP5, and GDF15 expression was significantly lower in DN samples than in control samples, but MARCKS expression was significantly higher ( Figure 5B ). Significantly, gene expression patterns were found to be consistent between the tubulointerstitium (GSE104954) and kidney biopsy (GSE142025) testing cohorts ( Figures 5C, D ).

LR generated the prediction model based on the multi-chip dataset and displayed it as a nomogram with a c-index of 0.994 ( Figure 5E ), indicating a high correlation degree. Five genes identified in the multi-chip dataset (AUC > 0.6) exhibit diagnostic efficacy in differentiating DN patients from healthy controls, as confirmed by the ROC curve analysis of diagnostic power of characteristic genes ( Figure 5F ). In addition, external validation of nomogram in the validation data set (GSE104954 and GSE142025) showed that the c-index was 0.941 and 0.63, respectively ( Supplementary Figures 2A, B ). Notably, compare to single-core genetic models, the model of prediction nomogram had the best diagnostic performance. According to our findings, the aforementioned ER stress-related characteristic genes possessed a high diagnostic power and had the potential to act as diagnostic biomarkers for DN.

Through the CIBERSORT method, we measured the characteristics of immunocytes between DN and control groups to further investigate the differential expression of immune components, and the cumulative histograms displayed the relative proportions of 22 immune cells ( Figure 6A ). The results indicated that activated Mast cells, activated Macrophages M2, activated Monocytes, activated NK cells, naive T cells CD4, and plasma cells constituted the majority. The correlation heat map for 22 immune cell types revealed a strong positive link between activated NK cells and Macrophages M1 and a significant negative correlation between activated Mast cells and resting Mast cells, Neutrophils and Macrophages M2, as well as Plasma cell and Macrophages M1 ( Figure 6B ). In addition, the violin plot regarding the difference in immune cell infiltration revealed significantly more B cells memory, T cells gamma delta, NK cells activated, Macrophages M2, Dendritic cells resting, and Mast cells resting in the experimental group than in the control group, whereas B cells naive, Macrophages M0, and Mast cells activated were significantly decreased ( Figure 6C ). Via analyzing immune cell correlations, we determined that specific genes may be implicated in the advancement of DN by controlling immune cells such as mast cells, neutrophils, and NK cells ( Figure 6D ).

To illustrate the ER stress-related patterns in DN, we designed a new consensus clustering approach to classify 90 DN samples based on the expression landscapes of five ER stress-related characteristic genes. The optimal number of subtypes was two, which was determined using a consensus matrix plot, a CDF plot, relative alterations in the area under the CDF curve, and consistent cluster score (> 0.9) ( Figures 7A–D ). We thus divided DN samples into two distinct subtypes, including subtype1 (n = 49) and subtype2 (n = 41). The PCA analysis demonstrated the striking distinction between the subtypes ( Figure 7E ).

To clarify the molecular differences between these subtypes, the differential expression of characteristic genes among different ER stress-related subtypes was first evaluated. Subtype1 exhibited higher expression of CDKN1B, EGR1, FKBP5, and GDF15, whereas subtype2 was characterized by higher MARCKS expression ( Figure 8A ). However, there was no substantial difference in immune cell infiltration between the two subtypes ( Figure 8B ). Then, the GSVA analysis was conduct to evaluate the differences of molecular mechanisms in the subtypes with different ER stress expression patterns. The functional enrichment results suggested that negative regulation of cell substrate adhesion, positive regulation of hematopoietic stem cell proliferation, negative regulation of cell matrix adhesion, and negative regulation of intracellular steroid hormone receptor signaling pathway were prominently upregulated in subtype2 ( Figure 8C ). In addition, pathways enrichment results revealed that TGF-β signaling pathway, WNT signaling pathway, focal adhesion, and ECM receptor interaction were significantly elevated in subtype2 ( Figure 8D ). We determined that subtype2 is more closely related with ER stress in DN.

Moreover, we performed a CMap (http://portals.broadinstitute.org/cmap/) analysis to find potential drug candidates for DN to reverse the altered expression of hub ERSRGs. Table 2 displays the top10 CMap compounds, with the protein synthesis inhibitor bruceantin having the highest negative connectivity score.