A Pyroptosis-Related Gene Signature for Predicting Survival in Glioblastoma

Background In this study, a prognostic model based on pyroptosis-related genes was established to predict overall survival (OS) in patients with glioblastoma (GBM). Methods The gene expression data and clinical information of GBM patients were obtained from The Cancer Genome Atlas (TCGA), and bioinformatics analysis of differentially expressed genes was performed. LASSO Cox regression model was used to construct a three-pyroptosis-related gene signature, and validation was performed using an experimental cohort. Results A total of three pyroptosis-related genes (CASP4, CASP9, and NOD2) were used to construct a survival prognostic model, and experimental validation was performed using an experimental cohort. Receiver operating characteristic (ROC) analysis was performed, and the area under the ROC curves (AUC) was 0.921, 0.840, and 0.905 at 1, 3, and 5 years, respectively. Functional analysis revealed that T-cell activation, regulation of T-cell activation, leukocyte cell-cell adhesion, and positive regulation of cell adhesion among other immune-related functions were enriched, and immune-related processes were different between the two risk groups. Conclusion In this study, a novel prognostic model based on three pyroptosis-related genes is constructed and used to predict the prognosis of GBM patients. The model can accurately and conveniently predict the 1-, 3-, and 5-year OS of GBM patients.


INTRODUCTION
Glioblastoma multiforme (GBM) is the most common and most aggressive human brain tumor. Increasing incidence of GBM is becoming evident in patients with advanced age (1). It is estimated that more than 15,000 people die each year from GBM in the USA (2). Formulating a treatment plan for GBM patients requires a multidisciplinary treatment approach (3). Despite several innovations in the treatment of GBM, it remains one of the most difficult and complex cancers to treat, and the median survival is approximately 1 year (4). The multifactorial etiology of GBM makes prognostic prediction challenging. Therefore, considering the limited treatment strategies for GBM, novel prognostic models should be developed, to accurately and conveniently predict the overall survival (OS) of GBM patients.
Pyroptosis is a new type of programmed cell death mediated by the gasdermin D protein and characterized by the release of inflammatory mediators (5). Recently, pyroptosis has become a research hotspot in cancer initiation and progression. Besides, it is reported to be closely related to gastric cancer, colorectal cancer, hepatocellular carcinoma, breast cancer, skin cancer, and malignant mesenchymal tumors (6)(7)(8)(9)(10)(11)(12). Recent studies on the relationship between pyroptosis and cancer are providing new research ideas for the prevention and treatment of cancer. However, whether pyroptosis-related genes are correlated with the prognosis of patients with GBM remains largely unknown.
In the present study, mRNA expression data and corresponding follow-up clinical information of GBM patients were obtained from the Cancer Genome Atlas (TCGA) database and GTEx web portal. A prognostic multigene signature was then constructed with pyroptosis-related differentially expressed genes (DEGs). Functional enrichment analysis was performed to explore the potential mechanisms associated with the identified DEGs.

