Your new experience awaits. Try the new design now and help us make it even better

ORIGINAL RESEARCH article

Front. Immunol., 05 November 2025

Sec. Inflammation

Volume 16 - 2025 | https://doi.org/10.3389/fimmu.2025.1684937

This article is part of the Research TopicAdvancing Immune Research in Chronic Liver Diseases Through New Approach MethodologiesView all 3 articles

Identification of glycolysis-related clusters and immune cell infiltration in hepatic fibrosis progression using machine learning models and experimental validation

Guanglin XiaoGuanglin XiaoZhiling DengZhiling DengKe QiuKe QiuAoyi LiAoyi LiXingyue YiXingyue YiHong Ren*Hong Ren*
  • Department of Infectious Diseases, Key Laboratory of Molecular Biology for Infectious Diseases, Ministry of Education, Institute for Viral Hepatitis, The Second Affiliated Hospital, Chongqing Medical University, Chongqing, China

Objectives: Although glycolytic reprogramming constitutes a fundamental driver of hepatic fibrosis (HF), its precise mechanistic contributions remain incompletely characterized. This investigation systematically identified molecular signatures of glycolysis-related genes (GRGs) in HF. We further developed a glycolytic activity-based model for HF risk stratification.

Methods: Integrated analysis of GEO datasets (GSE276114, GSE84044, GSE49541) identified differentially expressed genes (DEGs) associated with HF progression. Integrated weighted gene co-expression network analysis (WGCNA) with six machine learning algorithms to identify core GRGs genes associated with HF progression, and systematically characterized their biological functions and immunoregulatory roles through immune infiltration assessment, functional enrichment, consensus clustering, and single-cell differential state analysis. Glycolytic activity was evaluated in CCl4-induced fibrotic mice and TGF-β-stimulated LX-2 cells. Additionally, the expression of core GRGs was validated using immunohistochemical staining and RT-qPCR.

Results: Through the intersection of WGCNA, DEGs, and GRGs, machine learning identified six core GRGs: B3GNT3, CHST4, DCN, GPC3, SOX9, and VCAN. Based on the core GRGs, three GRG-based molecular subtypes were defined. Cluster C, with higher expression of the core GRGs, exhibited significantly enhanced immune infiltration, particularly of adaptive immune cells compared to Cluster A and B. Cluster C comprised a mixed landscape of T cells, mast cells, and pro-fibrogenic cells, distinct from the innate immune-dominant profiles of Clusters A and B. Both in vivo and in vitro analyses demonstrated enhanced glycolysis in fibrotic progression, accompanied by consistent upregulation of core GRGs.

Conclusions: Glycolytic reprogramming is a key pathogenic driver in HF progression and associated immune infiltration. Investigating this metabolic-immune dysregulation represents a promising therapeutic focus for progression of HF.

1 Introduction

Hepatic fibrosis (HF) initially functions as a self-limiting repair mechanism, where transient extracellular matrix deposition preserves tissue integrity during acute injury. This tightly regulated process becomes detrimental when chronic insults (e.g., MASLD) sustain fibrogenic activation, leading to disproportionate scar accumulation that disrupts hepatic architecture and function (1). This maladaptive remodeling drives cirrhosis, which accounts for 2.4% of global mortality and often necessitates liver transplantation due to vascular distortion and regenerative failure (2). Despite decades of research, current antifibrotic strategies remain limited to symptom management rather than targeting core regulators of progressive matrix deposition.

The Warburg effect, originally described in cancer, denotes a metabolic shift from oxidative phosphorylation to aerobic glycolysis. This adaptation allows cancer cells to rapidly generate biomass and energy to support proliferation. Subsequently, this metabolic reprogramming has been observed in fibrotic diseases across multiple organs, including the liver, kidneys, skin, and lungs (36). Hepatic stellate cells (HSCs) located in the space of Disse are the primary pro-fibrogenic cells in the liver. Upon activation, the quiescent HSCs transdifferentiate into proliferative, migratory, and contractile myofibroblasts (7). During HSC activation, key glycolytic enzymes, including HK2, PFK1, PKM2, GLUT1, and MCT4, are upregulated, while gluconeogenic enzymes such as PCK1 and FBP1 are downregulated (8, 9). This metabolic shift from oxidative phosphorylation to aerobic glycolysis supports the high energy demands of proliferating HSCs by enabling rapid ATP generation (10).

Glycolysis underpins the functional capacity of innate immune cells (such as inflammatory macrophages, DCs, and neutrophils) and is equally critical for innate-like and adaptive immune cells, particularly in the acquisition of an IL-17+ phenotype (11). Macrophage heterogeneity plays a central role in liver fibrogenesis (12). Resident Kupffer cells, upon activation, initiate and propagate inflammatory responses during liver injury. The transition of macrophages toward an inflammatory phenotype is regulated by metabolic reprogramming of glucose metabolism (13). Inflammatory macrophages engage glycolysis, the pentose phosphate pathway, and a TCA cycle configured for citrate production, alongside nitric oxide generation. Conversely, anti-inflammatory macrophages are dependent on fatty acid oxidation, oxidative phosphorylation, and the arginase pathway (14). In liver fibrosis, this metabolic shift is exploited through PKM2 activation, which augments glycolytic flux and drives macrophage reprogramming toward a profibrogenic phenotype, as demonstrated in mouse models (15).

Glycolysis in HF is under intensive study but far from fully elucidated. HF involves joint action from HSCs, immunocytes, and hepatocytes. However, glycolysis of HF and its impact on the infiltration and activation of immune cells has been scarcely studied. To address this, we combined transcriptomic data from multiple patient cohorts with comprehensive glycolysis-related genes (GRGs) sets to identify core GRGs consistently associated with HF progression. Utilizing machine learning algorithms, we developed a robust diagnostic model and validated its staging performance. Furthermore, we performed consensus clustering to define novel glycolysis-based molecular subtypes and characterized their distinct immunological profiles. Finally, we experimentally validated the relevance of core GRGs and glycolytic phenotypes in vitro and in animal models. Our work thus provides a holistic view of glycolytic reprogramming in HF, establishing a foundational resource for future mechanistic and therapeutic exploration.

2 Materials and methods

2.1 Data sources

Transcriptomic profiles and clinical staging data from patients with HF were obtained from three GEO datasets: GSE276114, GSE84044, and GSE49541. GSE276114 classified 177 specimens into histological fibrosis stages: 39 early fibrosis (F0-F2), 24 severe fibrosis (F3), and 114 cirrhosis (F4) cases (16). GSE84044 provided Scheuer Score-evaluated data encompassing 43 non-fibrotic livers, 20 stage 1, 33 stage 2, 18 stage 3, and 10 stage 4 fibrosis cases (17). GSE49541 contained 72 samples grouped as 40 stage 0–1 versus 32 stage 3–4 fibrosis specimens (18). Based on clinical consensus, stages F0-F2 and F3-F4 were designated mild and advanced fibrosis respectively. Additionally, 315 GRGs were retrieved from the MSigDB. The GSE136103 dataset included single-cell transcriptome data from 5 healthy patient liver tissue samples and 5 cirrhotic patient liver tissue samples (19).

To validate the generalizability of core GRGs in HF, we retrieved two independent GEO datasets: GSE130970 and GSE14323. The GSE130970 cohort comprised 23 healthy controls and 53 fibrosis patients stratified by disease severity (stage 1: n=28, stage 2: n=9, stage 3: n=14, stage 4: n=2) (20). GSE14323 provided transcriptomic profiles from liver specimens of individuals with cirrhosis attributed to hepatitis C virus (HCV) infection or alcohol consumption. Liver tissue samples were classified for the analysis as normal livers (n=19) and HCV cirrhosis (n=41). The characteristics of the samples are described in the Supplementary Tables 1, 2 (21). Additionally, gene expression data and corresponding survival information for hepatocellular carcinoma (HCC) cases were acquired from the TCGA.

