DLBCL specimens were obtained through biopsy in the First Affiliated Hospital of Zhengzhou University and frozen at −80°C for storage. All participants provided written informed consent for the use of their specimens in this study. Clinicopathological features of 45 DLBCL patients from the First Affiliated Hospital of Zhengzhou University are performed in Supplementary Table 1 . The study protocol was approved by the ethics committee of the First Affiliated Hospital of Zhengzhou University (ethics number 2021-KY-0835-001).

We systematically explored publicly available DLBCL gene expression datasets with corresponding clinical information of patients from the GEO (https://www.ncbi.nlm.nih.gov/geo/) and TCGA (https://portal.gdc.cancer.gov/) databases. For this study, we gathered a total of 2,335 patients with DLBCL from four cohorts, including GSE117556 (n = 928), GSE10846 (n = 420), GSE31312 (n = 498), and TCGA-NCICCR (n = 489) ( Table 1 ) and our clinical cohort (n = 45). Patients with incomplete transcriptomic data and clinical data were excluded. Ultimately, we included 928 patients from GSE117556, 414 from GSE10846, 470 from GSE31312, and 234 from TCGA-NCICCR. The GSE10846 and GSE31312 datasets used the GPL570 platform, while the GSE117556 dataset used the GPL14951 platform.

In this study, we used univariate Cox analysis to screen genes related to the OS of DLBCL patients in the GSE117556, GSE10846, and GSE31312 datasets. Genes with hazard ratio (HR) >1 and HR <1 were defined as the risk and protective genes. A p-value <0.01 was the cutoff point. The risk genes and protective genes shared by the three datasets were intersected and combined with Least Absolute Shrinkage and Selection Operator (LASSO) regression and multivariate Cox regression to build the final gene risk prediction model (22). We conducted univariate and multivariate Cox regression analyses using the R package “survival.” Another R package “glmnet” was used for the LASSO Cox regression analysis (23). The risk score (RS) of each sample was calculated by multivariate Cox regression analysis. The correlation analysis was based on the R package “corrplot,” and the forest plots for univariate and multivariate analyses were constructed by R package “forestplot.”

We employed two algorithms to assess immune infiltration in DLBCL. CIBERSORT (http://cibersort.stanford.edu/) algorithm was used to obtain the proportion of 22 immune cell types with a threshold of p < 0.05 (24). We applied Single-Sample Gene Set Enrichment Analysis (ssGSEA) to assess the infiltration level of 28 different immune cells in DLBCL expression profile data by the “GSVA” package (25). The R package “ggpubr” was used to visualize differences in the distributions of immune-related cells in the low- and high-risk patient groups from the overall cohort. ”^***,” “^**,” “^*,” and “ns” indicate p < 0.001, p < 0.01, p < 0.05, and not significant, respectively, for the Kruskal–Wallis test.

Total RNA was isolated from the specimens harvested using TRIzol reagent (Invitrogen Corporation, Carlsbad, CA, # A33250). The Prime Script RT reagent kit with genomic DNA eraser (TaKaRa, Tokyo, Japan, #RR037A) was used to synthesize complementary DNA (cDNA). Quantitative real-time polymerase chain reaction (qRT-PCR) analyses were detected by SYBR Green Master Mix (TaKaRa), and the primers for qRT-PCR analyses are listed in Supplementary Table 2 . The 2^−ΔCT method was utilized to calculate the relative mRNA expression of each gene.

Statistical analyses were performed with R (version 3.6.3). The Kaplan–Meier method was used to assess the differences in survival time, and the log-rank test was used to determine the statistical significance. Time-dependent receiver operating characteristic (ROC) curve analysis was used to measure the prognostic performance by comparing the areas under curves (AUC). The nomogram was plotted using the “rms” package (26). The difference in immune infiltration levels between high- and low-risk groups was calculated by Kruskal–Wallis test, and a p-value < 0.05 was statistically significant.

To identify the prognostic signature of DLBCL, we performed a multiple-step analysis ( Figure 1 ). We first screened the GEO database and selected three datasets for univariate Cox proportional hazards regression to identify candidate genes significantly related to OS. We conducted univariate Cox analysis in the GSE117556, GSE10846, and GSE31312 datasets with p < 0.01 as the cutoff value. In total, 1,426, 1,904, and 1,788 candidate protective genes (with hazard ratios (HR) <1) and 890, 2,958, and 2,525 candidate risk genes (with HR >1) were identified in GSE117556, GSE10846, and GSE31312, respectively. A total of 11 genes were candidate protective genes after intersecting the candidate protective genes ( Figure 2A ). Similarly, after matching the candidate risk genes identified in the three datasets, 13 common genes are retained ( Figure 2B ).

The detailed information on these genes is listed in Table 2 . Next, GSE117556 was assigned as the training dataset, and the LASSO was applied to screen the candidate genes, yielding eight genes ( Supplementary Figures 1A, B ). Ultimately, the HRs of the eight genes were acquired by conducting multivariate Cox regression analysis. The forest plot ( Supplementary Figure 1C ) showed that HK2, GAB1, GRPEL1, RCSD1, PLAC8L1, and RASAL1 were risk factors (HR >1), and CAPG and PDPN were protective factors (HR <1) for OS. Finally, the following risk score model was established: risk score = 0.271 × HK2 expression + 0.182 × GAB1 expression + 0.172 × RASAL1 expression − 0.254 × CAPG expression − 0.358 × PDPN expression + 0.362 × GRPEL1 expression + 0.370 × RCSD1 expression + 0.177 × PLAC8L1 expression. Each patient was assigned a risk score with the prognostic model. According to the median risk score, patients in the GSE117556 dataset were stratified into high-and low-risk groups. Of the eight-gene prognostic signature, HK2 have the most significant correlation with worse survival. Through analyzing the RT-PCR results deriving from our clinical specimens, We found HK2 had high expressed in high stage DLBCL ( Supplementary Figure 2A ). To further explore the function of HK2, we first identified genes that correlated with HK2 in the GSE117556 dataset ( Supplementary Figure 2B ). We then utilized positive and negative correlated genes to perform GO and KEGG enrichment analyses. Results revealed that positive correlated genes were involved in cell cycle ( Supplementary Figure 2C ) and negative correlated genes were involved in immune response ( Supplementary Figure 2D ). These results partly reflect genomic differences between high- and low-risk groups.

To further verify the predicting model, we analyzed the survival of high- and low-risk groups. Patients in the low-risk group demonstrated a longer survival time than those in the high-risk group ( Supplementary Figure 3 ). Consistently, the Kaplan–Meier curve indicated that patients in the high-risk group had significantly worse prognoses than low-risk patients (log-rank test p < 0.001) ( Figure 2C ). We further validated the risk score model in the testing dataset GSE 31312 (log-rank test p < 0.001) ( Supplementary Figure 4A ) and GSE 10864 (log-rank test p < 0.001) ( Supplementary Figure 4B ), TCGA-NCICCR dataset (log-rank test p = 0.0052) ( Figure 2D ), and our clinical specimens (log-rank test p = 0.024) ( Figure 2E ), the results were consistent with the training dataset findings. The predictive power of the risk score model was assessed by time-dependent ROC, which yielded good performance in the above datasets. In dataset GSE117556, the AUC values for 1-, 3-, and 5-year overall survival were 0.7, 0.74, and 0.89, respectively. In TCGA-NCICCR dataset, the AUC values for 1-, 3-, and 5-year overall survival predictions were 0.61, 0.60, and 0.62, respectively. Our clinical cohort yielded AUC values of 0.74, 0.79, and 0.78 for the 1-, 3-, and 5-year overall survival, respectively ( Figures 2F–H ). In dataset GSE31312, the AUC values for 1-, 3-, and 5-year overall survival were 0.64, 0.7, and 0.74, respectively. In dataset GSE10864, the AUC values for 1-year overall survival were 0.56 ( Supplementary Figures 4C, D ). The above AUC curves provided an objective validation of the clinical application value of our model.

We used GSE117556 as the training dataset to detect the correlation and interdependence between the eight risk genes, which proved that our prognostic signature could minimize data bias caused by gene collinearity ( Figure 3A ). Using the median risk score as the cutoff point, the training set was classified into low- and high-risk groups. Our results suggest that the risk score model has significant value for evaluating characteristic genes. The eight genes were differentially expressed between the two groups, indicating their role in contributing to the prognosis of DLBCL ( Figure 3B ). Analysis of validation dataset TCGA-NCICCR yielded consistent results ( Figures 3C, D ).

It has been established that DLBCL exhibits a significant heterogeneity in cell origin, clinical manifestations, gene expression profiles, and so on. To verify the effectiveness of the risk signature in the existing clinical subgroups, we conducted a series of subgroup analyses on dataset GSE117556. Stratification based on clinical characteristics such as COO, ECOG, MYC/BCL2 double expression, lactate dehydrogenase (LDH) value, age, stage, and gender was conducted. Among those subtypes, patients in the high-risk groups had worse survival outcomes than patients in the low-risk group (log-rank test p < 0.05; Figure 4 and Supplementary Figure 5 ). These findings validated that our risk model yielded good predictive performance after stratifying for different clinicopathological characteristics.

To further explore the potential survival mechanisms related to the risk score model, mRNA data from dataset GSE117556 were first used to detect the proportion of 22 immune cell types in each sample via the CIBERSORT algorithm. As shown in Figure 5A , the proportion of immune cells was significantly different between the high- and low-risk score groups. Compared with the low-risk group, the high-risk group exhibited increased B-cell infiltrations, with less-naive CD4 T cell, T follicular helper cell, M0 macrophages, M1 macrophages, and other proinflammatory cells. Consistent results were obtained when ssGSEA was applied ( Figure 5B ). In contrast with the high-risk group, tumor-infiltrating lymphocytes (TIL), antigen-presenting cells were significantly enriched in the low-risk groups. We further applied the above two algorithms to the validation dataset TCGA-NCICCR. Similar results were obtained from GSE117556 analysis results. Immune signatures between the high- and low-risk groups were different, and the low-score group was significantly infiltrated with proinflammatory immune cells ( Figures 5C, D ). In addition, we analyzed immune-related genes between the two groups and further approved the above findings ( Supplementary Figure 6 ).

To evaluate the performance of our risk score model on the prognosis of DLBCL patients, we integrated the clinicopathological characteristics with risk score signatures in different algorithms. As shown in Figure 6A , univariate Cox regression analysis demonstrated that the risk score model was a significant predictor of OS in patients with DLBCL (p < 0.0001, HR = 1.377), compared with other clinicopathological characteristics. Multivariate Cox regression analysis showed that the risk score model was an independent prognostic factor for poor prognosis (p < 0.0001, HR = 1.380) ( Figure 6B ). Compared with other indicators, the risk score was superior for predicting the patient prognosis (AUC = 0.820) ( Figure 6C ). Finally, we established a nomogram based on the above clinical parameters to predict patient prognosis quantitatively. Accordingly, our nomogram has huge prospects for clinical application for predicting the OS of individual DLBCL patients ( Figure 6D ).