Datasets
The corresponding clinical information and mRNAs-seq data for GBM were downloaded from the TCGA Data Portal (https:// portal.gdc.cancer.gov/repository). Normal tissues were matched TCGA adjacent tissue and GTEx normal tissue. A total of 33 pyroptosis-related genes from previous systematic reviews were extracted and are provided in Table S1 (13)(14)(15)(16)(17). A total of 62 GBM patients were recruited from The First Affiliated Hospital of Zhengzhou University as the validation cohort.
All validation cohort patients gave written informed consent, and ethical permission was obtained from The First Affiliated Hospital of Zhengzhou University.

Identification of Differentially Expressed Pyroptosis-Related Genes
Normalization of the read count values was performed using edger (R package). To identify the differentially expressed genes (DEGs) between tumor tissues and adjacent normal tissues, the "limma" R package was used with FDR <0.05 and |log2FC| ≥1 in the training cohort.

Construction of the Prognostic Signature
The prognostic values of DEGs were determined by univariate Cox analysis, and genes significantly related to OS of GBM patients were identified. To avoid the overfitting problem, least absolute shrinkage and selection operator (LASSO)-penalized Cox regression analysis was performed. The penalty parameter (l) adjustment was performed by tenfold crossvalidation based on minimum criteria ( Figure S1). The survival genes were selected and used to construct a prognostic model based on the multivariate Cox regression analysis results. Each GBM patient was assigned an individual risk score, and the risk score formula was defined based on the expression level of each gene and the regression coefficient derived from the multivariate Cox regression model. The prognostic signature as risk score = S n i=1 exp i * b i (where n, exp i , and b i represent the number of prognostic genes, the expression value, and the coefficient of gene i, respectively). Patients were divided into high-risk and low-risk groups based on the median values of the risk score. Principal Component Analysis (PCA) was carried out in R using the prcomp function from the stats package and visualized with the ggbiplot library. For the survival analysis of each gene, the "surv_cutpoint" function of the "survminer" R package was used to determine the optimal cutoff values of the risk scores. In addition, a time-dependent receiver operating characteristic (ROC) curve analysis was performed using the R package "survivalROC" to evaluate the discrimination ability of the gene signature.

Validation of the Prognostic Signature by Quantitative Real-Time PCR and Western Blot
To validate the prognostic model, tumor tissue and normal tissue were collected from 62 GBM patients from The First Affiliated Hospital of Zhengzhou University. Using the same risk score formula, the risk score of each patient in the validation cohort was calculated, and patients were classified into the high-or lowrisk groups based on the median-based cutoff values. Multivariate Cox analysis revealed that CASP9 was an independent prognostic factor. Ethical approval was obtained from the First Affiliated Hospital of Zhengzhou University Ethics Committee.

Quantitative Real-Time PCR
Total RNA was extracted from the target tissue samples and thoroughly ground in a mortar under liquid nitrogen. To lyse the cells, 1 ml of Trizol reagent (Life Technology, Grand Island, NY, USA) was added, and the sample was incubated for 15 min at room temperature on a shaker. To assess the mRNA expression level, the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, Lithuania) was used to synthesis the first-strand cDNA. Quantitative PCR was performed using Roche LightCycler ® 480 Real-Time PCR System with SYBR ® Green qPCR mix 2.0 kit. The primers used in this study were obtained from TsingKe biological technology (Nanjing, China), including CASP9 (forward 5′-CTGTCTACGGCACAGATGGAT-3′, reverse 5′-GGGACTCGTCTTCAGGGGAA-3′), b-actin (Forward: 5′-CGAGCACAGAGCCTCGCCTTTGCC-3′, Reverse: 5′-TGTCGACGACGAGCGCGGCGATAT-3′). The relative mRNA levels were calculated by the 2-DDCt method.

Western Blot
Western blot was performed to determine the protein expression level. Samples were isolated and lysed in RIPA buffer with protease inhibitors. Proteins (40 mg) were separated using SDS-PAGE with 10% acrylamide gels. Western blot analysis was performed using antibodies against mouse monoclonal antibody-anti-human CASP9 (9502) from Cell Signaling Technology, and mouse monoclonal antibody-anti-human bactin (sc-47778, Santa Cruz Biotechnology), followed by incubation with horseradish peroxidase (HRP)-coupled mouse secondary antibody (1:10,000, NA93, GE Healthcare). To confirm equal protein loading, the blots were reprobed with a b-actin antibody, and analysis of the data was performed using NIH ImageJ software.

Functional Enrichment Analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses of the DEGs were performed using the "clusterProfiler" R package and the Cytoscape plugin "ClueGO." The infiltrating score of 16 immune cells and the activity of 13 immune-related pathways were calculated using the single-sample gene set enrichment analysis (ssGSEA) method in the Gene Set Variation Analysis (GSVA) package of R software (18,19). Besides, correlation analysis between the expression levels of the three pyroptosisrelated DEGs and infiltrating immune cells was performed using ssGSEA.

Identification of Prognostic Pyroptosis-Related DEGs
A total of 130 normal tissue samples were downloaded from the TCGA and GTEx databases, while 169 GBM tumor tissue samples were obtained from the TCGA database. Tables 1, 2 show the clinicopathological data of the GBM patients included in the study. A total of 19 pyroptosis-related genes were found to be differentially expressed between the tumor and adjacent nontumor tissues. Differential gene expression in the two groups was represented in a heatmap ( Figure 1A, blue: low expression level; red: high expression level). The relationship between the genes is presented in Figure 1B. To build a prognostic model for GBM, the training dataset was used and LASSO Cox regression analysis was performed to identify stable markers from the survival-related genes. To find an optimal l, 10-fold crossvalidation based on the minimum criteria was employed, and the final model gave a minimum crossvalidation error. A total of three DEGs (CASP4, CASP9, and NOD2) were identified. Figure 2A illustrates the results of multivariate Cox regression.