2.2 Data preprocessing and differential gene expression analysis

The efficacy of batch effect correction was validated through two-dimensional principal component analysis (PCA) clustering. Following quality control, we generated a consolidated normalized gene expression dataset using the “limma” package (22). To improve the reliability of differentially expressed genes (DEGs), probe sets for which the adjusted P <0.05, and |log2FC| > 0.5 between mild and advanced fibrosis were defined as significant DEGs.

2.3 Functional analyses of DEGs and gene set enrichment analysis

Using the “ClusterProfiler” package (23), we conducted comprehensive GO and KEGG enrichment analyses. Statistical significance thresholds were established at a false discovery rate (FDR) < 0.05, with P-values adjusted via the Benjamini-Hochberg procedure for multiple test correction.

2.4 Weighted gene co-expression network analysis

WGCNA was implemented to delineate associations between gene modules and pathological traits, enabling identification of hub genes implicated in HF progression. We constructed a similarity matrix using Pearson correlation coefficients, subsequently transformed into adjacency and topological overlap matrices (TOM) through optimized soft-thresholding. Gene clustering was performed via dynamic tree-cutting algorithms with a minimum module size threshold of 50 genes. Modules were designated arbitrary color labels, while module eigengenes (MEs) represented their transcriptional signatures. We identified GRGs critically involved in HF pathogenesis by intersecting these candidate genes with DEGs and pre-defined GRGs.

2.5 Development and validation of a diagnostic signature for HF progression using machine learning algorithms

Six distinct machine learning algorithms—Random Forest (RF), Support Vector Machine (SVM), Generalized Linear Model (GLM), Gradient Boosting Machine (GBM), Neural Network (NNET), and Least Absolute Shrinkage and Selection Operator (LASSO)—were employed to assess gene significance. Using the “caret” package, models underwent training with 70% of data (stratified sampling via createDataPartition) and validation on 30% test sets. Model robustness was enhanced through 5-fold repeated cross-validation. Hyperparameter tuning was conducted using the default tuning grids of the caret package, while feature importance ranking was performed through a model-agnostic approach implemented in the “DALEX” package. Core GRGs were derived from intersecting the top-ranked features across all algorithms. Building on this, we established a logistic regression model using the “glm” package. This model subsequently served as the foundation for constructing a clinical nomogram via the “rms” package to stratify HF progression. Performance of the prognostic model was rigorously evaluated through calibration curves, decision curve analysis (DCA), and receiver operating characteristic (ROC) curves.

2.6 Gene set enrichment analysis

To elucidate signaling pathway disparities between mild and advanced fibrosis stages, GSEA was implemented (24). Annotated gene sets for disease-relevant pathways were curated from the MsigDB to establish the background reference. Gene sets that were significantly enriched were pinpointed using consistency scores (adjusted p value <0.05).

2.7 Gene set variation analysis

Employing an unsupervised, non-parametric methodology, GSVA quantifies pathway-level perturbations by transforming gene-centric expression data into functional enrichment scores (25). This approach characterizes biological activity through comprehensive scoring of predefined gene sets. For the current study, we retrieved pathway-specific gene collections from the MsigDB and computed enrichment profiles using the GSVA algorithm.

2.8 Consensus clustering analysis and differences between subtypes

Leveraging expression profiles of hepatic fibrosis-associated GRGs, we performed molecular subtyping of HF patients through unsupervised consensus clustering. This analysis was implemented in R using the “ConsensusClusterPlus” package (26) with the following parameters: 1,000 bootstrap iterations to ensure algorithmic stability, item resampling proportion of 0.8, feature resampling proportion of 1.0, and the k-means clustering algorithm as the foundational method.

2.9 Immune cell infiltration analysis

Immune cell infiltration was quantified using both CIBERSORT and ssGSEA. Spearman correlation analysis assessed relationships between immune cell proportions and core GRGs expression, while boxplots visualized infiltration differences of 28 immune cells across GRG-based subtypes.

2.10 Analysis of scRNA-seq data

Single-cell RNA sequencing data were processed and analyzed using the “Seurat” package, which facilitated normalization, dimensionality reduction, clustering, and visualization (27). Differentially expressed genes across clusters were identified using the FindAllMarkers function to support cell type classification. Cell identities were determined through a three-step annotation workflow: automated cell type prediction was performed with SingleR, followed by manual verification using well-established marker genes from literature, and finally validated through the CellMarker database (28, 29). Single-cell Phenotype-Associated Subpopulation identifier (scPAS) was applied to quantify the association between individual cells in scRNA-seq data and GRG-based molecular subtypes. Differential expression analysis between high- and low-glycolysis cell subpopulations was performed using the Wilcoxon rank-sum test via the FindAllMarkers function. Genes with log2FC > 1 and FDR < 0.05 in high-glycolysis cells were defined as a positive signature, while those with log2FC < –1 and FDR < 0.05 in low-glycolysis cells comprised the negative signature. GSVA was then used to compute signature scores for each sample, with a composite score derived by subtracting the negative signature score from the positive signature score (30).

2.11 CCl4-induced HF model and histological analysis

All experimental procedures involving animals were conducted in compliance with ethical standards approved by the Institutional Animal Care and Use Committee of Chongqing Medical University (Approval No. IACUC-SAHCQMU-2025-0116). To establish fibrosis models, mice received intraperitoneal injections of carbon tetrachloride (CCl4) dissolved in corn oil (1:4 v/v) at 1 ml/kg body weight. Administration occurred three times weekly with differential durations: 8 weeks for mild fibrosis induction and 12 weeks for advanced fibrosis induction. Liver sections underwent histological evaluation using standardized Hematoxylin and Eosin (H&E) and Masson’s trichrome staining protocols as described (31).

2.12 Immunohistochemistry

Deparaffinized and rehydrated sections underwent antigen retrieval in heated citrate buffer. After peroxidase blocking, tissues were incubated with primary antibody at 4 °C overnight, followed by a secondary antibody and DAB development using a standard IHC protocol.

2.13 Cell lines and culture

The human HSC line LX-2 was acquired from Ubigene Biosciences (China). Cells were maintained at 37 °C under 5% CO2 atmosphere with saturated humidity. Culture medium consisted of high-glucose DMEM supplemented with 10% FBS (ExCell Bio, China) and 1% penicillin/streptomycin.

2.14 EdU staining assay

LX-2 cells were plated in 6-well plates and maintained overnight. The cells were treated with TGF-β (10 ng/ml, Sinobiological) and/or the glycolytic inhibitor 2-DG (2mM, Selleck, USA). EdU incorporation assay was performed with an EdU assay kit (Beyotime, Shanghai, China), according to the manufacturer’s instructions (32). Finally, stain the nuclei with DAPI and observe them using a fluorescence microscope (Nikon, Tokyo, Japan).

2.15 Cell migration

Cell migratory capacity was quantified using 8-μm pore Transwell chambers (LABSELECT, China). LX-2 cells were plated in serum-free medium within the upper compartment, while the lower chamber contained 10% FBS. Following 24-hour incubation, traversed cells were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet (Beyotime, China). Invasive cells were manually enumerated across five randomly selected microscopic fields using an inverted fluorescence microscope (Nikon, Japan), with subsequent image acquisition and quantitative analysis.

2.16 Lactate and glucose measurement

