The transcriptome matrix was downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). In this study, GSE96804 and GSE30122 were downloaded from the GEO database, and the batch effect of the transcriptome matrix was removed via “SVA” R package. The GSE96804 contains 20 human glomeruli control samples (unaffected portion of tumor nephrectomies) and 41 human glomeruli DN samples ([HTA-2_0] Affymetrix Human Transcriptome Array 2.0 [transcript (gene) version]). The GSE30122 contains 26 glomerulus of control kidney samples, and 9 glomeruli of DN samples ([HG-U133A_2] Affymetrix Human Genome U133A 2.0 Array). Perl scripts were conducted to extract the gene expression matrix and the probes were annotated based on the platform annotation file.

The metabolism-related genes (MRGs) were obtained from the Molecular Signatures Database (MSigDB) (http://software.broadinstitute.org/gsea/msigdb). “C2.cp.kegg. v7.5.1.symbols.gmt” was used as the reference genes set and a total of 739 unique MRGs were extracted from the gene matrix using the Perl scripts. “limma” R package was conducted to identify the differential expression MRGs (DE-MRGs) in HC and RA group, and the threshold was set at |Fold Change| ≥ 1 and P-value < 0.05.

The WGCNA package was used for weighted gene co-expression network analysis (WGCNA) in R software. Firstly, DE-MRGs with variances greater than all quartile variances were selected to perform WGCNA analysis, because the results of network module analysis are easily affected by outlier samples, it was particularly important to remove outliers before constructing the network. In this study, the inter-array correlation (IAC) between chips was used to evaluate the distribution of microarray data and outliers with significantly lower mean IAC values will be removed. The co-expression network of genes conforms to no scale distribution, in another word, it follows a power law distribution. WGCNA selected the weighted coefficients to obtain the results that most accord with the scale-free network distribution. Finally, genes and gene modules were associated with clinical information to identify key genes with potential biological significance in the network. For any genes, gene significance (GS) relative to a certain dependent variable was defined as the correlation coefficient between its expression level and the level of the dependent variable. Pearson correlation coefficient was used for continuous dependent variables. For a Module, the definition of Module Significance (MS) relative to a certain dependent variable is the correlation coefficient between its characteristic gene and the level of the dependent variable. For any gene in a module, the module membership (MM) of this gene in the module is defined as the correlation coefficient between this gene and the characteristic gene of this module.

Multiple machine learning algorithms were utilized to identify the diagnostic feature biomarkers. The least absolute shrinkage and selection operator (LASSO) was a regression analysis algorithm and this algorithm obtained a refined model by constructing a penalty function: by finally determining that the coefficients of some indicators are zero, the LASSO algorithm achieves the purpose of reducing the set of indicators. Support vector machine (SVM) was a common discrimination method usually used for pattern recognition, classification and regression analysis. Moreover, it also was a supervised learning model in the field off machine learning. To avoid overfitting and achieve reliable accuracy, recursive feature elimination (RFE) algorithm was used to select the optimal genes from the training cohort. Therefore, SVM-RFE was performed to select the appropriate feature biomarkers. Finally, the overlapping genes of WGCNA, LASSO, RFE, and SVM-RFE was considered as diagnostic feature biomarker and the expression levels of candidate genes were further validated in the validation cohort.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to identify potential functional components and pathways using “clusterProfiler” R packages (12). Gene Set Enrichment Analysis (GSEA) was utilized to enrich the KEGG terms in the HC and DN group, with P < 0.05 was considered as statistically significant.

“CIBERSORT” R package was performed to investigate the immune infiltration landscape of HC and DN samples. The gene expression of immune cells based on “CIBERSORT R script v1.03” and the algorithm was run using the LM22 signature for 1000 permutations. Correlation analysis was performed to determine the relationships between immune cells via “Corrplot” R package and “ggplot2” R packages was utilized to determine the difference of immune cell between HC and DN samples.

The correlation of the 4 diagnostic feature biomarkers and immune cells were estimated via Spearman’s rank correlation algorithm. The correlation results were visualized in lollipop diagram via “ggplot2” R packages, and P-value < 0.05 were considered as statistically significant.

All statistical analyses were performed using R software (version 4.1.0) and Perl scripts. Correlation analyses between the two variables were performed using the Spearman’s rank correlation algorithm and P-value < 0.05 was considered significantly different. Differential functions were analyzed using the Wilcoxon rank-sum test between the two groups, and statistical significance was set at P-value < 0.05.

The workflow of this study was illustrated in the Figure 1 . A total of 20 HC samples and 41 DN samples were extracted from the GEO database (GSE30122, GSE96804). Under the threshold set at |Fold change| ≥ 1 and P-value < 0.05, a total of 314 down-regulated and 135 up-regulated DE-MRGs were identified ( Figure 2A ). Heatmap diagram displayed the top 30 up-regulated and down-regulated DE-MRGs in HC and DN group ( Figure 2B ). The principal component analysis (PCA) showed a clear separation between the HC and DN group based on the MRGs ( Figure 2C ).

To investigate the potential molecular mechanism of the DE-MRGs, the functional enrichment analysis was conducted. Biological process (BP) enrichment analysis showed that the DE-MRGs were significantly enriched in small molecule catabolic process, ribose phosphate metabolic process, and organic acid catabolic process ( Figure 3A ). Kyoto Encyclopedia of Genes and Genome (KEGG) enrichment analysis suggested that carbon metabolism, purine metabolism, and glycolysis/gluconeogenesis were enriched of the DE-MRGs ( Figure 3B ). GSEA results illustrated the top 5 KEGG signaling pathways in the HC and DN group. As shown in the Figures 3C, D , the calcium signaling pathway, GNRH signaling pathway, linoleic acid metabolism, long term depression and melanogenesis were enriched in the HC group, whereas arginine and proline metabolism, citrate cycle TCA cycle, fatty acid metabolism, propanoate metabolism, and valine leucine and isoleucine degradation were significantly enriched in the DN group.

Multiple machine learning algorithms were utilized to screen the diagnostic feature biomarkers in DN. A gene co-expression network was established according to weighted gene co-expression network analysis (WGCNA). The power of β = 6 (scale-free R² > 0.85) was selected as the soft-thresholding parameter to construct a scale-free network, and the correlation coefficient of the module eigengenes and disease characteristics were calculated ( Figure 4A ). The results of the module trait showed the correlation of the module eigengenes and the disease characteristics. Of note, the module blue was positively correlated with DN, and was selected for the subsequent analysis ( Figure 4B ). Based on the LASSO algorithm, 26 DE-MRGs were identified as diagnostic feature biomarkers for DN ( Figures 4C, D ). 16 DE-MRGs were identified as diagnostic biomarkers from the DE-MRGs using the SVM-RFE algorithm ( Figure 4E ). Additionally, 11 diagnostic feature biomarkers were selected from the DE-MRGs using the RFE algorithm ( Figure 4F ).

According to the four machine learning algorithms, 4 diagnostic feature biomarkers were identified ( Figure 5A ). The expression of the 4 diagnostic feature biomarkers in the training cohort showed that ADI1 and PTGS2 were expressed higher in the HC group, however, the DN group showed a higher expression of DGKH and POLR2B ( Figure 5B ). In validation cohort, the expression of ADI1 and PTGS2 were higher in the HC group, whereas the expression of DGKH and POLR2B were lower in the DN group compared to the HC group ( Figure 5C ). The ROC curve of the ADI1, PTGS2, DGKH, and POLR2B revealed a satisfactory diagnostic value for the 4 diagnostic feature biomarkers in the training cohort with AUC was 0.954, 0.994, 0.993, and 0.926, respectively ( Figure 5D ). The AUC of ADI1, PTGS2, DGKH, and POLR2B in the validation cohort was 0.893, 1, 1, and 0.996, respectively ( Figure 5E ). Overall, these above findings demonstrate a satisfactory diagnostic value of the 4 diagnostic feature biomarkers and could be applied to the clinical diagnosis of DN.

In addition to metabolic disorders, immune system disorders are another important aspect of diabetic nephropathy (13). To investigate the immune infiltration landscape of patients, CIBERSORT algorithm was performed and the fraction of 22-types immune cells were assessed ( Figure 6A ). The correlation analysis revealed a remarkable correlation in the 22-types immune cells ( Figure 6B ). T cells CD8 was positively correlated with T cells CD4 naïve, T cells CD4 memory resting, and T cells CD4 + memory; B cells memory was positively correlated with B cells naïve; neutrophils was negatively correlated with T cells CD4 + memory, B cells naïve, B cells memory, and Monocytes; T cells CD4 memory resting was negatively correlated with macrophages M2. As shown in Figure 6C , the violin diagram showed that the fraction of B cells naïve, B cells memory, T cells CD8, T cells CD4 naïve, T cells CD4 memory resting, Macrophages M0, and mast cells resting were significantly higher in the DN group than HC group. However, the HC group had a higher proportion of macrophages M2 and neutrophils.

The correlation analysis of the diagnostic biomarkers and immune infiltration landscape was further investigated. POLR2B was positively correlated with B cells naïve, B cells naïve, mast cells resting, T cells CD4 memory resting, macrophages M0, T cells CD4 memory activated, T cells gamma delta, and T cells CD8, but negatively correlated with macrophages M2 and neutrophils ( Figure 7A ). PTGS2 was positively correlated with NK cells activated, but negatively Macrophages M0 ( Figure 7B ). DGKH was positively correlated with B cells naïve, Mast cells resting, B cells memory, T cells CD4 memory resting, Macrophages M0, but negatively correlated with Macrophages M2 and Neutrophils ( Figure 7C ). In addition, ADI1 was positively correlated with macrophages M2, neutrophils, dendritic cells activated, but negatively correlated with macrophages M1, B cells memory, T cells CD4 memory resting, T cells CD4 memory activated, T cells gamma delta, T cells CD8, NK cells activated, monocytes, T cells follicular helper, macrophages M0, T cells CD4 naïve, and B cells naïve ( Figure 7D ). Taken together, these results demonstrate that the diagnostic feature biomarkers are correlated with immune infiltration landscape, providing a fresh insight for the future research in DN.