# A Novel Prognostic Pyroptosis-Related Gene Signature Correlates to Oxidative Stress and Immune-Related Features in Gliomas

Authors: Shulian Zeng; Wenjun Li; Hao Ouyang; Ying Xie; Xiaoqiang Feng; Li Huang
Journal: Oxidative Medicine and Cellular Longevity
 (2023)

---

## Abstract

Gliomas are highly invasive and aggressive tumors having the highest incidence rate of brain cancer. Identifying effective prognostic and potential therapeutic targets is necessitated. The relationship of pyroptosis, a form of programmed cellular death, with gliomas remains elusive. We constructed and validated a prognostic model for gliomas using pyroptosis-related genes. Differentially expressed pyroptosis-related genes were screened using the “limma” package. Based on LASSO-Cox regression, nine significant genes including CASP1, CASP3, CASP6, IL32, MKI67, MYD88, PRTN3, NOS1, and VIM were employed to construct a prognostic model in the TCGA cohort; the results were validated in the CGGA cohort. According to the median risk score, the patients were classified into two risk groups, namely, high- and low-risk groups. Patients at high risk had worse prognoses relative to those at low risk evidenced by the Kaplan-Meier curve analysis. The two groups exhibited differences in immune cell infiltration and TMB scores, with high immune checkpoint levels, TMB scores, and immune cell infiltration levels in the high-risk group. KEGG and GO analyses suggested enrichment in immune-related pathways. Furthermore, we found that the genes in our signature strongly correlated with oxidative stress-related pathways and the subgroups exhibited different ssGSEA scores. Some small molecules targeted the genes in the model, and we verified their drug sensitivities between the risk groups. The scRNA-seq dataset, GSE138794, was processed using the “Seurat” package to assess the level of risk gene expression in specific cell types. Finally, the MYD88 level was lowered in the U87 glioma cell line using si-RNA constructs. Cellular proliferation was impaired, and fewer pyroptosis-related cytokines were released upon exposure to LPS. In summary, we built a pyroptosis-related gene model that accurately classified glioma patients into high- and low-risk groups. The findings suggest that the signature may be an effective prognostic predictive tool for gliomas.


---

## Body

## 1. Introduction

Malignant gliomas are the most prevalent tumors in the central nervous system (CNS) which resulted in 17,000 deaths in 2017 in America, thus making it a serious public health concern [1]. GBM is highly invasive, and the estimated 5-year survival is less than 5%. LGG has a better prognosis but a higher recurrence rate [2, 3]. Standard treatments for gliomas include surgical resection, radiotherapy, chemotherapy, and immunotherapy; however, due to tumor heterogeneity, some patients show no response to these treatments [4]. Recently, several molecular markers implicated in the regulation of cancer cell proliferation and death have been utilized for prognostic prediction and pathological diagnosis of gliomas, including the mutations in isocitrate dehydrogenase (IDH), codeletion of the short and long arms of chromosomes 1 and 19, respectively (1p/19q), and O6-methylguanine-DNA methyltransferase (MGMT) methylation [5]. Clinical trials for several targeted therapies involving these molecular markers have been completed or are underway; however, only a few have been successful. Thus, the identification of effective diagnostic targets that can classify patients into different groups is necessary for individualized treatment.Pyroptosis is a form of programmed cellular death comprising the canonical or noncanonical type based on the caspase-1 (CASP-1) activating pathways [6]. It differs from apoptosis, necroptosis, or ferroptosis, in the case of the formation of inflammasome when cells are challenged. Subsequently, caspases are activated, which in turn cleave gasdermins, thereby causing cells to rupture and release their contents as well as inflammatory cytokines like IL-1β and IL-18 [7]. Pyroptosis is vital in the development of anti-infection and inflammatory responses [8]. The role of pyroptosis in tumors remains ambiguous. It may facilitate tumor growth, and metastasis, or participate in killing the cancer cells [9]. High levels of gasdermin B (GSDMB) correlated with a poor prognosis in breast cancer, where its isoform 2 promotes the progression, invasion, and metastasis of MCF7 cells [10]. However, in the melanoma mouse model, GSDMB cleaved by cytotoxic particles of perforin and granzyme A, triggers target cell pyroptosis, and promotes tumor clearance [11]. In gliomas, the levels of IL-1β and IL-18 increased, which are associated with lowered survival rates. Berberine targets ERK/CASP1 signaling, thus decreasing the secretion of IL-1β and IL-18 and inhibiting glioma cells [12]. Therefore, targeting pyroptosis is a promising strategy for tumor therapy. However, the underlying function of pyroptosis in gliomas remains elusive. It is important for patients to discover the link between pyroptosis and gliomas.Owing to the rapid developments in high-throughput sequencing technologies, genomic data have facilitated the understanding of heterogeneities and molecular characteristics of gliomas, which have helped improve the prognoses of these patients. Herein, we sourced public datasets comprising bulk and sc-RNA sequencing data to develop a prognostic model and nomogram using pyroptosis-related genes in an effect to predict the overall survival (OS) in glioma via comprehensive bioinformatic analysis and in vitro experiments. The results may bear implications for diagnosis and provide new insights into individualized treatment strategies for gliomas.


## 2. Methods and Materials

### 2.1. Data Sourcing from Public Databases