Lactate concentrations in LX-2 cell supernatant, mice serum, and mice liver tissue were quantified using lactate assay kits (A019-2-1, Nanjing Jiancheng Bioengineering Institute, China), strictly adhering to the manufacturer’s protocol. Glucose concentration in LX-2 cell supernatant was measured using a glucose assay kit (#S0201S, Beyotime, Shanghai, China).

2.17 2-NBDG glucose uptake

LX-2 cells were plated in 24-well plates and allowed to adhere overnight. Experimental groups were stimulated with TGF-β (10 ng/ml) for 24 hours. Following treatment, cellular glucose uptake was measured by incubating cells with 0.5 ml 2-NBDG working solution at 37 °C for 60 min per manufacturer’s instructions (Beyotime, China). After removal of supernatant and two PBS washes, nuclei counterstaining was performed. Fluorescence microscopy (Nikon, Japan) was employed to quantify glucose uptake.

2.18 Real-time reverse transcription PCR

Total RNA isolation employed Trizol reagent per manufacturer’s protocol (Invitrogen, USA). Subsequent cDNA synthesis utilized a Takara reverse transcription kit (Japan). Primer sequences are detailed in Supplementary Tables 3, 4.

2.19 Statistical analysis

The Wilcoxon test was utilized to complete the comparative analysis among the groups. Correlation analyses utilized Spearman’s rank-order method. All analyses were performed in R 4.4.1 and GraphPad Prism 8.0, with P < 0.05 defining statistical significance.

3 Results

3.1 DEGs between mild and advanced HF

To establish a robust integrated cohort, the three HF datasets (GSE276114, GSE84044, and GSE49541) were systematically combined. PCA revealed distinct clustering patterns among fibrotic cases both before and after normalization (Figure 1A). Furthermore, distribution characteristics were assessed using box plots, which demonstrated improved comparability after batch adjustment (Figure 1B). Comparative profiling revealed distinct transcriptional signatures between mild and advanced HF stages, yielding 517 significantly upregulated and 208 downregulated DEGs. These differential expression patterns were visualized through volcano plot and hierarchical clustering heatmap analyses (Figures 1C, D). Subsequent pathway enrichment analysis of hepatic fibrosis-associated DEGs identified critical biological mechanisms. GO analysis demonstrated significant enrichment (P<0.05) in fibrogenic processes, particularly elastic fiber assembly, collagen-activated signaling pathway, and complex of collagen trimers (Figure 1E). Complementary KEGG pathway analysis revealed prominent involvement in ECM-receptor interactions, galactose metabolism, and fructose and mannose metabolism (Figure 1F). Integrated GO and KEGG analyses indicate that hepatic fibrogenesis is associated with metabolic reprogramming in glycometabolism.

Figure 1
Panel A shows PCA plots with three distinct groups in different colors. Panel B depicts boxplots of two datasets. Panel C presents a volcano plot highlighting significantly upregulated and downregulated genes. Panel D contains a heatmap of gene expression with a color scale indicating control and treated groups. Panel E illustrates a GO enrichment analysis with a bubble plot, highlighting significant biological processes. Panel F shows a KEGG pathway enrichment analysis with another bubble plot, detailing pathways involved.

Figure 1. Merging the database and screening differential expressed genes (DEGs). (A) Principal Component Analysis (PCA) of samples from the three databases before and after data merging. (B) Sample distribution before and after homogenization of the datasets. (C, D) DEGs were screened using the limma package under the standard of adj.P.Val<0.05, and |log2FC| > 0.5. Upregulated genes were plotted in red and downregulated genes in green in a volcano plot (C). The heatmap (D) displays a group of genes differentially expressed in mild and advanced fibrosis. (E) GO annotation of the DEGs in association with annotated biological process (BP), cellular component (CC), and molecular function (MF). (F) Demonstration of KEGG enrichment results.

3.2 Identification of key genes related to HF progression using WGCNA

WGCNA identified HF-associated hub genes through integrated analysis of transcriptomic profiles and clinical metadata. Sample clustering based on Pearson’s correlation coefficients revealed distinct patient subgroups (Figure 2A). Optimization of the scale-free topology model determined a soft threshold power of 3, which maximized network interconnectivity while preserving gene similarity (Figure 2B). This approach delineated 11 co-expression modules, among which the blue module (976 genes) exhibited the strongest correlation with HF progression (r = 0.42, P = 3 × 10-17) (Figures 2C, D).

Figure 2
Panel A shows a sample dendrogram and trait heatmap, categorizing samples into control and treated groups with red blocks indicating treatment. Panel B features two graphs: scale independence and mean connectivity as power thresholds. Panel C displays a gene dendrogram with colored modules representing different gene groups. Panel D is a heatmap of module-trait relationships, showing correlations between gene modules labeled by color and traits like mild and advanced fibrosis, with correlation values and significance levels.

Figure 2. Construction and module analysis of WGCNA. (A) Sample dendrogram and trait heatmap. (B) Network topology analysis under different soft threshold powers. (C) Clustering Dendrogram, illustrating the hierarchical grouping of genes by topological overlap, with the specified module colors representing different gene clusters. (D) Correlation analysis for the relationship between different coexpression modules and clinical features.

3.3 Glycolysis underscores the core predictive DEGs in HF progression

Fifteen upregulated GRGs were identified through intersecting three gene sets: 517 upregulated DEGs, 976 key genes from the WGCNA module most relevant to HF progression, and 315 known GRGs (Figure 3A). These genes included B3GNT3, CHST4, DCN, GPC3, GPC4, HIF1A, HK1, HKDC1, LDHB, SLC2A1, SOX9, TFF3, TGFBI, TPBG, and VCAN. Heatmap analysis demonstrated pronounced upregulation of these 15 GRGs in advanced fibrosis compared to mild fibrosis (Figure 3B). Protein-protein interaction (PPI) network analysis of the 15 GRGs revealed robust connectivity among all genes except TPBG, TFF3, CHST4, and B3GNT3, with the transcription factor SOX9 emerging as a central hub genes (Figure 3C). Distinct expression patterns between advanced and mild HF further supported their pathological relevance (Figure 3D). Functional enrichment analysis via GO highlighted their association with glycolytic processes (e.g., “lactate metabolic process”) and fibrotic pathways, including “negative regulation of mesenchymal cell apoptotic process” (Figure 3E). KEGG pathway analysis underscored their involvement in critical pathways such as “Central carbon metabolism in cancer” and “Glycolysis/Gluconeogenesis” (Figure 3F), suggesting a mechanistic link between glycolysis dysregulation and fibrotic progression.

Figure 3
A series of scientific data visuals includes: A) A Venn diagram showing overlaps among WGCNA, Merage_up, and Glycolysis categories with intersection numbers. B) A heatmap of gene expressions for mild versus advanced fibrosis, highlighted by a blue to red gradient. C) A network diagram connecting proteins linked to fibrosis. D) Box plots comparing gene expression levels between fibrosis stages. E) A GO enrichment analysis bubble chart showing biological processes with varying significance levels. F) A KEGG pathway analysis chart illustrating metabolic pathways enriched in fibrosis.

Figure 3. Glycolysis of DEGs in mild and advanced fibrotic livers. (A) Venn plot displaying the intersection of upregulated DEGs, WGCNA and glycolysis-related genes (GRGs). (B) Heatmap displaying the expression patterns of the 15 GRGs. (C) PPI network of 15 overlapping genes. (D) Boxplot displaying the expression patterns of the 15 GRGs. (E, F) GO annotation and KEGG enrichment analyses of GRGs. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

3.4 Construction and validation of diagnostic model by machine learning

A diagnostic model was constructed using six machine learning algorithms: GBM, RF, GLM, SVM, LASSO, and NNET. Evaluation of the cumulative residual distribution plots and residual boxplots revealed that GBM, RF, GLM, SVM, and LASSO exhibited minimal residuals (Supplementary Figure S1), while all six models demonstrated robust discriminative performance for clinical outcomes (AUC > 0.8; Figure 4A). The top ten variables for each model, ranked by root-mean-square error (RMSE), were visualized in Figure 4B. To identify core diagnostic biomarkers, we extracted the intersection of the top ten most influential genes across all models, yielding six core GRGs: B3GNT3, CHST4, DCN, GPC3, SOX9, and VCAN. ROC analysis confirmed strong diagnostic utility for all genes except B3GNT3 (Figure 4C). A nomogram integrating these six core GRGs was subsequently developed using the “rms” package to stratify HF risk (Figure 4D). The calibration curve (Brier: 0.139) and DCA collectively demonstrated the model’s robust predictive performance (Figures 4E, F).