Construction and Validation of the Prognostic Model
The risk score was calculated using the following formula: 0.01204*expression level of CASP4+0.22295* expression level of NOD2-0.08962*expression level of CASP9. A median cutoff value was applied to stratify patients into a high-risk group (n = 76) and a low-risk group (n = 78) ( Figure 2B). The PCA showed that patients in the two risk groups were distributed in different directions ( Figure 2C). As shown in Figure 2D, patients in the high-risk score group were closely associated with a high risk of mortality, while patients in the low-risk group had a higher probability of survival. Kaplan-Meier analysis showed that the probability of survival was significantly higher in the low-risk group compared with patients in the high-risk group ( Figure 2E, p < 0.05). The time-dependent ROC curve was carried out to estimate the performance of the risk prediction model. The AUC of the prognostic risk assessment model for the three pyroptosis-related genes was 0.921, 0.840, and 0.905 at 1, 3, and 5 years, respectively ( Figure 2F). The calibration curve was close to the ideal curve, indicating that the model had a good prognostic effect ( Figure 2G). In the validation cohort, the risk score of each patient was calculated using the same formula and cutoff value. The overall survival of patients in the high-risk group was significantly worse compared with the low-risk group (p < 0.05) ( Figure 3A). Besides, the AUC of the three-gene signature was 0.869, 0.904, and 0.808 at 1, 3, and, 5 years,  respectively ( Figure 3B). In the validation group, the association between CASP9 and patient outcome was analyzed. The KM survival curve showed that low CASP9 gene expression patients had a significantly worse prognosis ( Figure 3C).

Expression Level of CASP9 in GBM Tissues
Quantitative real-time PCR (qRT-PCR) and Western blot analysis were performed to detect CASP9 expression. The results showed that CASP9 expression in GBM was significantly higher than in adjacent tissues ( Figures 4A, B).
Besides, protein expression of the three pyroptosis-related DEGs was upregulated in GBM tissues compared with normal tissues ( Figure 4C).

Independent Prognostic Value of the Risk Model
The prognostic significance of different clinical features in patients with GBM was assessed. Univariate Cox regression analysis was used to identify the univariate predictors, and the risk score was found to have a prognostic significance in OS (HR = 2.392; 95% CI, 1.861-3.075; Table 3). Even after adjusting for confounding factors, the risk score was found to be a prognostic factor (HR = 1.726; 95% CI, 1.302-2.287; Table 3) in GBM patients. The three-gene signature was experimentally validated as an independent prognostic factor ( Table 4). Figure 5 shows a heatmap of the pyroptosis-related DEGs and different distributions of patients' age and survival status between the low-and high-risk subgroups.

DEG-Based Tumor Classification
The relationship between the expression of pyroptosis-related DEGs and GBM subtypes was explored. Consensus clustering analysis was used to classify GBM patients into subtypes. By increasing the clustering variable (k) from 2 to 10, the highest intragroup correlation and low intergroup correlation were observed when k = 3 ( Figure 6A). The results showed that the 169 GBM patients were divided into three clusters based on the 19 DEGs. The gene expression profile was matched to the clinical data and presented in a heatmap ( Figure 6B). The OS in the three clusters was compared, and significant differences were reported ( Figure 6C).

Functional Analyses and Infiltrating Immune Cells
KEGG and GO function enrichment analyses were performed. GO was carried out in three functional ontologies: biological process (BP), cellular component (CC), and molecular function (MF). The DEGs were most enriched in biological processes, and specifically in immune-related processes such as T-cell activation, positive regulation of cell adhesion, positive regulation of cytokine production, and leukocyte cell−cell adhesion (p < 0.05; Figure 7A). In addition, KEGG pathway analyses indicated that the DEGs were highly enriched in cytokine-cytokine receptor interaction, PI3K−Akt signaling pathway, chemokine signaling pathway, hematopoietic cell lineage, rheumatoid arthritis, Staphylococcus aureus infection, etc. (p < 0.05; Figure 7B). Similarly, BP was performed using the  ClueGO network of the KEGG pathways and GO (Figures 7C, D).
To further explore the relationship between immune cell infiltrations and risk scores in the prognostic risk model, ssGSEA in R package GSVA was used to quantify the level of immune cell function or pathways in the cancer samples. Interestingly, the scores of most immune cell types (activated dendritic cells, macrophages, neutrophils, T helper cells, Th2 cells, tumor-infiltrating lymphocyte, and T regulatory cell) were significantly different between the two   Figure 8A; p < 0.05). Moreover, APC coinhibition, APC costimulation, cytokine-cytokine receptor, checkpoint, cytolytic activity, human leukocyte antigen, human leukocyte antigen, inflammation-promoting, parainflammation, T-cell coinhibition, T-cell costimulation, and type III FN response were significantly different between the two groups ( Figure 8B; p < 0.05). This finding was consistent with results from GO and KEGG analysis. ssGSEA was used to analyze the correlation between the various infiltrating immune cell types in GBM and the expression levels of the three DEGs, and the results showed that they were closely related to infiltrating immune cells ( Figures 8C-E and Table S2).

