In this study, we used both bulk RNA sequencing (RNA-seq) and single-cell RNA sequencing (scRNA-seq) data to explore gene expression profiles related to IPF. The bulk RNA-seq data were obtained from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/gds/) database (Goddard et al., 2000; Jorge et al., 1978). The GSE70866 dataset contains 196 samples, including 176 IPF samples and 20 normal samples. Gene expression profiling was carried out on cells from bronchoalveolar lavage (BAL) fluid collected during regular clinical procedures, with normal samples from healthy individuals undergoing bronchoscopy for reasons not related to IPF to ensure proper controls. The GSE218997 dataset consists of 137 samples from the lung tissues of male C57 mice. The subjects were divided into two age groups: young mice aged 8–12 weeks and aged mice at 21 months. Among these samples, 59 were assigned to the bleomycin (BLM)-induced fibrosis group, and the remaining 78 were used as control samples. Single-cell RNA-seq data on IPF and normal samples were also acquired from the GEO database (GSE128033). Lung tissues were collected during transplant surgeries in line with approved protocols, enabling a detailed examination of gene expression at the single-cell level (Morse et al., 2019). However, detailed patient characteristics like age, sex, and treatment status were not provided in these datasets and thus were not considered in our analysis. Specifically, the datasets did not indicate whether IPF patients were undergoing any treatments at the time of sample collection, so the potential effects of treatment on the studied pathways could not be evaluated. Our main focus was to identify gene expression patterns associated with IPF itself (Supplementary Table 1). GRGs were obtained from the Molecular Signatures Database v7.0 (MSigDB) (https://software.broadinstitute.org/gsea/msigdb/index.jsp) (Supplementary Table 2; Supplementary Data 1) (Male et al., 1987; Zacho et al., 1975; Wheeler, 1973).

The “limma” R package was used to analyze the expression differences of GRGs in IPF (Burch and Varkey, 1985). Differentially expressed genes were identified using thresholds of log fold-change (logFC) >1 and adjusted p-value (adj.P.Val) <0.05 to ensure robust selection of significant genes. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.kegg.jp/) pathway enrichment analyses were conducted using the “clusterProfiler” R package (Fenster and Downey, 2000; Vobecky et al., 1988). The threshold for each analysis was set at P < 0.05.

The least absolute shrinkage and selection operator (LASSO) regression was performed under 10-fold cross-validation for deeper disease-specific gene selection. The optimal λ generated a minimum cross-validation error (Friedman et al., 2010). We also applied the support vector machine-recursive feature elimination (SVM-RFE) model using the “SVM” R package. In this model, each feature’s score was sorted, and the next iteration was performed, followed by removing the minimum value until the best features were selected. Glycolysis-related differentially expressed genes (GRDEGs) for IPF were identified by overlapping the results of the above algorithms. Moreover, the diagnostic value of GRDEGs was assessed by calculating receiver operating characteristic (ROC) curves using the “pROC” R package (Robin et al., 2011).

