The public gene expression databases were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database and The Cancer Genome Atlas (TCGA) database (Project name: Lymphoid Neoplasm Diffuse Large B-cell Lymphoma, https://portal.gdc.cancer.gov/projects/TCGA-DLBC). We included GSE31312, GSE10846, GSE12195, GSE56315, and TCGA-DLBC in subsequent analyses. Those four GEO datasets were derived from the same platform (GPL570 Affymetrix Human Genome U133 Plus 2.0). We downloaded raw CEL files and adopted a robust multiarray average (RMA) method to normalize the data. GSE12195 and GSE56315 were used to determine the expression patterns of the pyroptosis-related genes (PRGs) between normal B cells and DLBCL. The results were demonstrated using a box diagram.

The GSE31312 dataset contains abundant clinical information, including the age, gender, survival information, Gene Expression Profiling (GEP), treatment regimen (RCHOP) and response, and the International Prognostic Index (IPI) as the main analysis objects. The GSE10846 dataset, which was also treated using the RCHOP regimen, was used as a validation set. We further processed GSE31312 and GSE10846 as follows: 1) patient sample data with an overall survival time of <30 days were excluded; and 2) for the GSE31312 dataset, samples missing the clinical information mentioned above were excluded.

We obtained information on the somatic mutation and copy number variation (CNV) of DLBCL from the TCGA-DLBC dataset. R package “maftools” was used to analyze the mutation frequency of PRGs and visualize the results using an oncoplot waterfall plot (Mayakonda et al., 2018). We calculated the CNV frequency of PRGs and displayed the results with a lollipop chart. The “RCircos” package of the R software package was used to visualize the location of those genes on chromosomes (Zhang et al., 2013).

We obtained 52 PRGs from “REACTOME PYROPTOSIS”, “GOBP PYROPTOSIS”, and previous studies (Man and Kanneganti, 2015; Wang and Yin, 2017; Karki and Kanneganti, 2019; Xia et al., 2019; Ye et al., 2021). “REACTOME PYROPTOSIS” and “GOBP PYROPTOSIS” were derived from MSigDB (Molecular Signatures Database v7.4). The R software package “ConsensusClusterPlus” was used to screen the distinct clusters based on the expression of the 52 PRGs (Wilkerson and Hayes, 2010). Resampling was performed 100 times to ensure classification reliability. We then investigated the relationship between clusters and clinical features.

To determine whether or not there were differences in biological processes among different clusters, we downloaded “c2.cp.kegg.v7.4.symbols.gmt” and “h.all.v7.4.symbols.gmt” from MSigDB. The pathways with adjusted p values < 0.05 were considered to have significant differences. The above analysis process was performed using the R software package “GSVA” (Hänzelmann et al., 2013). Finally, the significant pathways were presented as a heatmap.

We used the R software package “ESTIMATE” (Yoshihara et al., 2013) to assess the “immunescore” and “stromalscore” in DLBCL. Sngle-sample gene-set enrichment analysis (ssGSEA) was used to quantify estimation of infiltration abundance of different immune cells (Supplementary file). To avoid unnecessary interference, we excluded B cells and related immune cells. Through the ssGSEA algorithm, we obtained the score of each immune cell in each patient and standardized the findings by the Min-Max method. We collected other TME-related signatures through previous studies, including the 1) CD8 T-effector signature (Rosenberg et al., 2016); 2) pan-fibroblast TGFβ response signature (pan-F-TBRS) (Mariathasan et al., 2018); 3) antigen-presenting machinery (APM) signature (Şenbabaoğlu et al., 2016); 4) angiogenesis signature (Şenbabaoğlu et al., 2016); 5) three epithelial mesenchymal transition) (EMT) signatures (EMT1, EMT2, and EMT3) (Damrauer et al., 2014; Hedegaard et al., 2016; Hugo et al., 2016), and 6) two stromal gene signatures (stromal 1 and stromal 2). We then adopted the same process with those signatures (Lenz et al., 2008).

As the training dataset, the GSE31312 dataset was used to screen for PRGs associated with the prognosis of DLBCL patients. The expression of the PRGs was first determined using a univariate Cox analysis (Cox, 1972). Genes with p values < 0.01 were retained. We used a Lasso-penalized Cox regression analysis and Stepwise regression analysis to further screen for PRGs with the best predictive performance and build the risk score signature. The multivariate regression model was then internally validated using the bootstrapping method. After obtaining the regression coefficients for each enrolled gene, we calculated the risk score for each patient based on the expression of each gene using the following formula:

The risk score = (Exp _gene1 × coefficient _gene1) + (Exp _gene2 × coefficient _gene2) +···+ (Exp _gene7 × coefficient _gene7)

