Introduction

Non-Hodgkin lymphoma (NHL) ranks as the 10th most diagnosed and the 11th leading cause of cancer death in 20221. Diffuse large B-cell lymphoma (DLBCL), the most common subtype of NHL in adults worldwide, is a biologically and clinically heterogeneous disease. Although standard immunochemotherapy with R-CHOP can cure DLBCL in 60% of patients, up to 40% of patients remain refractory to treatment or relapse after remission, leading to significant morbidity and mortality2. Deficiency in tumor infiltrating lymphocyte (TIL) plays an important role in mediating drug resistance by altering the lymphoma microenvironment (LME) and offering protection from apoptosis3. DLBCL patients with higher TIL levels have a superior complete response rate after R-CHOP therapy and a better prognosis in overall survival (OS)4. Intriguingly, the degree of TIL is mainly determined by the immunogenicity of the DLBCL cells, which is partially induced by inflammatory changes in the LME5. Hence, identification of a molecular signature based on LME may effectively forecast the response of the immunochemotherapy.

Recent studies have identified a unique inflammatory cell death process named PANoptosis, which harbors pyroptosis, apoptosis and necroptosis, but cannot be accounted by any of them alone6. Intricate crosstalk occurs between above three programmed cell death (PCD) pathways during the executing of PANoptosis, even though it depends on a single inflammatory vesicle complex called PANoptosome7. Under some circumstances, interactive model of PANoptosis enables different cell death events to change from one mode to another8. This process is pivotal as PANoptosis triggering might conceivably provide compensatory inflammatory cell death events when apoptosis is hindered, thus overcoming drug resistance and priming anti-tumor immune responses in LME. Currently, growing evidence shows that prognostic models in compliance with PANoptosis have shown promise in other cancers9, yet there are no data assessing the PANoptosis-related genes on DLBCL. This gap underscores the need for a deeper understanding of immunogenic PANoptosis in DLBCL, which could provide potential opportunities to generate promising insights for DLBCL immunochemotherapy.

This study aims to fill this gap by investigating gene expression and immune landscape of DLBCL samples from the public databases. Robust analysis was carried out to identify critical PANoptosis-related genes, including nonnegative matrix factorization (NMF) clustering, least absolute shrinkage and selection operator (LASSO) Cox regression, multivariate Cox regression10. Using these insights, we constructed a PANoptosis-related gene prognostic index (PANGPI) to predict the prognosis and TIL landscape in DLBCL patients. We further analyzed the relationship of PANGPI with the immunotherapy response and related sensitive drugs11. Additionally, in silico docking and molecular dynamic (MD) simulation studies were presented to elucidate the mode of interaction of these predicted drugs12. Finally, the PANGPI showed satisfactory prediction efficiency in terms of prognostic risk stratifcation, TIL evaluation, and therapeutic response prediction for DLBCL.

Methods

Public data acquisition and preprocessing

RNA sequencing gene expression profile and corresponding clinical information of 48 DLBCL samples were obtained from the Cancer Genome Atlas (TCGA) (TCGA-DLBC cohort) (https://tcga-data.nci.nih.gov/tcga/). Clinicopathological information and microarray expression profiling data of 33 normal tissues and 807 DLBCL specimens in 4 DLBCL cohorts (GSE56315, GSE87371, GSE10846, GSE53786) were retrieved via the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo). A total of 34 PANoptosis-related genes (PANRGs) were acquired from the previous literature13, and the genes are listed in Supplementary Table S1.

Differentially expressed PANRGs and protein–protein interaction network

Differentially expressed PANRGs were identified by the “limma” R package14 and selected according to |fold change (FC) > 1.5 and false discovery rate (FDR) < 0.05. The protein–protein interactions (PPI) among PANRGs were constructed using the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/) with minimum required interaction score ≥ 0.615. Cytoscape software (version 3.7.2) was used for visualization of the PPI network16.

Unsupervised clustering of PANRGs