Figure 4
Panel A presents ROC curves comparing multiple models, with RF achieving the highest AUC of 0.873. Panel B displays feature importance across models like GLM, LASSO, and NNET. Panel C shows ROC curves with VCAN having the highest AUC of 0.792. Panel D depicts a scoring system linking gene expression to diagnostic possibility. Panel E plots observed frequency against predicted risk with a Brier score of 0.139. Panel F illustrates a decision curve for net benefit against high-risk thresholds for different models.

Figure 4. Construction and validation of diagnostic model by machine learning. (A) The ROC curves for GBM, RF, GLM, SVM, LASSO and NNET models and their corresponding AUC values. (B) The top ten variables in the RMSE ranking used for evaluating the feature importance of the models (GBM, RF, GLM, SVM, LASSO and NNET), and the significance contribution of each model to the input features was analyzed. (C) ROC curve and AUC values of core GRGs (B3GNT3, CHST4, DCN, GPC3, SOX9, and VCAN). (D) A nomogram of six core GRGs. (E) Calibration curve analysis of the prognostic model. (F) DCA.

External validation in the GSE14323 and GSE130970 datasets confirmed the robust discriminative power of these core GRGs. Furthermore, subgroup analyses stratified by etiology (steatohepatitis and viral hepatitis) consistently demonstrated their diagnostic efficacy across distinct fibrotic subtypes (Supplementary Figure S2).

3.5 Immune landscape analysis across HF progression

Immune infiltration profiling across HF patient samples characterized the relative abundance of 22 immune cell subtypes, visualized via stacked bar chart (Figure 5A). In addition, there was a close interplay between immune cells, with a significant negative correlation between M1 macrophages and M2 macrophages, and a strong positive correlation between activated dendritic cells and M2 macrophages (Figure 5B). In patients with advanced HF, the proportion of memory B cells, naive B cells, activated mast cells, plasma cells, CD4 memory activated T cells, M0 macrophages, M1 macrophages was found to be higher than those in mild HF. However, the percentage of activated dendritic cells, M2 macrophages, monocytes, resting NK cells, and neutrophils was lower in advanced HF samples (Figure 5C). All of the core GRGs had positive correlations with M0 macrophages, activated mast cells, plasma cells and activated T cells CD4 memory and negative correlation with M2 macrophages, resting mast cells, monocytes, resting NK cells and naive T cells CD4 (Figure 5D).

Figure 5
Panel A shows a stacked bar chart representing estimated cell type proportions in advanced and mild fibrosis. Panel B displays a correlation matrix heatmap of different cell types. Panel C is a box plot comparing cell type proportions between advanced and mild fibrosis groups, indicating significance. Panel D is another correlation matrix showing gene expression correlations with cell types, colored by correlation coefficient.

Figure 5. Immune infiltration analysis. (A) The bar chart displays immune cell infiltration results of 22 immune cells in each liver fibrosis patients. (B) The correlation matrix of immune cell proportions. (C) The group comparison chart illustrates differences in the abundance of immune cell infiltration between Cluster C and Cluster A and B. (D) Correlation analysis of the core GRGs with immune cells. *P < 0.05, **P < 0.01, ***P < 0.001. ns: no statistically significant (P > 0.05).

3.6 GSEA and GSVA analysis of core GRGs

GSEA analysis of six core GRGs revealed distinct pathway enrichments: B3GNT3 was predominantly enriched in integrin αβ-talin-vinculin signaling and translation initiation; CHST4 showed significant enrichment in translation initiation and aberrant TDP-43-mediated electron transfer in complex I; DCN primarily enriched Salmonella SopB-ARNO-ARF-G-actin signaling pathway and ARNO-ARF-G-actin signaling pathway; GPC3 and SOX9 both activated Shigella IpgD-ARNO-ARF-G-actin signaling pathway and ARNO-ARF-G-actin signaling pathway; VCAN was enriched in PI3K signaling pathway and RAS-PI3K signaling pathway (Supplementary Figure S3).

GSVA demonstrated that overexpression of all examined genes (B3GNT3, CHST4, DCN, GPC3, SOX9, and VCAN) consistently activated the epithelial-mesenchymal transition (EMT) pathway. Furthermore, angiogenesis was co-activated by CHST4, DCN, GPC3, SOX9, and VCAN, while hypoxia signaling was shared by B3GNT3 and CHST4. Notably, distinct pathways exhibited gene-specific activation: apical junction by B3GNT3, apoptosis by DCN, and TNF-α/NF-κB signaling by the GPC3, SOX9, and VCAN (Supplementary Figure S4).

3.7 Analysis of core GRGs based on the GSE136103 single-cell dataset

We further explored the expression of the core GRGs in HF patients through scRNA-seq data. After rigorous quality control screening, 51,972 high-quality cells were included in the subsequent analysis. The FindAllMarkers function and the Wilcoxon rank sum test were used to identify specific gene signatures for each cluster. 17 major cell clusters were identified, including exhausted CD8+ T cells (CD3D+CD8A+PDCD1+), Naive/central memory CD8+ T cells (CD3D+CD8A+LEF1+), NK cells (GNLY+NKG7+), B cells (MZB1+CD79A+), plasma cells (MZB1+CD38+), neutrophils (S100A8+VCAN+), conventional dendritic cells (CLEC9A+XCR1+), mast cells (TPSAB1+CPA3+), Liver sinusoidal endothelial cells (LSEC) (CLEC4G+OIT3+), Kupffer cells (ADGRE1+CD68+), endothelial cells(CLEC14A+PLVAP+), HSCs (ACTA2+COL1A1+), innate lymphoid cells (KLRF1+KLRC1+), myofibroblasts (COL1A1+COL3A1+), and vascular endothelial cells (VWF+ACKR1+), epithelial cells (EPCAM+ALB+), proliferating cells (MKI67+), plasmacytoid dendritic cells (LILRA4+LRRC26+) (Figures 6A, B). Analysis of the UMAP clustering and expression heatmap showed prominent expression of B3GNT3, CHST4, and SOX9 in epithelial cells; DCN and GPC3 in HSCs and myofibroblasts; and VCAN in neutrophils (Figure 6C).

Figure 6
Panel A shows a UMAP plot with data points representing Healthy and Cirrhotic groups in different colors. Panel B displays a UMAP plot with clusters colored by cell type, labeled such as epithelial and liver_b_endothelial_cells. Panel C consists of six UMAP plots, each labeled with gene expressions (B3GNT3, CHST4, DCN, GPC3, SOX9, VCAN) with a color gradient legend indicating expression levels.

Figure 6. Clustering of GSE136103 scRNA-seq data and identification of cell types. (A) UMAP visualization of different group. (B) UMAP visualization of different cell types. (C) UMAP visualization of core GRGs.

3.8 Phenotyping power of the core GRGs in HF

