The gene expression profile and clinicopathological characteristics of patients with DLBCL of GSE181063 (n = 1310), GSE10846 (n = 420), and GSE53786 (n = 119) datasets were downloaded from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). The GSE181063 profile was used as the training cohort and based on the Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip platform (Illumina, Bethesda, MD, United States). The independent external validation cohorts GSE10846 and GSE53786 were based on the Affymetrix Human Genome U133 Plus 2.0 Array platform (Affymetrix, Santa Clara, CA, United States). The inclusion criteria of DLBCL patients were as follows: 1) histologically confirmed DLBCL; 2) had detailed corresponding clinicopathological information; 3) initially treated with RCHOP or CHOP regimen; 4) overall survival (OS) time ≥30 days. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Glycolysis-related gene sets (REACTOME_GLYCOLYSIS, KEGG_GLYCOLYSIS_GLUCONEOGENESIS, HALLMARK_GLYCOLYSIS, BIOCARTA_GLYCOLYSIS_PATHWAY, BIOCARTA_FEEDER_PATHWAY, GOBP_POSITIVE_REGULATION_OF_GLYCOLYTIC_PROCESS, GOBP_NEGATIVE_REGULATION_OF_GLYCOLYTIC_PROCESS, GOBP_GLYCOLYTIC_PROCESS_THROUGH_FRUCTOSE_6_PHOSPHATE) were obtained from the Molecular Signatures Database (MSigDB v7.5.1, http://www.broad.mit.edu/gsea/msigdb/).

Cox and least absolute shrinkage and selection operator (LASSO) regression analyses were used in the training cohort to identify the glycolysis-related prognostic model. Firstly, we performed the univariate Cox regression analysis with the “survival” package in R to identify survival-associated glycolytic genes. Genes with p < 0.01 were selected. Then, LASSO regression was carried out to further select key glycolytic genes by using the “glmnet” package in R. Finally, the multivariate Cox analysis was performed to establish the glycolytic risk model. The risk model was constructed by each glycolytic gene expression value (G) and the corresponding coefficient (β). The r i s k   s c o r e = ∑ n = 1 j β j * G j . Patients with DLBCL were classified into low- and high-risk groups according to the median value of the risk score. Heatmap of the risk gene expression, risk score distribution plots, and survival status scatter plots of the low- and high-risk groups were established to assess the risk model. The predictive accuracy of the risk model was further evaluated by Kaplan-Meier (K-M) analyses and time-dependent receiver operating characteristic (ROC) curves by using the “survival”, “survminer”, and “timeROC” packages in R.

GSE181063 dataset was selected as the training cohort. The other two independent datasets, GSE510846 and GSE53786 were used as external validation cohorts. Additionally, samples in the GSE181063 dataset were randomly divided into two internal validation cohorts at a 1:1 ratio by using the “caret” package in R. Clinical features of the training and validation cohorts are shown in Table 1. The comparison of clinical characteristics between two internal validation cohorts is summarized in Supplementary Table S1. The performance of the risk model was assessed in all validation cohorts. p < 0.05 was considered as the significance threshold.

To explore whether the risk model could accurately predict prognosis in the subgroup of clinical parameters including age, gender, Ann Arbor stage, ECOG-PS, lactate dehydrogenase (LDH) level, COO (cell of origin), number of extranodal sites, and B symptoms, we compared the survival between different risk groups based on the stratification of clinical parameters. To further figure out the relationship between the risk model and clinical factors, we compared the risk score distribution in subgroups of the clinical parameters by using the “ggpubr” package in R.

Univariate and multivariate Cox regression analyses were applied in the training and validation (GSE10846) cohorts to determine whether the risk score and clinical parameters were independent prognostic factors in DLBCL patients. Next, we developed a prognostic nomogram to integrate independent prognostic factors with “rms”, “foreign”, and “survival” packages in R. According to the risk contribution to OS, prognostic factors were assigned score points. The total score points were calculated to predict the 1-, 3-, and 5-year prognosis of DLBCL patients. The Harrell’s concordance index (C-index) and calibration curves were performed to evaluate the consistency of the actual and predicted survival probabilities of patients by using the R packages “rms”, “foreign”, and “survival”. ROC curves were drawn and the area under the ROC curve (AUC) was calculated to assess the prognostic performance of nomogram and IPI by using the “survivalROC” package in R.

TME consists of tumor cells, immune cells, stromal cells, and non-cellular components (Bader et al., 2020). Immune scores and stromal scores were calculated in each sample using the “estimate” package in R software (Yoshihara et al., 2013). They represented the proportion of stromal and immune cells in TME and were incorporated as ESTIMATE scores, to further infer the tumor purity in each DLBCL patient. The survival was compared between different score levels with a t.test.

To further explore immune cells (encompass activated CD8 T cells, activated CD4 T cells, activated B cells, activated dendritic cells, natural killer cells, and macrophages) infiltration profile in DLBCL patients, the single-sample gene set enrichment analysis (ssGSEA) was performed based on gene expressions of immune cell-specific markers to quantify the enrichment level of each immune cell type with “GSVA” package in R (Bindea et al., 2013). The correlation between the immune cells infiltration and risk model was evaluated with Pearson correlation coefficient and displayed by using the “ggplot2” package in R.

Immune checkpoint inhibitors (ICIs) prevent checkpoint proteins from binding with the partner proteins and prompt T cells to kill cancer cells (Galluzzi et al., 2020). The immune checkpoint molecules, including PD-1, PD-L1, PD-L2, CTLA-4, LAG3, TIM-3, TIGIT, B7-H3, B7-H4, and CD47, were considered promising targets in B-cell lymphoma in previous studies (Ansell et al., 2009; Lee et al., 2017; Advani et al., 2018; Xu-Monette et al., 2018; Jiang et al., 2019; Hatic et al., 2021; Shimada et al., 2021; Lee et al., 2022). The association between immune checkpoint molecule expression and risk score were assessed with the Pearson correlation coefficient and presented using the “corrplot” package in R to explore the potential effect of our glycolytic risk model in tumor immunotherapy.

The immunohistochemical staining of the glycolytic risk genes in lymphoma patients and normal controls was obtained from the Human Protein Atlas (HPA) database (http://www.proteinatlas.org/) (Liu et al., 2022).

Details of DLBCL cell lines information were downloaded from Cancer Cell Line Encyclopedia (CCLE, https://portals.broadinstitute.org/ccle/data), including cell lines annotations, mRNA expression, and metabolites. The correlation analysis between gene expression and metabolites were performed. Metabolites with p-value < 0.05 and |cor| > 0.5 were selected in this study and shown in the heat map using “ggplot2” package in R software (Sui et al., 2022).

To predict the drug sensitivity based on the heterogenicity of tumor samples, the half-maximal inhibitory concentration (IC50) of DLBCL patients in high- and low-risk groups was estimated by using the “pRRophetic” package in R. The potential anti-cancer compounds were screened based on the drug response data in Genomics of Drug Sensitivity in Cancer (GDSC) database (Wang et al., 2021).

All data analyses were performed in R version 4.0.3. p < 0.05 was considered statistically significant.

The study included 947 eligible patients from the training cohort (GSE181063, n = 559) and two external validation cohorts (GSE10846, n = 305; GSE53786, n = 83) (Lenz et al., 2008; Scott et al., 2014; Lacy et al., 2020). Patient baseline characteristics are shown in Table 1. The training cohort was randomly allocated into two internal validation cohorts at a 1:1 ratio. There were no significant differences in clinical factors between the two groups (Chi-square test, p > 0.05). Results are shown in Supplementary Table S1.

We performed the univariate Cox regression analysis with the gene expression profile and survival information of DLBCL patients in the training cohort and identified 14 prognostic genes related to OS (p < 0.01), including eight protective genes (hazard ratio (HR) < 1.0) and six risk genes (HR > 1.0, Figure 1A). Lasso regression (Figures 1B, C) and multivariate analyses (Figure 1D) were further applied to screen prognostic genes. We established a prognostic signature composed of eight glycolysis-related genes in DLBCL, eventually. The glycolytic risk score for OS = (ADH1B × 0.17355739) + (ALDH2 × −0.121597937) + (ANGPTL4 × 0.301463871) + (BPGM × 0.396173449) + (CTH × 0.196458723) + (NUP98 × 0.655960453) + (PAM × −0.14314632) + (PLOD2 × −0.144167963). Among these prognostic genes, ADH1B (p < 0.001), ANGPTL4 (p = 0.030), BPGM (p < 0.001), NUP98 (p < 0.001), and PAM (p = 0.046) were independent risk factors for survival. Based on the risk score, 559 patients were assigned into high-risk (n = 279) and low-risk groups (n = 280) using the median risk score as the cutoff threshold. As revealed by the heatmap, ADH1B, ANGPTL4, BPGM, CTH, and NUP98 were significantly up-regulated while ALDH2, PAM, and PLOD2 were down-regulated in the high-risk group (all p < 0.001, Figure 2A). The risk score distribution and survival status of the individual patient indicated that the poor prognosis was tightly related to the high-risk score (Figure 2A). Patients in the high-risk group had a significantly worse survival than that in the low-risk group (p < 0.001, Figure 2B). The corresponding OS rates at 5-year were 57.2% and 75.7%, respectively. AUCs of the glycolytic risk model were 0.718, 0.695, and 0.688 at 1-, 3-, and 5-year, demonstrating its excellent predictive value (Figure 2C).

To validate the prediction performance of the glycolytic risk model in DLBCL patients, we assessed the gene signature in two external validation cohorts and two internal validation cohorts. Similarly, high-risk patients exhibited a significantly unfavorable prognosis, compared to low-risk patients (Figures 3A, C, E, G; p < 0.001 for internal validation cohorts and GSE10846; p < 0.05 for GSE53786). AUCs at 1-, 3-, 5-year ranged from 0.669 to 0.76 in two internal validation cohorts (Figures 3B, D). The range of AUC values in two external validation cohorts was between 0.62 and 0.70 (Figures 3F, H), which further confirmed the predictive efficacy of the risk model.

To assess the predictive performance of the glycolytic risk model in clinical subgroups, we conducted the clinical parameters stratification based on age (≤60 vs. >60), gender (male vs. female), ECOG-PS (<2 vs. ≥2), AnnAnbor Stage (I-II vs. III-IV), LDH (normal vs. elevated), COO (GCB vs. nonGCB), extranodal sites (<2 vs. ≥2), and B symptoms (no vs. yes) (Figures 4A–P). Overall, high-risk patients displayed poorer prognosis compared with low-risk patients regardless of age, gender, ECOG-PS, stage, LDH level, COO, number of extranodal sites, and with or without B symptoms, indicating our risk model was effective and stable.

The distribution of risk scores was analyzed among various clinical factors in the training cohort (Figures 5A–I). Age>60 (p = 0.0017), ECOG-PS≥2 (p < 0.001), advanced stage (p = 0.0038), elevated LDH (p < 0.001), nonGCB (p < 0.001), and with B symptoms (p = 0.019) subgroups contained more high-risk patients compared with age≤60, ECOG-PS<2, early stage, normal LDH, GCB, and without B symptoms subgroups.

By integrating the clinical factors and risk score, we constructed a nomogram to provide a quantitative method for predicting the 1-, 3-, and 5-year individual survival probability in DLBCL patients. Univariate Cox analysis showed that age, ECOG-PS, LDH, stage, rituximab, COO, and glycolytic risk score were OS-related factors in both training and validation cohorts (GSE10846) (Table 2, all p < 0.05). Multivariate Cox analysis revealed that the glycolytic risk score was an independent prognostic indicator in both cohorts (Table 2, all p ≤ 0.002). A higher total score in the nomogram indicated worse survival (Figure 6A). Calibration curves of the nomogram exhibited good agreement between the predicted value and the actual value (C-index = 0.73 for training cohort, Figure 6B; C- index = 0.73 for validation cohort, Figure 6C). In addition, AUCs of the nomogram at 1-, 3-, and 5-year showed its better predictive performance than IPI (Figures 6D–F for training cohort; Figures 6G–I for validation cohort), suggesting the glycolytic risk model improved the predictive value upon existing clinical predictors basis in DLBCL.

TME was represented by the stromal and immune scores (Liang et al., 2021). Patients in the high-risk group exhibited low immune scores, stromal scores, and ESTIMATE scores (Figures 7A–C; all p < 0.01) and high tumor purities (Figure 7D; p < 0.001), indicating the glycolytic risk model was associated with immunosuppression. The corresponding results of survival analysis (Figures 7E–G, all p ≤ 0.001) were coincident with the high-risk score-related poor prognosis. To explore the underlying mechanism, we analyzed the correlation between the risk model and infiltration levels of immune cells. As shown in Figure 7, the glycolytic risk score was negatively correlated with activated CD8 T cells (Figure 7I; r = −0.12, p = 0.006), activated dendritic cells (Figure 7L; r = −0.11, p = 0.007), natural killer cells (Figure 7M; r = −0.18, p < 0.001), and macrophages (Figure 7N; r = −0.09, p = 0.027), thus might explained why glycolytic risk score had a negative influence on clinical outcomes in DLBLC patients.

Encouraged by plenty of achievements in antitumor immunotherapy, we evaluated the correlations of our glycolytic risk model and the immune checkpoint molecules PD-1, PD-L1, PD-L2, CTLA-4, LAG3, TIM-3, TIGIT, B7-H3, B7-H4, and CD47 in DLBCL. As shown in Figure 8, the glycolytic gene signature was negatively associated with PD-L2 (r = −0.09, p = 0.039), CTLA4 (r = −0.18, p < 0.001), TIM-3 (r = −0.23, p < 0.001), TIGIT (r = −0.09, p = 0.028), B7-H3 (r = −0.09, p = 0.036).

To explore the protein expression of the eight glycolytic risk genes in patients with DLBLC, we compared the protein staining intensity of these genes by using the HPA database. ADH1B, CTH, and PLOD2 showed overexpression in the lymph node tissue of non-Hodgkin lymphoma patients compared to normal lymph node tissue. Unfortunately, the information of NUP98 was not available in the HPA database. No expression difference was found of ALDH2, ANGPTL4, BPGM, and PAM between lymphoma patients and control samples (Supplementary Figure S1).

Based on the GDSC database, we assessed the drug therapeutic responses between high- and low-risk groups. Of note, patients in the high-risk group were more sensitive to cytarabine, methotrexate, vinblastine, gemcitabine, lenalidomide, and rapamycin, while patients in the low-risk group were more sensitive to docetaxel, bortezomib, bleomycin and LFM.A13, which provided the reference directing clinical treatment (Supplementary Figure S3).

GSVA analysis was performed to explore the difference in biological signaling pathways activated/inhibited between high- and low-risk groups. Our results showed significant differences in biological processes between high- and low-risk patients. Notably, MYC targets V2, MYC targets V1, E2F targets, G2M checkpoint, and unfolded protein response signaling were the top enriched signatures in the high-risk group. In comparison, low-risk patients were more relevant to the apical surface, coagulation, epithelial-mesenchymal transition, angiogenesis, and estrogen response early signaling pathways (Supplementary Figure S4).