We downloaded the genomic data, copy number alteration, methylation and clinical data of glioma patients from the CGGA (http://www.cgga.org.cn/) and TCGA databases (https://portal.gdc.cancer.gov/). Additional gene-centric RMA-normalized gene expression profiles and drug response data of over 1000 cancer cell lines were accessed from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/downloads). Immune-associated data, including immune cells and immunophenoscores, were downloaded from TCIA (https://tcia.at/home). Thirty-three pyroptosis-related genes were defined from a previous publication and are provided in Supplementary Table S1 (Man and Kanneganti, 2015; Wang and Yin, 2017; Karki and Kanneganti, 2019; Xia et al., 2019).

We identified the optimal clustering number visualizing consensus matrix, tracking plot, and cumulative distribution function plot. In addition, a T-distributed stochastic neighbor embedding-based approach was used to validate the clustering in glioma patients. We calculated the enrichment scores for every sample using the GSVA R package.

We developed a pyroptosis-related prognostic signature based on the CGGA training cohort. Twenty differentially expressed genes with p < 0.05 were entered into LASSO Cox regression, which identified potential genes for the prognostic signature in the CGGA training cohort. Then, we calculated the risk score for each sample of the CGGA and TCGA validation cohorts using the obtained regression coefficient in the CGGA training cohort: risk score = coef1*gene₁ expression + coef₂* gene₂ expression + …coef_n*gene_n expression. The CGGA and TCGA samples were divided into a high-risk group and a low-risk group based on the median risk score. Receiver operating characteristic curves were plotted to evaluate the 1-year, 2-years, and 3-years sensitivity and specificity of the prognostic signature. We also established a prognostic nomogram to evaluate the clinical value of the prognostic signature. Calibration analysis of the prognostic predictive value of the nomogram was carried out.

Gene Ontology and KEGG pathway analyses were performed using the “clusterProfiler” package. We used single-sample gene set enrichment analysis (ssGSEA) to estimate the enrichment score of stem cell-like properties (RNAss, DNAss) and the TME (stromal score, immune score, and ESTIMATE score) in the TCGA cohort because the CGGA dataset did not provide such data. The immune-related cell and function scores were also calculated for each sample (downloaded from https://www.gsea-msigdb.org/).

Based on the risk groups in the TCGA cohort, we compared the somatic copy number alteration, mutation, and DNA methylation levels between the high-risk and low-risk groups using the “limma” R package.

To further explore the transcriptome regulation network of different risk groups, we used Cytoscape version 3.8.2 to establish a lncRNA-miRNA-mRNA regulatory network. We explore the correlation between small molecular drugs and the identified prognostic signature genes using Pearson correlation analysis |R|>0.25 and p > 0.05 were considered significant.

We further performed the Western blot, cell migration assays, cell scratchy assays, and clonogenic assays to verify the present finding. We selected the CASP8 to validate the molecular function because CASP8 showed significant differences between normal tissue and GBM or LGG, and the elevated expression is associated with poor prognosis.

The details of experiments process in vitro were supplied in Supplementary Material.

The log-rank test was used to compare the survival curves of Kaplan-Meier analysis. The hazard ratio (HR) and 95% confidence interval (CI) of each gene and clinical parameters were calculated when univariate and multivariate Cox regression were applied. All analyses were achieved using R software version 4.0. A two-sided p value < 0.05 was considered significant unless otherwise specified.

The flow chart of the data analysis is presented in Figure 1A. From two CGGA RNA-seq datasets, we obtained 1018 samples of gene expression data and further identified 30 pyroptosis-related genes based on MAD>0.5. The gene symbols and descriptions of the 30 pyroptosis-associated genes used for classification are listed in the Supplementary Table S1. We first explored the interactions among these genes using PPIs (Figure 1B), and the PPI network indicated that CASP8, CASP4, CASP1, NLRP3, NLRP1, and NLRC4 are hub genes. The correlation circle plot of the 30 genes is presented in Figure 1C (red: positive correlation; green: negative correlation). We identified the optimal k value as 2 by estimating the comprehensive correlation coefficient. Therefore, we divided the glioma samples into two different subclasses: cluster 1 and cluster 2. For the optimal k value (k = 2), the consensus matrix showed a relatively sharp and clear boundary, indicating stable and robust clustering (Figure 1D). To verify the subclass stability, we further performed t-sensitivity PCA and found that a two-dimensional t-sensitivity distribution supported subtype clustering (Figure 1E). The consensus clustering for each sample is listed in Supplementary Table S2. The Kaplan-Meier analysis indicates that the median survival time was significantly shorter in cluster 2 than in cluster 1 (MST: 1.87 vs. 6.92 years, p < 0.001, Figure 1F). This result indicated that the two subclasses had distinct prognostic patterns.

Two subclasses were obtained based on pyroptosis-related genes. To explore the pathway enrichment for the two subclasses, we performed GSVA by transforming the expression data from a gene-by-sample matrix to a gene set by two subclasses. Then, differential pathways were enriched in the two subclasses. Compared with cluster 1, the GSVA results indicated that cluster 2 had 182 kinds of significantly differential signaling pathways (Supplementary Table S3). The upregulated pathways were associated with immune-related pathways, such as autoimmune, allograft rejection, graft vs. host disease, primary immunodeficiency, antigen processing, and presentation. Some signaling pathways, such as the cytosolic DNA sensing pathway, NOD-like receptor signaling pathway, Toll-like receptor signaling pathway, and metabolism-related pathways, were also significantly enriched. The significantly downregulated pathway was long-term potentiation (Figure 2A).

We explored the correlation of subclasses with clinical characteristics (Figure 2B). Compared with patients in cluster 2 with a favorable prognosis, patients in cluster 1 tended to have GBM (p < 0.001), WHO grade IV (p < 0.001), a higher proportion of age >41 years, 1p19q non-codeletion status (p < 0.001), and IDH wildtype status (p < 0.001). Sex, PRS type and radiotherapy status were not associated with the molecular subclasses (p > 0.05). For the pyroptosis-related genes except CASP9, significant differential expression was observed in the two clusters. Among these differentially expressed genes, all genes were upregulated in cluster 1 and downregulated in cluster 2 (Figure 2B). We also compared the differences in pyroptosis-related genes in patients with different histologies, grades, IDH mutation statuses, and 1p19q statuses. Compared with the LGG group, the GBM group had one upregulated gene (AIM2) and 21 downregulated genes (Supplementary Figure S1A). Twenty-one DEGs were found for grade, and their expression increased with increasing WHO grade (p < 0.005, Supplementary Figure S1B). For IDH status, 25 DEGs were found (Supplementary Figure S1C). Thirty pyroptosis-related DEGs were found for 1p191 status (Supplementary Figure S1D).

We further performed differential expression analysis between cluster 1 and cluster 2. A total of 392 DEGs were found, 18 genes were upregulated, and 372 genes were downregulated in cluster 2 (Supplementary Table S4). GO and KEGG enrichment analyses were performed for all DEGs (Supplementary Tables S5, S6). A total of 874 differentially expressed functions were enriched, including 709 biological processes, 95 cellular components and 70 molecular functions. The top 30 enrichment results are presented in Supplementary Figure S2. Most of these functions were associated with immunity. In addition, 56 pathways were also identified in the KEGG analysis (Supplementary Figure S3), and the top five pathways were phagosome, Staphylococcus aureus infection, tuberculosis, complement and coagulation cascades, and human T-cell leukemia virus 1 infection.

To explore the tumor heterogeneity between the two subclasses, we investigated the immune cell and immune function differences. Compared with cluster 2, cluster 1 had higher aDC, CD8⁺ T cell, DC, iDC, macrophage, mast cell, neutrophil, NK cell, pDC, T helper cell, Tfh cell, Th2 cell, TIL, and Treg levels (all p < 0.001, Supplementary Figure S4A). Similarly, cluster 1 had higher immune function scores than cluster 2, including APC coinhibition, APC costimulation, CCR, checkpoint, cytolytic activity, HLA, inflammation promotion, MHC class I, parainflammation, T cell coinhibition, type I IFN response and type II IFN response (all p < 0.001, Supplementary Figure S4B).

Initially, we performed univariate Cox regression to identify the correlations of the 30 pyroptosis-related genes with OS (Supplementary Figure S5A) in the CGGA cohort. In total, 20 pyroptosis-related genes were identified as associated with the overall survival of glioma patients. The Kaplan-Meier plot indicated that high expression of CASP3, CASP4, CASP5, CASP6, CASP8, ELANE, GSMAD, IL6, NLRP3, NOD1, NOD2, PLCG1, PRKACA, PYCARD, and SCAF11 was associated with poorer OS in glioma. Using 20 prognostic pyroptosis-related genes, we developed a prognostic signature by performing LASSO regression in the CGGA training cohort (Supplementary Figures S5B,C). Fifteen of the 20 prognostic genes were used to develop the risk signature. We calculated the risk score for each sample using the regression coefficients of the 15 genes (Supplementary Table S7). Glioma patients with risk scores greater than the median value were divided into a high-risk group, and the others were divided into a low-risk group. Compared with the low-risk group, the high-risk group was more likely to have GBM (p < 0.001), a higher WHO grade (p < 0.001), recurrence (p < 0.001), older age (p < 0.001), IDH wildtype status (p < 0.001), 1p19q non-codeletion status (p < 0.001), and a history of chemotherapy (p < 0.001). The heatmap showed the association between the risk group and clinical parameters and differentially expressed genes of the high- and low-risk groups (Supplementary Figure S5D). Furthermore, we found that glioma patients belonging to cluster 1, patients with a poor prognosis, patients with GBM, patients with WHO grade IV patients with 1p19q non-codeletion status and patients with IDH wildtype status had higher risk scores (all p < 0.001, Supplementary Figure S6).

The Kaplan-Meier analysis showed that the high-risk group had a significantly poorer OS than the low-risk group (Figures 3A,B). Univariate Cox regression indicated that the risk score was positively associated with OS in glioma (HR = 3.105, 95% CI: 2.681–3.596, p < 0.001, Figure 3C). Multivariate Cox regression suggested that the risk score was an independent unfavorable prognostic predictor in glioma (HR = 1.685, 95% CI: 1.392–2.039, p < 0.001, Figure 3D). In addition, PRS type, tumor grade, and age were positively associated with OS. However, chemotherapy, wildtype IDH status, and 1p19q status were negatively associated with OS in the CGGA training cohort. The PCA plot indicated that patients in different risk groups were separated into obviously different clusters (Figure 3E). Time-dependent receiver operating characteristic analysis was performed to evaluate the predictability of the prognostic model. Our results showed that the AUCs at 1, 2, and 3 years were 0.717, 0.784, and 0.773 (Figure 3F), respectively. We further compared the OS status among different histology, IDH status, 1p19q codeletion status, and grade subgroups. The results showed that the OS of the high-risk group was still poorer than that of the low-risk group (Supplementary Figure S7, all p < 0.001).

To further validate the prognostic value of the pyroptosis-related gene model, we also calculated the risk score of glioma patients in the TCGA cohort using the regression coefficients of the CGGA cohort. The Kaplan-Meier analysis indicated a significant correlation of the high-risk group with worse OS than the low-risk group (Figures 4A–C). Univariate Cox regression showed that the risk score was significantly associated with OS in the TCGA cohort (HR = 2.084, 95% CI: 1.890–2.297, p < 0.001, Figure 4D). In multivariate Cox regression, the risk score was also an independent prognostic indicator (HR = 1.425, 95% CI: 1.247–1.629 p < 0.001, Figure 4E). The PCA plot validated the high- and low-risk distribution of all glioma patients based on the TCGA cohort. Furthermore, the AUCs of the risk score were 0.844 at 1 year, 0.863 at 2 years, and 0.874 at 3 years (Figure 4F).

To further evaluate the clinical prediction value of the prognostic signature, we constructed a prognostic nomogram model based on multivariate Cox regression analysis that included all clinical parameters in the CGGA cohort The calibration curves indicated that the clinical nomogram model could precisely predict the 1, 3, and 5-years OS of glioma patients (C-index = 0.799). The predictive accuracy of this nomogram was well validated in the TCGA cohort (C-index = 0.841, Figure 5).

We further explored the underlying biological functions that define the survival of glioma patients. We first performed DEG analysis between the high-risk and low-risk groups and then annotated the functions of the DEGs in terms of biological processes, cellular components, and molecular functions using GO enrichment and KEGG pathways. We identified 338 DEGs in the CGGA cohort (Supplementary Table S8) and 2600 DEGs in the TCGA cohort (Supplementary Table S9). The GO enrichment and KEGG pathway analyses indicated that the CGGA and TCGA cohorts shared some enrichment results, such as extracellular matrix organization, extracellular structure organization, immune response, ECM-receptor interaction, and cell adhesion molecules (Figures 6A–D).

We also explored the differences in immune cells and immune functions based on the risk score in the CGGA (Figures 6E,G) and TCGA datasets (Figures 6F,H). As shown in the box plots, the immune cell score showed a similar trend in the CGGA and TCGA datasets. All immune cell scores were significantly upregulated in the high-risk group. The immune function differences of the different risk groups were the same in the CGGA and TCGA datasets (all p < 0.001). All immune function scores were significantly upregulated in the high-risk group. Significant expression levels were also observed among different immune subtypes, which indicated that the glioma prognosis risk could be associated with immune status (Supplementary Figure S8). We also explored the correlation of the expression of target genes with cancer stem cell-like properties (RNAss, DNAss) and the TME (stromal score, immune score, and ESTIMATE score). We found that PCG1 was negatively associated with RNAss, the stromal score, the immune score, and the ESTIMATE score. SCAF11 was only negatively associated with DNAss. The rest of the genes showed positive correlations with RNAss, DNAss and the stromal, immune and ESTIMATE scores (Supplementary Figure S9).

Molecular alterations of pyroptosis-related genes were also evaluated based on histology in the TCGA dataset. NLRP2, NLRP7, and PLCG1 were the only gene alterations in LGG, and NLRP3, NLRP7, NLRP2, SCAF11, NOD1, PLCG1, NLRP1, and CASP1 were gene alterations in GBM. All gene alterations were within 2% (Figure 7). The somatic copy number alteration analysis indicated significant differences among the pyroptosis-related genes. Among these genes, the copy variation number was significantly increased in GPX4, NLRP7, NLRP2, CASP3, CASP6, IL1B, CASP8, IL6, AIM2, NLRP4, NLRP3, PRKACA, ELANE, SCAF11, CASP9, NOD1, and PLCG1 and was significantly decreased in GSDMB, GSDMD, NLRP1, CASP9, TIRAP, CASP1, CASP4, NOD2, CASP5, PYCARD, GSMDC, GSMDA, and IL18 in the high-risk group. The DNA methylation levels of the pyroptosis-related genes were also compared. The results showed that the overall DNA methylation levels were significantly decreased in the high-risk group and increased in the low-risk group.

A ceRNA network was constructed based on the differentially expressed mRNAs, lncRNAs and miRNAs between the high-risk and low-risk groups in the TCGA dataset. We identified 763 downregulated mRNAs, 1176 upregulated mRNAs, 116 downregulated lncRNAs, 132 upregulated lncRNAs (Supplementary Table S10), 47 downregulated miRNAs and 71 upregulated miRNAs (Supplementary Table S11). Finally, 39 mRNAs (28 upregulated and 11 downregulated), 26 lncRNAs (15 upregulated and 15 downregulated) and 14 miRNAs (13 upregulated and 1 downregulated) were included in the ceRNA network (Figure 8). The Kaplan-Meier curves suggested that 13 lncRNAs (positive correlation: AC025211.1, AC068643.1, GDNF-AS1, and LINC00519; negative correlation: ADH1L1-AS2, CRNDE, FAM181A-AS1, HOTAIRM1, MCF2L-AS1, MIR210HG, NEAT1, SLC6A1, and SNHG9; Supplementary Figure S10), 41 mRNAs (Supplementary Table S12, Supplementary Figure S11) and 8 miRNAs (miR-21, miR-155, miR-200a, miR-216a, miR-221, miR-222, miR-429, and miR-503; Supplementary Figure S12) were associated with OS in glioma patients.

To identify potential target drugs, we performed correlations of the identified prognostic signature genes with drugs. We identified 257 pairs of significant gene-drug correlations (Supplementary Table S13). There were 9 pairs with correlation coefficients >0.5 or < −0.5.

ELANE-hydroxyurea, ELANE-cyclophosphamide, CASP3-nelarabine, NOD2-imiquimod, NLRP3-rebimastat, ELANE-ABT-199, ELANE-imexon, and NOD2-isotretinoin showed drug sensitivity. PRKACA-cobimetinib showed drug resistance (Figure 9).

We selected the CASP8 to validate the molecular function because: First, the significantly high expression of CASP8 is found between glioma and normal control (Supplementary Figure S16). Second, CASP8 is included the prognosis model based on pyroptosis-relate genes (Supplementary Table S7). Third, the elevated expression of CASP8 is more significantly associated with poor overall survival (Supplementary Table S7). We firstly detected the expression of CASP8 in glioma cell lines using the Western blot analysis, and found CASP8 is the most highly expressed in LN299 cell. We built the CASP8-si LN229, H4 and U87 cells of glioma. The qPCR indicated mRNA level of CASP8 is significantly down-regulated in U87 and LN229 cells. Furthermore, the silence of CASP8 expression inhibited the cell migration ability (Figure 10). The clonogenic assay also showed that the number of clonogenicity of U87 and LN229 cells were significantly suppressed after knockout of CASP8. These results suggested that CASP8 promotes the progression of glioma cells.