Leveraging core GRG expression profiles, we stratified 373 multi-database samples into three distinct molecular phenotypes (Figure 7A). PCA validated clear separation among these subtypes (Figure 7B). Clinically, advanced fibrosis cases predominantly accumulated in Clusters C versus Cluster A and B (Figure 7C). A gradient expression pattern of core GRGs was observed across the clusters, which correlated with low (Cluster A), intermediate (Cluster B), and high (Cluster C) glycolytic levels (Figure 7D). To functionally characterize the differences among GRG-based molecular subtypes identified in our analysis, we performed GSEA comparing the transcriptional profiles of Cluster A and B vs Cluster C. The most notably enriched gene sets in Cluster C were Epithelial-Mesenchymal Transition (EMT), Hypoxia, and TNFA Signaling via NFKB (p < 0.001) (Figure 7E). ssGSEA was performed to clarify the characteristics of the three GRG-based molecular subtypes in the immune microenvironment. The proportions of activated B cells, activated CD4 T cells, activated CD8 T cells, central memory CD4 T cells, effector memory CD8 T cells, eosinophils, mast cells, MDSCs, memory B cells, NK cells, NKT cells, plasmacytoid dendritic cells, Th1 cells, and Th2 cells were positively correlated with glycolytic levels. In contrast, the proportions of macrophages and monocytes showed a negative correlation (Figure 7F).

Figure 7
A multi-panel scientific figure includes: A) a consensus matrix heatmap showing clustering; B) a principal component analysis (PCA) scatter plot differentiating three clusters; C) a heatmap illustrating gene expression patterns across clusters; D) box plots depicting expression levels of specific genes across clusters; E) an enrichment plot for gene sets highlighting significant pathways; F) box plots comparing immune infiltration levels among differing types; G) UMAP plots illustrating risk score and cluster distribution; H) box plots showing cell type proportions across clusters A and B.

Figure 7. Phenotype clustering by the expression of core GRGs. (A) Consensus clustering on liver fibrosis samples based on the six core GRGs. (B) PCA of the sample distribution across different phenotypes. (C) Heatmap showing the association between gene expression and different phenotypes plotted. (D) Expression distinction of core GRGs across different phenotypes. (E) GSEA analysis based on the canonical pathways gene sets. (F) Box plot showing immune infiltration differences among GRG-based molecular subtypes by ssGSEA. (G) The UMAP visualization of the GRG subtype-calculated risk scores and the UMAP visualization of the GRG subtype-selected cells. (H) GRG subtype-selected cells with the corresponding bar plots shows the detailed constitutions in each cell type. *P < 0.05, **P < 0.01, ***P < 0.001.

Utilizing 373 bulk tissue samples as reference, scPAS analysis identified 3,511 cells associated with Cluster C and 6,791 cells associated with Clusters A and B (Figure 7G). Notably, Cluster C was predominantly enriched with immune cell populations, including exhausted CD8+ T cells, naive T cells, mast cells, and plasma cells. Additionally, this cluster contained substantial contributions from epithelial cells, myofibroblasts and HSCs. In contrast, Clusters A and B were primarily composed of NK cells, Kupffer cells, and dendritic cells (Figure 7H).

3.9 Analysis of the core GRGs in HCC

HF poses significant diagnostic challenges and frequently progresses to malignancy, with many patients presenting at advanced HCC stages. Analysis of TCGA-LIHC data revealed elevated expression of SOX9, DCN, GPC3, and B3GNT3 in tumor tissues versus non-tumor counterparts (Figure 8A). Notably, survival analysis identified CHST4 (log-rank P = 0.040), SOX9 (log-rank P = 0.018), and VCAN (log-rank P = 0.037) as prognostic biomarkers significantly correlated with overall survival in HCC patients (Figure 8B), implicating these genes as potential oncogenic drivers during fibrotic malignant transformation.

Figure 8
Figure A contains six box plots comparing gene expression levels between tumor and normal tissue for CHST4, SOX9, DCN, GPC3, VCAN, and B3GNT3. Figure B presents survival analysis with Kaplan-Meier curves for disease-free and overall survival across different gene expression levels: CHST4, SOX9, DCN, GPC3, VCAN, and B3GNT3. Each plot shows blue and red lines for low and high expression groups, with log-rank p-values and hazard ratios detailed.

Figure 8. Expression and prognostic values of core GRGs in patients with HCC. (A) Expressions of core GRGs in the TCGA-LIHC patients. (B) Overall survival of HCC patients with high or low expression of core GRGs.

3.10 Glycolysis and core GRGs expression CCl4-induced HF mice model

Mice models of HF with varying severity were induced by intraperitoneal injection of CCl4 for 8 or 12 weeks (Figure 9A). The results of HE, Masson and IHC staining showed that compared with the oil controls, the CCl4 groups had significantly increased collagen deposition (Figure 9B). As HF progresses, plasma levels of alanine aminotransferase (ALT) and aspartate aminotransferase (AST) increased (Figure 9C). Given that lactate is the end product of glycolysis, we measured lactate levels in both plasma and liver tissues of mice. Our results demonstrated a concomitant elevation in lactate levels with the progression of HF (Figure 9D). IHC staining for the core GRGs demonstrated a clear and significant upregulation of these proteins in fibrotic liver tissues compared to oil controls (Figure 9E). RT-PCR analyses demonstrated compared to oil controls, CCl4-induced fibrotic mice exhibited significant elevation of B3gnt3, Chst4, Dcn, Gpc3, Sox9, and Vcan (all p<0.05) (Figure 9F).

Figure 9
Experimental study results on mice liver fibrosis depicting various analyses. Panel A outlines the experimental setup with oil and CCl₄ injections. Panel B shows stained liver tissue sections with H&E, Masson's trichrome, and α-SMA staining at different time points. Panels C and D present bar graphs of ALT, AST, and lactate measurements in IU/L and mmol/L. Panel E displays liver sections with immunostaining for B3GNT3, CHST4, DCN, GPC3, VCAN, and SOX9. Panel F features a graph of relative mRNA levels for α-SMA, B3GNT3, CHST4, DCN, GPC3, SOX9, and VCAN with statistical significance marked.

Figure 9. Glycolysis and core GRGs expression in vivo model of liver fibrosis in mice. (A) Scheme of protocol. (B) Paraffin sections were stained with H&E, Masso’s staining, and anti-α-SMA. (C) Plasma ALT and AST levels. (D) Lactate levels in plasma and liver tissue. (E) Core GRGs were examined by immunohistochemistry. (F) Core GRGs were examined by qPCR. Mice, n = 3/group. Data are presented as mean ± SD. *p < 0.05; **p < 0.01; ***p < 0.001.

3.11 Glycolysis and core GRGs expression in LX-2 cells

Given that HSC activation is a central event in HF, we further investigated alterations in glycolysis and the expression of core GRGs during LX-2 cell activation. Following TGF-β stimulation, LX-2 cells exhibited increased glucose uptake (Figure 10A), reduced supernatant glucose levels, and elevated lactate production (Figure 10B). Treatment with the glycolytic inhibitor 2-DG suppressed LX-2 cell migration (Figure 10C), proliferation (Figure 10D), glucose consumption and lactate production (Figure 10E). Finally, we examined the expression of core GRGs (Figure 10F). TGF-β stimulation significantly upregulated the expression of B3GNT3,GPC3, and SOX9 (p < 0.05).

Figure 10
Panel A displays glucose uptake activity in control and TGF-β treated cells, with higher positive cells in TGF-β. Panel B shows bar graphs of extracellular lactate and glucose levels, with significant increases in TGF-β. Panel C includes images of cell migration and bar charts indicating increased migration and EDU positivity in TGF-β treated cells. Panel D presents EDU and DAPI staining, with more EDU-positive cells in TGF-β conditions. Panel E shows bar graphs of extracellular lactate and glucose with significant differences across control, TGF-β, and TGF-β+2-DG conditions. Panel F displays relative mRNA levels of various markers, showing significant changes upon TGF-β treatment.

Figure 10. Glycolysis and core GRGs expression in TGF-β-stimulated LX-2 cells. (A, B) Glycolytic alterations in LX-2 cells after 48-hour stimulation with TGF-β (10 ng/ml). (A) Representative images of glucose uptake with quantifications; (B) Glucose consumption and lactate production in LX-2 cell supernatants. (C-E) Effects of the glycolytic inhibitor 2-DG on LX-2 cell proliferation, migration, and glycolysis. (C) Representative images of transwell migration assay with quantifications; (D) Representative images of EdU assay (proliferation) with quantifications; (E) Glucose consumption and lactate production; (F) Core GRGs were examined by qPCR. n = 3/group. Data are presented as mean ± SD. *p < 0.05; **p < 0.01; ***p < 0.001.

