Defining the expression of PRGs in GBM
A total of 42 differential expressed PRGs in GBM were recognized through the comparison of the differential expressions of 49 PRGs in 120 normal subjects and 155 GBM from the TCGA and GTEx database (P < 0.05). The expression levels of 13 genes (NLRP2, NLRP7, TNF, IL1B, IL6, NLRP1, IL1A, PRKACA, GSDMB, CHMP3, CHMP4B, CHMP2B, and CHMP7) in GBM were downregulated relative to those in the brain tissues of normal subjects, whereas the expression levels of 29 genes (PLCG1, CASP9, CHMP6, SCAF11, CHMP4A, CHMP2A, HMGB1, IRF2, NLRC4, GZMB, NOD1, TP63, BAK1, CASP8, BAX, IRF1, CASP3, PYCARD, IL18, GSDMD, AIM2, GSDMA, NOD2, GZMA, CASP5, CASP6, TP53, CASP1, and CASP4) were upregulated (Figure 1A). To further explore the interactions of these differentially expressed PRGs, we carried out PPI analysis. As shown in Figure 1B, the minimum required interaction score was set at 0.4, and the degree of each node in the network was expressed and arranged according to color. The top 10 genes in the term of degree were TNF, IL1B, IL18, CASP8, IL6, PYCARD, AIM2, NLRC4, AIM2, and NLRC4. The correlation network results of the differentially expressed PRGs are shown in Figure 1C. Most differentially expressed PRGs were positively related, whereas negative correlations between GZMA and CHMP2A, between IL6 and PLCG1, and between GSDMB and CHMP4A were found.
Analysis of the somatic mutation characteristics of PRGs
Somatic mutation frequency of 49 PRGs in GBM was analyzed. The results demonstrated that missense mutation was the most common mutation classification, and SNP was the most common variant type. Single-nucleotide variation mainly occurred in the form of C > T (Figure 2A). In the GBM samples, the top 10 mutation genes related to pyroptosis were NLRP3, NLRP7, NLRP2, NOD1, CASP1, NLRP1, PLCG1, GZMB, NOD2, and CHMP4C (Figure 2B).
Functional enrichment analysis of the PRGs
GO and KEGG databases were used in analyzing the related pathways of differentially expressed PRGs and clarifying their functions. As shown in Figure 3A, the differentially expressed PRGs were mainly enriched in the positive regulation of cytokine production, regulation of interleukin-1 production, positive regulation of cysteine-type endopeptidase activity, and pyroptosis biological processes. In addition, the genes were mainly related to multiple cell components, such as inflammasome complex, ESCRT complex, multivesicular body, late endosome membrane, pore complex, and immunological synapse. Moreover, molecular function analysis verified that the PRGs were correlated to endopeptidase activity, cytokine receptor binding, protease binding, cysteine-type peptidase activity, peptidase activator activity, and CARD domain binding. KEGG pathway enrichment analysis demonstrated that the PRGs were closely related to NOD-like receptor signaling pathway, necroptosis, Salmonella infection, legionellosis, lipid and atherosclerosis, and pathogenic Escherichia coli infection (Figure 3B and C).
Construction of the prognostic PRG model
According to the survival-related genes screened through univariate Cox regression analysis, 34 differentially expressed PRGs (CHMP4A, GSDMA, CHMP3, AIM2, CHMP2B, CASP5, NLRC4, BAX, IL18, IRF1, IL1B, CASP8, NLRP1, IL6, IL1A, CASP1, CHMP6, TP63, PRKACA, NLRP2, PYCARD, GZMA, CHMP7, CHMP2A, CHMP4B, CASP6, GSDMD, CASP3, NOD2, NOD1, IRF2, CASP4, GZMB, and NLRP7) had HRs of >1, indicating that they were related to increased risk. The other genes (CASP9, BAK1, BAK1, HMGB1, GSDMB, PLCG1, SCAF11, and TNF) had HRs of <1 and were protective factors (Figure 4A).
To evaluate the prognosis value of these differentially expressed PRGs, a Lasso-Cox regression model was built with “glmnet” R package and used in further screening. According to the minimum penalty parameter (λ), five of the 42 differentially expressed PRGs were retained, namely, CASP3, NLRP2, TP63, GZMB, and CASP9 (Figure 4B and C). A prognosis risk score formula was established according to weighted regression coefficient and expression levels based on multivariate Cox regression analysis: Risk score = 0.0046 × (expression value of CASP3) + 0.0060 × (expression value of NLRP2) + (0.0227) × (expression value of TP63) + 0.1804 × (expression value of GZMB) - (0.0481) × (expression value of CASP9). According to the median score calculated in the prognosis risk score formula, 152 patients with GBM were divided into low- and high-risk groups. PCA and t-SNE analysis demonstrated that patients with different risks can be divided into two types (Figure 4D and E). The risk score distribution and survival states of the patients are shown in Figure 4F and G. As risk score increased, the risk of death in the patients increased, whereas survival time was shortened. According to the Kaplan–Meier survival curve, the low-risk group had more obvious survival advantages (P = 0.009, Figure 4H). The area under curve values of the 1-, 3-, and 5-year ROC curves were 0.675, 0.669, and 0.723, respectively (Figure 4I), indicating that the pyroptosis model can predict the 1-, 3-, and 5-year survival rates of patients with GBM.
Independent prognostic analysis and nomogram construction based on risk model
Whether risk score from the gene characteristic model can be used as an independent prognosis biomarker of GBM was determined using univariate and multivariate Cox regression analysis. Univariate and multivariate Cox regression analyses demonstrated that the risk score was an independent factor that influenced the prognosis of patients with GBM (HR = 6.917, 95% CI: 2.760–17.330 and HR = 6.168, 95% CI: 2.467–15.420; Figure 5A and B). In addition, the clinical feature heatmap in Figure 5C showed that THE gender and age of patients showed no significant differences between the low- and high-risk groups (P > 0.05).
To discuss THE importance of these five prognosis PRGs to THE prognosis prediction of patients with GBM, we constructed a nomogram to visualize the Cox regression model as the survival probability of patients. The performance of the model was evaluated using a calibration curve. Compared with the ideal model, the predicted total 3-year survival rate of patients with GBM based on prognosis PRGs was slightly different from the practical value. All the five prognosis PRGs had good prognosis prediction ability for GBM (Figure 5D and E).
The risk score was closely related to the expression levels of the five prognosis PRGs. It was an important factor for predicting the prognosis results of patients with GBM and had good prediction ability. The five prognosis PRGs might participate in the progress of GBM.
Correlation between the expression of prognostic PRGs and immunocyte levels and immunoregulators in GBM
Using the TIMER database, we detected the correlation of five prognosis PRGs with the levels of immunocyte infiltration in GBM. As shown in Figure 6A, CASP3 was significantly associated with CD8+ T cell (r = 0.109; P = 2.65e-02), macrophage (r = 0.084; P = 8.72e-02), neutrophil (r = 0.155; P = 1.50e-03), and dendritic cell (r = 0.193; P = 6.89e-05). NLRP2 was significantly associated with CD8+ T cell (r = −0.197; P = 4.96e-05). TP63 was significantly associated with CD4+ T cell (r = −0.143; P = 3.34e-03), macrophage (r = −0.167; P = 5.93e-04), and neutrophil (r = −0.196; P = 5.65e-05). GZMB was significantly associated with CD8+ T cell (r = −0.203; P = 3.00e05) and CD4+ T cell (r = −0.174; P = 350e-04). CASP9 was significantly associated with CD4+ T cell (r = 0.199; P = 4.22e-05) and macrophage (r = 0.174; P = 3.48e-04). The above results showed a significant correlation between prognosis PRGs and tumor immunocyte infiltration.
To further address the influences of prognostic PRGs on tumor immune response, the correlation between the expression of prognostic PRGs and immunoregulators was calculated according to the TISIDB database. CASP3 and CASP9 were negatively correlated with immunoinhibitor, immunostimulator and MHC molecule in GBM, whereas NLRP2, TP63, and GZMB were positively related to them (Figure 6B). Hence, the expression levels of the prognosis PRGs were closely related to the abundance of immunocyte infiltration and immunoregulators.
Immunocyte infiltration between the high-risk and low-risk groups
We used the CIBERSORT algorithm to analyze the abundance of immunocyte infiltration between the high-risk and low-risk groups in patients with GBM. The high-risk group presented higher distribution levels of monocytes, activated dendritic cells, activated mast cells, and eosinophils and lower distribution levels of resting NK cells and resting mast cells (Figure 7A). According to the correlations among 22 types of immunocytes, activated dendritic cells were positively related with naive CD4 T cells (r = 0.47), whereas M0 macrophages were negatively related with regulatory T cells (r = 0.42). The resting NK cells were negatively correlated with activated NK cells (r = −0.69), monocytes, and M2 macrophages were negatively correlated with M0 macrophages (r = −0.67; r = −0.6; Figure 7B).
The differential expression of immunocyte infiltration and immune-related function score between the high- and low-risk groups was analyzed by ssGSEA algorithm. Significant differences in B cells, CD8+ T cells, dendritic cells, macrophages, neutrophils, NK cells, T helper cells, tumor-infiltrating lymphocytes, and regulatory T cells were observed between the groups (Figure 7C). Moreover, significant differences in the scores of APC co-inhibition, APC co-stimulation, CCR, check point, cytolytic activity, HLA, inflammation-promoting process, parainflammation, T-cell co inhibition, T-cell co stimulation, type I FN response, and type II FN response were found between the groups (Figure 7D).
Significant difference in immunocyte infiltration was found between the groups. Therefore, monocytes, dendritic cells, eosinophils, NK cells, and mast cells might be potential immunocytes that play important roles and participated in the progress of GBM.
Construction of the network of mRNA–miRNA–lncRNA
To elaborate the potential molecular mechanism of the five pyroptosis-related mRNAs (CASP3, NLRP2, TP63, GZMB, and CASP9) which had prognosis values in GBM, an mRNA–miRNA–lncRNA co-expression network was built. miRNA targets that combined with these five prognostic PRGs were predicted reversely by using mirTarBase and Starbase. A total of 514 miRNA targets of CASP3, NLRP2, TP63, and CASP9 were obtained (Figure 8A). The degree values among miRNAs nodes were calculated using the CytoHubba plug-in of Cytoscape software, and the top 10 highly connected miRNAs were determined (hsa-miR-342-3p, hsa-miR-582-5p, hsa-miR-3163, hsa-miR-320a, hsa-miR-320c, hsa-miR-320d, hsa-miR-320b, hsa-miR-526a, hsa-miR-519b-5p, and hsa-miR-519c-5p; Figure 8B). The OncoLnc database determined a miRNA, hsa-miR-519c-5p, and GBM patients with high expression levels of hsa-miR-519c-5p. These patients had relatively high total survival rates (P = 0.0337; Figure 8C). Based on the above miRNA, the upstream lncRNA targets were explored for the construction of a miRNA–lncRNA axis. As shown in Figure 8D, for the intersection lncRNAs of Lncbase and Starbase, KCNQ1OT1, GABPB1-AS1, ENTPD1-AS1, XIST, LINC01018, and LINC00662 were determined as targets. According to survival analysis and differential expression of the lncRNA targets, only GABPB1-AS1 significantly decreased the survival probability of patients with GBM (P = 0.05; Figure 8E) and it was upregulated significantly in the GBM samples (|log2FC| < 1; P < 0.01; Figure 9F). These results indicated that the regulatory axis of lncRNA GABPB1-AS1/hsa-miR-519c-5p/CASP3/TP63 plays an important role in the expression of GBM.
The GSEA results showed that CASP3 and TP63 were enriched on ECM-related signaling pathways (e.g. ECM receptor interaction, adherens junction, and TGFβ signaling pathways) and immune-related signaling pathways (e.g. Toll-like receptor and NOD-like receptor signaling pathways; Figure 9). This result showed that the two pyroptosis-related mRNAs are closely related to the growth, metastasis, and diffusion of GBM, and they play important roles in the immunoregulation of GBM TME.