Gene regulatory network analysis comprises two main components: miRNAs and transcription factors (TFs). Then, the intersection between differential genes was identified in the DEG-miRNA and TF-DEG regulatory networks using NetworkAnalyst (https://www.networkanalyst.ca/home.xhtml) (Xia et al., 2015). TarBase (https://www.microrna.gr/tarbase) and miRTarBase (https://mirtarbase.cuhk.edu.cn/) were used to explore DEG-miRNA interaction networks (Karagkouni et al., 2018). The JASPAR database (https://software.broadinstitute.org/gsea/msigdb/index.jsp) was used to analyze TF-DEG interaction networks (Fornes et al., 2020; Vergoulis et al., 2012; Huang et al., 2020). The DGIdb database (https://www.dgidb.org/) was used to identify potential therapeutic agents for IPF (Miyadera and Iwai, 1964).

Candidate GRDEGs were used in the univariate Cox and LASSO regressions to construct the prognostic signature using the “glmnet” R package. The formula risk score = (each gene’s expression × corresponding coefficient) was used to classify IPF patients. Patients with a risk score above the median were placed in the high-risk group, and the rest were placed in the low-risk group. The “timeROC” R package was used to construct 1-, 2- and 3-year ROC curves analyses, and the results were quantified by the area under the ROC curve (AUC). Meanwhile, the Kaplan–Meier survival analysis was conducted using ‘‘survival” and ‘‘survminer” R packages. The universal datasets were divided into train and test sets, keeping the ratio to about 1:1. The above analyses were conducted for the whole, train, and test sets.

The GSEA was conducted using the “GSEABase” R package (Subramanian et al., 2005). According to the log₂[Fold Change (FC)] value from the differential analysis, genes were ranked from high to low to be defined as the test gene set. Then, the KEGG signaling pathway set was used to detect the possible relevance between hub genes and the gene set.

This analysis was complemented using the “GSVA” R package (Hänzelmann et al., 2013). The KEGG pathway set was used to conduct the GSVA of hub genes. The “limma” R package was used to identify the GSVA score differences between high and low-expression samples, with P < 0.05 as the cut-off. The pathway was classified as active if t > 0 in the high expression group and if t < 0 in the low expression group.

CIBERSORT (Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts) is a computational method for estimating the composition of complex tissues from gene expression profiles (Newman et al., 2019). In this study, we applied CIBERSORT to analyze the proportion of 22 types of infiltrating immune cells in the GSE70866 dataset, which comprises gene expression profiles derived from BAL fluid samples of both IPF patients and healthy controls. Additionally, Pearson correlation coefficients were calculated to quantify the association between the relative abundance of immune cells and the expression levels of hub genes.

We obtained the raw data from the GEO database (GSE128033) and processed it using the “Seurat” R package (Hao et al., 2021). The data preprocessing steps included normalization, scaling, and cell clustering, resulting in the identification of four predominant cell types. To ensure data quality, we applied specific criteria: nFeature_RNA >200, nFeature_RNA <5500, percent.mt <10, nCount_RNA >1000, nCount_RNA <35,000 (Ilicic et al., 2016). Single cells were extracted while excluding doublets and dead cells based on these criteria. Next, we performed principal component analysis (PCA) on the highly variable genes and reduced the dimensions of the data. The resulting clusters were identified using the functions “FindNeighbors,” “FindClusters,” and ‘runTSNE'.

We initially identified 14 glycolysis-related genes using LASSO and SVM-RFE. LASSO regression then narrowed the list to seven genes with prognostic significance, minimal redundancy, and reduced multicollinearity. To validate these genes in vivo, we used a bleomycin (BLM)-induced pulmonary fibrosis model in 24 male C57 mice (8–12 weeks old), randomly assigned to either the PBS group (n = 12) or the BLM group (n = 12; 2.5 U/kg, intratracheal). After 21 days, pulmonary fibrosis was confirmed by histology (HE and Masson staining) and Western blot (Col1a1, Fibronectin). We subsequently measured the expression of the seven selected genes by qPCR (primer sequences in Supplementary Table 3). All animal procedures were approved by the Ethics Committee of Wuhan University (No. WP20220442).

R software (version 4.0.3) was used for statistical analyses and visualization of results. ALL group’s difference calculated by Mann-Whitney U test. A P < 0.05 was considered statistically significant. Significance correlation coefficients were defined by absolute values >0.2 and P < 0.05.

The analysis of GRGs in the GSE70866 dataset revealed significant expression differences between IPF patients and healthy controls, including 63 upregulated genes and 55 downregulated genes. Key glycolytic enzymes, such as ENO1 , were significantly upregulated in IPF (P = 0.028), indicating enhanced glycolytic activity, while genes involved in oxidative metabolism, such as DLD and ADH1C , were markedly downregulated (P < 0.05), suggesting impaired energy homeostasis (Figure 1A; Supplementary Data 2). GO analysis further highlighted their localization to key cellular compartments, including the lysosomal lumen and vascular lumen, emphasizing roles in cellular transport and energy regulation (Figure 1B). Functional enrichment analysis demonstrated that differentially expressed GRGs were involved in critical processes, including energy generation, carbohydrate metabolism, and oxidative stress regulation, with significant enrichment in biological pathways such as glycolysis/gluconeogenesis, HIF-1 signaling, and carbon metabolism. These findings indicate a hypoxia-driven metabolic reprogramming in IPF (Figure 1C). Collectively, these results suggest that metabolic alterations in GRGs contribute to fibroblast activation, extracellular matrix remodeling, and disease progression, offering potential targets for therapeutic intervention in IPF.

To identify robust biomarkers for IPF, a combination of LASSO regression and SVM-RFE algorithms was applied to optimize feature selection. LASSO regression analysis, guided by 10-fold cross-validation, selected 37 candidate genes with minimal binomial deviance, as shown in Figure 2A. The corresponding coefficient profiles (Figure 2B) highlight the effect of regularization, where increased λ progressively reduced the number of retained variables. SVM-RFE further refined the selection to 19 genes, achieving the highest cross-validation accuracy (0.969) and the lowest error rate (0.0308), as depicted in Figures 2C, D. The overlap between the LASSO and SVM-RFE methods was visualized in a Venn diagram (Figure 2E), revealing 14 genes shared by both models. These shared genes, including PFKM, AURKA, CXCR4, and SDC1, were prioritized for their diagnostic potential. ROC curve analysis (Figure 2F) demonstrated the diagnostic performance of each gene, with AUC values ranging from 0.704 (PFKM) to 0.887 (VCAN), underscoring their utility in distinguishing IPF patients from controls. Overall, this integrated computational approach identified a robust panel of 14 genes with high diagnostic accuracy, providing a basis for further experimental validation and mechanistic exploration in IPF pathogenesis.

To identify potential therapeutic agents targeting the 14 key genes (PFKM, PKLR, FBP2, AGRN, ARTN, AURKA, CXCR4, DEPDC1, MERTK, PLOD2, SDC1, TPBG, VCAN, and PGP) implicated in IPF, we used the DGIdb. Through Cytoscape for visualization, the analysis showed a comprehensive network of drug-gene interactions, providing practical insights for targeted therapeutic strategies (Figure 3). AURKA was related to a variety of inhibitors, including Tozasertib, Alisertib, and Barasertib, which are known to control mitosis and cellular proliferation (Dickson et al., 2016; Wang-Bishop et al., 2019). CXCR4, a crucial chemokine receptor, was connected to antagonists like Plerixafor, Mavorixafor, and CTCE-9908, which might suppress fibrosis by adjusting immune cell migration and inflammation (Hassan et al., 2011; Kwong et al., 2009). MERTK was recognized as a target of inhibitors such as Paclitaxel and Cisplatin, indicating its function in apoptosis regulation and tissue remodeling (Choi et al., 2020; Kang et al., 2015). Other notable interactions included VCAN, targeted by agents such as Cyclosporine and Balixafortide, which could influence extracellular matrix dynamics (Luo et al., 2023; Choocheep et al., 2010), and SDC1, which was linked to drugs like Heparin and Indatuximab, emphasizing its involvement in cell adhesion and matrix interactions (Schönfeld et al., 2018; Floer et al., 2010). TPBG was targeted by antibody-based therapies such as Naptumomab, underlining its potential in immune modulation (Sun et al., 2024; Hu et al., 2015). Additionally, metabolic regulators FBP2 and PGP were associated with chemotherapeutic agents such as Etoposide and Docetaxel, suggesting their roles in stress response pathways (Kemper et al., 2004; Roy et al., 2014). These discoveries emphasize the therapeutic significance of the identified genes and their druggable targets, providing a basis for precision medicine approaches in IPF.

To elucidate the regulatory mechanisms of the 14 key genes involved in IPF, we used NetworkAnalyst to construct miRNA-gene and TF-gene interaction networks. These analyses show the complex layers of transcriptional and post-transcriptional regulation that might contribute to the pathogenesis of IPF. The miRNA network found many regulatory miRNAs interacting with the target genes, with VCAN, AURKA, and TPBG emerging as central hubs. Key miRNAs (Figure Figure4A), such as hsa-mir-125a-5p, hsa-mir-146a-5p, and hsa-mir-335-5p, were predicted to target multiple genes, indicating their roles as important post-transcriptional regulators (Cheng et al., 2024; Zhou et al., 2024; Li et al., 2014). These miRNAs probably affect processes crucial to IPF, including extracellular matrix remodeling, immune modulation, and cellular proliferation. The TF network revealed major transcriptional regulators such as TP53, STAT1, and GATA2, which are connected to multiple target genes like CXCR4, VCAN, and AURKA (Figure 4B) (Dauch et al., 2016; Feng et al., 2020; Rugge et al., 2010). Additional TFs like PPARG, RELA, and NRF1 are linked to immune responses, oxidative stress pathways, and metabolism, emphasizing their role in coordinating the gene dysregulation in IPF (Koudstaal and Wijsenbeek, 2023; Tian et al., 2023). Genes like CXCR4 and AURKA are identified as key targets for these regulatory TFs, highlighting their functional significance in IPF pathophysiology (Smadja et al., 2014; Jaffar et al., 2020). These findings offer a detailed regulatory framework for the key genes related to IPF, emphasizing the interaction between miRNAs and TFs in regulating their expression. The networks identify potential upstream regulators, providing valuable insights into disease mechanisms and potential therapeutic targets.

We constructed a survival model using LASSO regression based on 7 genes (PFKM, FBP2, ARTN, AURKA, DEPDC1, MERTK, and SDC1) selected from an initial pool of 14, integrating survival data from IPF patients (Figures 5A, B). The results demonstrate the model’s predictive accuracy and effective risk stratification. In the overall cohort, the model demonstrated predictive performance, with AUC values of 0.831, 0.824, and 0.793 at 1, 2, and 3 years, respectively (Figure 5C). Kaplan-Meier survival analysis revealed significant survival differences between high- and low-risk groups (P < 0.001, Figure 5F). After a 1:1 split into training and testing cohorts, the model performed well in the training set, achieving AUC values of 0.878, 0.845, and 0.873 at 1, 2, and 3 years, respectively (Figure 5D). Survival analysis confirmed the model’s ability to stratify patients effectively (P < 0.001, Figure 5G). Validation in the testing cohort yielded consistent AUC values of 0.785, 0.819, and 0.749 at 1, 2, and 3 years, respectively (Figure 5E). Kaplan-Meier analysis further demonstrated significant survival differences between risk groups (P < 0.001, Figure 5H). Risk scores, derived from the survival model, stratified patients into high- and low-risk groups. Risk scores correlated positively with mortality, as visualized in the scatterplots and gene expression heatmaps. In the overall cohort (Figure 6A), as well as the training (Figure 6B) and testing subsets (Figure 6C), higher risk scores were associated with a greater number of deaths. The heatmaps further demonstrated the differential expression of the 7 genes across risk groups, reinforcing their contribution to the model.

To further explore the potential function of marker genes to distinguish IPF from normal samples, we conducted a single-GSEA-KEGG pathway analysis. The top six pathways enriched for each marker gene were presented in Figures 7A–G. ARTN was enriched in cytokine-cytokine receptor interaction, lysine degradation, neuroactive ligand-receptor, glycan biosynthesis, olfactory transduction, and tight junction. AURKA was enriched in the cell cycle, DNA replication, nucleotide excision repair, olfactory transduction, P53 signaling pathway, and the ribosome. FBP2 was also related to the P53 signaling pathway, olfactory transduction, and ribosome. Besides, FBP2 was enriched in other pathways, including adherens junction, colorectal cancer, and the proteasome. DEPDC1 was enriched in DNA replication, homologous recombination, neuroactive ligand-receptor interaction, olfactory transduction, P53 signaling pathway, and the ribosome. MERTK was enriched in bladder cancer, cardiac muscle contraction, cytokine-cytokine receptor interaction, Fc gamma R-mediated phagocytosis, hematopoietic cell lineage, and olfactory transduction. PFKM was enriched for cardiac muscle contraction, glycosaminoglycan biosynthesis, chondroitin sulfa-, maturity-onset diabetes of the young, neuroactive ligand-receptor interaction, olfactory transduction, and ubiquitin-mediated proteolysis. Finally, SDC1 was enriched in bladder cancer, neuroactive ligand-receptor interaction, nod-like receptor signaling pathway, olfactory transduction, pathogenic Escherichia coli infection, and sphingolipid metabolism.

Moreover, the GSVA showed that the low ARTN expression might induce IPF by activating ECM receptor interaction, ascorbate and aldarate metabolism, and maturity-onset diabetes of the young, while ARTN overexpression activated mismatch repair, DNA replication, and primary bile acid biosynthesis (Figure 8A). Moreover, AURK upregulation activated the taurine and hypotaurine metabolism and olfactory transduction, while Its downregulation activated DNA replication and homologous recombination (Figure 8B). High DEPDC1 expression activated neuroactive ligand-receptor interaction and phenylalanine metabolism, while low DEPDC1 expression activated DNA replication (Figure 8C). Riboflavin metabolism, steroid hormone biosynthesis, and linoleic acid metabolism were enriched in the high SDC1 expression group, while NOD-like receptor signaling, ECM receptor interaction, glycosaminoglycan biosynthesis chondroitin sulfate, and glycosaminoglycan biosynthesis keratan sulfate were enriched in the low-expression group (Figure 8D). High MERTK expression was related to riboflavin metabolism, while its downregulation was associated with glycosaminoglycan biosynthesis heparan sulfate- (Figure 8E). Notably, for FBP2, various pathways related to IPF pathogenesis were enriched, including limonene and pinene degradation, aminoacyl tRNA biosynthesis, and biosynthesis of unsaturated fatty acids. Low FBP2 expression was related to taurine and hypotaurine metabolism and olfactory transduction (Figure 8F). Furthermore, PFKM, whose expression was limited in IPF, was more closely related to one carbon pool by folate, glycosylphosphatidylinositol GPI anchor biosynthesis, and the mammal circadian rhythm. Noteworthy, high PFKM expression activated the pathways such as taste transduction, maturity-onset diabetes of the young, neuroactive ligand-receptor interaction, and allograft rejection that might induce IPF (Figure 8G).

Then, we evaluated the relationship between the infiltration of immune cells in IPF and normal samples (Figure 9A). Nine of the 22 immune cell types differed between the two groups: naive B cells, resting CD4 memory T cells, activated CD4 memory T cells, monocytes, M1 macrophages, M2 macrophages, activated dendritic cells, activated mast cells, and neutrophils. Next, we assessed the relationship between immune cells and risk scores using Pearson correlation analysis (Figures 9B–I). Naive B cells (r = −0.26, P = 0.00057), M0 macrophages (r = −0.2, P = 0.0093), resting mast cells (r = 0.39, P = 8.4e-08), CD4 memory T cells (r = −0.15, P = 0.049), and NK memory T cells (r = 0.4, P = 7e-08) were negatively correlated with the risk score. Meanwhile, activated mast cells, neutrophils, and activated NK cells were positively associated with risk scores. Pearson’s correlation analysis also revealed the relationships between immune cells and hub genes. MERTK had a positive correlation with activated mast cells and a negative correlation with resting mast cells. SDC1 had a positive correlation with activated NK cells and negative correlations with resting NK cells. AURKA was positively correlated with M0 macrophages and negatively correlated with gamma delta T cells. DEPDC1 was positively correlated with resting CD4 memory T cells. PFKM was negatively correlated with activated CD4 memory T cells (Figure 9J).

To validate the expression trends and diagnostic value of the seven selected genes in IPF, we analyzed an external dataset (GSE218997). Differential expression analysis revealed that Artn, Aurka, Mertk, Pfkm, and Sdcl were significantly upregulated in the lung tissue of mice induced with bleomycin compared to saline controls (P < 0.001, Figures 10A, B, E–G), while Fbp2 also showed increased expression with moderate statistical significance (P = 0.0031, Figure 10D). Depdc1 exhibited a marginally higher expression trend in bleomycin groups, though this did not reach statistical significance (P = 0.053, Figure 10C). These results align with the findings from our primary dataset, supporting the robustness of these genes' expression patterns. ROC analysis further confirmed the diagnostic value of these genes. Artn (AUC: 0.926, Figure10H), Sdcl (AUC: 0.961, Figure10N), and Pfkm (AUC: 0.889, Figure10M) demonstrated excellent diagnostic accuracy, while Aurka (AUC: 0.766, Figure10I) and Mertk (AUC: 0.780, Figure10L) showed good performance. Fbp2 and Depdc1 had lower diagnostic utility, with AUCs of 0.648 and 0.597, respectively(Figures 10K, J).These findings further validate the clinical relevance of the survival model and its utility for IPF diagnosis and prognosis.

Single-cell RNA sequencing analysis was performed to investigate the cell-type-specific expression patterns of the selected genes in IPF. The t-SNE plot (Figure 11A) identified major cell populations, including epithelial cells, fibroblasts, macrophages, and endothelial cells, providing a comprehensive cellular landscape. Gene expression mapping (Figures 11B–F) revealed distinct cellular localization for each gene. DEPDC1 (Figure 11B) and AURKA (Figure 11C) showed expression in macrophages, suggesting their roles in inflammatory and immune responses. SDC1 (Figure 11D) was predominantly expressed in epithelial cells, indicating its involvement in maintaining epithelial integrity and extracellular matrix interactions. PFKM (Figure 11E) and MERTK (Figure 11F) were primarily expressed in fibroblasts and macrophages, implicating them in metabolic pathways and tissue remodeling processes. Violin plots (Figures 11G–K) further quantified gene expression across cell types, confirming these patterns. DEPDC1 and AURKA exhibited significantly higher expression in macrophages, reinforcing their roles in immune regulation. SDC1 was highly expressed in epithelial cells, while PFKM and MERTK were enriched in fibroblasts and macrophages, aligning with their respective roles in metabolism and reparative functions. These findings provide a detailed characterization of the cellular distribution of these genes, highlighting their contributions to IPF pathophysiology and their potential as cell-type-specific therapeutic targets.

In order to further verify the results of cell experiments in vitro and simulate the pathological process of patients with pulmonary fibrosis as much as possible, we successfully constructed a bleomycin-induced pulmonary fibrosis mouse model for in-vivo experiments. The results of HE, Masson staining and Western blot showed that compared with the control group, the bleomycin group had significantly increased lung tissue inflammation and collagen deposition (Figures 12A, B). And Body weight changes showed that compared with the control group, the bleomycin group manifested a significant weight loss (Figure 12C). RT-qPCR results indicated that the glycolysis-related genes Artn, Aurka, Sdc1, Mertk, and Fbp2 were significantly upregulated in bleomycin-induced mouse lung tissue, whereas Pfkm was significantly downregulated. These findings are in agreement with our previous datasets, suggesting that these six differentially expressed glycolysis-related genes may play an important role in the development of pulmonary fibrosis. In contrast, the expression level of Depdc1 was markedly elevated compared to the earlier datasets, which may be attributable to differences in the timing or grouping criteria for bleomycin intervention in the animal models (Figures 12D–J).