A flowchart is shown to demonstrate the procedure and results of our study ( Figure 1 ). RNA-seq data of 570 DLBCL patients were obtained from The Cancer Genome Atlas (TCGA) (48 from TCGA-DLBC, 41 from CTSP-DLBCL1, and 481 from NCICCR-DLBCL). Among the 566 DLBCL patients with overall survival (OS) data, 321 (56.71%) were male and 245 (43.29%) were female. The age of the patients ranged from 14 to 92 years (median, 62 years) at initial diagnosis. Other clinical characteristics, including follow-up period, Ann Arbor stages, lactate dehydrogenase (LDH) ratio, Eastern Cooperative Oncology Group (ECOG) performance status, and number of extranodal sites, are documented in Table 1 and Supplementary Material 1 . To remove the batch effect among these three projects, we utilized the ComBat-seq function to transform the raw count data using the sva R package. Then, principal component analysis was performed to show that there was no obvious batch effect among the samples ( Figure 2A ). After excluding samples with unrecorded IPI scores or with an IPI score crossing risk groups (such as 1–5 or 3–4), 458 DLBCL patients were divided into low-risk (n = 118), intermediate-risk (n = 221; 106 at low-intermediate risk, and 92 at high-intermediate risk), and high-risk groups (n = 109) ( Table 1 ). Consistent with previous publications, patients in the low-IPI risk group had a much longer OS ( Figure 2B ).

As shown by gene set variation analysis (GSVA), ECM receptor interaction, CD28-dependent PI3K/AKT signaling, IL-6-type cytokine receptor ligand interaction, and other immunologic signaling pathways were significantly enriched in the low-risk group ( Figure 2C and Supplementary Figures S1A, B ). A total of 4,651 genes (633 upregulated and 4,018 downregulated) were significantly differentially expressed between the IPI high- and low-risk groups ( Supplementary Figures S1C, D and Supplementary Material 2 ). By intersecting with the immunologic signature gene sets (combining 20,837 genes from ImmuneSigDB and Immport, Supplementary Material 3 ), 1,927 immune differentially expressed genes (DEGs) (254 upregulated and 1,673 downregulated) were identified for further analysis ( Figure 2D and Supplementary Figures S1E, F ).

We applied variance-stability-transformed (VST) expression data via DESeq2 as the input data for WGCNA, including 13,329 genes with the top 25% variance among all samples ( Figure 3A ). All clinical characteristics were enrolled as trait variables, and the best β value in the co-expression network was calculated to be 9 ( Supplementary Figure S2B ). The distance threshold for merging modules was set to be 0.30, so as to construct a reasonable number of merged modules ( Supplementary Figures S2A, D ). As shown in the module–trait relationship, eight modules were significantly correlated with the IPI group ( Figure 3B ), and a high correlation (p < 0.0001) between gene significance of IPI risk groups and gene module membership was found in the genes of three modules (brown, pink, dark red) ( Figure 3C and Supplementary Figures S2E, G ). By intersecting 4,106 genes from the top three IPI-correlating modules with 1,927 immune DEGs, a total of 371 genes were identified as immune-related genes associated with IPI, which were used for prognostic risk model construction ( Supplementary Figure S2H ). Several gene sets, such as ECM organization, ECM receptor interaction, PI3K/AKT signaling, and integrin cell surface interaction were enriched in the top three IPI-correlating modules ( Supplementary Figure S2I and Supplementary Material 4 ).

