The study’s flow chart is shown in Figure 1. We conducted a systematic search of the Gene Expression Omnibus database (GEO) (http://www.ncbi.nlm.nih.gov/geo/) using the terms: “Homo sapiens” and “NAFLD”. This study included datasets that met the following eligibility criteria: 1) the type of study was expression profiling by array; 2) the dataset contained liver tissue samples from NAFLD patients and controls; 3) the sample size was greater than 15; 4) the raw data or array gene expression profiles were available in the GEO. Finally, five datasets were identified from the GEO database. GSE89632 (Microarray, platform GPL14951) (Arendt et al., 2015), GSE48452 (Microarray, platform GPL11532) (Ahrens et al., 2013), GSE66676 (Microarray, platform GPL6244) (Xanthakos et al., 2015), GSE63067 (Microarray, platform GPL570) (Frades et al., 2015) and GSE164760 (Microarray, platform GPL13667) (Pinyol et al., 2021)were included in this study.

First, we combined the datasets GSE89632, GSE48452, and GSE66676, which had 72 normal samples and 104 NAFLD samples. Then, the merged gene expression datasets were normalized by using the “sva” package (Leek et al., 2012). DEGs between the NAFLD and control groups were discovered using the “limma” package (Colaprico et al., 2016). p-values lower than 0.05 were regarded as statistically significant. The datasets GSE63067, including 14 normal samples and 32 NAFLD samples, and GSE164760 including 6 normal samples and 74 NAFLD samples, were used to be the validation set. In the Liu et al. research, ten genes were connected to disulfidptosis (Liu et al., 2023). GYS1, NDUFS1, OXSM, LRPPRC, NDUFA11, NUBPL, NCKAP1, RPN1, SLC3A2, and SLC7A11 are among these genes (Liu et al., 2023).

To determine the relative abundances of 22 immune cells in each sample, the CIBERSORT method was used (Supplementary Table S1) (Zhao et al., 2023). The sum of the 22 immune cells proportions in each sample was 1, and p < 0.05 represented a significant correlation. We assessed the correlations between the proportional percentage of immune cells and the expression of the DRGs using the spearman correlation coefficient. Histograms, heat maps, and box plots were plotted using the “ggplot2” R packages (version 0.92).

Based on DE-DRGs, unsupervised clustering analysis of NAFLD patients was performed by using the R package “ConsensusClusterPlus” (version 2.60) (Wilkerson and Hayes, 2010). Cumulative distribution function (CDF) curves, consensus matrix, and consistent cluster score were used to estimate the optimal cluster number.

The “c2. cp.Kegg.symbols” and “c5. go.symbols” files from the MSigDB database were first downloaded for this investigation. The enrichment differences between GO and KEGG pathways were then clarified using a non-parametric unsupervised gene set variation analysis (GSVA) approach by the R package of “GSVA” (version 2.11) (Hänzelmann et al., 2013). p-value lower than 0.05 was the cutoff value for the statistical significance term.

The “clusterProfiler” package (version 3.16.1) was used to perform a Gene Set Enrichment Analysis (GSEA), which was then used to determine whether KEGG pathways had an abundance of relevant genes (Kumar and Futschik, 2007).

Utilizing the “WGCNA” R package (version 1,70.3), WGCNA analysis was carried out on the training set to look into the connection between gene networks and diseases (Langfelder and Horvath, 2008). To ensure the accuracy of the study, we first grouped the samples and eliminated outliers. A soft threshold from 1 to 20 was used for topology calculation to determine the optimal soft threshold. When the minimum module size was set to 100, a “dynamic tree cutting” algorithm was used to group genes with similar patterns into modules. Finally, Pearson correlation analysis was performed to calculate the correlation between modules and traits. Based on the correlation between modules and clinical features, the most relevant module to the disease is selected as the key module.

Through the use of the “caret” R package (version 6.0.91), we constructed four machine learning models by using data from two distinct DRGs clusters. These models were the random forest model (RF), the support vector machine model (SVM), the generalized linear model (GLM), and the eXtreme Gradient Boosting (XGB) (Nelder and Wedderburn, 1972; Gold and Sollich, 2003; Chen and Guestrin, 2016; Rigatti, 2017). The distinct clusters were considered as the response variable, and the cluster-specific DEGs were selected as explanatory variables. 104 NAFLD samples were randomly classified into training and testing cohort. All these machine learning models were using default parameters and evaluated by 5-fold cross-validation. The caret package (Long and Yang, 2020) was used to automatically adjust the parameters of these models. The aforementioned four machine learning models are analyzed by using the “DALEX” R package (version 2.4.0) to illustrate the relevance of features across these machine learning models as well as the distribution of residuals. The “pROC” R package (version 1.18.0) was used to visualize the area under the ROC curve. As a result, we determined the top machine-learning models for NAFLD and selected significant prognostic genes linked to NAFLD for additional research. The diagnostic utility of this diagnostic model was subsequently validated using ROC curve analysis on the GSE63067 and GSE164760 datasets.

Using the R package “rms” (version 6.2.0) and the screened model genes, we built a diagnostic nomogram model and then used calibration curves and DCA to evaluate the nomogram model’s predictive ability.

Twenty-four 6-week-old male C57BL/6J mice were kept in conventional housing (room temperature: 23°C ± 2°C; 12-h light/dark cycle) with unrestricted access to food and water. Twelve mice each were randomly assigned to the normal chow (NC) and high-fat diet (HFD) groups after 1 week of acclimatization, and each group received either the high-fat diet (60 kcal% fat; d12492, medicine, Jiangsu, China) or standard laboratory food. After 12 weeks of HFD feeding, a mouse model with non-alcoholic fatty liver was created (Dusabimana et al., 2021; Qian et al., 2021). All mice were anesthetized using 2% isoflurane after 12 weeks. After mice were sacrificed by cervical dislocation, blood samples and liver tissues were collected. All studies were approved by the Professional Committee for Animal Protection of Harbin Medical University (2022-DWSYLLCZ-20).

For histological analysis, formalin-fixed mouse liver tissues were processed and 5-μm-thick paraffin sections were cut and stained with hematoxylin-eosin (H&E). Triglyceride (TG) content in the liver was measured by a commercial kit (A110-1-1; Jiancheng, Nanjing, China) according to the manufacturer’s instructions.

TRIzol reagent (Invitrogen, Carlsbad, CA, USA) was used to extract total RNA from homogenized tissues. Then, 2*SYBR Green qPCR (Vazyme, Nanjing, China) was used to analyze gene expression after 1 ug of total RNA was reverse transcribed using PrimeScript reverse transcriptase (Takara, Kusatsu, Japan). The expression of β-actin was taken as the control. In Supplementary Table S2, the primer sequences are presented.

R software (version 4.2.0) was used to conduct bioinformatics analyses. The software GraphPad Prism version 9.0 was used to statistically analyze and visualize the data from animal experiments. Using an unpaired Student’s t-test, the means of two groups of normally distributed variables were compared. Data are presented as the mean ± standard deviation. p < 0.05 were considered significantly different.

First, we merged and normalized three GEO datasets: GSE89632, GSE48452, and GSE66676 (Supplementary Figures S1, S2). Subsequently, we investigated the expression patterns of the 10 DRGs in control and NAFLD samples (Figure 2A), and the results showed that NCKAP1 was significantly highly expressed in NAFLD samples, while SLC3A2 was lowly expressed in NAFLD samples (Figure 2B). Also, the chromosomal positions of the 10 DRGs were visualized (Figure 2C).

Previous research indicates a direct connection between NAFLD and the immunological microenvironment (Huby and Gautier, 2022). To investigate the variations in the immune microenvironment between NAFLD patients and normal samples, we applied the CIBERSORT algorithm. As seen in Figure 2D, compared to normal samples, NAFLD samples had reduced proportions of monocytes, activated dendritic cells, activated mast cells, eosinophils, and neutrophils, as well as higher proportions of M1 macrophages, M2 macrophages, resting dendritic cells, and resting mast cells. This raises the possibility that immune system changes may be the primary cause of NAFLD. Additionally, Figure 2E’s correlation analysis revealed that NCKAP1 was considerably negatively connected with memory B cells, M2 macrophages, and regulatory T cells, and significantly positively correlated with activated dendritic cells. While SLC3A2 was strongly inversely connected with resting mast cells and δγT cells, it was considered positively correlated with neutrophils and monocytes. This data suggests that the altered immune milieu of NAFLD may be significantly influenced by NCKAP1 and SLC3A2, which act as disulfidptosis molecules in NAFLD patients.

To elucidate the expression patterns associated with DRGs in NAFLD, we grouped the NAFLD samples based on DE-DRGs using a consensus clustering algorithm. When k = 2, examination of the Delta area and the CDF value showed that the results of clustering was relatively stable (Figure 3A, B). The highest consensus values (all over 0.9) were observed for each subtype when k = 2 (Figure 3C). The consensus matrix showed that each sample in the cluster exhibits strong correlation when k = 2, when the samples in the two subtypes are most stable (Figure 3D). Therefore, k = 2 was the best choice. To verify the differences between Cluster 1 and Cluster 2 samples, PCA analysis showed significant differences between these subtypes (Figure 3E). Finally, based on the expression of DE-DRGs, we categorized the 104 NAFLD samples into two different subtypes, including Cluster 1 (n = 69) and Cluster 2 (n = 35).

To investigate the molecular differences between clusters, we completely analyzed the DRGs expression differences between Cluster 1 and Cluster 2. We found that the two clusters had distinct DRGs expression landscapes (Figure 4A). Meanwhile, we analyzed the differences in DRGs between different DRGs clusters and found that SLC3A2 and NDUFA11 were upregulated in Cluster 2 and NCKAP1 was downregulated in Cluster 2 (Figure 4B). We then examined the variations in immune cells and their immune activities across various clusters to further study the differences in immune microenvironment features between various DRGs clusters. The findings demonstrated that the immunological microenvironment between Cluster 1 and Cluster 2 was different. (Figure 4C). Cluster 1 had relatively high levels of resting mast cells and γδ T cells; whereas monocytes had higher levels of abundance in Cluster 2 (Figure 4D).

We also performed a GSVA analysis to explore the functional differences between the two clusters. Functional enrichment results showed that non-homologous end-joining, oocyte meiosis, propanoate metabolism, and regulation of autophagy pathways were enhanced in Cluster 2, while the Adipocytokine signaling pathway, JAK-STAT signaling pathway, and amino acid metabolic pathway were upregulated in Cluster 1 (Figure 4E). The results showed that Cluster 2 was closely associated with type III interferon production, asymmetric cell division, and regulation of protein localization to the cilium. However, the neutrophil apoptotic process, positive regulation of dephosphorylation, defense response to the bacterium, and response to lipids were enriched in Cluster 1 (Figure 4F). Therefore, we hypothesize that Cluster 1 may be involved in various immunoregulatory and energy metabolic responses.

In addition, we analyzed key gene modules closely associated with disulfidptosis clusters using the WGCNA algorithm. All samples were clustered in the dataset and no samples were excluded (Supplementary Figures S3, S4). The value of the optimal soft power for the training set according to the WGCNA method was 14 (Figure 5A), while for the samples related to the disulfidptosis clusters was 6 (Figure 5B). And based on the co-expression network with the optimal soft power, 4 modules were identified in the training set, including MEblue, MEbrown, MEgrey, and MEturquise (Figure 5E) and visualized by hierarchical clustering (Figure 5C). In the disulfidptosis clusters, 9 modules including MEblack, MEblue, MEbrown, MEgreen, MEgrey, MEpink, MEred, MEturquise and MEyellow were identified (Figure 5F) and visualized by hierarchical clustering (Figure 5D). The blue module showed the strongest association with NAFLD. The blue module had the strongest correlation with the disulfidptosis clusters. By analyzing the interaction between module-associated genes of disulfidptosis clusters and module-associated genes of NAFLD and non-NAFLD people, a total of 67 cluster-specific DEGs were discovered (Figure 5G; Supplementary Table S3).

Based on the expression profiles of 67 cluster-specific DEGs, we constructed four machine-learning models to further uncover subtype-specific genes with high diagnostic values. The outcomes demonstrated that the residuals for the SVM and RF machine learning models were reasonably low (Figures 6A, B). By computing the receiver operating characteristic (ROC) curves, we also assessed the diagnostic performance of the four machine learning algorithms in the testing cohort. The root means square error (RMSE) was then used to order the top 10 significant feature variables in each model (Figure 6C). The SVM machine learning model, which had an AUC of 0.836 compared to RF’s 0.773, XGB’s 0.760, and GLM’s 0.552, had the highest AUC (Figure 6D). Combining these findings, the SVM model was found to be the most effective at distinguishing patients with various clusters among the four machine-learning models. We chose the top five predicted genes (MMP9, DDO, SLC45A3, FRK, and TMEM19) from the SVM model for additional investigation. DDO, FRK, and TMEM19 were three of the five genes with AUC values larger than 0.7 (Figure 6E), and these three genes were considerably elevated in NAFLD in comparison to controls (Figure 6F). We also compared their expression differences between disulfidptosis-associated clusters and the result showed that the expression of all three genes was greater in Cluster 1 than in Cluster 2, but the expression of FRK and TMEM19 was significantly different between Clusters 1 and 2, while the expression of DDO was not significantly different between disulfidptosis-associated clusters (Supplementary Figure S5).

To better predict the risk of patient morbidity, we constructed a nomogram model based on the three genes in the SVM model (Figure 7A). The results of calibration curves showed that the predictive ability of the nomogram model was accurate (Figure 7B). In addition, patients with NAFLD could benefit from column line graphs as shown in the DCA (Figure 7C). After that, we tested our 3-gene prediction model on two independent liver tissue datasets to validate its accuracy. The ROC curves demonstrated that the performance of the 3-gene prediction model was satisfactory, with an AUC value of 1.000 in the GSE63067 dataset (Figure 7D), and 0.909 in the GSE164760 dataset (Figure 7E). This indicates that our diagnostic model is equally effective in differentiating NAFLD patients from normal individuals.

To further investigate the potential role of the three biomarkers in NAFLD, we performed GSEA on each biomarker in the training set. The results of GSEA showed that the “cytokine receptor interaction” pathway and the “JAK-STAT signaling pathway” were enriched in the groups with low expression of DDO, FRK, and TMEM19 (Figures 8A–C).

Next, we analyzed whether or not there was a connection between the expression of three diagnostic genes and infiltrated immune cells in the training set. The results showed that the DDO gene was positively correlated with resting mast cells, δγT cells, M2 macrophages, and CD8T cells; and negatively correlated with activated mast cells, monocytes, activated NK cells, neutrophils, naive B cells, naive CD4 T cells and activated dendritic cells (Figure 8D). And FRK gene was positively correlated with resting mast cells, M1 macrophages, resting dendritic cells, and δγT cells; and negatively correlated with activated mast cells, neutrophils, monocytes, and naive B cells (Figure 8E). TMEM19 gene was positively correlated with naive CD4 T cells, resting dendritic cells, M1 macrophages, andδγT cells; and it was negatively correlated with activated mast cells, neutrophils, monocytes, naive B cells, activated dendritic cells, and Plasma cells (Figure 8F). In general, the expression of these genes may be related to the amount of infiltration of various immune cells, which suggests that these critical diagnostic genes may be engaged in immune control in the pathogenesis of NAFLD.

We assessed the expression levels of these biomarkers in a 12-week HFD-fed mice model to further confirm the validity of the biomarkers discovered in NAFLD. The HFD group had severe hepatic steatosis with sporadic inflammation, as seen by H&E staining of liver tissue sections (Figures 9A, B). In comparison to the NC group, the levels of hepatic TG were considerably higher in the HFD group (Figure 9C). Finally, we looked at three model genes’ expression in the livers of HFD and NC groups. We discovered that the HFD group had significantly higher levels of DDO, FRK, and TMEM19 expression than the NC group (Figures 9D–F).