4 Discussion

This investigation established the prognostic relevance of GRGs in HF progression through integrated analysis of GEO datasets. Initial screening prioritized 15 candidate marker genes (B3GNT3, CHST4, DCN, GPC3, GPC4, HIF1A, HK1, HKDC1, LDHB, SLC2A1, SOX9, TFF3, TGFBI, TPBG, and VCAN). Subsequent validation via six machine learning algorithms consistently identified six core GRGs (B3GNT3, CHST4, DCN, GPC3, SOX9, and VCAN) as robust predictors of HF progression. Immunological landscape analysis further uncovered significant correlations between GRG expression and immune cell infiltration within the fibrotic microenvironment. These findings were experimentally confirmed in a CCl4-induced murine HF model, where all six core GRGs demonstrated marked upregulation in fibrotic liver tissues (P<0.05). Collectively, our results delineate a glycolytic regulatory axis driving fibrogenesis and provide mechanistic insights for therapeutic targeting.

Metabolic reprogramming toward enhanced glycolysis represents a pivotal mechanism in fibrogenesis, facilitating fibroblast activation and pathological extracellular matrix deposition (33). Our findings align with prior studies demonstrating SOX9 upregulation in fibrotic tissues, where it promotes myofibroblast differentiation as a key pathological feature of fibrosis (34). Notably, SOX9 modulates disease progression in various pathologies by regulating glycolysis (35, 36). Our findings support this, as we observed a significant SOX9 dysregulation in HF, suggesting a crucial role of glycolysis in the fibrotic process. We additionally identified VCAN as an orchestrator of fibroblast migration/proliferation and collagen deposition (37). DCN, a minor chondroitin-dermatan sulfate proteoglycan in normal liver ECM, becomes progressively upregulated during fibrogenesis and serves as a structural component in cirrhotic ECM scaffolding (38). GPC3 is a heparan sulfate proteoglycan. Although its expression is elevated in HF, numerous studies have demonstrated its utility as a biomarker for HCC (39, 40). The relationship between glycosyltransferase B3GNT3 and carbohydrate sulfotransferase CHST4 in HF remains unreported; however, both enzymes exhibit upregulated expression in cancers (41, 42). Since elevated glycolytic activity is similarly observed in both cancer and HF, the upregulation of these genes in HF may provide a plausible explanation for this shared metabolic alteration.

The distinct immune landscapes associated with glycolysis in liver fibrosis, revealed by bulk and single-cell analyses, reflect complementary biological processes within the chronic inflammatory microenvironment. The broad, positive correlation between glycolytic levels and numerous activated and memory lymphocyte populations in bulk tissue signifies a systemic state of immune engagement. This pattern is characteristic of persistent antigenic stimulation, where adaptive immune cells undergo activation, clonal expansion, and differentiation (43). The concurrent rise in myeloid and innate cells like MDSCs and pDCs further illustrates a coordinated, multi-lineage immune response (44). Previous studies have demonstrated that glycolytic signatures, including LDHA, promote the accumulation and immunosuppressive function of MDSCs. In tumors, glycolytic metabolism orchestrates a molecular network involving the AMPK-ULK1 axis, autophagy, and the transcription factor CEBPB to sustain this MDSC-mediated immunosuppression (45). Single-cell analysis refines this view by identifying the specific metabolic states of key cellular players. The high glycolytic flux in CD8+ exhausted T cells represents a critical adaptation to chronic activation (46). Similarly, the elevated glycolysis in plasma cells is a fundamental requirement for their role as antibody factories, meeting the substantial biosynthetic demands of high-rate protein secretion (47). The convergence of both datasets on mast cells confirms their active participation in the glycolytic milieu, likely fueling their rapid degranulation and cytokine production (48). The finding that naive T cells also display a glycolytic phenotype suggests that the inflammatory microenvironment can impose metabolic reprogramming even on quiescent cells, potentially priming them for future activation (49).

This study has several limitations. First, the functional relationship between glycolytic reprogramming and fibrosis remains primarily correlative in our study. Future work necessitates genetic perturbation of key GRGs in relevant cell models to establish direct causality. Second, although we validated core GRGs in mouse liver, their therapeutic potential remains unexplored. Subsequent studies will employ liver-specific AAV or Cre-loxP systems to modulate these genes in both preventive and therapeutic regimens, with parallel assessment of fibrosis and glycolytic flux. Finally, transcriptome-based immune infiltration algorithms require complementary validation via experimental approaches such as immunohistochemistry and mass cytometry performed directly on liver tissues.

5 Conclusions

This study holds significant translational potential. Elevated expression of B3GNT3, CHST4, DCN, GPC3, SOX9, and VCAN may represent therapeutic targets for HF through modulation of glycolytic pathways. The glycolysis-derived model robustly distinguished fibrotic progression stages, enabling early identification of high-risk patients and supporting personalized therapeutic stratification. Collectively, these findings expand the mechanistic understanding of metabolic dysregulation in fibrosis and provide actionable insights for clinical intervention.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: GEO database, https://www.ncbi.nlm.nih.gov/geo/.

Ethics statement

The animal study was approved by Institutional Animal Care and Use Committee of Chongqing Medical University (Approval No. IACUC-SAHCQMU-2025-0116). The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

GX: Conceptualization, Data curation, Formal analysis, Methodology, Visualization, Writing – original draft, Writing – review & editing. ZD: Formal analysis, Investigation, Methodology, Writing – review & editing. KQ: Data curation, Formal analysis, Methodology, Writing – review & editing. AL: Supervision, Visualization, Writing – review & editing. XY: Methodology, Software, Writing – review & editing. HR: Conceptualization, Funding acquisition, Investigation, Project administration, Resources, Supervision, Validation, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research and/or publication of this article. This work was supported by the Remarkable Innovation-Clinical Research Project and the Joint Project of Pinnacle Disciplinary Group of The Second Affiliated Hospital of Chongqing Medical University, the first batch of key disciplines on public health in Chongqing, Health Commission of Chongqing, China.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1684937/full#supplementary-material

References

1. Lee YA, Wallace MC, and Friedman SL. Pathobiology of liver fibrosis: A translational success story. Gut. (2015) 64:830–41. doi: 10.1136/gutjnl-2014-306842

PubMed Abstract | Crossref Full Text | Google Scholar

2. Moon AM, Singal AG, and Tapper EB. Contemporary epidemiology of chronic liver disease and cirrhosis. Clin Gastroenterol Hepatol. (2020) 18:2650–66. doi: 10.1016/j.cgh.2019.07.060

PubMed Abstract | Crossref Full Text | Google Scholar

3. Peng H, Wang Q, Lou T, Qin J, Jung S, Shetty V, et al. Myokine mediated muscle-kidney crosstalk suppresses metabolic reprogramming and fibrosis in damaged kidneys. Nat Commun. (2017) 8:1493. doi: 10.1038/s41467-017-01646-6

PubMed Abstract | Crossref Full Text | Google Scholar

4. Gilgenkrantz H, Mallat A, Moreau R, and Lotersztajn S. Targeting cell-intrinsic metabolism for antifibrotic therapy. J Hepatol. (2021) 74:1442–54. doi: 10.1016/j.jhep.2021.02.012

PubMed Abstract | Crossref Full Text | Google Scholar

5. Zhao X, Psarianos P, Ghoraie LS, Yip K, Goldstein D, Gilbert R, et al. Metabolic regulation of dermal fibroblasts contributes to skin extracellular matrix homeostasis and fibrosis. Nat Metab. (2019) 1:147–57. doi: 10.1038/s42255-018-0008-5