DISCUSSION
GBM is a highly invasive tumor, characterized by infiltration, aggressiveness, and resistance to treatment. The incidence of GBM is 3.2 per 100,000 population (20), and relapse has been reported in GBM patients treated with surgery alone. Despite the continuous improvements in surgical resection, chemotherapy, and radiation therapies, GBM remains difficult to treat and the overall mortality rates remain high (21). Thus, there is an urgent need to develop new prognostic models (22). Pyroptosis is a new lytic, and proinflammatory type of programmed cell death that is mainly involved in proinflammatory events (23). It is mediated by cysteine aspartic acid-specific protein kinases 1, 4, 5, and 11, which mainly depends on the poreforming activity of the Gasdermin protein family. At present, pyroptosis is reported to participate in the occurrence and development of various diseases, especially its dual role in promoting and inhibiting tumor formation and tumor microenvironment (24,25). However, the prognostic value and mechanism of pyroptosisrelated genes in GBM remain to be investigated. Thus, in this study, the expression of 33 pyroptosis-related genes in GBM tissues and their associations with OS was systematically explored. The results showed that 32 pyroptosis-related genes were differentially expressed between GBM and adjacent normal tissues in the derivation cohort. Moreover, a novel pyroptosisrelated gene signature was constructed and validated and found to have good accuracy for predicting the survival of GBM  patients. A total of three pyroptosis-related genes (CASP4, CASP9, and NOD2) were identified and included in the prognostic model. Studies have demonstrated that dysregulation of CASP4 is implicated in the initiation and development of many human diseases, including malignancies (26,27). Abnormal expression of CASP4 in clear cell renal cell carcinomas has been reported as an important prognostic factor affecting prognosis (27). CASP9 is involved in many cellular processes, and its main role is an initiator caspase. Besides, it is an important therapy target for many diseases related to pyroptosis (28). A combination of CASP9 overexpression with radiation has shown promise in enhancing the efficacy of tumor treatment (29). This is consistent with the current study findings. The role of NOD2 in inducing innate and adaptive immunity is complicated and thought to play a decisive role in maintaining microbial tolerance at the intestinal barrier (30). More experiments and clinical studies are, however, needed to establish the role of NOD2 in GBM.
The three pyroptosis-related genes were used to build a prognostic model with a good model fit, and high prediction power (0.921, 0.840, and 0.905 at 1, 3, and 5 years, respectively). The risk score of each patient was calculated using the prognostic model and was found to be an independent predictor of OS for GBM patients. These findings can contribute to the development of accurate and sensitive diagnostic and prognostic biomarkers for GBM. The three pyroptosis-related genes were further experimentally validated by qRT-PCR and Western blotting, thus providing strong evidence for the above conclusion.
The GO and KEGG enrichment analysis revealed significant enrichment of genes in immune-related pathways. In addition, the scores of most immune cell types (activated dendritic cell, macrophages, neutrophils, T helper cells, Th2 cells, tumorinfiltrating lymphocyte, and T regulatory cell) and APC coinhibition, APC costimulation, cytokine-cytokine receptor, checkpoint, cytolytic activity, human leukocyte antigen, human leukocyte antigen, inflammation-promoting, parainflammation, T-cell coinhibition, T-cell costimulation, and type III FN response were found to be statistically different between the two groups.
In summary, a novel prognostic model based on three pyroptosis-related genes was constructed and used to predict the prognosis of GBM patients. This model proved to be significantly associated with OS and can provide insights into the prediction of GBM prognosis.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Independent Ethics Committee of The First Affiliated Hospital of Zhengzhou University. The patients/ participants provided their written informed consent to participate in this study.