A pyroptosis gene-based prognostic model for predicting survival in low-grade glioma

Data collection

Gene expression profiles and corresponding clinical information from 529 LGG patients were collected from The Cancer Genome Atlas database (TCGA, https://portal.gdc.cancer.gov/). We excluded 53 samples with very short survival times (OS<60 days) and used the remaining 476 samples for the subsequent analysis. Due to the lack of relevant data on normal brain tissue in the TCGA database, data from 1,152 normal brain tissue samples were downloaded from the Genotype–Tissue Expression database (GTEx, https://commonfund.nih.gov/gtex) and served as a normal control group. In addition, the gene expression profiles and corresponding clinical information for 50 LGG patients in the GSE43378 dataset from the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/) were used for independent validation.

Acquiring PGs

We searched for “Pyroptosis” in the Gene Set Enrichment Analysis (GSEA, http://www.gsea-msigdb.org/gsea/index.jsp) database and downloaded gene sets associated with PGs, which allowed us to obtain 27 PGs. We combined this list with the list of 51 PGs in a previous article (Zhang et al., 2021) to yield a total of 51 PGs.

Identification of DEPGs

Tumor tissue samples from the TCGA database and normal tissue samples from the GTEx database were integrated for subsequent research. Thereafter, we used the R package “limma” to perform variation analysis with the integrated gene expression. P < 0.05 were the essential screening criteria for DEPGs. The DEPGs were then visualized by using the “ggplot2” package with parameters of P < 0.05 for the volcano map and the “pheatmap” package with cut P = 0.05 for the heatmap.

Functional enrichment analysis

To further explore the DEPGs, the R package “clusterProfiler” was used for gene function enrichment analysis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were used to reveal the functions of the DEPGs and their related pathways. GO enrichment analysis was used to explore the biological processes (BP), molecular functions (MF) and cellular components (CC) of DEPGs. KEGG pathway analysis was used to analyze the signaling pathways of the DEPGs.

Construction of a PG risk score model

To further explore the prognosis-related genes in the DEPGs, we used the R package “survival” to perform univariate Cox regression analysis on the TCGA training set (OS ≥ 60 days) to identify prognosis-related genes. In addition, P < 0.05 was chosen as the filter to select significantly associated genes for subsequent model construction. Next, the mRNAs identified in the Cox regression results were subjected to LASSO regression dimensionality reduction, mainly relying on the R package “glmnet”. To construct a further accurate regression model, lambda screening with cross-validation was conducted.

Finally, we used multivariate Cox regression to identify the most relevant hub genes for LGG prognosis in DEPGs and constructed a predictive prognostic genetic risk score model. The sum of the products of the expression of the individual genes and the corresponding coefficients was used to calculate the risk score. Kaplan–Meier (KM) curves were constructed using the R packages “tidyverse” and “survminer” to determine the difference in OS between the high- and low-risk groups. ROC curves were constructed using the R package “timeROC” to assess the predictive accuracy of the genetic risk score model.

Analysis of genetic changes

Gene mutation analysis of nine HPGs was performed using the cBioPortal database (http://www.cbioportal.org). The cBioPortal database can be used to analyze the type and frequency of alterations in genes, including copy number amplifications, deep deletions and missense mutations of unknown significance.

Construction and validation of a prognostic nomogram for OS

Univariate and multivariate Cox analyses were performed on clinicopathological parameters (including age, grade and stage) and risk scores to assess whether the risk score model was independent of each parameter. Independent prognostic factors identified by multivariate Cox regression were combined to construct a nomogram using the R package “rms” and the R package “survivor” to predict 3-, 5- and 7-year survival. Moreover, calibration curves were used to validate their predictive values.

Immune cell correlation analysis

Immune cell infiltration tables for TCGA tumors were downloaded through TIMER2.0 (http://timer.comp-genomics.org/). Risk scores and immune cell correlations were analyzed using the R packages “limma”, “scales”, “ggplot2”, “ggtext”, and “ggpubr”. Immunological differences between the high- and low-risk groups were compared based on the results of the analysis.

External validation of the risk score model

To validate the predictive prognostic value of the risk score model, we used the GSE43378 dataset from the GEO database as an external validation group. Risk scores were calculated using the formula developed using the training group. KM curves were plotted to compare the OS between the two groups, and ROC curves were plotted to assess the stability of the model.

Establishment of the protein–protein interaction (PPI) network

To understand the interactions between DEPGs, we constructed a PPI network with the open-source software platform STRING (https://string-db.org/) to build the network model. Gene scores >0.4 and P < 0.05 were the essential screening criteria for the network model.

Tumor and paracancerous tissue collection from LGG patients

Our study was approved by the Clinical Research Ethics Committee of The First Affiliated Hospital of Wenzhou Medical University (permission: 2016-076), and written informed consent was obtained from all patients. The study protocol was carried out according to the principles of the Helsinki Declaration. Ten LGG patients (Grade II, n = 5; Grade III, n = 5) who had not received preoperative radiotherapy or chemotherapy underwent surgery between 2019 and 2021 at the First Affiliated Hospital of Wenzhou Medical University. The samples in this study were from five males and five females, aged between 26 and 62 years (mean, 41.5 years). The clinical characteristics of the patients are shown in Table 1.

Detection of the gene expression levels of nine HPGs in tumors and paracancerous tissues

Total RNA was extracted from the tumor tissue and paracancerous tissue using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, United States) according to the manufacturer’s instructions. Single-stranded cDNA was synthesized from 1 mg of total RNA using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara Biotechnology Co., Ltd., Dalian, China). Each cDNA sample (2 µl) was amplified in SYBR Green Real-time PCR Master Mix (final volume, 20 µL, YEASEN, China, #11203ES03), and the mRNA expression of the hub genes in the amplified samples was detected by a 7500 PCR system (Thermo Fisher Scientific). Thermal cycling conditions were as follows: 95 °C for 5 min, followed by 40 cycles of 95 °C for 20 s and 60 °C for 30 s. qPCR assays were conducted in triplicate in a 10 mL reaction volume for each sample. The relative expression of mRNA was calculated by the 2- ΔCt method. The primers used are shown in Table 2.

Statistical analysis

All analyses in this study were performed using R version 4.1.1 (R Core Team, 2021) and GraphPad Prism (version 8.3; GraphPad Software, La Jolla, CA, USA). The data between the two groups were compared with Student’s t test. Unless otherwise specified, P < 0.05 was considered significant.

Clinical characteristics of the patient sample

The analysis process of our study is shown in Fig. 1. A total of 476 LGG patients in the TCGA dataset met the criteria for the training group. Fifty LGG patients in the GSE43378 dataset were used for the validation group. The detailed clinical characteristics of the training group are shown in Table 3.

Identification of DEPGs in LGG

We screened 45 DEPGs between tumor samples (TCGA) and normal samples (GTEx) by comparing the expression levels of 51 PGs. Among them, 25 genes (CHMP6, NOD2, NLRC4, CHMP4B, CASP3, BAX, TIRAP, CASP9, BAK1, NLRP6, IL18, NLRP3, CASP5, GSDMC, TNF, IL1A, CASP1, AIM2, PYCARD, CHMP7, CHMP2B, CASP6, IL1B, TP53 and IRF2) were upregulated, while the remaining 20 genes (CASP4, NLRP1, TP63, CHMP4C, SCAF11, CYCS, HMGB1, NOD1, PRKACA, CHMP2A, GPX4, CASP8, NLRP2, ELANE, GSDMD, PLCG1, GZMB, IL6, IRF1 and GSDMB) were downregulated, as illustrated by heatmap (Fig. 2A) and volcano map (Fig. 2B). In order to further explore the signal pathway and interrelationship of 45 DEPGs, KEGG analysis showed that 45 DEPGs were mainly enriched in the NOD—like receptor signaling pathway, Salmonella infection and Necroptosis (Fig. 2C).

In summary, 45 differentially expressed pyroptosis genes were selected for subsequent analysis.

Construction of a risk score model to predict patient prognosis

To further screen for DEPGs associated with LGG prognosis, we analyzed the relationship between the expression levels of 45 DEPGs and the OS of LGG patients in the training group by univariate Cox regression. The results showed that the expression levels of 29 prognostic genes were significantly associated with OS in 476 LGG patients. Among them, five prognostic genes (HR <1) were protective genes, while the remaining 24 prognostic genes (HR > 1) were associated with poorer outcome in LGG patients (P < 0.05) (Fig. S1A). Subsequently, we reduced the number of candidate genes to 13 by LASSO analysis to obtain a better fitting model (Figs. S1B and S1C). Finally, the results of multivariate Cox regression analysis showed that nine hub genes associated with PGs, namely, CASP4, TP63, TIRAP, CASP9, IL18, TNF, IL1A, PLCG1, and TP53, could be used to construct a risk score model. The prognosis-related risk score was calculated as follows: risk score = (0.476 × Expression of CASP4) + (0.509 × Expression of TP63) + (−0.399 × Expression of TIRAP) + (−0.618 × Expression of CASP9) + (0.279 × Expression of IL18) + (−0.341 × Expression of TNF) + (0. 563 × Expression of IL1A) + (1.093 × Expression of PLCG1) + (0. 356 × Expression of TP53). As shown in Fig. 3A, the 476 LGG patients were divided into high-risk and low-risk groups based on the median risk score obtained from the formula. Survival analysis showed that patients in the low-risk group had a longer survival time and fewer deaths than those in the high-risk group (Fig. 3B). The expression of the nine HPGs in the high-risk and low-risk groups is shown in Fig. 3C. The expression levels of TIRAP, CASP9 and TNF were increased in the low-risk group compared to the high-risk group. In contrast, the expression levels of the remaining six HPGs were decreased in the low-risk group. The results of survival analysis showed that patients with high expression in TIRAP, CASP9 and TNF are significantly higher than the OS of patients with low expression (Figs. 3D–3F). The remaining genes show the opposite results (Figs. 3G–3L).

These results showed that we could predict the prognosis of LGG patients by calculating their risk score.

HPG expression level analysis

To further explore the expression levels of the nine HPGs, we analyzed data from 529 LGG patients in the TCGA dataset as well as 1,152 normal brain tissues in the GTEx. Comparing the expression levels of the nine HPGs between tumor samples in the TCGA dataset and normal samples in the GTEx dataset, we found that although CASP4, TP63, and PLCG1 were considered risk genes, the expression of these genes was elevated in normal brain tissue compared to the tumor group (Figs. 4A–4C) (P < 0.001). In contrast, the expression of TIRAP, CASP9, IL18, TNF, IL1A, and TP53 was reduced in normal brain tissues (Figs. 4D–4I) (P < 0.001). This conclusion was also confirmed by RT-qPCR results, as shown in Figs. 4J–4R. By analyzing the differential expression of RNA levels of these nine hub genes in the tumors of LGG patients and paired paracancerous tissues, we also found that compared with the tumor tissues, the expression of paracancerous tissues in LGG patients was significantly increased by 940.93 ± 235.41% at CASP4 (p < 0.001), 76.05 ±  21.42% at PLCG1 (p < 0.001), and 196.50 ± 29.66% at TP63 (p < 0.001) (Figs. 4J–4L). However, the expression of paracancerous tissues in LGG patients was significantly reduced by 271.33 ± 52.07% at CASP9 (p < 0.001), 108.95 ± 20.35% at IL1A (p < 0.001), and 1,292.90 ± 527.72% at IL18 (p < 0.001), 2,496.74 ±  519.85% at TIRAP (p < 0.001), 9,446.02 ± 1,240.85% at TNF (p < 0.001), 6,240.67 ± 1387.22% at TP53 (p < 0.001), compared with the tumor tissues (Figs. 4M–4R).

The above results showed that there was a significant difference in the expression of nine HPGS between the paracancerous tissues and tumor tissues, suggesting that these genes might be mutated in tumor tissues.

Analysis of the structural changes in HPGs

To further explore gene mutations, we analyzed nine HPGs in 283 LGG patients from the cBioPortal database. The results showed that 156 of 283 LGG patients (55%) exhibited significant changes in genetic sequence, with TP53 (52%) being the most frequently altered gene, followed by TIRAP (2.8%), TP63 (2.5%), PLCG1 (1.8%), CASP4 (1.4%), IL18 (1.4%), TNF (0.7%), CASP9 (0.4%), and IL1A (0.4%) (Fig. 5A). Interestingly, TP53 is the most commonly mutated gene in almost all tumors, especially in astrocytoma (Ohgaki & Kleihues, 2007; Synoradzki et al., 2021). In addition, we further analyzed the types of changes in the genes. As shown in Fig. 5B, 136 patients (48.06%) had mutations, 15 patients (5.3%) had multiple alterations, three patients (1.06%) had amplifications, and two patients (0.71%) had large-scale deletions. The domains altered within the hub genes also varied. In the TP63 gene, the R343L mutation appeared in the P53 DNA-binding domain (P53) structural domain, and the R487C and G635Vfs*69 mutations appeared outside the structural domain (Fig. 5C). In the IL18 gene, the L180F mutation appeared in the interleukin-1/18 (IL1/18) structural domain, and the T22M mutation appeared outside the conformational domain (Fig. 5D). In the PLCG1 gene, the E1163K and E1163_E1164dup mutations appeared in the C2 domain (C2), and the Y754C mutation appeared outside the conformational domain (Fig. 5E). In the TNF gene, there was only one mutation, a G224RE mutation in the structural domain of TNF (Fig. 5F). In addition, there were a large number of mutations in the TP53 gene. The most common mutations were R273C, R273H and R273L of the P53 DNA binding domain (P53), and mutations in the P53 tetramerization motif (P53_tetramer) and a small number of mutations outside the structural domain were also observed. In addition, the CNV status of nine HPGs has been shown in Fig. S2.

These results demonstrated that nine HPGs are mutated in LGG patients, and among these HPGs, TP53 presents the highest mutation rate.

Validation of the risk score model for survival prediction

On the basis of the median risk score, LGG patients in GSE43378 (n = 50) were divided into high-risk and low-risk groups. Meanwhile, to further confirm the validity of the model for predicting OS in LGG patients, we plotted the Kaplan −Meier curve and ROC curve associated with the model. The results of survival analysis showed that the OS of the high-risk group was significantly lower than that of the low-risk group in both the training and validation groups (Figs. 6A, 6E). In addition, ROC curve analysis was used to assess the accuracy of the model predictions. The results showed that the AUC values in the training group for 1-, 2- and 3-year survival were 0.871, 0.856 and 0.858, respectively (Figs. 6B–6D). The AUC values in the validation group for 1-, 2- and 3-year survival were 0.653, 0.764 and 0.738, respectively (Figs. 6F–6H).

Therefore, the results indicated that the model has good specificity and sensitivity for predicting LGG prognosis.

Evaluation of the risk score model as an independent prognostic factor for LGG

To verify whether the genetic risk score obtained by the model could be used as an independent prognostic factor, we performed univariate and multivariate Cox analyses on clinical factors such as risk score, age, tumor grade, and gender in the training set. The results of the univariate Cox analysis showed that risk score (P < 0.001), age (P < 0.001), and tumor grade (P < 0.001) were significantly associated with OS in LGG patients (Fig. 7A). In addition, the results of the multivariate Cox analysis also showed that risk score (P < 0.001), age (P < 0.001) and tumor grade (P < 0.001) were independent prognostic factors in patients with LGG (Fig. 7B). Thus, the risk score model is a reliable prognostic indicator independent of other clinicopathological parameters. To further clarify the independence of the risk score model, we performed stratified analysis for each indicator. We divided the patients into two age groups using 65 years as the basis of division and divided the patients into male and female groups according to sex. We also divided the patients into G2 and G3 groups (WHO grade II and III) according to pathological classification. The results showed that the survival curve was not affected by age, sex, or tumor grade, which further confirmed that the risk score model could be used to predict the prognosis of LGG patients (Fig. 7C).

Based on the results obtained from the above analysis, we constructed a nomogram incorporating age, sex, tumor stage, and risk score to predict 3-, 5- and 7-year OS in LGG patients (Fig. 8A). In addition, 3-, 5- and 7-year survival were predicted with good accuracy using the nomogram calibration curves (Figs. 8B–8D). Thus, our model can be used to better predict the prognosis of LGG patients.

HPG functional enrichment analysis

To further clarify the correlations between the nine HPGs in the model and signal path, gene enrichment analysis was used to further understand the potential functions of the HPGs. We investigated the molecular processes involved in LGG progression in depth by GO enrichment analysis and KEGG pathway analysis. The GO enrichment analysis showed that the nine HPGs were mainly related to cellular response to biotic stimulus, intrinsic apoptotic signaling pathway and response to lipopolysaccharide among the BP category terms. Among CC terms, the HPGs were significantly enriched in ruffle membrane, leading edge membrane and ruffle. Among MF terms, the HPGs were mainly enriched in cytokine receptor binding, cytokine activity and receptor ligand activity involved in the apoptotic process. Together, these findings showed that the nine HPGs are mainly enriched in immune-related signaling pathways (Fig. 9A). In addition, KEGG analysis showed that nine HPGs were mainly enriched in the NF-kappa B signaling pathway and inflammatory bowel disease (Fig. 9B).

These findings indicated that the nine HPGs play important roles in immunity and participate in immune factor regulation.

Results of immune infiltration analysis

To further clarify the correlation between the nine HPGs in the model and immunity, as shown in Fig. 10, we further analyzed the association of the risk score model and immune cells. We found that the infiltration levels of B cells (Fig. 10A, Cor = 0.329 and P = 1.565e−13), CD8 ⁺ T cells (Fig. 10B, Cor = 0.407 and P = 1.851e−20), CD4 ⁺ T cells (Fig. 10C, Cor = 0.255 and P = 1.569e−08), Macrophages (Fig. 10D, Cor = 0.372 and P = 3.977e−17), Dendritic cells (Fig. 10E, Cor = 0.406 and P = 2.531e- 20) and Neutrophils (Fig. 10F, Cor = 0.412 and P = 6.28e−21) were positively correlated with risk values.

The results showed that the expression of immune cells in the high-risk group was increased relative to that in the low-risk group among LGG patients.