In the training cohort, of which 563 patients with matched RNA-seq data and overall survival follow up data, 93 out of 371 IPI-related immune genes were significantly correlated with OS in the univariate Cox regression analysis ( Table 2 and Supplementary Material 5 ). Next, we applied Lasso-penalized Cox regression to identify the optimal number of genes (n = 10) for the risk score model ( Figure 4A and Supplementary Figure S3A ). As a result of multivariate Cox regression analysis with variable selection by Akaike information criterion, eight genes were selected to construct the most optimal IPI-IPM ( Figure 4B ). Risk score = (expression level of CMBL * [0.360] + expression level of TLCD3B * [-0.350] + expression level of SYNDIG1 * [-0.247] + expression level of ESM1 * [-0.238] + expression level of EPHA3 * [-0.163] + expression level of HUNK * [-0.156] + expression level of PTX3 * [0.138] + expression level of IL12A * [0.111]). As shown in the time-dependent receiver operating characteristic (ROC) curve, areas under the curve (AUC) were 0.703, 0.738, and 0.733 for 1, 3, and 5 years, respectively ( Figure 4C ). As shown in Supplementary Figure S3B , there was a significant difference in OS and PFS between high- and low-risk groups, taking the median of IPI-IPM risk scores as the cutoff value. According to the ROC curve of the median survival time ( Supplementary Figure S3C ), we defined a cutoff value of 0.982 and then divided all patients into high- and low-risk groups according to their risk scores. As shown in Figure 4D , Kaplan–Meier survival analysis showed worse OS of patients in the high-risk score group (log-rank p = 3.13e-14). Similar results showing that high-risk score patients had worse PFS (log-rank p = 6.95e-14) are shown in Figures 4E, F . The C-index for OS of the risk score model was 0.732 (95% CI: 0.684–0.779, p = 1.38e-21), and that for PFS was 0.728 (95% CI: 0.679–0.776, p = 3.11e-20). These results show that IPI-IPM has a good capacity for OS and PFS prediction. The distribution of the risk score, survival status, and the eight-gene expression between the high- and low-risk score groups is shown in Supplementary Figures S3D, E and Figure 4G . In addition, we conducted Spearman’s rank-order analysis to examine the correlation between IPI risk group and the eight-gene expression. As shown in Supplementary Figure S3F , the expression of almost all genes except PTX3 was significantly correlated with IPI risk groups. Then, we tested the difference of the eight-gene expression between high- and low-IPI risk groups. There existed significantly difference of gene expression in seven genes except PTX3.

Moreover, 335 patients with available clinicopathologic parameters were enrolled in the multivariate Cox regression analysis, presenting the risk scores and age as independent prognostic factors of OS ( Figure 4H ). Risk scores, along with age, Ann Arbor clinical stage, LDH ratio, ECOG performance status, and number of extranodal sites, were integrated to construct a nomogram model ( Figure 5A ). As shown in the decision curve analysis (DCA), the nomogram and the risk score from IPI-IPM showed a relatively high net benefit ( Figure 5B ). Moreover, the bias-corrected lines for the nomogram were close to the ideal line in the 1-, 3-, and 5-year and median survival time periods ( Figure 5C ). The C-index for the nomogram was 0.790 (95% CI: 0.736-0.843, p = 2.38e-26). Altogether, these results suggest that the nomogram has excellent capacity and consistency for OS prediction in the training cohort. Moreover, we collected patients with matched IPI-IPM risk score and available IPI information. A total of 445 patients with OS data and 386 patients with PFS data were enrolled to compare the difference of IPI and IPI-IPM in prognostic predictability. As shown in Figures 5D, E , the AUCs of time-dependent ROC for IPI-IPM were almost larger than those for the IPI risk group, indicating fairly equivalent prognostic predictability between IPI and IPI-IPM. In addition, the C-indices for OS of IPI-IPM and IPI were 0.749 (p = 3.78e-23) and 0.756 (p = 1.72e-13), and the C-indices for PFS were 0.739 (p = 6.36e-21) and 0.738 (p = 8.40e-11), respectively.

Three independent cohorts from different clinical centers (the cohort of Staudt et al., GSE10846, n = 414; the cohort of Dubois et al., GES87341, n = 221; the cohort of Sha et al., GSE117556, n = 928) were enrolled for further validation of IPI-IPM. The risk score for each patient was calculated, and all patients were divided into the high- and low-risk groups. As for the cohort of Staudt et al., the AUCs were 0.619, 0.603, and 0.601 for 1, 3, and 5 years, respectively ( Figure 5F ). Kaplan–Meier survival analysis showed significantly shorter OS of patients in the high-risk group (log-rank p = 5.30e-05, Figure 5G and Supplementary Figures S3G, H ). Moreover, patients in the high-risk group were shown to have a shorter OS, regardless of the treatment regimens they received. Similar results of large AUC, shorter OS, and PFS of high-risk patients were also shown in the other two cohorts ( Figures 5H–K and Supplementary Figures S3F, G ). Taking the results of the training and testing cohorts together, the IPI-IPM and the nomogram combined risk score with relevant clinical characteristics (age, Ann Arbor clinical stage, LDH ratio, ECOG performance status, and number of extranodal sites) was an excellent model for predicting short-term or long-term OS in DLBCL patients, which may guide therapeutic strategy decisions and long-term prognosis.