Patients were then divided into high and low-risk groups according to the median risk score. Differences in the overall survival (OS) of patients in different risk groups were assessed using Kaplan-Meier curves (Ranstam and Cook, 2017) and AUC of the ROC curve (Kamarudin et al., 2017). The regression coefficients obtained from the training dataset were then applied to the GSE10846 testing dataset with complete clinical information to calculate the patient’s risk score for external validation.

To determine whether or not the gene signature has independent prognostic value and to obtain clinically independent prognostic factors associated with the DLBCL prognosis, we performed univariate and multivariate Cox regression analyses on prognostic gene signature and clinicopathological parameters, including the age, gender, GEP, treatment regimen (RCHOP) and response, and IPI, in the training dataset. Factors with p < 0.05 in the univariate analysis were selected for a further multivariate Cox regression analysis. The above-incorporated independent parameters were used to construct a prognostic nomogram by stepwise Cox regression to predict the OS in patients with DLBCL at 2, 4, 6, and 8 years. The AUC-ROC curve, Harrell’s concordance index (Harrell et al., 1982) and a calibration plot comparing the predicted and observed OS were used to evaluate the performance of the nomogram.

An alluvial diagram was developed showing the changes in pyroptosis-related clusters, pyroptosis risk score, IPI score, and treatment response using the R software package “ggalluvial” (Brunson and Read, 2020) and “ggplot2” (Villanueva and Chen, 2019). We performed the R function “cor” to calculate the correlation coefficients between the TME (including infiltrating immune cells and other TME signatures), mRNA expression of IL-1β and IL-18 and pyroptosis risk score, and the correlation matrix was drawn using the “corrplot” package (Wei et al., 2017).

We downloaded the target databases from TCGA and GEO databases and performed preliminary processing of the data, after which we explored the landscape of the genetic and expression variation of PRGs in DLBCL (Figure 1A). We then performed consensus clustering according to the expression of 52 PRGs and compared the differences in the survival, biological processes and TME among different clusters (Figure 1B). Subsequently, we constructed a prognostic signature based on the gene expression and survival information of the 421 patients in the training set. The time-ROC curve was used to evaluate the stability, and external verification was performed to prove the reliability of the signature. We then integrated the clinical features, constructed a predictive nomogram and evaluated the model with calibration curves (Figures 1C,D).

First, we summarized the incidence of CNV and somatic mutations of the 52 PRGs in DLBCL based on the TCGA-DLBC dataset. Twelve of 37 patient samples had somatic mutations, the most common form being missense mutations. At the gene level, a total of 12 genes were mutated, with TP53 showing the highest frequency (up to 14%), followed by NLRP9; the remaining 10 genes, including ZBP1 and TNF, showed roughly 3% frequency (Figure 2A). To our surprise, GSDMD did not show any mutations in DLBCL samples. By estimating the frequency of CNV, we found that PRGs showed prevalent CNV alterations. However, PLCG1, NAIP, GZMB, ELANE, ELANE, and APIP were neither amplified nor deleted (Figure 2C).

The locations of the CNV alteration of PRGs on chromosomes are shown in Figure 2B. To determine whether or not CNV alteration of PRGs affects the mRNA expression and to ensure the stability of the results, we verified these findings in the GSE12195 and GSE56315 datasets at the same time. We first excluded genes with inconsistent expression trends in the two GEO datasets and then presented the remaining genes in the form of boxplots (Figures 2D,E). Most genes were upregulated in DLBCL, and we also found that among the upregulated PRGs, PYCARD, IRF1, GZMA, GSDMD, GSDMC, GPX4, CHMP2A, CASP5, CASP1, IL-18, and BAK1 showed an amplified CNV status, suggesting a significant relationship between the PRG expression and CNV. After the discovery that the expression of a large proportion of PRGs was imbalanced in normal B cells and DLBCL samples, we analyzed the association of these genes with the prognosis by a univariate Cox regression analysis. As shown in Supplementary Figure S1A, there was a significant correlation between the PRGs and prognosis. These findings suggested that genetic and expression variation of PRGs is common in DLBCL, so pyroptosis may play an important role in the occurrence and development of DLBCL.

We clustered GSE31312 based on the expression of 52 PRGs quantities by the R package “ConsensusClusterPlus” and ultimately identified 3 completely different clusters, including 98 samples in cluster A, 229 samples in cluster B, and 94 samples in cluster C (Supplementary Figures S1B,C). We subsequently performed a survival analysis of the three clusters and found that they had distinct prognoses, with cluster A having the best prognosis and cluster B having a more favorable long-term prognosis than cluster C (Figure 3A). The ssGSEA method was used to establish a pyroptosis signature score for the three clusters, with results showing that cluster B had the highest score, while there was no significant difference between the scores of A and C (Figure 3B).