The HTSeq-FPKM RNA-seq data of normal (n=5), GBM, and LGG (n=698) brain tissues were sourced (https://portal.gdc.cancer.gov/) along with the corresponding clinical data from TCGA [13]. Another dataset (mRNAseq_325 and mRNA_693, n=1018) was downloaded from CGGA (http://www.cgga.org.cn/) as the test group [14]. The “SVA” package in R was used to remove batch effects between the different matrices. We assessed GeneCards (https://www.genecards.org/) using the keyword, “pyroptosis,” and obtained 155 pyroptosis-related protein-coded genes (Supplementary Table S1) [15].


### 2.2. Construction of the Pyroptosis-Related Gene Signature

First, a univariate COX regression analysis was conducted to screen all the prognostic pyroptosis-related genes related to OS; those withp values < 0.01 were screened, following which the “limma” package was employed to select the pyroptosis-related differentially expressed genes (DEGs) between normal (n=5) and tumor (n=698) tissues from TCGA transcriptome data. DEGs with log2FC≥1 and falsediscoveryrateFDR<0.05 were used for further analysis. The overlapping genes between DEGs and prognostic pyroptosis-related genes were visualized using a Venn diagram (https://bioinformatics.psb.ugent.be/webtools/Venn/). Using the “glmnet” package, a least absolute shrinkage and selection operator- (LASSO-) COX regression for the overlapping genes was performed to minimize overfitting. Thus, a pyroptosis-related signature was obtained. The coefficient of each gene in the signature was computed to estimate the risk score of patients using the following formula:
(1)Riskscore=∑i=1nCoefβi∗ExpXi,where Coefβi is the coefficient of each pyroptosis-related gene obtained from analysis of LASSO-COX regression, βi, and the level of gene expression is denoted by ExpXi. Using the median risk score, two risk groups of patients were stratified, namely, low- and high-risk groups, in both TCGA and CGGA cohorts. Subsequently, principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were performed using the “Rtsne” package. Survival analysis of each patient based on the genes using “survival” and “surviminer” packages was performed. The “timeROC” package was utilized for evaluating the predictive power of the constructed model. A predictive nomogram was used to assess the OS using the “rms” package.


### 2.3. Construction of the Protein-Protein Interaction (PPI) Network

STRING (https://cn.string-db.org/) comprises predicted and known PPIs [16]. These interactions comprise direct (physical) and indirect (functional) associations based on knowledge transfer between organisms, interactions aggregated from primary databases, and computational predictions. The pyroptosis-related DEGs were used to acquire the relationship between every two proteins. The cytoHubba plug-in in Cytoscape v3.8.2 was downloaded, and its algorithm of maximal clique centrality (MCC) was widely applied to search top genes in the constructed PPI network [17, 18]. In this study, the top 5 hubba nodes ranked by MCC were screened out for further analysis.


### 2.4. Functional Enrichment and Analysis of Tumor Mutational Burden (TMB)

The “clusterProfiler” package was employed for Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses for the above-mentioned significant DEGs between the risk groups. Copy number variation data for glioma patients were also obtained from TCGA and analyzed using the “maftools” package. Spearman’s analysis was used to detect the correlation between TMB scores and risk scores.


### 2.5. Single-Sample Gene Set Enrichment Analysis (ssGSEA) for Immune-Related Features and Oxidative Stress-Related Pathways

ssGSEA is an extension of the GSEA method that calculates the enrichment score for each sample in the context of specific gene sets [19]. Thirteen immune-related pathways, sixteen immune cell types, and fourteen oxidative stress-related pathways were assessed by ssGSEA function in the “gsva” package.


### 2.6. Analysis of Gene Expression and Drug Sensitivity

CellMiner, a query tool, and database facilitate the study and integration of pharmacological and molecular data from NCI-60 cancerous cell lines [20]. The transcriptome data were downloaded. Spearman’s analysis was performed to determine the relationship between drug sensitivity and gene expression with cor>0.3 and p value < 0.01 as the criteria of statistical significance.


### 2.7. Validation of Model Genes by HPA and scRNA-seq Dataset

To verify the protein expression of these 9 genes, the Human Protein Atlas (HPA) (https://www.proteinatlas.org/) database was queried to obtain immunohistochemistry (IHC) images of the cerebral cortex and glioma tissues [21]. GSE138794 (https://www.ncbi.nlm.nih) was sourced from the GEO database [22]. The “Seurat” package was utilized for processing the scRNA-seq dataset. Each gene was expressed in at least 10 cells with each cell exhibiting the expression of a minimum of 200 genes. Cells having >10% mitochondrial-expressed genes were excluded. The fastMNN algorithm was used to correct batch effects, and the top 2000 highly variable genes were identified using the “FindVariableFeatures” function. t-SNE was executed for reducing the dimension of the single cells. The function of “FindAllMarkers” in the “Seurat” package with logFC>0.5 and adjusted p value < 0.05 as threshold criteria was employed to screen marker genes. Subsequently, we sourced the Human Primary Cell Atlas data [23, 24] and applied the “SingleR” package to automatically annotate the identified cell clusters, and t-SNE plots were drawn for their visualization.


### 2.8. Cell Culture and Transfection

The human normal microglial cell line, HMC3, and glioma cell line, U87 (BNCC, Wuhan, China), were grown at 37°C and 5% CO2 in DMEM containing 10% fetal bovine serum (FBS, Gibco). siRNA constructs for targeting MYD88 and si-NC were synthesized by GenePharma (Shanghai, China). The sequences of si-MYD88 and si-NC were 5′-CCGGCAACUGGAGACACAATT-3′ and 5′-UUCUCCGAACGUGUCACGUTT-3′, respectively. Lipofectamine 2000 (Invitrogen) was used to transfect the siRNA constructs into the U87 cells. And a blank control (BC) that only contained transfection components was performed, too. Cells were harvested for subsequent experiments 48-hour posttransfection.


### 2.9. Real-Time Quantitative Polymerase Chain Reaction (RT-qPCR)

The EZBioscience kit (Roseville, USA) was used for total RNA extraction, and its concentration and purity were assessed using a NanoDrop (Thermo Fisher). Complementary DNA (cDNA) was synthesized following the protocol provided by the manufacturer (HiScript III RT SuperMix for qPCR (+gDNA wiper), Vazyme, Nanjing, China). qPCR was conducted using the ChamQ Universal SYBR qPCR kit (Vazyme, Nanjing, China) on the 7500 Real-Time PCR System (Thermo Fisher, USA). The relative mRNA levels were computed by the 2-ΔΔCT methodology with GAPDH as an internal reference. Sangon Biotech (Shanghai, China) synthesized the primers for MYD88 and GAPDH. The MYD88 primer sequences were as follows: 5′-CTGTGTCCGCACGTTCAAGA-3′ (reverse) and 5′-GGCTGCTCTCAACATGCGA-3′ (forward). The GAPDH primer sequences were as follows: 5′-GAAGATGGTGATGGGATT TC-3′ (reverse) and 5′-GAAGGTGAAGGTCGGAGTC-3′ (forward).


### 2.10. Cell Counting Kit-8 (CCK-8) Assay

Cellular proliferation after transfection was measured by the CCK-8 assay (GLPBIO, Montclair, CA, USA). Briefly, at a density of 1000 cells/well, cells were grown in 96-well plates containing DMEM with 10% FBS, following which 10μl CCK-8 solution was added per well at the indicated time points. These cells were incubated for 2 h, and the absorbance values at 450 nm were detected on a microplate spectrophotometer (Thermo Fisher, USA) to determine cellular proliferation.


### 2.11. Wound Healing Assay

In a 6-well plate, U87 cells were allowed to attain 90% confluency, following which, a sterile 200μl pipette tip was used to form a wound on the cells. Deciduous cells were washed off with phosphate-buffered saline. The cells were incubated in DMEM without serum. Images were captured at 100x (Leica, Germany) at indicated time points.


### 2.12. Enzyme-Linked Immunosorbent Assay (ELISA)

LPS (10μg/ml) was used to induce the secretion of IL-18 and IL-1β in si-MYD88 or si-NC-treated U87 cells. Cell supernatant was collected and used for ELISA (Boster, Wuhan, China, Cat. No. EK0864-IL-18 and Cat. No. EK0932-IL-1β). Briefly, 10 μl cell supernatant was added to coated plates and incubated with ELISA detection antibodies following the kit protocols.


### 2.13. Statistical Analysis

All statistical analyses were performed on GraphPad Prism (8.1) and R software (4.0.1). Log-rank tests were conducted to analyze the significance of the Kaplan-Meier survival differences between groups. ROC curve analysis was used to predict OS using the “pROC” package. Hazard ratios and the coefficients of independent prognostic factors were estimated by univariate and multivariate COX regression analyses. The Wilcoxon test was performed to assess immune-related functions and immune cell infiltration levels between risk groups. Two or more groups comprising continuous variables were compared by Student’st-test or one-way ANOVA, respectively (in terms of risk scores and gene expression). Data are presented as mean±SD with p<0.05 representing a statistically significant result.


## 2.1. Data Sourcing from Public Databases

The HTSeq-FPKM RNA-seq data of normal (n=5), GBM, and LGG (n=698) brain tissues were sourced (https://portal.gdc.cancer.gov/) along with the corresponding clinical data from TCGA [13]. Another dataset (mRNAseq_325 and mRNA_693, n=1018) was downloaded from CGGA (http://www.cgga.org.cn/) as the test group [14]. The “SVA” package in R was used to remove batch effects between the different matrices. We assessed GeneCards (https://www.genecards.org/) using the keyword, “pyroptosis,” and obtained 155 pyroptosis-related protein-coded genes (Supplementary Table S1) [15].


## 2.2. Construction of the Pyroptosis-Related Gene Signature

First, a univariate COX regression analysis was conducted to screen all the prognostic pyroptosis-related genes related to OS; those withp values < 0.01 were screened, following which the “limma” package was employed to select the pyroptosis-related differentially expressed genes (DEGs) between normal (n=5) and tumor (n=698) tissues from TCGA transcriptome data. DEGs with log2FC≥1 and falsediscoveryrateFDR<0.05 were used for further analysis. The overlapping genes between DEGs and prognostic pyroptosis-related genes were visualized using a Venn diagram (https://bioinformatics.psb.ugent.be/webtools/Venn/). Using the “glmnet” package, a least absolute shrinkage and selection operator- (LASSO-) COX regression for the overlapping genes was performed to minimize overfitting. Thus, a pyroptosis-related signature was obtained. The coefficient of each gene in the signature was computed to estimate the risk score of patients using the following formula:
(1)Riskscore=∑i=1nCoefβi∗ExpXi,where Coefβi is the coefficient of each pyroptosis-related gene obtained from analysis of LASSO-COX regression, βi, and the level of gene expression is denoted by ExpXi. Using the median risk score, two risk groups of patients were stratified, namely, low- and high-risk groups, in both TCGA and CGGA cohorts. Subsequently, principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were performed using the “Rtsne” package. Survival analysis of each patient based on the genes using “survival” and “surviminer” packages was performed. The “timeROC” package was utilized for evaluating the predictive power of the constructed model. A predictive nomogram was used to assess the OS using the “rms” package.


## 2.3. Construction of the Protein-Protein Interaction (PPI) Network

STRING (https://cn.string-db.org/) comprises predicted and known PPIs [16]. These interactions comprise direct (physical) and indirect (functional) associations based on knowledge transfer between organisms, interactions aggregated from primary databases, and computational predictions. The pyroptosis-related DEGs were used to acquire the relationship between every two proteins. The cytoHubba plug-in in Cytoscape v3.8.2 was downloaded, and its algorithm of maximal clique centrality (MCC) was widely applied to search top genes in the constructed PPI network [17, 18]. In this study, the top 5 hubba nodes ranked by MCC were screened out for further analysis.


## 2.4. Functional Enrichment and Analysis of Tumor Mutational Burden (TMB)

The “clusterProfiler” package was employed for Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses for the above-mentioned significant DEGs between the risk groups. Copy number variation data for glioma patients were also obtained from TCGA and analyzed using the “maftools” package. Spearman’s analysis was used to detect the correlation between TMB scores and risk scores.


## 2.5. Single-Sample Gene Set Enrichment Analysis (ssGSEA) for Immune-Related Features and Oxidative Stress-Related Pathways

ssGSEA is an extension of the GSEA method that calculates the enrichment score for each sample in the context of specific gene sets [19]. Thirteen immune-related pathways, sixteen immune cell types, and fourteen oxidative stress-related pathways were assessed by ssGSEA function in the “gsva” package.


## 2.6. Analysis of Gene Expression and Drug Sensitivity

CellMiner, a query tool, and database facilitate the study and integration of pharmacological and molecular data from NCI-60 cancerous cell lines [20]. The transcriptome data were downloaded. Spearman’s analysis was performed to determine the relationship between drug sensitivity and gene expression with cor>0.3 and p value < 0.01 as the criteria of statistical significance.


## 2.7. Validation of Model Genes by HPA and scRNA-seq Dataset

To verify the protein expression of these 9 genes, the Human Protein Atlas (HPA) (https://www.proteinatlas.org/) database was queried to obtain immunohistochemistry (IHC) images of the cerebral cortex and glioma tissues [21]. GSE138794 (https://www.ncbi.nlm.nih) was sourced from the GEO database [22]. The “Seurat” package was utilized for processing the scRNA-seq dataset. Each gene was expressed in at least 10 cells with each cell exhibiting the expression of a minimum of 200 genes. Cells having >10% mitochondrial-expressed genes were excluded. The fastMNN algorithm was used to correct batch effects, and the top 2000 highly variable genes were identified using the “FindVariableFeatures” function. t-SNE was executed for reducing the dimension of the single cells. The function of “FindAllMarkers” in the “Seurat” package with logFC>0.5 and adjusted p value < 0.05 as threshold criteria was employed to screen marker genes. Subsequently, we sourced the Human Primary Cell Atlas data [23, 24] and applied the “SingleR” package to automatically annotate the identified cell clusters, and t-SNE plots were drawn for their visualization.


## 2.8. Cell Culture and Transfection

The human normal microglial cell line, HMC3, and glioma cell line, U87 (BNCC, Wuhan, China), were grown at 37°C and 5% CO2 in DMEM containing 10% fetal bovine serum (FBS, Gibco). siRNA constructs for targeting MYD88 and si-NC were synthesized by GenePharma (Shanghai, China). The sequences of si-MYD88 and si-NC were 5′-CCGGCAACUGGAGACACAATT-3′ and 5′-UUCUCCGAACGUGUCACGUTT-3′, respectively. Lipofectamine 2000 (Invitrogen) was used to transfect the siRNA constructs into the U87 cells. And a blank control (BC) that only contained transfection components was performed, too. Cells were harvested for subsequent experiments 48-hour posttransfection.


## 2.9. Real-Time Quantitative Polymerase Chain Reaction (RT-qPCR)

The EZBioscience kit (Roseville, USA) was used for total RNA extraction, and its concentration and purity were assessed using a NanoDrop (Thermo Fisher). Complementary DNA (cDNA) was synthesized following the protocol provided by the manufacturer (HiScript III RT SuperMix for qPCR (+gDNA wiper), Vazyme, Nanjing, China). qPCR was conducted using the ChamQ Universal SYBR qPCR kit (Vazyme, Nanjing, China) on the 7500 Real-Time PCR System (Thermo Fisher, USA). The relative mRNA levels were computed by the 2-ΔΔCT methodology with GAPDH as an internal reference. Sangon Biotech (Shanghai, China) synthesized the primers for MYD88 and GAPDH. The MYD88 primer sequences were as follows: 5′-CTGTGTCCGCACGTTCAAGA-3′ (reverse) and 5′-GGCTGCTCTCAACATGCGA-3′ (forward). The GAPDH primer sequences were as follows: 5′-GAAGATGGTGATGGGATT TC-3′ (reverse) and 5′-GAAGGTGAAGGTCGGAGTC-3′ (forward).


## 2.10. Cell Counting Kit-8 (CCK-8) Assay

Cellular proliferation after transfection was measured by the CCK-8 assay (GLPBIO, Montclair, CA, USA). Briefly, at a density of 1000 cells/well, cells were grown in 96-well plates containing DMEM with 10% FBS, following which 10μl CCK-8 solution was added per well at the indicated time points. These cells were incubated for 2 h, and the absorbance values at 450 nm were detected on a microplate spectrophotometer (Thermo Fisher, USA) to determine cellular proliferation.


## 2.11. Wound Healing Assay

In a 6-well plate, U87 cells were allowed to attain 90% confluency, following which, a sterile 200μl pipette tip was used to form a wound on the cells. Deciduous cells were washed off with phosphate-buffered saline. The cells were incubated in DMEM without serum. Images were captured at 100x (Leica, Germany) at indicated time points.


## 2.12. Enzyme-Linked Immunosorbent Assay (ELISA)

LPS (10μg/ml) was used to induce the secretion of IL-18 and IL-1β in si-MYD88 or si-NC-treated U87 cells. Cell supernatant was collected and used for ELISA (Boster, Wuhan, China, Cat. No. EK0864-IL-18 and Cat. No. EK0932-IL-1β). Briefly, 10 μl cell supernatant was added to coated plates and incubated with ELISA detection antibodies following the kit protocols.


## 2.13. Statistical Analysis

All statistical analyses were performed on GraphPad Prism (8.1) and R software (4.0.1). Log-rank tests were conducted to analyze the significance of the Kaplan-Meier survival differences between groups. ROC curve analysis was used to predict OS using the “pROC” package. Hazard ratios and the coefficients of independent prognostic factors were estimated by univariate and multivariate COX regression analyses. The Wilcoxon test was performed to assess immune-related functions and immune cell infiltration levels between risk groups. Two or more groups comprising continuous variables were compared by Student’st-test or one-way ANOVA, respectively (in terms of risk scores and gene expression). Data are presented as mean±SD with p<0.05 representing a statistically significant result.


## 3. Results

### 3.1. Identification of the Candidate Prognostic Pyroptosis-Related DEGs in the TCGA Cohort

We screened 25 pyroptosis-related DEGs between tumor and normal tissues and 128 genes that were correlated with OS by univariate COX regression analysis. Finally, 21 overlapping pyroptosis-related prognostic genes were identified (Figure1(a)). Figure 1(b) depicts the findings of the univariate COX regression (all p<0.001). APOE, MST1, and NOS1 were protective factors, while the remaining were unfavorable prognostic factors. CASP1, CASP3, EGFR, TP53, and MYD88 were the top 5 hubs in the PPI network (Figure 1(c)). Figure 1(d) shows the relationship among the genes in another way. For example, there is a negative correlation between TP53 and APOE, while MST1 is positively correlated with APOE. It implied the distinct interacting patterns of complicated pyroptosis-related genes.Figure 1
Identification of the candidate prognostic pyroptosis-related DEGs in the TCGA cohort. (a) The Venn diagram for identifying prognostic pyroptosis-related differentially expressed genes between tumor and adjacent normal tissues. (b) Forest plot based on univariate Cox regression showing the hazard ratio andP values of candidate pyroptosis-related genes (all P<0.001); 3 and 18 genes correlate with good and poor prognoses, respectively. (c) Protein-protein interaction (PPI) network constructed using STRING. (d) The correlation network shows the relationship among candidate prognostic pyroptosis-related genes. Red line: positive correlation; blue line: negative correlation. Deeper colors indicate stronger relevance.
(a)(b)(c)(d)


### 3.2. Generation of a Pyroptosis-Related Gene Model in the TCGA Cohort

For identifying prognostic markers among the 21 pyroptosis-related prognostic genes, we performed LASSO regression analysis and obtained an optimalλ value (Figures 2(a) and 2(b)). Thus, 9 genes including CASP1, CASP3, CASP6, IL32, MKI67, MYD88, PRTN3, NOS1, and VIM were identified and used subsequently for constructing a prognostic model. Survival analysis as shown in Supplementary Figure S1 identified 8 genes that were significantly correlated with a poor prognosis, and only NOS1 with low expression in gliomas was associated with a good prognosis. We derived the risk scores as follows: 0.4654×MYD88+0.421×CASP6+0.2282×CASP3+0.1407×PRTN3+0.1217×MKI67+0.0859×IL32+0.0729×VIM+0.0673×CASP1–0.0375×NOS1. Patients were then stratified into high- (n=333 patients) and low-risk groups (n=334 patients) according to the median value of the risk score (Figure 2(c)). Discrete patient distribution was confirmed by t-SNE and PCA results (Figures 2(d) and 2(f)). The survival distribution for each patient was ranked from left to right according to their risk score (Figure 2(e)). The OS status of those in the high-risk group was remarkably poorer than patients in the low-risk group. Consistently, the Kaplan-Meier curve suggested that patients at high risk had a worse prognosis (Figure 2(g), p<0.001). We plotted time-dependent ROC curves to assess the predictive performance of the risk score for patients’ OS. Notably, the pyroptosis-gene signature showed potent prognostic validity to predict survival status, and the area under the curve (AUC) was 0.866, 0.895, and 0.847 for 1-, 3-, and 5-year OS, respectively (Figure 2(h)). Moreover, the expression of all genes in the model increased with the risk score except for NOS1 (Figure 2(i)). Thus, we successfully built a pyroptosis-related gene signature for the TCGA cohort.Figure 2
Construction of the pyroptosis-related gene signature in TCGA cohort. (a) Cross-validation for tuning parameters in the proportional hazards model. (b) The least absolute shrinkage and selection operator (LASSO) coefficients of 21 genes related to gliomas. (c) The Kaplan-Meier curves for overall survival (OS) in the high- and low-risk groups. (d) The receiver operating characteristic (ROC) curve analysis for verification of predictive performance of the risk model for 1-,3-, and 5-year OS. (e, f) Distribution of patients’ risk scores and survival status. (g, h) PCA andt-SNE plot. (i) The heat map showed the levels of gene expression in the high- and low-risk groups.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


### 3.3. Verification of the Prognostic Signature in the CGGA Cohort

The CGGA cohort including two datasets was used for verifying the robustness of the signature constructed in the TCGA cohort. The patients were also stratified into two risk groups based on the median risk score in the CGGA cohort (Figure3(a)). Similarly, t-SNE and PCA results for the CGGA cohort confirmed that patients were discretely distributed into two risk groups (Figures 3(c) and 3(d)). Likewise, relative to the patients in the low-risk group, those in the high-risk group exhibited a poorer prognosis with shorter OS (Figures 3(b) and 3(e)). The AUC for 1-, 3-, and 5-year OS was 0.686, 0.753, and 0.778, respectively (Figure 3(f)). With an increase in the risk score, the expression of all genes in the signature increased, except that of NOS1 in the CGGA cohort. Collectively, the results from the CGGA cohort were consistent with those of the TCGA cohort indicating that the 9-gene pyroptosis-related signature correlated positively with gliomas.Figure 3
Prognostic validation of the nine-gene signature in the CGGA cohort. (a, b) Distribution of patients’ risk scores and survival status. (c) The Kaplan-Meier curves for overall survival (OS) in the high- and low-risk groups. (d) ROC curve analysis for verification of predictive performance of the risk model for 1-, 3-, and 5-year OS. (e, f) PCA andt-SNE plot. (g) The heat map shows the levels of gene expression in the high- and low-risk groups.
(a)(b)(c)(d)(e)(f)(g)


### 3.4. Validation of Clinical Independence of the Pyroptosis-Related Prognostic Signature in the CGGA and TCGA Cohorts and Construction of a Predictive Nomogram for the Latter

The findings from univariate COX regression suggested a remarkable association between risk score and OS in both CGGA and TCGA cohorts (CGGA cohort:HR=1.676, 95% CI=1.557-1.804, p<0.001; TCGA cohort: HR=3.150, 95% CI=2.741-3.620, p<0.001) (Figures 4(a) and 4(c)). The risk score remained an independent predictive factor for OS after correcting for confounders, as evidenced by multivariate COX regression analysis results (CGGA cohort: HR=1.267, 95% CI=1.164-1.378, p<0.001; TCGA cohort: HR=2.054, 95% CI=1.654-2.550, p<0.001) (Figures 4(b) and 4(d)).Figure 4
Validation of the clinical independence of the signature in TCGA and CGGA cohorts. (a, b) Univariate and multivariate COX regression analyses in TCGA cohort. (c, d) Univariate and multivariate COX regression analyses in the CGGA cohort.
(a)(b)(c)(d)We developed a nomogram incorporating clinical characteristics of radio status, age, grade, and risk score in the TCGA cohort based on the findings of the multivariate COX regression analyses. The concordance index (C-index) was 0.862, suggesting its relatively high performance for clinical diagnosis (Figure5(a)), wherein a higher risk score predicted a worse prognosis. The nomogram had high accuracy for predicting 1-, 3-, and 5-year OS (AUC=0.901, 0.924, and 0.879, respectively) (Figures 5(b)–5(d)). With excellent performance, the calibration curves demonstrated relatively good fits for predicting OS (Figures 5(e)–5(g)). Taken together, the nomogram may serve as a useful quantitative tool for predicting prognoses of gliomas.Figure 5
Construction of the predictive nomogram for the TCGA cohort. (a) Significant factors obtained from multivariate COX regression analysis in the TCGA cohort were used to construct a nomogram for predicting 1-, 3-, and 5-year survival. (b–d) ROC curve analysis of the nomogram for predicting 1-, 3-, and 5-year survival. (e–g) Calibration curve analysis for the accuracy of nomogram in predicting 1-, 3-, and 5-year survival.∗∗∗p<0.001.
(a)(b)(c)(d)(e)(f)(g)


### 3.5. Functional Annotation and Immune-Related Features of the Pyroptosis-Related Gene Signature

After constructing the pyroptosis-related model, we assessed differential pathways and biological functions between the risk groups in CGGA and TCGA cohorts. A total of 548 DEGs were screened in the TCGA cohort (38 and 510 genes were down- and upregulated, respectively) and 377 DEGs in the CGGA cohort (39 and 338 genes were down- and upregulated, respectively). GO functional enrichment results for the two cohorts suggested that the DEGs were primarily enriched in immune-associated biological processes, cellular components, and molecular functions, including adaptive immune response based on antigen-binding, somatic recombination of immune receptors built from immunoglobulin superfamily domains, humoral immune response, and MHC II protein complex (Figures6(a) and 6(c)). Consistently, the enriched KEGG pathways comprised infectious immune pathways, including the Epstein-Barr virus infection, influenza A, tuberculosis, human T cell leukemia virus 1 infection, and Th17 differentiation (Figures 6(b) and 6(d)). Therefore, the immune microenvironment was remarkably associated with the pyroptosis-related signature.Figure 6
Functional enrichment analyses in TCGA and CGGA cohorts. (a, b) Results of GO and KEGG analyses in TCGA cohort. (c, d) Results of GO and KEGG analyses in the CGGA cohort.
(a)(b)(c)(d)As several immune-related pathways were enriched, thereby highlighting the significance of the immune environment in tumors, we examined the differences in the immune landscape between the risk groups. Immune-related functions, immune checkpoints, and infiltrated immune cell types were analyzed. The results of ssGSEA demonstrated that immune-related functional pathways were remarkably enriched in the high-risk group relative to the low-risk group; similar trends were observed for most immune cells experiencing infiltration and checkpoints (Figures7(a), 7(b), and 7(e)). Consistent results were obtained in the CGGA cohort (Figures 7(c), 7(d), and 7(f)). Collectively, the risk groups showed correlation and obvious differences in immune-related functional pathways, levels of infiltrated immune cells, and immune checkpoints.Figure 7
The risk groups differ in immune-related functions, immune cell infiltration levels, and immune checkpoint expression levels. The scores of immune-related functions and immune cell infiltration in TCGA (a, b) and CGGA (c, d) cohorts. Immune checkpoint expression in TCGA (e) and CGGA (f) cohorts.∗p<0.05, ∗∗p<0.01, and ∗∗∗p<0.001. ns: no significance.
(a)(b)(c)(d)(e)(f)


### 3.6. The Relationship between Oxidative Stress and Pyroptosis Signature

Recently, studies found that cell pyroptosis was associated with oxidative stress in triple-negative breast cancer and non-small-cell lung cancer [25, 26]. However, whether pyroptosis is related to oxidative stress in glioma remains unknown. We acquired fourteen oxidative stress-related pathways from MsigDB (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) by the keyword of “oxidative stress” [27, 28]. The GO terms ID and their corresponding standard names are listed in Supplementary Table S2. Correlation analysis showed that the nine genes in our model were correlated to most of the specific oxidative stress-related pathways. We found that the model genes were almost consistently negatively associated with the pathways of GOBP_OXIDATIVE_DEMETHYLATION, GOBP_REGULATION_OF_TRANSCRIPTION_FROM_RNA_POLYMERASE_II_PROMOTER_IN_RESPONSE_TO_OXIDATIVE_STRESS, GOBP_OXIDATIVE_RNA_DEMETHYLATION, and GOBP_OXIDATIVE_DNA_DEMETHYLATION. On the other hand, the genes, except for NOS1 which was mainly negatively related to the oxidative stress pathway, were strongly positively correlated with the pathways of GOBP_CELLULAR_RESPONSE_TO_OXIDATIVE_STRESS, GOBP_INTRINSIC_APOPTOTIC_SIGNALING_PATHWAY_IN_RESPONSE_TO_OXIDATIVE_STRESS, and GOBP_RESPONSE_TO_OXIDATIVE_STRESS (Figure 8(a)). Moreover, we compared the ssGSEA scores of oxidative stress-related pathways between the high- and low-risk groups. It displayed that the ssGSEA scores of eight GO terms (GOBP_INTRINSIC_APOPTOTIC_SIGNALING_PATHWAY_IN_RESPONSE_TO_OXIDATIVE_STRESS, GOBP_NEGATIVE_REGULATION_OF_CELLULAR_RESPONSE_TO_OXIDATIVE_STRESS, GOBP_NEGATIVE_REGULATION_OF_RESPONSE_TO_OXIDATIVE_STRESS, GOBP_OXIDATIVE_PHOSPHORYLATION, GOBP_REGULATION_OF_OXIDATIVE_PHOSPHORYLATION, GOBP_REGULATION_OF_OXIDATIVE_STRESS_INDUCED_CELL_DEATH, GOBP_REGULATION_OF_OXIDATIVE_STRESS_INDUCED_INTRINSIC_APOPTOTIC_SIGNALING_PATHWAY, and GOBP_RESPONSE_TO_OXIDATIVE_STRESS) were higher in the high-risk group, and four pathways (GOBP_OXIDATIVE_DEMETHYLATION, GOBP_OXIDATIVE_DNA_DEMETHYLATION, GOBP_OXIDATIVE_RNA_DEMETHYLATION, and GOBP_REGULATION_OF_TRANSCRIPTION_FROM_RNA_POLYMERASE_II_PROMOTER_IN_RESPONSE_TO_OXIDATIVE_STRESS) were enriched in the low-risk group (Figure 8(b)). These results were generally the same in the CGGA cohort (Figure S2). It suggested that our novel model of pyroptosis was closely related to oxidative stress pathways.Figure 8
The pyroptosis-related gene signature was correlated with oxidative stress pathways in the TCGA cohort. (a) The relationship between the model gene and oxidative stress pathways. (b) The differences in oxidative stress pathway scores between the high- and low-risk groups.∗p<0.05, ∗∗p<0.01, and ∗∗∗p<0.001. ns: no significance.
(a)(b)


### 3.7. Difference between TMB and Drug Sensitivities in the Two Subgroups

Tumor mutational burden is one of the characteristics and biomarkers for gliomas [29]. We analyzed the correlation of risk scores and TMB scores; the results showed that the risk score correlated positively with TMB and the high-risk group with higher TMB scores than another group (Figures 9(a) and 9(b)). And the top 20 genes with the highest mutational frequencies are shown in Figures 9(c) and 9(d). Most patients carried mutations in the low- and high-risk groups (87.77% and 95.95%, respectively), and missense mutations were the most frequent. TP53 is an important tumor inhibitor gene whose mutational frequency was the highest in the low-risk group (40%), followed by IDH1 (31%). The corresponding mutational frequencies were 43% and 89%, respectively, in the high-risk group.Figure 9
Alteration of TMB between the risk groups in the TCGA cohort. Correlation of TMB and risk score (a). Differences in TMB between the risk groups (b). Top 20 DEGs in glioma between the risk groups (c, d).
(a)(b)(c)(d)To examine the link between drug sensitivities and pyroptosis-related genes, the CellMiner database based on a suite of genomic and pharmacological information including more than 20,000 compounds tested in NCI-60 cell lines was queried. We assessed 75 compounds under clinical trials and 188 FDA-approved drugs. The top 16 most relevant correlations are shown in Figure10. The 9 genes were associated with chemotherapeutic drugs (all p<0.01 and cor>0.3). For instance, the expression of PRTN3 correlated positively with the IC50 value of imexon, ABT-199, cyclophosphamide, and nandrolone phenpr indicating that tumor cells were more resistant to elevated gene expression levels. Conversely, enhanced expression of IL32 was associated with heightened cancer cell-drug sensitivity to tanespimycin, dolastatin 10, tyrothricin, bafetinib, and vinblastine. The differences in the IC50 values of 16 drugs between the risk groups were assessed. Seven drugs (tanespimycin, dolastatin, tyrothricin, staurosporine, geldanamycin analog, bafetinib, and vinblastine) correlated significantly with the risk scores, and except for the staurosporine, all drugs had enhanced sensitivity in the low-risk group (Supplementary Figure S3). Together, our model may provide insights for drugs targeting the pyroptosis-related genes and may bear implications for developing personalized treatment strategies for glioma patients.Figure 10
The correlation between model gene expression and drugs. The top 16 most relevant correlations were visualized: (a) imexon, (b) ABT-199, (c) tanespimycin, (d) cyclophosphamide, (e) dolastatin, (f) nelarabine, (g) tyrothricin, (h) nandrolone phenpr, (i) SR16157, (j) staurosporine, (k) hydroxyurea, (l) geldanamycin analog, (m) bafetinib, (n) E-7820, (o) vinblastine, and (p) chelerythrine.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)


### 3.8. Validation of the Gene Signature in GBM Tissues by HPA and scRNA-Seq

We then verified the gene signature in GBM tissues using IHC images from the HPA database. The expression levels of CASP1, CASP3, CASP6, IL32, MKI67, MYD88, and VIM were higher in tumor tissues, while those of PRTN3 and NOS1 were lower (Figure11). Furthermore, the distribution of these genes in glioma tumors was assessed in GSE138794. Eight samples of gliomas were integrated using “fastMNN” and analyzed using the “Seurat” package. The top 2000 variate genes and features are shown in Supplementary Figure S4. Twenty-two clusters were clarified, and six cell types, including astrocytes (cluster 6, 5.88%), endothelial cells (cluster 18, 1.12%), macrophages (clusters 4 and 9, 12.21%), monocytes (cluster 8, 4.49%), oligodendrocytes (cluster 13, 3.31%), and malignants classified into 14 clusters with the highest proportion of up to 72.99% (Figures 12(a) and 12(b) and Supplementary Figure S5), were screened. The distribution of 8 genes in different cell types is shown in the feature plot (Figure 12(c)) and violin plot (Supplementary Figure S6). PRTN3 was not detected in this dataset. CASP1 was mainly present in clusters 4, 8, and 9, comprising macrophages and monocytes. CASP3 was located in clusters 2, 7, 15, and 17 which were malignant. CASP6 was detected in all clusters at low levels. IL32 was present in cluster 18 comprising endothelial cells. MKI67 showed a remarkably high expression level in cluster 7, indicating that the malignants in cluster 7 have enhanced growth abilities. VIM was distributed in all clusters at high levels but primarily located in malignants. MYD88 was primarily present in clusters 1 (malignants) and 4 (monocytes). NOS1 was rarely found, consistent with the results of IHC. Collectively, we confirmed the levels of protein expression of the gene signature in specific cells by scRNA. Thus, pyroptosis-related genes distributed across different cell types may influence glioma progression via the microenvironment.Figure 11
IHC results for the 9 genes were validated in the HPA database: (a) CASP1, (b) CASP3, (c) CASP6, (d) IL32, (e) MKI67, (f) MYD88, (g) NOS1, (h) PRTN3, and (i) VIM3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 12
Single-cell RNA-seq results for the model genes were validated in GSE138794. (a) Thet-SNE plot showed 22 clusters in GSE138794. (b) The t-SNE plot displayed the annotation of cells in the 22 clusters. (c) The scatter plot showed the expression patterns of eight model genes in different cell types.
(a)(b)(c)


### 3.9. Knocking Down MYD88 Impaired Cell Proliferation and Reduced the Release of Pyroptosis-Related Cytokines

MYD88 was selected to verify the reliability of our model in vitro because of the following reasons: (i) MYD88 was a risk factor upregulated significantly in gliomas relative to normal tissues and correlated with a poor prognosis, (ii) MYD88 was one of the top 5 hub genes in the PPI network, (iii) the coefficient of MYD88 was the maximum in our signature, and (iv) the involvement of MYD88 has been previously reported in LPS-induced pyroptosis [30]. We detected the expression of MYD88 in the normal cell line (HMC-3) versus the glioma cell line (U87) by qRT-PCR. MYD88 was upregulated in U87 cells (1.00±0.07 vs. 2.29±0.33, p=0.003, Figure 13(a)). si-RNA was used to successfully reduce the expression of MYD88 in the U87 cell line (1.00±0.16 vs. 0.56±0.09, p=0.014, Figure 13(a)). Cellular proliferation (72 h and 96 h) was inhibited following treatment with si-MYD88 relative to si-NC-transfected cells (Figure 13(b)). The wound healing assay showed that cellular migration reduced significantly in si-MYD88-treated cells (Figure 13(c)). MYD88 is known to be indispensable in LPS-induced cell pyroptosis via TLR4/MYD88/NF-κB/NLRP3 or TLR-4/MYD88/PI3K/AKT signaling [30–34]. Therefore, we treated si-MYD88- and si-NC-transfected U87 cells with LPS (10 μg/ml) for 24 hours and detected the levels of IL-1β and IL-18 to verify cellular pyroptosis. A simultaneous reduction in the levels of IL-1β (22.63±5.73 vs. 13.36±0.23, p=0.049) and IL-18 (92.10±6.07 vs. 51.86±16.25, p=0.016) was observed in the si-MYD88 group relative to the si-NC group (Figure 13(d)), indicating that downregulation of MYD88 in U87 cells attenuated LPS-induced cell pyroptosis. Thus, MYD88 may exert an oncogenic role, and its activation is a potential treatment target for glioma.Figure 13
Knockdown of MYD88 impaired cell proliferation and reduced release of pyroptosis-related cytokines. (a) The level of MYD88 expression in HMC3 and U87 cell lines (left), and the si-RNA construct effectively attenuated MYD88 in the U87 cell line. (b) Proliferation curves for the U87 cell line after treatment with si-RNA. (c) Wound healing between the si-NC and si-MYD88 groups at 0 h and 48 h. (d) Release of pyroptosis-related cytokines, IL-1β and IL-18, in the si-NC and si-MYD88 groups challenged with LPS. ∗p<0.05;∗∗p<0.01; ∗∗∗p<0.001.
(a)(b)(c)(d)


## 3.1. Identification of the Candidate Prognostic Pyroptosis-Related DEGs in the TCGA Cohort

We screened 25 pyroptosis-related DEGs between tumor and normal tissues and 128 genes that were correlated with OS by univariate COX regression analysis. Finally, 21 overlapping pyroptosis-related prognostic genes were identified (Figure1(a)). Figure 1(b) depicts the findings of the univariate COX regression (all p<0.001). APOE, MST1, and NOS1 were protective factors, while the remaining were unfavorable prognostic factors. CASP1, CASP3, EGFR, TP53, and MYD88 were the top 5 hubs in the PPI network (Figure 1(c)). Figure 1(d) shows the relationship among the genes in another way. For example, there is a negative correlation between TP53 and APOE, while MST1 is positively correlated with APOE. It implied the distinct interacting patterns of complicated pyroptosis-related genes.Figure 1
Identification of the candidate prognostic pyroptosis-related DEGs in the TCGA cohort. (a) The Venn diagram for identifying prognostic pyroptosis-related differentially expressed genes between tumor and adjacent normal tissues. (b) Forest plot based on univariate Cox regression showing the hazard ratio andP values of candidate pyroptosis-related genes (all P<0.001); 3 and 18 genes correlate with good and poor prognoses, respectively. (c) Protein-protein interaction (PPI) network constructed using STRING. (d) The correlation network shows the relationship among candidate prognostic pyroptosis-related genes. Red line: positive correlation; blue line: negative correlation. Deeper colors indicate stronger relevance.
(a)(b)(c)(d)


## 3.2. Generation of a Pyroptosis-Related Gene Model in the TCGA Cohort

For identifying prognostic markers among the 21 pyroptosis-related prognostic genes, we performed LASSO regression analysis and obtained an optimalλ value (Figures 2(a) and 2(b)). Thus, 9 genes including CASP1, CASP3, CASP6, IL32, MKI67, MYD88, PRTN3, NOS1, and VIM were identified and used subsequently for constructing a prognostic model. Survival analysis as shown in Supplementary Figure S1 identified 8 genes that were significantly correlated with a poor prognosis, and only NOS1 with low expression in gliomas was associated with a good prognosis. We derived the risk scores as follows: 0.4654×MYD88+0.421×CASP6+0.2282×CASP3+0.1407×PRTN3+0.1217×MKI67+0.0859×IL32+0.0729×VIM+0.0673×CASP1–0.0375×NOS1. Patients were then stratified into high- (n=333 patients) and low-risk groups (n=334 patients) according to the median value of the risk score (Figure 2(c)). Discrete patient distribution was confirmed by t-SNE and PCA results (Figures 2(d) and 2(f)). The survival distribution for each patient was ranked from left to right according to their risk score (Figure 2(e)). The OS status of those in the high-risk group was remarkably poorer than patients in the low-risk group. Consistently, the Kaplan-Meier curve suggested that patients at high risk had a worse prognosis (Figure 2(g), p<0.001). We plotted time-dependent ROC curves to assess the predictive performance of the risk score for patients’ OS. Notably, the pyroptosis-gene signature showed potent prognostic validity to predict survival status, and the area under the curve (AUC) was 0.866, 0.895, and 0.847 for 1-, 3-, and 5-year OS, respectively (Figure 2(h)). Moreover, the expression of all genes in the model increased with the risk score except for NOS1 (Figure 2(i)). Thus, we successfully built a pyroptosis-related gene signature for the TCGA cohort.Figure 2
Construction of the pyroptosis-related gene signature in TCGA cohort. (a) Cross-validation for tuning parameters in the proportional hazards model. (b) The least absolute shrinkage and selection operator (LASSO) coefficients of 21 genes related to gliomas. (c) The Kaplan-Meier curves for overall survival (OS) in the high- and low-risk groups. (d) The receiver operating characteristic (ROC) curve analysis for verification of predictive performance of the risk model for 1-,3-, and 5-year OS. (e, f) Distribution of patients’ risk scores and survival status. (g, h) PCA andt-SNE plot. (i) The heat map showed the levels of gene expression in the high- and low-risk groups.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 3.3. Verification of the Prognostic Signature in the CGGA Cohort

The CGGA cohort including two datasets was used for verifying the robustness of the signature constructed in the TCGA cohort. The patients were also stratified into two risk groups based on the median risk score in the CGGA cohort (Figure3(a)). Similarly, t-SNE and PCA results for the CGGA cohort confirmed that patients were discretely distributed into two risk groups (Figures 3(c) and 3(d)). Likewise, relative to the patients in the low-risk group, those in the high-risk group exhibited a poorer prognosis with shorter OS (Figures 3(b) and 3(e)). The AUC for 1-, 3-, and 5-year OS was 0.686, 0.753, and 0.778, respectively (Figure 3(f)). With an increase in the risk score, the expression of all genes in the signature increased, except that of NOS1 in the CGGA cohort. Collectively, the results from the CGGA cohort were consistent with those of the TCGA cohort indicating that the 9-gene pyroptosis-related signature correlated positively with gliomas.Figure 3
Prognostic validation of the nine-gene signature in the CGGA cohort. (a, b) Distribution of patients’ risk scores and survival status. (c) The Kaplan-Meier curves for overall survival (OS) in the high- and low-risk groups. (d) ROC curve analysis for verification of predictive performance of the risk model for 1-, 3-, and 5-year OS. (e, f) PCA andt-SNE plot. (g) The heat map shows the levels of gene expression in the high- and low-risk groups.
(a)(b)(c)(d)(e)(f)(g)


## 3.4. Validation of Clinical Independence of the Pyroptosis-Related Prognostic Signature in the CGGA and TCGA Cohorts and Construction of a Predictive Nomogram for the Latter

The findings from univariate COX regression suggested a remarkable association between risk score and OS in both CGGA and TCGA cohorts (CGGA cohort:HR=1.676, 95% CI=1.557-1.804, p<0.001; TCGA cohort: HR=3.150, 95% CI=2.741-3.620, p<0.001) (Figures 4(a) and 4(c)). The risk score remained an independent predictive factor for OS after correcting for confounders, as evidenced by multivariate COX regression analysis results (CGGA cohort: HR=1.267, 95% CI=1.164-1.378, p<0.001; TCGA cohort: HR=2.054, 95% CI=1.654-2.550, p<0.001) (Figures 4(b) and 4(d)).Figure 4
Validation of the clinical independence of the signature in TCGA and CGGA cohorts. (a, b) Univariate and multivariate COX regression analyses in TCGA cohort. (c, d) Univariate and multivariate COX regression analyses in the CGGA cohort.
(a)(b)(c)(d)We developed a nomogram incorporating clinical characteristics of radio status, age, grade, and risk score in the TCGA cohort based on the findings of the multivariate COX regression analyses. The concordance index (C-index) was 0.862, suggesting its relatively high performance for clinical diagnosis (Figure5(a)), wherein a higher risk score predicted a worse prognosis. The nomogram had high accuracy for predicting 1-, 3-, and 5-year OS (AUC=0.901, 0.924, and 0.879, respectively) (Figures 5(b)–5(d)). With excellent performance, the calibration curves demonstrated relatively good fits for predicting OS (Figures 5(e)–5(g)). Taken together, the nomogram may serve as a useful quantitative tool for predicting prognoses of gliomas.Figure 5
Construction of the predictive nomogram for the TCGA cohort. (a) Significant factors obtained from multivariate COX regression analysis in the TCGA cohort were used to construct a nomogram for predicting 1-, 3-, and 5-year survival. (b–d) ROC curve analysis of the nomogram for predicting 1-, 3-, and 5-year survival. (e–g) Calibration curve analysis for the accuracy of nomogram in predicting 1-, 3-, and 5-year survival.∗∗∗p<0.001.
(a)(b)(c)(d)(e)(f)(g)


## 3.5. Functional Annotation and Immune-Related Features of the Pyroptosis-Related Gene Signature

After constructing the pyroptosis-related model, we assessed differential pathways and biological functions between the risk groups in CGGA and TCGA cohorts. A total of 548 DEGs were screened in the TCGA cohort (38 and 510 genes were down- and upregulated, respectively) and 377 DEGs in the CGGA cohort (39 and 338 genes were down- and upregulated, respectively). GO functional enrichment results for the two cohorts suggested that the DEGs were primarily enriched in immune-associated biological processes, cellular components, and molecular functions, including adaptive immune response based on antigen-binding, somatic recombination of immune receptors built from immunoglobulin superfamily domains, humoral immune response, and MHC II protein complex (Figures6(a) and 6(c)). Consistently, the enriched KEGG pathways comprised infectious immune pathways, including the Epstein-Barr virus infection, influenza A, tuberculosis, human T cell leukemia virus 1 infection, and Th17 differentiation (Figures 6(b) and 6(d)). Therefore, the immune microenvironment was remarkably associated with the pyroptosis-related signature.Figure 6
Functional enrichment analyses in TCGA and CGGA cohorts. (a, b) Results of GO and KEGG analyses in TCGA cohort. (c, d) Results of GO and KEGG analyses in the CGGA cohort.
(a)(b)(c)(d)As several immune-related pathways were enriched, thereby highlighting the significance of the immune environment in tumors, we examined the differences in the immune landscape between the risk groups. Immune-related functions, immune checkpoints, and infiltrated immune cell types were analyzed. The results of ssGSEA demonstrated that immune-related functional pathways were remarkably enriched in the high-risk group relative to the low-risk group; similar trends were observed for most immune cells experiencing infiltration and checkpoints (Figures7(a), 7(b), and 7(e)). Consistent results were obtained in the CGGA cohort (Figures 7(c), 7(d), and 7(f)). Collectively, the risk groups showed correlation and obvious differences in immune-related functional pathways, levels of infiltrated immune cells, and immune checkpoints.Figure 7
The risk groups differ in immune-related functions, immune cell infiltration levels, and immune checkpoint expression levels. The scores of immune-related functions and immune cell infiltration in TCGA (a, b) and CGGA (c, d) cohorts. Immune checkpoint expression in TCGA (e) and CGGA (f) cohorts.∗p<0.05, ∗∗p<0.01, and ∗∗∗p<0.001. ns: no significance.
(a)(b)(c)(d)(e)(f)


## 3.6. The Relationship between Oxidative Stress and Pyroptosis Signature

Recently, studies found that cell pyroptosis was associated with oxidative stress in triple-negative breast cancer and non-small-cell lung cancer [25, 26]. However, whether pyroptosis is related to oxidative stress in glioma remains unknown. We acquired fourteen oxidative stress-related pathways from MsigDB (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) by the keyword of “oxidative stress” [27, 28]. The GO terms ID and their corresponding standard names are listed in Supplementary Table S2. Correlation analysis showed that the nine genes in our model were correlated to most of the specific oxidative stress-related pathways. We found that the model genes were almost consistently negatively associated with the pathways of GOBP_OXIDATIVE_DEMETHYLATION, GOBP_REGULATION_OF_TRANSCRIPTION_FROM_RNA_POLYMERASE_II_PROMOTER_IN_RESPONSE_TO_OXIDATIVE_STRESS, GOBP_OXIDATIVE_RNA_DEMETHYLATION, and GOBP_OXIDATIVE_DNA_DEMETHYLATION. On the other hand, the genes, except for NOS1 which was mainly negatively related to the oxidative stress pathway, were strongly positively correlated with the pathways of GOBP_CELLULAR_RESPONSE_TO_OXIDATIVE_STRESS, GOBP_INTRINSIC_APOPTOTIC_SIGNALING_PATHWAY_IN_RESPONSE_TO_OXIDATIVE_STRESS, and GOBP_RESPONSE_TO_OXIDATIVE_STRESS (Figure 8(a)). Moreover, we compared the ssGSEA scores of oxidative stress-related pathways between the high- and low-risk groups. It displayed that the ssGSEA scores of eight GO terms (GOBP_INTRINSIC_APOPTOTIC_SIGNALING_PATHWAY_IN_RESPONSE_TO_OXIDATIVE_STRESS, GOBP_NEGATIVE_REGULATION_OF_CELLULAR_RESPONSE_TO_OXIDATIVE_STRESS, GOBP_NEGATIVE_REGULATION_OF_RESPONSE_TO_OXIDATIVE_STRESS, GOBP_OXIDATIVE_PHOSPHORYLATION, GOBP_REGULATION_OF_OXIDATIVE_PHOSPHORYLATION, GOBP_REGULATION_OF_OXIDATIVE_STRESS_INDUCED_CELL_DEATH, GOBP_REGULATION_OF_OXIDATIVE_STRESS_INDUCED_INTRINSIC_APOPTOTIC_SIGNALING_PATHWAY, and GOBP_RESPONSE_TO_OXIDATIVE_STRESS) were higher in the high-risk group, and four pathways (GOBP_OXIDATIVE_DEMETHYLATION, GOBP_OXIDATIVE_DNA_DEMETHYLATION, GOBP_OXIDATIVE_RNA_DEMETHYLATION, and GOBP_REGULATION_OF_TRANSCRIPTION_FROM_RNA_POLYMERASE_II_PROMOTER_IN_RESPONSE_TO_OXIDATIVE_STRESS) were enriched in the low-risk group (Figure 8(b)). These results were generally the same in the CGGA cohort (Figure S2). It suggested that our novel model of pyroptosis was closely related to oxidative stress pathways.Figure 8
The pyroptosis-related gene signature was correlated with oxidative stress pathways in the TCGA cohort. (a) The relationship between the model gene and oxidative stress pathways. (b) The differences in oxidative stress pathway scores between the high- and low-risk groups.∗p<0.05, ∗∗p<0.01, and ∗∗∗p<0.001. ns: no significance.
(a)(b)


## 3.7. Difference between TMB and Drug Sensitivities in the Two Subgroups

Tumor mutational burden is one of the characteristics and biomarkers for gliomas [29]. We analyzed the correlation of risk scores and TMB scores; the results showed that the risk score correlated positively with TMB and the high-risk group with higher TMB scores than another group (Figures 9(a) and 9(b)). And the top 20 genes with the highest mutational frequencies are shown in Figures 9(c) and 9(d). Most patients carried mutations in the low- and high-risk groups (87.77% and 95.95%, respectively), and missense mutations were the most frequent. TP53 is an important tumor inhibitor gene whose mutational frequency was the highest in the low-risk group (40%), followed by IDH1 (31%). The corresponding mutational frequencies were 43% and 89%, respectively, in the high-risk group.Figure 9
Alteration of TMB between the risk groups in the TCGA cohort. Correlation of TMB and risk score (a). Differences in TMB between the risk groups (b). Top 20 DEGs in glioma between the risk groups (c, d).
(a)(b)(c)(d)To examine the link between drug sensitivities and pyroptosis-related genes, the CellMiner database based on a suite of genomic and pharmacological information including more than 20,000 compounds tested in NCI-60 cell lines was queried. We assessed 75 compounds under clinical trials and 188 FDA-approved drugs. The top 16 most relevant correlations are shown in Figure10. The 9 genes were associated with chemotherapeutic drugs (all p<0.01 and cor>0.3). For instance, the expression of PRTN3 correlated positively with the IC50 value of imexon, ABT-199, cyclophosphamide, and nandrolone phenpr indicating that tumor cells were more resistant to elevated gene expression levels. Conversely, enhanced expression of IL32 was associated with heightened cancer cell-drug sensitivity to tanespimycin, dolastatin 10, tyrothricin, bafetinib, and vinblastine. The differences in the IC50 values of 16 drugs between the risk groups were assessed. Seven drugs (tanespimycin, dolastatin, tyrothricin, staurosporine, geldanamycin analog, bafetinib, and vinblastine) correlated significantly with the risk scores, and except for the staurosporine, all drugs had enhanced sensitivity in the low-risk group (Supplementary Figure S3). Together, our model may provide insights for drugs targeting the pyroptosis-related genes and may bear implications for developing personalized treatment strategies for glioma patients.Figure 10
The correlation between model gene expression and drugs. The top 16 most relevant correlations were visualized: (a) imexon, (b) ABT-199, (c) tanespimycin, (d) cyclophosphamide, (e) dolastatin, (f) nelarabine, (g) tyrothricin, (h) nandrolone phenpr, (i) SR16157, (j) staurosporine, (k) hydroxyurea, (l) geldanamycin analog, (m) bafetinib, (n) E-7820, (o) vinblastine, and (p) chelerythrine.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)


## 3.8. Validation of the Gene Signature in GBM Tissues by HPA and scRNA-Seq

We then verified the gene signature in GBM tissues using IHC images from the HPA database. The expression levels of CASP1, CASP3, CASP6, IL32, MKI67, MYD88, and VIM were higher in tumor tissues, while those of PRTN3 and NOS1 were lower (Figure11). Furthermore, the distribution of these genes in glioma tumors was assessed in GSE138794. Eight samples of gliomas were integrated using “fastMNN” and analyzed using the “Seurat” package. The top 2000 variate genes and features are shown in Supplementary Figure S4. Twenty-two clusters were clarified, and six cell types, including astrocytes (cluster 6, 5.88%), endothelial cells (cluster 18, 1.12%), macrophages (clusters 4 and 9, 12.21%), monocytes (cluster 8, 4.49%), oligodendrocytes (cluster 13, 3.31%), and malignants classified into 14 clusters with the highest proportion of up to 72.99% (Figures 12(a) and 12(b) and Supplementary Figure S5), were screened. The distribution of 8 genes in different cell types is shown in the feature plot (Figure 12(c)) and violin plot (Supplementary Figure S6). PRTN3 was not detected in this dataset. CASP1 was mainly present in clusters 4, 8, and 9, comprising macrophages and monocytes. CASP3 was located in clusters 2, 7, 15, and 17 which were malignant. CASP6 was detected in all clusters at low levels. IL32 was present in cluster 18 comprising endothelial cells. MKI67 showed a remarkably high expression level in cluster 7, indicating that the malignants in cluster 7 have enhanced growth abilities. VIM was distributed in all clusters at high levels but primarily located in malignants. MYD88 was primarily present in clusters 1 (malignants) and 4 (monocytes). NOS1 was rarely found, consistent with the results of IHC. Collectively, we confirmed the levels of protein expression of the gene signature in specific cells by scRNA. Thus, pyroptosis-related genes distributed across different cell types may influence glioma progression via the microenvironment.Figure 11
IHC results for the 9 genes were validated in the HPA database: (a) CASP1, (b) CASP3, (c) CASP6, (d) IL32, (e) MKI67, (f) MYD88, (g) NOS1, (h) PRTN3, and (i) VIM3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 12
Single-cell RNA-seq results for the model genes were validated in GSE138794. (a) Thet-SNE plot showed 22 clusters in GSE138794. (b) The t-SNE plot displayed the annotation of cells in the 22 clusters. (c) The scatter plot showed the expression patterns of eight model genes in different cell types.
(a)(b)(c)


## 3.9. Knocking Down MYD88 Impaired Cell Proliferation and Reduced the Release of Pyroptosis-Related Cytokines

MYD88 was selected to verify the reliability of our model in vitro because of the following reasons: (i) MYD88 was a risk factor upregulated significantly in gliomas relative to normal tissues and correlated with a poor prognosis, (ii) MYD88 was one of the top 5 hub genes in the PPI network, (iii) the coefficient of MYD88 was the maximum in our signature, and (iv) the involvement of MYD88 has been previously reported in LPS-induced pyroptosis [30]. We detected the expression of MYD88 in the normal cell line (HMC-3) versus the glioma cell line (U87) by qRT-PCR. MYD88 was upregulated in U87 cells (1.00±0.07 vs. 2.29±0.33, p=0.003, Figure 13(a)). si-RNA was used to successfully reduce the expression of MYD88 in the U87 cell line (1.00±0.16 vs. 0.56±0.09, p=0.014, Figure 13(a)). Cellular proliferation (72 h and 96 h) was inhibited following treatment with si-MYD88 relative to si-NC-transfected cells (Figure 13(b)). The wound healing assay showed that cellular migration reduced significantly in si-MYD88-treated cells (Figure 13(c)). MYD88 is known to be indispensable in LPS-induced cell pyroptosis via TLR4/MYD88/NF-κB/NLRP3 or TLR-4/MYD88/PI3K/AKT signaling [30–34]. Therefore, we treated si-MYD88- and si-NC-transfected U87 cells with LPS (10 μg/ml) for 24 hours and detected the levels of IL-1β and IL-18 to verify cellular pyroptosis. A simultaneous reduction in the levels of IL-1β (22.63±5.73 vs. 13.36±0.23, p=0.049) and IL-18 (92.10±6.07 vs. 51.86±16.25, p=0.016) was observed in the si-MYD88 group relative to the si-NC group (Figure 13(d)), indicating that downregulation of MYD88 in U87 cells attenuated LPS-induced cell pyroptosis. Thus, MYD88 may exert an oncogenic role, and its activation is a potential treatment target for glioma.Figure 13
Knockdown of MYD88 impaired cell proliferation and reduced release of pyroptosis-related cytokines. (a) The level of MYD88 expression in HMC3 and U87 cell lines (left), and the si-RNA construct effectively attenuated MYD88 in the U87 cell line. (b) Proliferation curves for the U87 cell line after treatment with si-RNA. (c) Wound healing between the si-NC and si-MYD88 groups at 0 h and 48 h. (d) Release of pyroptosis-related cytokines, IL-1β and IL-18, in the si-NC and si-MYD88 groups challenged with LPS. ∗p<0.05;∗∗p<0.01; ∗∗∗p<0.001.
(a)(b)(c)(d)


## 4. Discussion

Glioma is the most prevalent brain cancer with high invasiveness and proliferation which contribute to a high mortality rate [35]. Tumor heterogeneity and resistance to chemotherapy induce a high rate of recurrence in gliomas [4, 36]. Recently, the induction of tumor cell death via a programmed pathway has been suggested as a promising strategy for cancer treatment [37]. Pyroptosis is a newly identified proinflammatory pathway of programmed cell death. The canonical pathway is dependent on activated CASP1 which causes the secretion and maturation of IL-18 and IL-1β. It also cleaves the C- and N-terminal domains of GSDMD. The N-terminal domain is assembled and forms 10-15 nm pores on the cell membrane, thereby promoting cell rupture. The noncanonical pyroptosis pathway relies on the activation of CASP-4/5/11 [6, 38]. Pyroptosis is implicated in several tumors, including colon cancer [39], lung adenocarcinoma [40], and hepatocellular carcinoma [41]. Although tumor progression can be suppressed by cell pyroptosis without infection, a proinflammatory microenvironment with IL-1β, IL-18, and cellular components may facilitate tumor growth [12, 42]. Assessing the role of pyroptosis in tumors is beneficial for understanding their progression and developing targeted drugs. Therefore, we focused on the pyroptosis-related genes in gliomas and constructed a prognostic signature for predicting the OS of these patients.First, we identified 21 pyroptosis-related prognostic DEGs between gliomas and normal tissues. LASSO-Cox regression was performed to select optimal genes for constructing the prognostic model. Finally, nine genes (CASP1, CASP3, CASP6, IL32, MYD88, MKI67, NOS1, PRTN3, and VIM) comprised the signature. Subsequently, the risk scores of each patient were computed, and the cohort was stratified according to the median value into low- and high-risk groups. The Kaplan-Meier curves demonstrated that patients at high risk had a worse prognosis relative to those at low risk. The AUC for 1-, 3-, and 5-year OS was 0.866, 0.895, and 0.847, respectively, indicating that the signature was reliable in terms of OS prediction for gliomas. Similar results were noted in the validation cohort suggesting that the signature was stable and reliable. A nomogram including clinical factors and risk score was built comprising independent factors, as confirmed by multivariate and univariate COX regression analysis results. ROC and calibration curves suggested that the nomogram was excellent in predicting 1-, 3-, and 5-year OS. Several similar models comprising pyroptosis-related genes have been proposed [43–45]; however, the analysis criteria and source of pyroptosis-related genes were different, and our signature showed a similar ability in predicting the OS for glioma. The signature comprising pyroptosis-related genes was reliable and may provide a clinical tool for predicting OS.Among the nine model genes, all except for NOS1 were risk factors. They were not only upregulated in glioma tissues but also increased with the risk scores analyzed by bulk sequencing and IHC data. Furthermore, we assessed the gene expression in the sc-RNA dataset. CASP-1 activated by inflammasome is vital in the canonical pyroptosis pathway [6]. CASP-1 was upregulated in the glioma tissues and linked to a poor prognosis. scRNA-seq results showed that CASP-1 localized primarily to the macrophages and monocytes. Pyroptotic activities are prominent in macrophages and monocyte-derived dendritic cells [46]. Thus, we speculated that CASP1 may promote glioma progression through pyroptosis in macrophages and monocytes. The expression of CASP-3 was higher in gliomas relative to normal tissues, and it was a risk factor, mainly present in “malignants.” These findings were consistent with the results of Zarnescu et al. and Bodey et al. [47, 48]. The activation of CASP3 can induce tumor cell pyroptosis via the GSDME-dependent pathway in colon cancer following treatment with lobaplatin [39]. Therefore, triggering CASP-3 in tumor cells using specific drugs is a promising treatment strategy for gliomas. CASP-6 participates in cellular apoptosis and necroptosis [49]. Inhibition of CASP-6 may be promising for various neurodegenerative disorders [50]. CASP-6, a risk gene, was overexpressed in our signature. Downregulation of CASP-6 may be a potential strategy for the treatment of gliomas. IL32 is a proinflammatory cytokine that can contribute to the onset and progression of different tumors, such as lung adenocarcinoma [51], breast cancer [52], and gastric cancer [53]. Herein, IL32 was localized primarily to endothelial cells and upregulated in glioma tissues. It may favor glioma growth by facilitating an inflammatory microenvironment for the tumor. We searched for FDA-approved drugs and compounds undergoing clinical trials that targeted the genes in our signature using the CellMiner database and performed a correlation analysis. Among the top 16 drugs, tanespimycin, dolastatin 10, tyrothricin, geldanamycin analog, bafetinib, vinblastine, and staurosporine all targeting IL-32 showed differential drug sensitivities between the risk groups. Stratifying patients into risk groups based on the signature may facilitate better treatment outcomes, and further studies are needed to validate its utility. Furthermore, the blood-brain barrier (BBB) is the obstacle that must be overcome in the treatment of malignant gliomas and should be taken into account when using the drugs mentioned above. Only small molecules (<500 Da and 400 nm) and lipophilic elements can pass through the BBB passively, and the others cross the BBB via pinocytosis, receptor- or carrier-mediated transcytosis, and solute-carrier-protein mechanisms [54, 55]. Tanespimycin [56], geldanamycin analog [57], bafetinib [58], and vinblastine [59] as small lipophilic molecules were capable of penetrating the blood-brain barrier. As for the toxic compound staurosporine, it may disrupt the endothelial barrier and increase the permeability of barrier function [60]. Therefore, the five drugs listed above could cross the blood-brain barrier and become potential glioma treatments. MKI67 is a marker of cell proliferation that highly correlates with glioma grade [61]. Pyroptosis reduces cell growth, thus attenuating the expression of MIK67. MKI67 was upregulated in cluster 7 comprising the malignants. It indicated that cluster 7 had a high proliferative ability, which in turn led to glioma progression. VIM is a key member of the intermediate filament family of proteins overexpressed in several tumors and acts as a marker of epithelial-mesenchymal transition (EMT). It also contributes to drug resistance [62]. In glioma patients, elevated VIM levels are correlated with a poor prognosis and temozolomide resistance [63]. Similarly, VIM was expressed at a high level in glioma tissues and majorly localized to the malignant cell types in GSE138794. Thus, targeting VIM may restrain EMT and resistance to chemotherapy, ultimately inhibiting glioma progression. In the peripheral nervous system and brain, NOS1 exhibits the properties of a neurotransmitter. It is involved in neurotoxicity associated with neurodegenerative disorders and stroke, including gliomas [64]. NOS1 was downregulated in gliomas and correlated with a good prognosis. NOS1 was negligibly detected in the sc-RNA-seq results. Modulating the expression of NO is thus a promising therapeutic target. PRTN3 is upregulated in pancreatic cancer [65], vulvar squamous cell carcinoma [66], and clear cell renal carcinoma [67]. It is a prognostic marker that contributes to cancer progression, invasion, and metastasis. An important paralog of this gene is ELANE which has been included in another pyroptosis-related gene signature for glioma [67, 68], and thus, PRTN3 is putatively associated with pyroptosis. Like the case in other various cancers, PRTN3 was downregulated in glioma and was a risk factor in our signature. Guo et al. have recently reported that MYD88 is a prognostic gene for gliomas and is highly associated with macrophage infiltration [69]. MYD88 is implicated in adenovirus-induced cancer stem cells. In vivo, MYD88 deficiency attenuates glioma expansion via the p38-MAPK pathway [70]. Furthermore, MYD88 participates in LPS which induces cell pyroptosis via the TLR4/MYD88/NF-κB/NLRP3 or TLR-4/MYD88/PI3K/AKT pathway; deletion or inhibition of MYD88 attenuates the cell response [30, 31]. Therefore, MYD88 is indispensable in cellular pyroptosis. To some extent, these results can explain our findings. We performed in vitro assays and lowered the expression level of MYD88 in the U87 glioma cell line. The proliferation and migration of U87 cells were inhibited. Furthermore, si-MYD88-treated cells secreted less pyroptosis-related cytokines, IL-1β and IL-18, when challenged with LPS. Our results further suggested that MYD88 may play an oncogenic role in glioma and activated MYD88 could inhibit glioma progression via cell pyroptosis.The immune microenvironment including T cells, macrophages, monocytes, and neutrophils can affect glioma prognosis to a large extent [71]. Our pyroptosis-related signature was remarkably associated with immune functions, immune cell infiltration, and immune checkpoints. Pyroptosis in cancer can attract immune cells into the tumor [72]. The abundances of CD8 + T cells, NK cells, and Th1 cells were higher, and those of Tregs and Th2 cells increased simultaneously. Similar results were obtained for immune functions, where the scores of APC and T cell inhibition/stimulation were higher in the high-risk group. These findings suggested that immune suppression is the main reason for worse prognoses among patients at high risk. Immunotherapy based on the blockade of the immune checkpoints is a promising strategy for glioma treatment [73]. Targeting the PD-1/PD-L1 pathway has been extensively investigated in immunotherapy. Zhao et al. report that less than 10% of glioma patients show long-term responses to the anti-PD-1 immunotherapy, but the responders have longer survival relative to the nonresponders [74]. However, markers for predicting immunotherapeutic outcomes are scarce. Although glioma patients exhibit higher expression of immune checkpoints, their further identification can improve individualized immunotherapeutic options [75]. We compared the immune checkpoint levels of LAG-3, CTLA-4, and PD-1 between the low- and high-risk groups. The high-risk group showed high expression of most immune checkpoints. It not only suggests that the immune microenvironment in the high-risk group is highly immunosuppressive but also indicates a greater likelihood of patients at high risk benefiting from immunotherapy. TMB is also a prognostic predictor and indicator of immunotherapeutic efficacy [76]. We compared the scores of TMB between the risk groups. The high-risk group had elevated mutational levels, and most patients herein carried IDH-1 mutations. Although IDH-1 mutation is a marker for favorable prognostic outcomes in glioma, the patients at high risk had worse survival, and a highly immunosuppressive microenvironment may contribute to the phenomenon. Taken together, our pyroptosis-related signature may provide insights into personalized immunotherapeutic options.Oxidative stress has been reported to be involved in favoring cancer progressions, such as hepatocellular carcinoma [77], breast cancer [78], gastric cancer [79], and gliomas [80]. Reactive oxygen species (ROS) brought on by oxidative stress compromise DNA, lipids, and protein stability and may be the main factor contributing to the development of malignancies [81]. Recently, some studies discovered cross-talk between pyroptosis and oxidative stress in cancer. In triple-negative breast cancer, cadmium treatment of cell lines raised the production of ROS, which then activated the NLRP3 pathway and promoted cell pyroptosis [26]. In non-small-cell lung cancer, ROS generation and NLRP3 inflammasome activation were also facilitated by lncRNA-XIST knockdown [25]. However, the relationship between oxidative stress and pyroptosis in glioma is still unclear. In our study, we obtained fourteen oxidative stress-related pathways from MsigDB. The nine genes in our model showed a close correlation with most of the specific pathways. It seems that the genes are negatively associated with the epigenetic regulation related to oxidative stress, including GO:0070989 (OXIDATIVE DEMETHYLATION), GO:0035511 (OXIDATIVE DNA DEMETHYLATION), GO:0035513 (OXIDATIVE RNA DEMETHYLATION), and GO:0043619 (REGULATION OF TRANSCRIPTION FROM RNA POLYMERASE II PROMOTER IN RESPONSE TO OXIDATIVE STRESS). Meanwhile, the ssGSEA scores of the four mentioned pathways were higher in the low-risk group. It further implied that oxidative stress-induced epigenetic regulation may be contributed to an unfavorable prognosis in gliomas. In addition, the genes were positive with GO:0034599 (CELLULAR RESPONSE TO OXIDATIVE_STRESS), GO:0008631 (INTRINSIC APOPTOTIC SIGNALING PATHWAY IN RESPONSE TO OXIDATIVE STRESS), and GO:0006979 (RESPONSE TO OXIDATIVE STRESS). Interestingly, the three pathways were enriched in the high-risk group, which indicated that the pyroptosis-related genes may participate in the pathways of response to oxidative stress; the underlying mechanism needs further experiments. In a word, our novel pyroptosis-related signature was strongly associated with oxidative stress pathways and showed different enrichment in the two subgroups.Undoubtedly, some inevitable limitations exist in this study. First, we sourced public data for the construction of the model. Similar predictive powers were obtained in TCGA and CGGA datasets, but we did not include an internal cohort to verify the effectiveness of the signature, making it less reliable to some extent. From the translational aspect in clinical settings, we plan to include glioma patients in the future. Second, the prediction of drug sensitivities should be experimentally validated. Lastly, though we only conducted in vitro experiments to reveal the oncogenic role of MYD88, other pyroptosis-related genes in this model also need validation both in vitro and in vivo.


## 5. Conclusion

In conclusion, we constructed a pyroptosis-related gene signature for predicting survival and drug sensitivities in gliomas via comprehensive bioinformatic analysis methods and experimental tests. Furthermore, the risk groups classified using the signature exhibited differences in immune functions, infiltration levels of immune cells, immune checkpoint expression, oxidative stress-related pathways, and TMB scores. These findings may have a predictive value and clinical translation for glioma therapy.



---

*Source: 4256116-2023-02-01.xml*