CIBERSORT was applied to analyze the infiltrating abundance of various immune cell types in the different IPI-IPM subgroups ( Figure 9A and Supplementary Figure S6A ). Activated memory CD4⁺ T cells and resting NK cells showed high infiltration in the high-risk group, whereas memory B cells, CD8⁺ T cells, follicular helper T cells, regulatory T cells (Tregs), and non-activated macrophages were more abundant in the low-risk group ( Figure 9B ). In addition, the MCPCounter and xCell algorithms were applied to show that myeloid dendritic cells (mDCs) and common lymphoid progenitors showed high infiltration in the high-risk group, whereas hematopoietic stem cells, cancer-associated fibroblasts (CAFs), and T cells, especially CD8⁺ T cells, were more abundant in the low-risk group ( Figure 9C and Supplementary Figure S6B ).

In addition, DLBCL subtypes varied between the high- and low-risk groups. As shown in Figure 10A , more patients with higher IPI levels were classified into the IPI-IPM high-risk group, with similar results for age, Ann Arbor clinical stage, LDH ratio, ECOG performance status, and the number of extranodal sites ( Supplementary Figure S6C ). The high-risk group contained a higher proportion of activated B-cell-like DLBCLs, whereas the low-risk group was more enriched in germinal center B-cell-like DLBCLs (p = 1.34e-14). Schmitz et al. (4) identified four prominent genetic subtypes in DLBCL with different responses to immunochemotherapy: MCD (the co-occurrence of MYD88^L265P and CD79B mutations), BN2 (BCL6 fusions and NOTCH2 mutations), N1 (NOTCH1 mutations), and EZB (EZH2 mutations and BCL2 translocations), and uneven distributions of these four subtypes were found between the high- and low-risk groups. For example, the poor-prognostic MCD subtype represented 40% of the high-risk group and 12% of the low-risk group. In addition, the correlation analysis showed that IPI-IPM risk scores were positively correlated with the expression of BCL-2 and MYC and negatively correlated with BCL-6 expression ( Figure 10B ). Taken together, these data indicate that IPI-IPM provides additional orthogonal information from previous lymphoma classifications.

To further explore the interaction between DLBCL cells and the microenvironment, Kotlov et al. (11) defined the LME into four major transcriptionally defined categories with distinct biological properties and clinical behavior, including GC, MS, IN, and DP forms. Similarly, we utilized an unsupervised clustering method to assign the samples into four groups by using expression data from the 22 functional gene expression signature (F^GES) sets ( Figure 10C ). GSVA enrichment scores were calculated to demonstrate distinct TME characteristics among the four biological patterns ( Figure 10D ). Consistent with the CIBERSORT results, the DP LME categories represented 35% of the high-risk group whereas the GC-like and MS LME categories were more enriched in the low-risk group ( Figures 10E, F and Supplementary Material 8 ).

To further understand the effects of the risk score on drug response, 184 IPI-IPM-associated immune genes were mapped into the connectivity map database (17). As shown in Figure 11A , 14 genes (ADORA1, ADRA2A, CACNA1C, DBH, DNM1, ELN, ENPP1, EPHA1, IL12A, MMP9, PLAU, S1PR2, SLC1A2, and SV2A) were associated with 122 inhibitors, involving 54 mechanisms of action. The details are documented in Supplementary Material 9 .