Given that the three clusters had different prognoses under the same treatment regimen, we speculated that the clusters were biologically functionally different, and the “GSVA” enrichment analysis validated our hypothesis. Compared with cluster B, cluster A showed enrichment of “KEGG HEDGEHOG SIGNALING”, “HALLMARK KRAS SIGNALING Down” and “KEGG CALCIUM SIGNALING”, among others. For cluster B, a large number of pathways associated with nutrient metabolism were enriched, including “HALLMARK ADIPOGENESIS”, “HALLMARK FATTY ACID METABOLISM”, “KEGG PENTOSE PHOSPHATE”, “HALLMARK MTORC1 SIGNALING”, “HALLMARK/KEGG UNFOLDED PROTEIN RESPONSE”, and “KEGG/HALLMARK OXIDATIVE PHOSPHORYLATION”. In addition, cluster B also showed significant enrichment of the DNA damage response (DDR) pathway and its related pathways, such as “KEGG MISMATCH REPAIR”, “HALLMARK DNA REPAIR”, “KEGG NUCLEOTIDE EXCISION REPAIR”, and “HALLMARK G2M_CHECKPOINT”. The enrichment of cancer-related pathways, such as “HALLMARK TGF BETA SIGNALING” and “HALLMARK MYC TARGETS” (Figure 3C), in cluster B was similarly observed in the cohort of B versus C. In contrast, cluster C showed enrichment of the “KEGG ECM RECEPTOR INTERACTION” pathway and downregulation of the “KEGG B CELL RECEPTOR SIGNALING” pathway (Figure 3D). This suggests that not only the biological processes but also the immune microenvironment may differ among the clusters.

Next, we analyzed the clinical features of the clusters. Unexpectedly, cluster C, which had the worst survival, showed a higher proportion of the GCB type, whereas cluster B, which was intermediate between A and C, possessed the highest proportion of the non-GCB type (Figure 4A). Regarding the IPI score, cluster A, which had the best survival, showed the highest proportion in the low-risk group (IPI score: 0–1). Compared with cluster B, cluster C accounts for more in low IPI scores (Figure 4B). Given that the dataset GSE31312 as the RCHOP treatment cohort, we then explored the response of three clusters to the treatment regimen. These findings seemed to explain the difference in the survival among the three clusters, such cluster A showing the highest complete remission (CR) rate, cluster B the next highest, and cluster C the worst response to RCHOP treatment (Figure 4C).

We then verified several previously reported genes associated with resistance to the RCHOP regimen and found no significant difference in the expression of “ADAM12” or “ABCB1” among the three clusters, but “ABCG2” was significantly overexpressed in cluster C (Figure 4D) (Ohsawa et al., 2005; Kim et al., 2009; Yin et al., 2017).

GSVA enrichment analysis suggested that there were differences in the immune microenvironment status among the three clusters. Through an ssGSEA analysis, we found that cluster B had strong adaptive immune activation, including CD8⁺ and CD4⁺ T cells, and it had the highest Treg cell infiltration among the three clusters. Regarding innate immune cells, such as NK T cells, NK cells, and MDSC cells, there were no significant differences among the three clusters (Figure 5A).

We then used the “ESTIMATE” package to score the immunity and stromal of the three clusters. As shown in Figure 5B, cluster B had the highest immune score, which is consistent with the results of the ssGSEA, proving that the TME of cluster B had greater immune cell infiltration than that of the other clusters. Incidentally, we also evaluated the differences antigen-presenting machinery (APM) among the three clusters, which has been shown to be significantly associated with the T cell infiltration score (TIS) (Şenbabaoğlu et al., 2016); the results remained consistent with the former, with cluster B being significantly higher than clusters A and C, while no significant differences were observed between the remaining two clusters (Figure 5C). However, patients in cluster B did not derive a consequential survival benefit. Tumors with the immune-excluded phenotype show infiltration of a large number of immune cells, but due to stromal entrapment, these immune cells are only distributed close to tumor cells and are unable to penetrate the stroma to exert their effect. We therefore also determined the matrix score of the three clusters. A significant difference was noted, which was consistent with the prognosis trend. Cluster A’s stromal scores were smaller than those of clusters B and C, whereas no significant differences were observed between clusters B and C (Figure 5D).

One of the characteristics of the immune-excluded phenotype TME is active angiogenesis (Hegde et al., 2016), so we evaluated the angiogenesis pathways using the ssGSEA algorithm and compared the findings among the three clusters. We found that cluster B had a significantly higher angiogenesis score than cluster A but no significant difference from cluster C was noted (Figure 5E). We also found that clusters B and C had higher levels of hypoxia than cluster A (Figure 5F). IL-1β and IL-18 are the two most important inflammatory mediators of pyroptosis, and we found that IL-1β was highly expressed in cluster C, while IL-18 was highly expressed in cluster B (Figure 5G). Therefore, we speculate that the immune cell activities in cluster B were suppressed by inhibitory cellular, molecular, and other components of the TME and thus were unable to exert their pro-survival effects.

