Based on the public and well-known large-scale database, GEO (https://www.ncbi.nlm.nih.gov/geo/), we searched and collected several data sets including patients with DN (Edgar et al., 2002). The inclusion criteria for data sets were as follows: 1) Organism: Homo sapines, 2) Tissue: Glomerulus, 3) datasets containing corresponding normal (CN) samples and comprehensive introduction for each sample. Moreover, the detailed introduction of the GEO data sets included in our research was exhibited in Table 1. In addition, the necroptosis-related regulators (NRGs) were downloaded from the public website GeneCards (https://www.genecards.org/), in which the selection criteria was the correlation score >1 and searched with the key word “necroptosis” (Safran et al., 2010).

To investigate the expression difference of genes in DN and explore the expression characteristic at the mRNA level, we utilized differential expression analysis via the “limma” package in R software, which was utilized to obtain the differentially expressed genes (DEGs) in our training set between CN and DN samples, respectively (Ritchie et al., 2015). The p < 0.05 were used to define the threshold of DEGs. To visualize the results of differential expression analysis, the “ggplot2” package in R software was used to create a volcano plot and a heat map.

The WGCNA was performed to screen the significant gene clusters correlated with DN (Langfelder et al., 2008). The “gplots” package in R software was utilized for hierarchical cluster analysis, and abnormal samples or values were deleted. The gene correlation between samples was calculated via WGCNA algorithm. Then we selected the best soft threshold power and established a standard non-proportional network. The WGCNA model or network was related to the characteristics of external samples. Via a dynamic tree-cutting strategy, various different modules named by different colors were constructed through the hierarchical clustering of genes.

Several machine learning methods were used to investigate the diagnostic efficiency of the regulators. Random Forest (RF) can evaluate the significance of each feature in classification problems (Ishwaran and Kogalur, 2010; Wang and Zhou, 2017). Support vector machine recursive feature elimination (SVM-RFE), analgorithm that adds or removes features to obtain the optimal combination variable that maximizes model performance for specific feature variables. We use the “rfe” function in the R package “caret” for feature recursive elimination, set functions = caretFuncs, and obtain the best model. Moreover, least absolute shrinkage and selection operator-cox (Lasso-cox) regression analysis is an algorithm to evaluate the influence for binary outcome via merging some variables using the R package “glmnet” (Tang et al., 2017).

According to the expression level of the selected significant necroptosis-regulators, we divided all the 41 DN samples in our training set into different necroptosis-related subtypes via the “ConsensusClusterPlus” package in R software (Wilkerson and Hayes, 2010). Then, with the most appropriate k-value, two different subtypes were identified.

Moreover, to exhibit the heterogeneity between the above subtypes, the PCA algorithm was performed (Diaz-Papkovich et al., 2021). And then, the necroptosis score was calculated with the results of PCA algorithm to explore the biological background of the two subtypes: Necroptosis Score = Σ _i (PCA1 _i + PCA2 _i ).

We utilized the ClueGO plug-in in the Cytoscape software with the threshold (p < 0.05) to conduct gene ontology (GO) enrichment analysis (Bindea et al., 2009). And the Molecular complex detection (MCODE) plug-in in cytoscape was used to extract several significant models in the gene-pathway interacted network of the results of ClueGO. Moreover, the R software package “clusterProfiler” was selected to perform the KEGG and Reactome enrichment analysis (Yu et al., 2012).

In addition, the gene set enrichment analysis (GSEA) was performed to comprehensively explore the biological background between the two necroptosis-related subtypes of DN. With the standard gene set of Hallmark, in the analysis, the size of the gene set is limited to 5–5,000 genes. Results with a p-value less than 0.05 were considered statistically significant. (Hacisalihoglu et al., 2016).

Similarly and firstly, the differential expression analysis was performed between the two different necroptosis-related subtypes to obtain the DEGs (Hao et al., 2023). To select several more prominent genes to construct the diagnostic model, the criteria for identifying DEGs was the absolute value of log₂ fold change (|log₂ FC|) > 1 and p < 0.05.

Then, the Search Tool for the Retrieval of Interacting Genes (STRING, https://cn.string-db.org/) was used to construct the protein-protein interaction (PPI) network of our candidate genes and the cypscape software was utilized to visualized the PPI network (Shannon et al., 2003; Szklarczyk et al., 2023). Furthermore, we used the cyttoHubba plug-in in cytoscape software to screen the significant features from our candidate genes (Chin et al., 2014). Then, four machine learning methods: SVM-RFE, RF, univariate Cox regression analysis and Lasso-Cox regression analysis were conducted to further construct our diagnostic model. Finally, we performed multivariate logistics regression to establish our diagnostic model: Diagnostic Score = Σ _i (Coefficient _i * Expression level of feature _i ).

Single sample gene set enrichment analysis (ssGSEA) is a novel type of gene set variation analysis which is a method of unsupervised clustering based on a specific gene set to evaluate the score of each sample. In our study, ssGSEA was used to calculate the infiltration score of 28 different scores in DN. We collected genes from previous studies related to 28 different types of immune cells and predicted the infiltration abundance of immune cells based on these genes (Supplementary Table S1) (Barbie et al., 2009; Charoentong et al., 2017). Moreover, based on our expression matrix, we used the “IOBR” package in R software to select the ESTIMATE method to calculate the ESTIMATE Score of each sample in our training set (Yoshihara et al., 2013; Zeng et al., 2021). The CYT Score was calculated via averaging the expression value of GZMA and PRF1 in each sample to evaluate the inflammatory status. Furthermore, Mantel correlation analysis was performed to explore the correlated relationship between the significant necroptosis-related regulators and the infiltration score of each immune cell (Sunagawa et al., 2015).

The Study: Human Diabetic Kidney: Wilson et al., PNAS 2019. (*Note: The CFH + cluster was renamed to parietal epithelial cells) in the public website: Humphreyslab (http://humphreyslab.com/SingleCell/) was included in our research to explore the expression trend of significant necroptosis-regulators in DN (Wilson et al., 2019).

Male 8-week-old diabetic db/db (BKS-Leprem2Cd479/Nju; n = 8) and their nondiabetic wild type (C57BL/Ksj; n = 8) mice were purchased from the Model Animal Research Center of Nanjing University (Nanjing, China). All mice were housed under specific pathogen-free conditions with controlled temperature and humidity (22°C ± 2°C, 50% ± 5% RH) and a standard 12-h light/dark cycle. All experimental procedures were approved by the Institutional Animal Care and Use Committee at Nanchang University (NO. 0064257) and performed in accordance with the guidelines for the ethical treatment of laboratory animals. All mice were randomly divided into two groups: control (CN, n = 8) group and diabetic nephropathy (DN, n = 8) group. The successful establishment of the diabetic model was defined as random blood glucose levels ≥16.7 mmol/L. At 20 weeks of age, the mice were anesthetized with isoflurane and subsequently euthanized. Postmortem, the body weight, kidney weight, and tibial length were measured. The kidneys were then excised and divided into sections for further analysis.

The total RNA of serum was extracted from fresh glomerulus samples utilizing Trizol reagent (Ambion, Singapore). The quantity and purity of the total RNA were determined using the Nanodrop®ND1000 (TIANJIN). The total RNA was reverse transcribed into complementary DNA (cDNA) employing the Script cDNA synthesis kit (TianGEN). Quantitative PCR reaction was performed on an EDC-810 Real-Time PCR Detection Instrument (Eastwin life, China), using SYBR Premix ExTaq kit (Takara, China), All genes were analyzed Relatively and calibrated to the expression of control groups. The 2^{^−ΔΔCt} method was employed to calculate the relative RNA expression values. The primer pairs used for the amplification of target genes are listed in Supplementary Table S2. All the above experiments were strictly carried out in accordance with the instructions.

Serum blood urea nitrogen (BUN) levels were quantified in murine samples using a commercial colorimetric assay kit (DEIABL-M3) per the manufacturer’s protocol. An enzymatic creatinine reagent kit (DEIABL-M4S) was employed to determine serum creatinine concentrations in the sourced murine serum samples. Urine specimens harvested from the mice underwent analysis of microalbumin content by enzyme-linked immunosorbent assay utilizing a microalbumin-specific antibody (ab108792). Murine urine creatinine concentrations were evaluated by a commercial colorimetric assay kit (LS-F13025) following the provided manufacturer guidelines.

The data were shown as mean ± SD. R software (version 4.2.2) and its related software packages were used to process and analyze data. The unpaired t-test was utilized to analysis the statistical difference between two different groups. p < 0.05 were considered statistically significant. And the “RMS” package in R software was used to visualize Nomograms. The receiver operating characteristic (ROC) curves were visualized via the Sangerbox website (http://vip.sangerbox.com/home.html).

Firstly, to exhibit the remarkable characteristic of necroptosis in the transcriptomic level in DN. Initially, we conducted a differential expression analysis comparing the NC group and DN group. Genes with a p-value less than 0.05 were identified as differentially expressed genes (DEGs), as depicted in Figure 2A. A total of 10,779 DEGs were identified, with 5,159 genes being upregulated and 5,620 genes being downregulated. The detailed expression status of these DEGs is displayed in Figure 2B, providing a comprehensive overview of their differential expression patterns. Then, the location in chromosomes of the top 50 downregulated and upregulated genes were listed in Figure 2C with the corresponding p-value and normalized expression value. The primary NRGs in DN were obtained via intersecting the upregulated DEGs, downregulated DEGs, and NRGs in Genecards. Finally, a total of 45 necroptosis-related DEGs (NRDEGs) were identified (Figure 2D).

Furthermore, the WGCNA was performed to extract several gene modules significantly correlated with DN (Figures 2E–G). We intersected the above NRDEGs and the genes in the significant modules of WGCNA to ensure several NRDEGs prominently correlated with DN (Figure 2H). As a result, a total of 19 NRDEGs were identified. Moreover, a Spearman correlation analysis was performed among the above 19 NRDEGs and almost all regulators indicated a remarkably positive correlation (Figure 2I).

According to the above general exploration of the role of necroptosis-related regulators in DN, we subsequently excavated several significant NRDEGs in DN to evaluate the crucial characteristic of necroptosis in DN. To begin with, we performed the Lasso regression analysis to obtain the key features in distinguishing DN and NC (Figures 3A, B). Moreover, another ML method, SVM-RFE was performed to extract the features with prominent diagnostic efficiency. As a result, we obtained 10 NRDEGs (Figures 3C, D). Via the cross validation of a 10×fold, a total of 12 NRDEGs were identified. Furthermore, RF was utilized to assess the importance of 19 NRDEGs in recognizing NC and DN patients. With the importance >25 as the criteria, 10 NRDEGs were selected (Figure 3E). Finally, after intersecting results of the above 3 different ML methods, a total of 8 NRDEGs (CFLAR, FMR1, GSDMD, IKBKB, MAP3K7, NFKBIA, PTGES3, and SFTPA1) were finally regarded as the significant NRDEGs (Figure 3F). The co-expression network of the eight significant NRDEGs was construction with co-expression of 61.91%, physical interactions of 28.79%, shared protein domains of 7.12%, and pathway of 2.17%, wherein the enrichment analysis revealed that the 8 regulators was mainly associated with the inflammatory and immunal pathways, such as “pattern recognition receptor signaling pathways” “toll-like receptor signaling pathways”, and “I-kappaB kinase/NF-kappaB signaling” (Figure 3G).

Moreover, to better demonstrate the diagnostic characteristic of the above 8 regulators, the ROC curves were exhibited, the area under curve (AUC) of which were 0.96 0f IKBKB, 0.95 of GSDMD, 0.91 of PTGES3, 0.94 of SFTRA1, 0.91 of CFLAR, 0.88 of FMR1, 0.90 of MAP3K7 and 0.81 of NFKBIA, which emphasized the powerful diagnostic efficiency (Figures 3H–O). Meanwhile, the heat map exhibited the expression level of the above 8 features between CN and DN (Figure 3P). Furthermore, utilizing the expression levels of the eight features in GSE96804, we developed a nomogram to visualize the predictive ability of our model for diagnosing DN patients (Figure 3Q). Additionally, the calibration curve was generated to assess the accuracy and calibration of our model (Figure 3R).

According to the expression value of the above 8 significant NRDEGs (IKBKB, GSDMD, PTGES3, SFTRA1, CFLAR, FMR1, MAP3K7, and NFKBIA) in DN samples, all the 41 DN samples in our training set were classified into two various necroptosis-related subtypes (C1 and C2) with the most compatible K-value (K = 2) using a consensus cluster analysis (Figures 4A–C). Moreover, the PCA-gram exhibited an obvious differentiation between the above two DN subtypes (Figure 4D). In addition, the detailed expression level of all our 8 crucial NRDEGs between the two DN subtypes were exhibited in the heat map, which emphasized the heterogeneity of our two subtypes (Figure 4E). Meanwhile, via the PCA algorithm, we calculated the necroptosis-score to objectively recognize the two necroptosis-related subtypes, in which subtype C1 was prominently lower than subtype C2 (p < 0.05, Figure 4F).

According to the above analysis, the significant difference between the two necroptosis-related subtypes has been exhibited. Thus, we further investigated the biological differential characteristic between the two subtypes. First, we performed differential expression analysis to obtain the DEGs between the two subtypes (C2 vs. C1). As a result, 121 DEGs were identified based on the criteria of p < 0.05 and |log₂FC| > 1 (Figure 4G). Similarly, differential expression analysis was performed in the whole training set to determine DEGs between CN and DN group, which results in a total of 616 DEGs, including 290 upregulated and 326 downregulated. Moreover, to select some candidate features with potential diagnostic efficiency, we intersected the DEGs between two necroptosis-related subtypes and the DEGs between CN and DN group (Figure 4H).

Subsequently, to better exhibit the biological characteristic between the two subtypes, we performed GO and KEGG enrichment analysis for the 58 candidate features. Via ClueGO plug-in in Cytoscape software, the GO enrichment analysis revealed that pathways related to metabolism and immunity played a crucial role in necroptosis-related subtypes, especially “hemoglobin complex” with the proportion of 38.46% and “complement activation, alternative pathway” with the proportion of 23.08% (Supplementary Figure S1A). Moreover, the MCODE plug-in was utilized for the extract of some significant functional clusters and these clusters were mainly related to metabolism and immunity (Supplementary Figures S1B–D). Meanwhile, the results of KEGG and Reactome enrichment analysis proved the crucial role of immune landscape in the two necroptosis-related subtypes (Supplementary Figures S1E, F). Moreover, the GSEA enrichment analysis was also performed to further explore the functional backgrounds between the two subtypes, the results of which emphasized the role of “DNA Repair”, “TGF-β Signaling”, and “Heme Metabolism” (Supplementary Figures S1G–I).

According to the crucial role of significant NRDEGs in DN and the lack of effective diagnostic markers or models, we exploited a novel necroptosis-mediated diagnostic model for DN patients. Via the above rounds of selection, we finally ensured 58 genes for the construction of our model. Firstly, all the 58 features were imported into the STRING website to establish a PPI network (Figure 5A). Then via the cyttoHubba plug-in in Cytoscape software, we utilized 7 algorithms to obtain the top 15 features with highest interaction score in each algorithm (MCC, MNC, Degree, EPC, Closeness, Radiality, Stress). After intersecting the 10 features in each algorithm, a total of 5 features were identified in all the 7 algorithms, which were selected for the subsequent analysis (Figure 5B).

Furthermore, we utilized 4 different ML methods to select several significant features. Originally, one classical ML methods, SVM-RFE, was performed to screen features with more accurate diagnostic efficiency. As a result, all the 5 features were regarded as regulators with significant diagnostic efficiency (Figures 5C, D). Then, we utilized an additional ML method (RF) to calculate the importance of all the 5 features in recognizing the DN and CN samples. The lollipop chart revealed the importance in detail and 4 features with importance >4 were regarded significant features (Figures 5E, F). Moreover, to ensure the accuracy of our selection, we further performed another ML method, lasso-cox regression analysis to extract the features, which can better distinguish DN and CN patients. Five features were considered to play a key role in the diagnosis of DN (Figures 5G, H). In addition, the univariate-cox regression analysis was also conducted to select the significant features with the criteria: p < 0.05 (Table 2). Finally, the 4 features (COL6A3, FBLN5, LUM, and VCAN) intersected in the whole 4 ML methods were regarded as the source elements of the construction of our diagnostic model (Figure 5I).

Subsequently, for establishing our diagnostic model for DN, the logistic regression analysis was performed: Diagnostic Score = (12.1605262153362 × COL6A3) + (-3.03974790703008 × FBLN5) + (2.42256407679845 × LUM) + (-3.13831314409003 × VCAN). And the calibration curve revealed that there was no potential difference between our model and the ideal situation (p = 0.985, Figure 5J). In our training set (GSE96804) we performed unpaired t-test to inspect the statistical significance between CN and DN group based on the formula of the diagnostic score described above. As a result, our necroptosis-mediated diagnostic model revealed a statistically significant difference in the training set (p < 0.0001, Figure 5K). Meanwhile, the AUC of our model in the training set was 0.97, which emphasized the powerful diagnostic efficiency (Figure 5L). Additionally, to strengthen the reliability of our diagnostic model, we chose two external verification sets (GSE30528 and GSE30122) to validate our model. Interestingly, in GSE30528 and GSE30122, the statistical significance between CN and DN group both occurred (Figures 5M, O, Q). And the AUC was 0.86 the glomerulus tissue in GSE30122, 0.83 the tubules tissue in GSE 30122 and 0.86 in GSE30528, which indicated the accuracy of the diagnostic efficiency of our model (Figures 5N, P, R).

Based on the above functional exploration of the two necroptosis-related clusters, the immune microencironment was considered to keep a critical role in DN and corresponding clusters. Two different algorithms, ESTIMATE and ssGSEA, were utilized to comprehensively explore the immune microenvironment of DN. Via the ESTIMATE algorithm, the ESTIMATE Score was lower in C1 than that in C2, which revealed there were potentially different immune response and infiltration between the above two necroptosis-related subtypes of DN (Figure 6A). Then, via ssGSEA, we visualized the superficial landscape of various immune cells infiltrations in all samples of our training set, which exhibited the prominent difference of immune microenvironment between CN and two different DN clusters (Figure 6B). And between DN and CN, abundant immune cells exhibited significant infiltration status (Figure 6C). Moreover and intriguingly, we concuded the Mantel correlation analysis to compute the correlation coefficient between the expression level of our 8 significant NRDEGs and the immune infiltration score. As a result, the necroptosis score and expression of 8 NRDEGs in C1 was more correlated with the immune microenvironment than in C2, further proving the result of ESTIMATE (Figures 6D, E).

According to the above investigation of the immune landscape in two different clusters, we ensure the diversified immune microenvironment between C1 and C2. And based on the significant role of immune checkpoints (ICKs) in the regulation of immune response, we further explore the role of ICKs in DN. The unpaired t-test was used to detect the expression of 14 classical ICKs between C1 and C2. Interestingly, a total of 8 ICKs exhibited a significantly different expression in C2 than C1 (). Moreover, the CYT Score calculated via the average of the expression of GZMA and PRF1 also exhibited the different inflammatory status between the two necroptosis-related subtypes of DN (Figure 6T). Finally, the Sankey plot demonstrated the potential relationship between the two necroptosis-related subtypes, the immune microenvironment and the inflammatory status (Figure 6U).

To better describe the expression pattern of our 8 significant NRDEGs in DN, we selected a data set including 12 different cell clusters (PCT, CD-ICA, PEC, DCT, DCT/CT, CD-PC, CD-ICB, MES, LEUK, ENDO, LOH, PODO, n = 23,980). The general landscape of the above 12 cell clusters in different samples was exhibited in Figure 7A. Interestingly, the expression level of our 8 significant NRDEGs (CFLAR, FMR1, GSDMD, IKBKB, MAP3K7, NFKBIA, PTGES3, and SFTPA1) was different in the different cells between the DN and CN group, which revealed that our 8 significant NRDEGs played a crucial and various role in the characteristic of DN (Figures 7B–Q).