PubMed Abstract | Crossref Full Text | Google Scholar

6. Xie N, Tan Z, Banerjee S, Cui H, Ge J, Liu RM, et al. Glycolytic reprogramming in myofibroblast differentiation and lung fibrosis. Am J Respir Crit Care Med. (2015) 192:1462–74. doi: 10.1164/rccm.201504-0780OC

PubMed Abstract | Crossref Full Text | Google Scholar

7. Wang F, Chen L, Kong D, Zhang X, Xia S, Liang B, et al. Canonical wnt signaling promotes hsc glycolysis and liver fibrosis through an ldh-a/hif-1α Transcriptional complex. Hepatol (Blt Md). (2024) 79:606–23. doi: 10.1097/hep.0000000000000569

PubMed Abstract | Crossref Full Text | Google Scholar

8. Hou W and Syn WK. Role of metabolism in hepatic stellate cell activation and fibrogenesis. Front Cell Dev Biol. (2018) 6:150. doi: 10.3389/fcell.2018.00150

PubMed Abstract | Crossref Full Text | Google Scholar

9. Yin KL, Li M, Song PP, Duan YX, Ye WT, Tang W, et al. Unraveling the emerging niche role of hepatic stellate cell-derived exosomes in liver diseases. J Clin Trans Hepatol. (2023) 11:441–51. doi: 10.14218/jcth.2022.00326

PubMed Abstract | Crossref Full Text | Google Scholar

10. Liberti MV and Locasale JW. The warburg effect: how does it benefit cancer cells? Trends Biochem Sci. (2016) 41:211–8. doi: 10.1016/j.tibs.2015.12.001

PubMed Abstract | Crossref Full Text | Google Scholar

11. Gilgenkrantz H, Sayegh RA, and Lotersztajn S. Immunoregulation of liver fibrosis: new opportunities for antifibrotic therapy. Annu Rev Pharmacol Toxicol. (2025) 65:281–99. doi: 10.1146/annurev-pharmtox-020524-012013

PubMed Abstract | Crossref Full Text | Google Scholar

12. Lefere S and Tacke F. Macrophages in obesity and non-alcoholic fatty liver disease: crosstalk with metabolism. JHEP Rep. (2019) 1:30–43. doi: 10.1016/j.jhepr.2019.02.004

PubMed Abstract | Crossref Full Text | Google Scholar

13. Voss K, Hong HS, Bader JE, Sugiura A, Lyssiotis CA, and Rathmell JC. A guide to interrogating immunometabolism. Nat Rev Immunol. (2021) 21:637–52. doi: 10.1038/s41577-021-00529-8

PubMed Abstract | Crossref Full Text | Google Scholar

14. Caputa G, Castoldi A, and Pearce EJ. Metabolic adaptations of tissue-resident immune cells. Nat Immunol. (2019) 20:793–801. doi: 10.1038/s41590-019-0407-0

PubMed Abstract | Crossref Full Text | Google Scholar

15. Rao J, Wang H, Ni M, Wang Z, Wang Z, Wei S, et al. Fstl1 promotes liver fibrosis by reprogramming macrophage function through modulating the intracellular function of pkm2. Gut. (2022) 71:2539–50. doi: 10.1136/gutjnl-2021-325150

PubMed Abstract | Crossref Full Text | Google Scholar

16. Yang H, Atak D, Yuan M, Li M, Altay O, Demirtas E, et al. Integrative proteo-transcriptomic characterization of advanced fibrosis in chronic liver disease across etiologies. Cell Rep Med. (2025) 6:101935. doi: 10.1016/j.xcrm.2025.101935

PubMed Abstract | Crossref Full Text | Google Scholar

17. Wang M, Gong Q, Zhang J, Chen L, Zhang Z, Lu L, et al. Characterization of gene expression profiles in hbv-related liver fibrosis patients and identification of itgbl1 as a key regulator of fibrogenesis. Sci Rep. (2017) 7:43446. doi: 10.1038/srep43446

PubMed Abstract | Crossref Full Text | Google Scholar

18. Moylan CA, Pang H, Dellinger A, Suzuki A, Garrett ME, Guy CD, et al. Hepatic gene expression profiles differentiate presymptomatic patients with mild versus severe nonalcoholic fatty liver disease. Hepatol (Blt Md). (2014) 59:471–82. doi: 10.1002/hep.26661

PubMed Abstract | Crossref Full Text | Google Scholar

19. Ramachandran P, Dobie R, Wilson-Kanamori JR, Dora EF, Henderson BEP, Luu NT, et al. Resolving the fibrotic niche of human liver cirrhosis at single-cell level. Nature. (2019) 575:512–8. doi: 10.1038/s41586-019-1631-3

PubMed Abstract | Crossref Full Text | Google Scholar

20. Hoang SA, Oseini A, Feaver RE, Cole BK, Asgharpour A, Vincent R, et al. Gene expression predicts histological severity and reveals distinct molecular profiles of nonalcoholic fatty liver disease. Sci Rep. (2019) 9:12541. doi: 10.1038/s41598-019-48746-5

PubMed Abstract | Crossref Full Text | Google Scholar

21. Mas VR, Maluf DG, Archer KJ, Yanek K, Kong X, Kulik L, et al. Genes involved in viral carcinogenesis and tumor initiation in hepatitis C virus-induced hepatocellular carcinoma. Mol Med (Camb Mass). (2009) 15:85–94. doi: 10.2119/molmed.2008.00110

PubMed Abstract | Crossref Full Text | Google Scholar

22. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for rna-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | Crossref Full Text | Google Scholar

23. Yu G, Wang LG, Han Y, and He QY. Clusterprofiler: an R package for comparing biological themes among gene clusters. Omics. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | Crossref Full Text | Google Scholar

24. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | Crossref Full Text | Google Scholar

25. Zhang J, Gu J, Guo S, Huang W, Zheng Y, Wang X, et al. Establishing and validating a pathway prognostic signature in pancreatic cancer based on mirna and mrna sets using gsva. Aging. (2020) 12:22840–58. doi: 10.18632/aging.103965

PubMed Abstract | Crossref Full Text | Google Scholar

26. Seiler M, Huang CC, Szalma S, and Bhanot G. Consensuscluster: A software tool for unsupervised cluster discovery in numerical data. Omics. (2010) 14:109–13. doi: 10.1089/omi.2009.0083

PubMed Abstract | Crossref Full Text | Google Scholar

27. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd, et al. Comprehensive integration of single-cell data. Cell. (2019) 177:1888–902.e21. doi: 10.1016/j.cell.2019.05.031

PubMed Abstract | Crossref Full Text | Google Scholar

28. Zhang C, Li J, Cheng Y, Meng F, Song JW, Fan X, et al. Single-cell rna sequencing reveals intrahepatic and peripheral immune characteristics related to disease phases in hbv-infected patients. Gut. (2023) 72:153–67. doi: 10.1136/gutjnl-2021-325915

PubMed Abstract | Crossref Full Text | Google Scholar

29. Zhang X, Lan Y, Xu J, Quan F, Zhao E, Deng C, et al. Cellmarker: A manually curated resource of cell markers in human and mouse. Nucleic Acids Res. (2019) 47:D721–d8. doi: 10.1093/nar/gky900

PubMed Abstract | Crossref Full Text | Google Scholar

30. Xie A, Wang H, Zhao J, Wang Z, Xu J, and Xu Y. Scpas: single-cell phenotype-associated subpopulation identifier. Briefings Bioinf. (2024) 26. doi: 10.1093/bib/bbae655

PubMed Abstract | Crossref Full Text | Google Scholar

31. Liao X, Liu L, Du J, Fan H, Yu Y, Luo Y, et al. Mettl3 inhibits liver fibrosis via rbx1 stability and txnip-mediated ferroptosis. Hepatol Commun. (2025) 9:e0752. doi: 10.1097/hc9.0000000000000752