Unsupervised NMF clustering was applied to separate DLBCL patients into distinct molecular subtypes based on the PANRGs via “NMF” R package (https://cran.r-project.org/package=NMF,version 0.26).

Construction of prognostic risk model based on PANRGs

The univariate Cox regression was performed to identify prognosis-related PANRGs, and LASSO Cox regression was performed to construct the prognostic risk model using the “glmnet” package (version 4.1.7) in R17. We used the following formula to calculate the risk score for each DLBCL patient in the cohort:

$$ {\text{Risk}}\;{\text{score}} = \sum\nolimits_{{{\text{i}} = 1}}^{{\text{n}}} {{\text{Coeffficient}}_{{\text{i}}} \times {\text{expression}}_{{\text{i}}} } $$

In this formula, “Coefficient” is the coefficient from LASSO-Cox regression, “Expression” is the expression of PANRGs and “n” is the number of genes we included in the prognostic model.

Survival analysis

R packages “survival” (https://CRAN.R-project.org/package=survival, version 3.5–5) and “survminer” (https://CRAN.R-project.org/package=survminer, version 0.4.9) were used to analyze DLBCL patient survival and prognosis in high and low risk groups. The survival curve was plotted using the Kaplan–Meier method, and the log-rank test was used to assess statistical significance. The multivariate Cox proportional hazards regression analysis was carried out to determine the PANGPI as an independent prognostic factor in DLBCL.

Nomogram plot

For prediction of survival and recurrence at 1, 3 and 5 years in individual DLBCL patients, the “rms” R package (https://CRAN.R-project.org/package=rms) was utilized to establish a nomogram plot. To evaluate the fit between the actual survival and nomogram-predicted survival, the calibration curve was generated.

Somatic mutation and functional annotation

The “maftools” R package was utilized to evaluate and visualize the frequency and somatic interactions of the differentially mutated genes18. The Gene Set Enrichment Analysis (GSEA) was used to analyze the association between the risk score and hallmarks. The GSEA was applied to recognize the most important feature between the PANGPI high and low risk groups. The genetic subtypes of DLBCL patients from TCGA-DLBC cohort were classified by LymphGen 2.0 online website (https://llmpp.nih.gov/lymphgen/index.php)19.

Immune infiltration analysis

The tumor immune microenvironment of patients with DLBCL was quantified using single-sample GSEA (ssGSEA) and xCell algorithms20, (https://github.com/dviraran/xCell). The tumor immune dysfunction and exclusion (TIDE) scores (http://tide.dfci.harvard.edu/) was used to predict the correlation between PANGPI and immune checkpoint inhibitor responses21. The “easier” package was used to calculate cytolytic activity (CYT) scores22 and T cell inflamed gene expression profile (GEP) signature23.

Drug sensitivity prediction analysis

The Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) was utilized to locate the drugs with the specific target potential as anti-cancer candidates based on the correlation between the PANGPI and drug response24. Anti-tumor drug response in the PANGPI high and low risk groups was predicted by half maximal inhibitory concentration (IC50) by “oncoPredict” R package25.

Molecular docking analysis

The protein structures of 9 PANoptosis core targets26 were downloaded from PDB (Protein Database) (https://www.rcsb.org). The chemical structures of ligands in 2D were drawn by using Chemdraw software and were converted into 3D. Molecular docking of predicated drugs with the core targets was performed using the AutoDock Vina 1.1.2, and visualized by Pymol and Discovery Studio. If the binding free energy was less than − 5 kcal mol−1, it indicated that the two had good binding affinity27.

Molecular dynamics simulation

Molecular dynamics simulation was conducted for the optimal core protein ligand obtained by molecular docking by Gromacs 2020.6 and Amber 20 software28. The complex was solvated in TIP3P water model with the proper number of Na+ and Cl- counterions to mimic a physiological ionic concentration. A leap frog algorithm was used to equilibrate NVT and NPT for 100 ps steps29.

Ethical approval

TCGA, GEO and GDSC datasets used human genomic data deposited in public repositories, so ethics approval is not applicable to these datasets.

Statistical analysis

Wilcoxon test was employed to calculate the statistical difference in two groups. Pearson correlation analysis was performed to imply the relationship between the different variates. P < 0.05 was considered as statistical different in this study.

Results

Identification of differentially expressed PANRGs

Differential gene expression analysis comparing tumor samples to normal counterparts within the GEO cohort (GSE56315) was conducted. A total of 9302 DEGs were identified in DLBCL samples compared with normal counterparts, including 3194 upregulated genes and 6108 downregulated genes (Fig. 1A,B). After overlapping, 10 DEGs (ZBP1, RIPK1, MLKL, TLR2, TLR3, TLR4, RAVER1, CASP1, MEFV, TUFM) were identified between 9302 DEGs and 34 PANRGs (Supplementary Table S2). Among them, MEFV, ZBP1, RAVER1, RIPK1 were downregulated in DLBCL samples, while TLR3, MLKL, CASP1, TUFM, TLR4, TLR2 were upregulated (Fig. 1C). GO analysis indicated that these genes were mainly involved in the positive regulation of inflammatory response (GO:0,050,729) (Fig. 1D) (Supplementary Table S3).

Fig. 1
figure 1

Volcano plot (A) and heatmap (B) of the 9302 DEGs in DLBCL cohorts of the GEO database (GSE56315). (C) Expression levels of 10 differently expressed PANRGs in tumor and normal tissues. (D) Go analysis of the PANRGs. **p < 0.01 and ***p < 0.001. Heatmap was generated with R software (v. 4.3.1; https://cran.r-project.org/web/packages/pheatmap/index.html).

Full size image

Establishment of molecular subtype based on the PANRGs

Unsupervised NMF consensus clustering was conducted on the GSE87371 cohort on the basis of 10 differently expressed PANRGs. Following the prompts of silhouette, consensus, as well as cophenetic(Fig. 2A), 221 samples were eventually split into 2 clusters (Fig. 2B). Survival analysis also revealed that there was an obvious difference in the OS of DLBCL patients between cluster C1 and C2 (Fig. 2C). The PPI network of 10 differentially expressed PANRGs was fashioned in STRING and MLKL, TLR2, TLR3, TLR4, CASP1, MEFV, ZBP1 and RIPK1 were determined as 8 hub genes of network (Fig. 2D). The differential expression patterns of 8 PANRGs between two clusters were visualized in heatmap, among them, MLKL, TLR4, CASP1, TLR3, TLR2 were significantly upregulated in cluster C1 (Fig. 2E).

Fig. 2
figure 2

Unsupervised clustering analysis based on PANRGs expression spectrum. (A) NMF rank survey with optimal number of clusters: rank = 2. (B) Consensusmap defining two different subtypes of PANRGs expression patterns. (C) Kaplan–Meier curves of overall survival for the two PANRGs-related clusters. PPI network (D) and heat map (E) of 8 differently expressed PANRGs. Heatmap was generated with R software (v. 4.3.1; https://cran.r-project.org/web/packages/pheatmap/index.html).

Full size image

Development of a PANRGs‑related prognostic risk signature

To clarify the prognostic value of the 8 hub genes, we applied univariate Cox regression analysis using OS data from the GEO dataset (GSE10846). We found 3 genes were still associated with the prognosis of patients with DLBCL, including MLKL, TLR4, MEFV (P < 0.05) (Fig. 3A). Notably, the Kaplan–Meier analysis indicated that MLKL or MEFV high expression patients showed poorer survival outcomes, while TLR4 high expression patients exhibited higher survival probability (Fig. 3B–D). Mantel tests indicate that there was a significant correlation between MLKL and MEFV (Fig. 3E).

Fig. 3
figure 3

Development of a PANRGs-related prognostic risk signature. (A) The results of univariate Cox analysis in the training dataset. (B–D) Kaplan–Meier analysis for patients with differentially expressed prognostic gene (MLKL, TLR4 and MEFV). (E) Mantel test analysis of MLKL, TLR4 and MEFV. Risk score curves, risk state dot plots and prognostic gene expression heat maps at the training (F) and validation (G) datasets. Heatmap was generated with R software (v. 4.3.1; https://cran.r-project.org/web/packages/pheatmap/index.html).

Full size image

According to the λ-value of every gene discovered based on LASSO Cox regression, a PANRGs‑related prognostic risk model was then constructed, and can be expressed as follows: PANGPI = 0.3110775* express value of MLKL + (− 0.2826979)* express value of TLR4 + 0.1107089* express value of MEFV. Then DLBCL patients were classified as the high- and low-risk groups according to the median value of risk scores in the training (GSE10846) and validation set (GSE53786). Consequently, the dead samples were mainly centralized in high-risk group with a low expression of TLR4, and a high expression of MLKL and MEFV (Fig. 3F,G).

Building a predictive nomogram

The results of multivariate Cox analysis of preoperative clinical parameters showed that age, stage, and PANGPI were independent prognostic factors. Among them, PANGPI was associated with unfavorable prognosis (HR = 2.87, 95%CI 1.94–4.26, P < 0.001) (Fig. 4A). The PANGPI high-risk group was associated with significantly reduced OS in both training and validation set (Fig. 4B,C). These 3 independent prognostic factors were employed to establish a nomogram model to predict the future probability of a specific outcome event of DLBCL patients. Each variable corresponds to a specific point, through drawing a straight line on the point axis, and then the sum of the values corresponding to these points can predict the probability of OS in 1, 3 and 5 years (Fig. 4D). The calibration plots were drawn to evaluate the consistency of the nomogram. In this study, the calibration graphs visually demonstrated excellent agreement between the predicted and actual outcomes at 1, 3, and 5 years (Fig. 4E).

Fig. 4
figure 4

Nomogram construction and prediction accuracy assessment. (A) Forest plot showing the multivariate Cox regression model of PANGPI and other clinical phenotypes. Kaplan–Meier analysis of PANGPI in training (B) and validation (C) datasets. (D) Nomogram construction based on the PANGPI prognostic signature and clinical characteristics. (E) Calibration curves comparing the accuracy of nomogram predicted OS with the actual OS.

Full size image

Molecular characteristics of patients with different PANGPI

We characterize the landscape of somatic mutations in the PANGPI high- and low-risk groups using data from TCGA-DLBC cohort (Supplementary Table S4). Although IGHV2-70, IGLL5 and IGLV3-1 were the most frequent mutations, the relative frequency varied between two groups (Fig. 5A,B). Differences in biological function were assessed by GESA, NOD-like receptor, Notch, p53, Toll-like receptor and VEGF signaling pathway differed significantly between two groups (Fig. 5C) (Supplementary Table S5). In addition, maftools analysis results showed that 4 genes mutated more frequently in DLBCL patients with high-PANGPI risk, including MUC16, CD79B, IRF4 and ZNF208, and there was extensive co-mutation between IRF4 and ZNF208, CD79B (Fig. 5D,E). Notably, genetic subtyping based on LymphGen 2.0 classifier show that there was a total of 5 MCD patients and 4 of them were distributed in the high-PANGPI-risk group (Supplementary Table S6). Moreover, many key inflammatory genes were differentially expressed between two groups (Fig. 5F) (Supplementary Table S7), among them MEFV was also confirmed upregulated in the high-PANGPI group.

Fig. 5
figure 5

Molecular characteristics and functional pathways in different PANoptosis patterns. (A, B) Visual summary showing common genetic alterations in the high and low-PANGPI groups. (C) GSEA plot showing the top 5 pathways in the high-risk group compared to the low-risk group. (D) Forest plot of gene mutations in the patients. (E) Interaction effect of genes mutating differentially in patients. (F) Heatmap displaying the inflammatory response genes expression between high and low-PANGPI groups. *p < 0.05. Heatmap was generated with R software (v. 4.3.1; https://cran.r-project.org/web/packages/pheatmap/index.html).

Full size image

Depiction of immune landscape on account of PANoptosis signature

Due to the significance of the LME in tumor development and outcome, we then conducted the immune infiltration estimations of DLBCL samples (GSE10846) using xCell and ssGSEA algorithms. We observed that diverse TILs were conspicuously distributed in high-risk DLBCL patients according to xCell algorithm, for instance, Th1 cells, CD4+ naive T cells and class switched memory B cells. In contrast, Th2 cells and CD8+ naive T cells were enriched in the low PANGPI patients (Fig. 6A). Moreover, results of ssGSEA algorithm revealed that infiltrating level of effector and central memory CD8+ T cells was positively correlated with PANGPI risk scores, effector memory CD4+ T cells and NK cells were negatively correlated with PANGPI risk scores (Fig. 6B).

Fig. 6
figure 6

Depiction of immune landscape on account of PANoptosis signature. Analysis of immune activity between the two risk groups using xCELL (A) and ssGSEA (B). Comparisons of the T cell-inflamed GEP (C), CYT (D), TIDE (E), T cell dysfunction (F), and T cell exclusion (G) scores between the high- and low-risk groups. (H) Comparison of immune checkpoint molecule expression between the high- and low-risk groups. **p < 0.01, ***p < 0.001 and ****p < 0.0001.

Full size image

The potential of immunotherapy response was subsequently assessed by quantifying the cytolytic activity and inflammatory states of the TILs. We found that the degree of CYT and T cell-inflamed GEP were higher in PANGPI high-risk group (Fig. 6C,D). However, the TIDE scores and T cell dysfunction scores were higher in PANGPI high-risk group, while T cell exclusion scores were higher in PANGPI low-risk group (Fig. 6E–G). Compared to the PANGPI low-risk group, DLBCL patients in PANGPI high-risk group exhibited a significant upregulation in the expression levels of immune checkpoint biomarkers, including CTLA4, CD274, LAG3, PDCD1 and CD80 (Fig. 6H).

Chemotherapeutic drug sensitivity analysis

To predict potential drugs for DLBCL patients in PANGPI high-risk group, we explored the correlation between the IC50 and PANGPI scores in tumor cell lines using the GDSC database (Fig. 7A,B) (Supplementary Table S8, Supplementary Table S9). After intersecting results in GDSC1 and GDSC2 datasets, we found that 7 drugs had lower IC50 values in the high-risk group, including PF-4708671, BMS-536924, Mirin, Gefitinib, Axitinib, Navitoclax and GSK269962Ab (Fig. 7C), indicating increased sensitivity. Notably, BMS-536924 (IGF-1R inhibitor), Gefitinib (EGFR inhibitor) and Navitoclax (Bcl-2 inhibitor) have an IC50 value lower than 10 µM (P < 0.01) (Fig. 7D,E). Among them, BMS-536924’activity was significantly correlated with the expression of PANGPI gene MEFV (Fig. 7A,C) and MLKL (Fig. 7B,C). The IC50 values of other agents were displayed in Fig. 7F–I.

Fig. 7
figure 7

Identification of candidate agents with higher drug sensitivity in high-risk group patients. The bubble chart showed the correlation coefficient between PANGPI, MEFV, MLKL, TLR4 and the estimate IC50 of candidate drugs in GDSC1 (A) and GDSC2 (B) datasets. (C) The intersecting results in GDSC1 and GDSC2 datasets. The estimated IC50 of 3 candidate drugs were compared between high- and low-risk groups in GDSC1 (D–F) and GDSC2 (G–I) datasets.

Full size image

Pharmaceutical analysis based on up-regulated PANGPI gene

To further verify the strong correlation between the up-regulated PANGPI gene and chemotherapeutic drug BMS-536924, we performed molecular docking to determine the binding modes of BMS-536924 with 9 PANoptosis biomarkers including Apoptosis (Caspase-3, Caspase-7 and Caspase-8), Pyroptosis (Caspase-1, GSDMD and GSDME) and Necroptosis (MLKL, RIPK1 and RIPK3). The results showed that all of them could successfully bind to BMS-536924 (docking score in the range of − 5.9 to − 9.4 kcal/mol), which provided more evidence for its modulation effect of PANoptosis. Among the docked protein, MLKL has the second highest binding affinity (− 8.5 kcal/mol) with 17 intermolecular bonds and 2 hydrogen bonds (Fig. 8A–R) (Supplementary Table S10).

Fig. 8
figure 8

Molecular docking of predicated drugs with the core targets of PANoptosis. (A–F) The 3D and 2D docking patterns of Apoptosis targets. (G–L) The 3D and 2D docking patterns of Pyroptosis targets. (M–R) The 3D and 2D docking patterns of Necroptosis targets.

Full size image

In order to enhance the reliability of the obtained results and confirm the stability of the system, a 100 ns MD simulation of the MLKL protein and BMS-536924 complex was assessed. The root-mean-square deviation (RMSD) plot indicated that the system gained stability after 20 ns. The average radius of gyration (Rg) valves (1.25–1.85 nm) and the resulting solvent accessible surface area (SASA) signify the stability of the ligands with MLKL. The root-mean-square fluctuation (RMSF) data showed that the fluctuations of most amino acids are less than 0.45 nm (Fig. 9A–D) (Supplementary Table S11). During MD simulation, BMS-536924 has hydrogen bonds occupancy with Lys230, Gly349, Glu351, Gln356, Thr357 and Thr363, and has more intermolecular bonds after simulation (Fig. 9E,F). Based on the free binding energy calculation (Fig. 9G–I) (Supplementary Table S12), the predicted results in this study were mostly energetically favorable.

Fig. 9
figure 9

Molecular dynamics simulation of BMS-536924 with the MLKL protein. (A) Protein–ligand complex RMSD. (B) Protein–ligand complex Rg. (C) Protein–ligand complex SASA. (D) Protein–ligand complex RMSF. (E) Hydrogen bond dynamics between BMS-536924 and MLKL protein. (F) Docking patterns of BMS-536924 and MLKL protein after simulation. (G) Graphical representation of binding free energies of protein–ligand complexes. (H) Binding modes of the BMS-536924 and MLKL protein. (I) Per-residue energy decomposition of the key residues.

Full size image

Discussion

DLBCL is one of the most common aggressive malignancies globally, and the OS of DLBCL patients remains suboptimal. Due to its diverse characteristics, the prognosis of DLBCL patients may no longer be accurately predicted using traditional scoring systems30. This highlights the urgent need to develop novel biomarkers that can help in precise stratification and individualized therapy of relapsed/refractory DLBCL patients. Transcriptomic analysis of the DLBCL microenvironment reveals cancer cell-extrinsic biology, by which the inflammatory subtype has been identified as one of the major LME categories31. PANoptosis, a recently defined inflammatory PCD, is currently the subject of intense research32, however the association between PANoptosis and DLBCL is poorly studied. We aimed to link PANoptosis to the prognosis of DLBCL patients by uncovering the role of PANRGs in LME and therapeutic response.

In our study, MLKL, TLR4 and MEFV were finally pinpointed to be associated with OS in DLBCL, and a robust risk model was successfully constructed based on these three prognosis related PANRGs. Mixed lineage kinase ___domain-like protein (MLKL) emerged as executioner of necroptosis, and its oligomerization was observed in releasing of pro-inflammatory intracellular contents during PANoptosis33. The priming of Toll-like receptors (TLRs) accelerates immune cell inflammatory responses, among them, TLR4-dependent mediated signals are pivotal for the activation of PANoptosome34. The Mediterranean fever (MEFV) gene encodes pyrin, an innate immune sensor, which drives the assembly of PANoptosome and cell death35. In conjunction with our study, above three PANRGs exhibit a significant association with patient prognosis, molecular characteristics and immune landscape.

DLBCL patients were classified into high- and low-risk groups based on the PANGPI scores, and notably, patients with high PANGPI scores exhibited a worse prognosis in both the external validation and the training cohorts. A nomogram model was subsequently established according to clinical characteristics as well as the risk model, to accurately predict the survival time of patients. The calibration plots revealed that the nomogram predicted survival rates were very close to the actual survival rates. From a clinical point of view, the PANoptosis signature has the potential to revolutionize treatment paradigms of DLBCL, identify patients with a poor prognosis and provide precise interventions.

Our study indicated that DLBCL patients with low- and high-risk groups exhibited significant variations in somatic mutations. The mutation frequency of MUC16, CD79B, IRF4 and ZNF208 gene was higher in the high-PANGPI group, and there was extensive co-mutation between IRF4 and CD79B. Coincidently, it has been reported that IRF4-rearranged DLBCL patients shared recurrent CD79B mutations with subsequent activation of the nuclear factor kappa-B (NF-κB) pathway that led to inflammatory responses36. Besides, there were more MCD genotype (including MYD88L265P and CD79B mutations) patients in the high-PANGPI group, suggesting a shared pathogenesis between PANoptosis and DLBCL. GSEA analysis depicted that the pathways of pattern recognition receptors (PRRs) such as NOD-like receptor (NLRs) and TLRs might be activated in the high-risk group. Recent genetic and biochemical analyses found that an NLR network regulated by TLR signaling act as central components of PANoptosome which drove inflammatory cell death37.

TILs encompass various immune cell such as B-cells, T-cell and natural killer cells with diverse phenotypic and functional properties, that contribute to “hot” or “T-cell inflamed” tumors38. However, the majority of DLBCL are either “cold” or “non-inflamed”39, thus offering the potential for PANoptosis to reprogram the LME. In our researtch, Th1 cells, CD4+ T cells, CD8+ T cells and B cells were positively correlated with PANGPI risk scores, while NK cells and Th2 cells were negatively correlated with PANGPI risk scores. We further measured the cytolytic activity and inflammatory states of TILs and found that the degree of CYT and T cell-inflamed GEP were higher in the high-PANGPI group. However, this immune cytolytic activity may not be directly associated with CD8+T cells, since the TIDE dysfunction scores was higher in the high-PANGPI group. It is noteworthy that CD4+ T cells also exert tumor antigen-specific cytolytic effect against DLBCL cells40. Furthermore, Th1 derived cytotoxic CD4+ T cells can initiate an indirect inflammatory cell death of remote cancer cells, by which rise the tumor immunogenicity41. Strikingly, we found that the expression of immune checkpoint biomarkers interrelated with immune cytolytic activity in the high-PANGPI group, which may be attributed to the inciting of other functional effectors by inadequately primed TILs42, or the inflammatory conditions caused by PANoptosis43.

A growing number of anticancer drugs inducing PANoptosis to overcome tumor resistance or circumvent immune suppression represent a prospective avenue for DLBCL treatments44,45. Due to the worse outcome of high-PANGPI patients, we explored candidate agents with negative correlation between their IC50 and PANGPI scores. Through statistical analysis, BMS-536924, Gefitinib and Navitoclax were obtained. Among them, the IC50 of BMS-536924 was negatively correlated with the up-regulated PANGPI gene as well, suggesting its therapeutic benefits may be involved in PANoptosis modulation. BMS-536924, a small molecule inhibitor of the insulin-like growth factor 1 receptor (IGF1R), has been experimentally confirmed of anticancer activity46. The results of in silico pharmacology revealed that BMS-536924 has a good binding affinity and stability to PANoptosis biomarkers, which may be applied to achieve drug repurposing for DLBCL.

Inevitably, our research is subject to several limitations. Firstly, this study relies on publicly available retrospective database, thus the clinical applicability of the PANGPI scoring system should be further validated in the prospective real-world studies. Additionally, we lack adequate single cell sequencing and spatial transcriptome data to support our analysis of the somatic mutation and immune landscape of DLBCL. Finally, our in-silico drug sensitivity results should be further confirmed in vivo or in vitro experiments to explore its potential of combination therapies for challenging cases of DLBCL.

Conclusion

To our best knowledge, it is the first study to evaluate the clinical, molecular, and immune characteristics associated with PANoptosis in DLBCL. A PANoptosis related gene signature was identified, DLBCL patients with the higher PANGPI-score may have a dismal survival but a potentially better cytolytic activity and inflammatory states of the TILs. Overall, this article provides a preliminary idea which may shed light on future investigations into the regulatory mechanisms of PANoptosis in DLBCL for prognostication and therapeutic decision-making.