The clinical development of cancer immunotherapy and the advances in genomic analysis have validated the important role of the TME in predicting cancer response to ICB therapy (13, 18, 19). We investigated the expression of several inhibitory immune checkpoints between the high- and low-risk groups by using the normalized gene expression data (VST transformed). As shown in Figure 11B , the expression of PD-L1, BTLA, and SIGLEC7 was significantly upregulated in the high-risk group. Expressions of other inhibitory immune checkpoints such as B7-H3 were found downregulated in the high-risk group ( Supplementary Figure S6D ).

Various biomarkers were reported to predict the response to immunotherapy, including tumor mutation burden and expression of immune checkpoints, such as PD-L1. We examined the value of the IPI-IPM to predict the response of patients to ICB therapy based on a publicly accessible dataset, the IMvigor210 Cohort (20). Using an optimal cutoff point calculated through the “surv_cutpoint” function (a minimal proportion of observations per group was set to 30%) of the survminer R package for group assignment, we observed that patients with high IPI-IPM risk scores had shorter OS after anti-PD-L1 treatment (log-rank p = 0.029, Figure 11C ).

Data acquisition of the present study is fully under the TCGA publication guidelines (https://www.cancer.gov/tcga). Gene expression data (RNA-seq) of 570 samples, masked somatic mutation data of 37 samples, and clinical follow-up data with clinicopathological characteristics of 566 patients of DLBCL projects (TCGA-DLBC, CTSP-DLBCL1, NCICCR-DLBCL) were obtained from the Cancer Genome Atlas (TCGA) by using the TCGAbiolinks (46) R package in R software (version 4.0.2, https://www.r-project.org). Clinical information, gene expression subtype, and genetic subtype of the DLBCL patients were supplemented by referring to open-access supplementary files of the GDC DLBCL publication (4). Matched gene expression data and survival follow-up data can be obtained in 563 samples. Matched gene expression data and survival follow-up data with available IPI data can be obtained in 445 samples ( Supplementary Table 1 ). Gene expression data (Microarray) with matched clinical information of 414, 221, and 928 DLBCL patients were obtained from GSE10846, GSE87371, and GSE117556, respectively (12, 47–49), by accessing the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo). The immunologic gene lists used in the present study were downloaded from the ImmPort (www.immport.org/shared/home) and ImmuneSigDB (50) through MSigDB Collections (www.gsea-msigdb.org/gsea/msigdb). For the RNA-Seq data, HTSeq-count data were downloaded. The Combat-Seq function of the Sva (51) R package was used to remove batch effects among different projects. The principal component analysis was performed and visualized to examine the batch effect (52). For each gene, the effective gene length was extracted by using the EDASeq (53) R package and TPM (Transcripts Per Kilobase Million) gene expression data were calculated through the count2TPM function utilized in the IOBR (54) R package by using the corresponding effective gene length. A Homo sapiens GRCh38 annotation file (Ensembl 103) downloaded from Ensembl (55) was used for gene symbol annotation. The DESeq2 (56) R package was applied for the normalization of RNA-seq count data, and variance stabilizing transformation (VST) data were used for downstream analysis. DLBCL microenvironmental (LME) signatures and LME subtype categorizing method were referred to the publication of Kotlov et al. on Cancer Discovery (11).

Samples with available IPI information were categorized to high-, intermediate-, and low-risk groups according to criteria of IPI, and the DESeq2 (56) R package was applied to identify differentially expressed genes (DEGs) between high- and low-risk groups. The differential expression was defined with a fold change of threshold at 1.5 and a false discovery rate (FDR) value < 0.05.

The clusterProfiler (57) R package was used for both overrepresentation analysis and pre-ranked gene set enrichment analysis (GSEA). Analysis of Gene Ontology (GO) (58), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (59), and Reactome pathway (60) were contained in the present study. An adjusted p value<0.05 was considered statistically significance. The threshold for GSEA was set at the p-value < 0.05, FDR < 0.05, and | normalized enrichment score (NES) | > 1.0. The non-parametric gene set variation analysis was further performed with the GSVA (61) package of R.

Weighted gene co-expression network analysis (WGCNA) is commonly used for analyzing high-throughput gene expression data with different characteristics, so as to mine gene co-expression networks and intramodular hub genes based on pairwise correlations in genomic applications. In the present study, we applied the WGCNA (62, 63) R package to analyze key gene clusters that were most relevant to IPI risk groups in DLBCL samples.

IPI-based immune-related genes were selected to construct the prognostic risk model. The training cohort (563 patients with matched normalized RNA-seq data and survival data from TCGA) was used for the construction of IPI-IPM, and three testing cohorts (GSE10846, GES87341, and GSE117556) were used for validation of the prognostic risk model. The Survival R package was used to analyze the correlation between the expression of objective gene sets and DLBCL patients’ overall survival (OS). Univariate Cox regression analysis was performed to screen genes, of which the expression was associated with OS with a p value < 0.05. Lasso (least absolute shrinkage and selection operator) regression analysis was applied for variable selection and regularization to enhance the prediction accuracy and interpretability by using the glmnet (64) R package. Multivariate Cox regression analysis was then carried out to select the optimal genes, based on the method of the Akaike information criterion (AIC) (65). For each sample, the risk score equals the sum of the normalized expression of each gene multiplying its corresponding regression coefficient. Time-dependent ROC curves were plotted by using the survivalROC (66) R package. Five hundred sixty-three DLBCL patients in the training cohort were divided to low- and high-risk score groups according to the optimal cutoff value with largest AUC in the receiver operating characteristic (ROC) curve of the median survival time. Then, Kaplan–Meier survival analysis and time-dependent ROC curve analysis were performed to evaluate the prognostic significance and accuracy of IPI-IPM (67). Besides, Harrell’s concordance index (C-index) was calculated by using the survcomp (68) R package. Univariate and multivariate Cox regression analyses were performed on the risk score and all available clinicopathologic parameters, including age, gender, Ann Arbor clinical stage, LDH ratio, ECOG performance status, and number of extranodal site. Then, we utilized the rms (69) R packages which set up a prognostic nomogram for OS probability assessment by enrolling all the independent prognostic factors. The discriminative efficacy, consistency, and clinical judgment utility of the nomogram score was evaluated by time-dependent calibration plots and decision curve analysis (DCA) (70) using the rms and rmda (71) R package.

DEGs between high- and low-risk score groups were analyzed following the fold change of a threshold at 1.5 and FDR value < 0.05. The gene expression of samples between IPI-IPM subgroups were analyzed with the t-distributed stochastic neighbor embedding (t-SNE) method by using the Rtsne (72) package of R and then visualized on the 3D map with the scatterplot3d (73) package of R. Somatic mutations of IPI-IPM subgroups were analyzed by using the Maftools (74) R package. The intersection of DEGs and immune-related gene set was used to construct a protein–protein interaction (PPI) network based on the STRING (75) database. Cytoscape (76), plugin MCODE (77), and cytoHubba (78) were utilized to identify top 10 degree genes in the network. The TRRUST (79) database was browsed to explore the curated transcriptional regulatory networks. Gene expression data (TPM) of 570 DLBCL samples were imported into CIBERSORT (80), MCPCounter (81), and xCell (82) to calculate the score to estimate the proportion of TME cells including immune and stromal cells.

An unsupervised clustering algorithm was applied to analyze the gene expression pattern of LME signatures in 563 DLBCL samples. By using the ConsensusClusterplus (83) R package, we performed the k-means clustering algorithm with 1,000 repetitions to ensure the stability. The Pearson distance matrix calculated from the clustering was then imported to Cytoscape for the visualization of distribution of samples corresponding to LME signature clustering.

All statistical data were analyzed in the R software (version 4.0.3). A Wilcoxon test was applied to compare continuous variables between two groups of sample data. A Kruskal–Wallis test was applied to compare continuous variables among three or more groups of sample data. A test is considered with statistical significance at two-sided p < 0.05. Pearson correlation was used to test the correlation between two sets of continuous data, and an absolute Pearson correlation coefficient larger than 0.3 was considered to be correlated. A Pearson correlation is considered with statistical significance at FDR < 0.05. Spearman correlation was used to test the correlation between gene expression and IPI groups. We used ggplot2, ggstatsplot, and ggpubr R packages (84, 85) for data analysis and visualization.