Identification of Glycolysis-Related Genes in MAFLD and Their Immune Infiltration Implications: A Multi-Omics Analysis with Experimental Validation

Figure 1. Flowchart of the bioinformatics analysis and experimental validation.

Figure 2. Identification and functional analysis of glycolysis-related differentially expressed genes (DEGs). (A) Volcano plot displaying DEGs in liver tissues from MASH patients compared to healthy controls. (B) Heatmap showing the expression profiles of glycolysis-related DEGs in healthy controls and MAFLD patients. (C) Venn diagram illustrating the intersection of DEGs and glycolysis-related genes. (D) GSEA of the glycolysis pathway. (E) KEGG pathway enrichment analysis of glycolysis-related DEGs. (F) GO pathway enrichment analysis of glycolysis-related DEGs.

Figure 3. Identification of key modules associated with MAFLD and glycolysis-related key genes. (A) Determination of the soft-thresholding power in WGCNA. (B) Cluster dendrogram of the WGCNA analysis. (C) Module–trait relationship heatmap. (D) Venn diagram showing the overlap between key modules and glycolysis-related DEGs. (E,F) Scatter plots showing the correlation between GS with MM in the brown and yellow modules, respectively.

Figure 4. Analysis of the immune microenvironment and immune-related functions. (A) Stacked histogram showing changes in immune cell proportions. (B) Comparison of immune cell infiltration between the control and MAFLD groups. (C) Correlation between infiltrating immune cells and key DEGs. *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001; ns: not significant.

Figure 5. The glycolytic metabolic pathway in the characteristics of single-cell transcriptomics. (A) UMAP plot illustrating cell type distribution based on canonical marker genes. (B) Heatmap displaying the top two marker genes for each identified cell cluster. (C) Dot plot demonstrating significant differences in glycolysis-related metabolic pathway activity between control and MAFLD groups. (D) UMAP plot showing glycolysis/gluconeogenesis pathway activity across various cell types in MAFLD samples. (E) Box plot depicting glycolysis/gluconeogenesis activity levels among different cell types. (F) Circle plot indicating the number of intercellular interactions between immune cells and hepatocytes in high-glycolysis versus low-glycolysis groups. ****, p < 0.0001.

Figure 6. Spatial co-localization analysis of the key genes and monocyte-derived macrophage markers. (A) Identified spatial clusters in representative liver sections from MASH and control mice: PV (perivenous zone), PP (periportal zone), Mid (middle zone), and MoM (monocyte-derived macrophage-enriched cluster). (B) Violin plots showing the expression levels of markers and the key genes across different zones. Glul, Cyp2e1, and Oat serve as markers for the PP zone; Cyp2f2 for the PV zone; and Ccr2 and Lyz2 for the MoM region. (C,D) Co-localization of CDK1 and SOX9 with the MoM region marker gene Ccr2. Note: Panels (C,D) show microscope images obtained from GEO database lacking scale or magnification information.

Figure 7. Identification of optimal feature genes among key genes using machine learning to construct the TF–mRNA–miRNA regulatory mulberry plot. (A) LASSO regression analysis identified four genes with the lowest binomial deviance as candidate feature genes for MAFLD diagnosis. (B) SVM-RFE selected the top four genes based on the highest diagnostic accuracy. (C) Random Forest error-rate curves across 0–500 trees (black: overall error; red and green: class-specific errors for the MASH and Control groups, respectively), and gene importance ranked by MeanDecreaseGini. (D) Gene importance scores from the XGBoost algorithm. (E) Overlap of feature genes Identified by four ML algorithms (F) Construction of the TF–mRNA–miRNA regulatory mulberry plot (H. sapiens). (G) mRNA–miRNA interaction network (M. musculus). (H) The association between genes and microRNAs.

Figure 8. Diagnostic value of optimal feature genes and validation of key genes. (A) Nomogram construction using optimal feature genes in dataset GSE213621. (B) Volcano plot of the external validation dataset GSE126848, highlighting differential expression of key genes. (C,D) Expression levels of optimal feature genes in GSE213621 and the four key genes in GSE126848, respectively. ****, p < 0.0001. (E) Expression levels of SOX9 in GSE33814. *, p < 0.05; ****, p < 0.0001. (F,G) ROC curve of key genes in datasets GSE213621 and GSE126848.

Figure 9. Experimental validation of key gene expression in a mouse MASH model. (A) The successful establishment of the MASH mouse model confirmed by HE and Oil Red O staining (n = 4). WT + MCS: Wild-type mice on MCS diet; WT + MCD: Wild-type mice on MCD diet. (B) Volcano plot showing differential gene expression between the MASH and control groups, including key genes. (C) GSEA reveals enrichment of the glycolysis pathway in MASH model mice. (D) Bar chart showing expression differences of CDK1, SOX9, and HKDC1 between MASH and control groups. **, p < 0.01. (E) Heatmap showing differential expression of key genes between MASH and control groups. (F) Protein expression levels of SOX9 in MASH and control groups.