PubMed Abstract | Crossref Full Text | Google Scholar

32. Jiang Y, Han S, Cheng W, Wang Z, and Wu A. Nfat1-regulated il6 signalling contributes to aggressive phenotypes of glioma. Cell Comm Signaling. (2017) 15:54. doi: 10.1186/s12964-017-0210-1

PubMed Abstract | Crossref Full Text | Google Scholar

33. Baudo G, Wu S, Massaro M, Liu H, Lee H, Zhang A, et al. Polymer-functionalized mitochondrial transplantation to fibroblasts counteracts a pro-fibrotic phenotype. Int J Mol Sci. (2023) 24:10913. doi: 10.3390/ijms241310913

PubMed Abstract | Crossref Full Text | Google Scholar

34. Trogisch FA, Abouissa A, Keles M, Birke A, Fuhrmann M, Dittrich GM, et al. Endothelial cells drive organ fibrosis in mice by inducing expression of the transcription factor sox9. Sci Trans Med. (2024) 16:eabq4581. doi: 10.1126/scitranslmed.abq4581

PubMed Abstract | Crossref Full Text | Google Scholar

35. Jain L, Jardim CA, Yulo R, Bolam SM, Monk AP, Munro JT, et al. Phenotype and energy metabolism differ between osteoarthritic chondrocytes from male compared to female patients: implications for sexual dimorphism in osteoarthritis development? Oa Cart. (2024) 32:1084–96. doi: 10.1016/j.joca.2023.09.013

PubMed Abstract | Crossref Full Text | Google Scholar

36. Yan F, Teng Y, Li X, Zhong Y, Li C, Yan F, et al. Hypoxia promotes non-small cell lung cancer cell stemness, migration, and invasion via promoting glycolysis by lactylation of sox9. Cancer Biol Ther. (2024) 25:2304161. doi: 10.1080/15384047.2024.2304161

PubMed Abstract | Crossref Full Text | Google Scholar

37. Bukong TN, Maurice SB, Chahal B, Schaeffer DF, and Winwood PJ. Versican: A novel modulator of hepatic fibrosis. Lab Invest. (2016) 96:361–74. doi: 10.1038/labinvest.2015.152

PubMed Abstract | Crossref Full Text | Google Scholar

38. Jármay K, Gallai M, Karácsony G, Ozsvár Z, Schaff Z, Lonovics J, et al. Decorin and actin expression and distribution in patients with chronic hepatitis C following interferon-alfa-2b treatment. J Hepatol. (2000) 32:993–1002. doi: 10.1016/s0168-8278(00)80104-x

PubMed Abstract | Crossref Full Text | Google Scholar

39. Llovet JM, Chen Y, Wurmbach E, Roayaie S, Fiel MI, Schwartz M, et al. A molecular signature to discriminate dysplastic nodules from early hepatocellular carcinoma in hcv cirrhosis. Gastroenterology. (2006) 131:1758–67. doi: 10.1053/j.gastro.2006.09.014

PubMed Abstract | Crossref Full Text | Google Scholar

40. Di Tommaso L, Franchi G, Park YN, Fiamengo B, Destro A, Morenghi E, et al. Diagnostic value of hsp70, glypican 3, and glutamine synthetase in hepatocellular nodules in cirrhosis. Hepatol (Blt Md). (2007) 45:725–34. doi: 10.1002/hep.21531

PubMed Abstract | Crossref Full Text | Google Scholar

41. Li CW, Lim SO, Chung EM, Kim YS, Park AH, Yao J, et al. Eradication of triple-negative breast cancer cells by targeting glycosylated pd-L1. Cancer Cell. (2018) 33:187–201.e10. doi: 10.1016/j.ccell.2018.01.009

PubMed Abstract | Crossref Full Text | Google Scholar

42. Jinawath N, Chamgramol Y, Furukawa Y, Obama K, Tsunoda T, Sripa B, et al. Comparison of gene expression profiles between opisthorchis viverrini and non-opisthorchis viverrini associated human intrahepatic cholangiocarcinoma. Hepatol (Blt Md). (2006) 44:1025–38. doi: 10.1002/hep.21330

PubMed Abstract | Crossref Full Text | Google Scholar

43. Hoytema Van Konijnenburg DP, Reis BS, Pedicord VA, Farache J, Victora GD, and Mucida D. Intestinal epithelial and intraepithelial T cell crosstalk mediates a dynamic response to infection. Cell. (2017) 171:783–94.e13. doi: 10.1016/j.cell.2017.08.046

PubMed Abstract | Crossref Full Text | Google Scholar

44. Thwe PM, Pelgrom LR, Cooper R, Beauchamp S, Reisz JA, D’Alessandro A, et al. Cell-intrinsic glycogen metabolism supports early glycolytic reprogramming required for dendritic cell immune responses. Cell Metab. (2017) 26:558–67.e5. doi: 10.1016/j.cmet.2017.08.012

PubMed Abstract | Crossref Full Text | Google Scholar

45. Li W, Tanikawa T, Kryczek I, Xia H, Li G, Wu K, et al. Aerobic glycolysis controls myeloid-derived suppressor cells and tumor immunity via a specific cebpb isoform in triple-negative breast cancer. Cell Metab. (2018) 28:87–103.e6. doi: 10.1016/j.cmet.2018.04.022

PubMed Abstract | Crossref Full Text | Google Scholar

46. Chen R, Zheng Y, Fei C, Ye J, and Fei H. Machine learning developed a cd8(+) exhausted T cells signature for predicting prognosis, immune infiltration and drug sensitivity in ovarian cancer. Sci Rep. (2024) 14:5794. doi: 10.1038/s41598-024-55919-4

PubMed Abstract | Crossref Full Text | Google Scholar

47. Lam WY, Becker AM, Kennerly KM, Wong R, Curtis JD, Llufrio EM, et al. Mitochondrial pyruvate import promotes long-term survival of antibody-secreting plasma cells. Immunity. (2016) 45:60–73. doi: 10.1016/j.immuni.2016.06.011

PubMed Abstract | Crossref Full Text | Google Scholar

48. Caslin HL, Taruselli MT, Haque T, Pondicherry N, Baldwin EA, Barnstein BO, et al. Inhibiting glycolysis and atp production attenuates il-33-mediated mast cell function and peritonitis. Front Immunol. (2018) 9:3026. doi: 10.3389/fimmu.2018.03026

PubMed Abstract | Crossref Full Text | Google Scholar

49. Ricciardi S, Manfrini N, Alfieri R, Calamita P, Crosti MC, Gallo S, et al. The translational machinery of human cd4(+) T cells is poised for activation and controls the switch from quiescence to metabolic remodeling. Cell Metab. (2018) 28:895–906.e5. doi: 10.1016/j.cmet.2018.08.009

PubMed Abstract | Crossref Full Text | Google Scholar

Keywords: glycolytic reprogramming, hepatic fibrosis, WGCNA, machine learning, immune cell infiltration

Citation: Xiao G, Deng Z, Qiu K, Li A, Yi X and Ren H (2025) Identification of glycolysis-related clusters and immune cell infiltration in hepatic fibrosis progression using machine learning models and experimental validation. Front. Immunol. 16:1684937. doi: 10.3389/fimmu.2025.1684937

Received: 13 August 2025; Accepted: 23 October 2025;
Published: 05 November 2025.

Edited by:

Ana C. Maretti-Mira, University of Southern California, United States

Reviewed by:

Jin-hua He, Guangzhou Panyu Central Hospital, China
Jiawei Yang, Zunyi Medical University, China

Copyright © 2025 Xiao, Deng, Qiu, Li, Yi and Ren. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Hong Ren, cmVuaG9uZzA1MzFAY3FtdS5lZHUuY24=

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.