Based on the above findings, we found that the three pyroptosis-related clusters have significantly different TME situations. Cluster B, with pyroptosis-related inflammation dominated by IL-18, closely resembles the immune-excluded phenotype, while cluster C is more consistent with an immune-desert phenotype, with pyroptosis-related inflammation dominated by IL-1β. Cluster A seems characterized by a low pyroptosis-related inflammatory phenotype.

After filtering and sorting the data, a total of 421 patients in the GSE31312 dataset met the criteria for inclusion in the survival analysis. A univariate Cox regression analysis showed that a total of 18 genes were significantly associated with the prognosis of DLBCL patients (p < 0.05). The results of the univariate analysis of these genes are shown in Table 1. A prognostic signature comprising seven genes—PRKACA, PLCG1, NLRP9, CASP6, CASP4, BAK1, and AIM2—was developed by a Lasso-penalized Cox analysis and stepwise Cox analysis. Based on the hazard ratio (HR), the downregulated PLCG1, CASP6, CASP4, and AIM2 were considered tumor suppressors, whereas the upregulated PRKACA, NLRP9, and BAK1 were regarded as oncogenes (Figure 6A).

We used the following formula to calculate the risk score:

[(Exp PRKACA × (0.70221) + (Exp PLCG1 × (−0.81499) + Exp NLRP9 × (1.21405) + Exp CASP6× (−0.88874) + Exp CASP4 × (−0.61049) + Exp BAK1 × (0.67189) + Exp AIM2 × (−0.21797)]

The median risk score was then used as the critical value for risk grouping. Patients from the training dataset were stratified into two groups. The Kaplan-Meier survival curves showed a significantly better OS in the low-risk group than in the high-risk group (p < 0.0001) (Figures 6B–D). The time-dependent ROC curve was used to determine the prognostic value of the signature. The AUCs for 2-, 4-, 6-, and 8-years OS predictions were 0.71, 0.71, 0.75, and 0.78, respectively, which were all higher than 0.7, proving that the risk score signature had favorable prognostic ability in the training dataset (Figure 6E).

Internal validation of the signature using the bootstrapping method showed that it was reasonably stable (Supplementary Figures S2A–D). The predictive performance of the seven-gene signature was externally validated in the GSE10846 dataset. Kaplan Meier curves showed that the OS remained significantly different between risk groups in the validation dataset (Figure 6F). The prognosis of the high-risk group was significantly worse than that of the low-risk group (p < 0.05). External validation showed that the seven-gene signature performed well in predicting the OS in patients with DLBCL.

We combined the risk status with the age, gender, GEP, treatment regimen (RCHOP) and response, and IPI score as candidate predictors and created a prognostic nomogram predicting the OS of DLBCL patients based on the LASSO regression analysis and stepwise Cox regression analysis. The risk status, age, GEP, response, and IPI score were considered to be the final parameters in the nomogram (Figures 7A,B). The AUCs of the OS for 2, 4, 6, and 8 years were 0.84, 0.86, 0.86, and 0.81, respectively, (Figure 7C). The C-index of the nomogram was 0.833 (95% confidence interval [CI]: 0.816–0.85), and the calibration curves revealed that the predicted the OS accorded with the observed OS (Figure 7D). Taken together, these results show that the nomogram has an excellent ability to predict the survival time in DLBCL patients.

The pyroptosis risk score gradually increased among the three clusters of A, B, and C, which was consistent with the difference in the prognosis between the clusters (Figure 8A). The low-risk group had a higher pyroptosis signature score than the high-risk group (Figure 8B). We further analyzed the relationship between the pyroptosis risk group and IPI score, response to treatment and PRG-related clusters and visualized them using an alluvial diagram (Figure 8C). It was evident that the low-risk group had lower IPI scores and better treatment responses than the high-risk group. The above analyses suggested that the low-risk group was significantly associated with more favorable clinical features for the survival than the high-risk group.

To evaluate the correlation between the TME and pyroptosis in DLBCL, we included the pyroptosis risk score, pyroptosis signature score and TME components, including immune cells, APM, CD8⁺ T cell effort, and stroma-related pathways, in the correlation matrix for the comprehensive analysis. Infiltrating immune cells, such as activated CD8⁺ T cells and NK cells, showed a negative correlation with the pyroptosis risk score, while stromal pathways, such as pan-F-TBRS, EMT1-3 and stromal 2, were not correlated (Figure 8D). Notably, the prognostically favorable stromal 1 signature was also negatively correlated with the pyroptosis risk score. The pyroptosis signature score showed an extremely strong correlation with the TME, associated with either adaptive or innate immunity, and stromal components.