Pyroptosis-related Prognosis Model, Immunocyte Infiltration Characterization, and Competing Endogenous RNA Network of Glioblastoma

doi:10.21203/rs.3.rs-1095246/v1

Download PDF

Primary research

Pyroptosis-related Prognosis Model, Immunocyte Infiltration Characterization, and Competing Endogenous RNA Network of Glioblastoma

https://doi.org/10.21203/rs.3.rs-1095246/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background: Glioblastoma (GBM) has a high incidence rate, invasive growth, and easy recurrence, and the current therapeutic effect is less than satisfying. Pyroptosis plays an important role in morbidity and progress of GBM. Meanwhile, the tumor microenvironment (TME) is involved in the progress and treatment tolerance of GBM. In the present study, we analyzed prognosis model, immunocyte infiltration characterization, and competing endogenous RNA (ceRNA) network of GBM on the basis of pyroptosis-related genes (PRGs).

Methods: The transcriptome and clinical data of 155 patients with GBM and 120 normal subjects were obtained from The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx). Lasso (Least absolute shrinkage and selection operator) Cox expression analysis was used in predicting prognostic markers, and its predictive ability was tested using a nomogram. A prognostic risk score formula was constructed, and CIBERSORT, ssGSEA algorithm, Tumor IMmune Estimation Resource (TIMER), and TISIDB database were used in evaluating the immunocyte infiltration characterization and tumor immune response of differential risk samples. A ceRNA network was constructed with Starbase, mirtarbase, and lncbase, and the mechanism of this regulatory axis was explored using Gene Set Enrichment Analysis (GSEA).

Results: Five PRGs (CASP3, NLRP2, TP63, GZMB, and CASP9) were identified as the independent prognostic biomarkers of GBM. Prognostic risk score formula analysis showed that the low-risk group had obvious survival advantage compared with the high-risk group, and significant differences in immunocyte infiltration and immune related function score were found. In addition, a ceRNA network of messenger RNA (CASP3, TP63)–microRNA (hsa-miR-519c-5p)–long noncoding RNA (GABPB1-AS1) was established. GSEA analysis showed that the regulatory axis played a considerable role in the extracellular matrix (ECM) and immune inflammatory response.

Conclusions: Pyroptosis and TME-related independent prognostic markers were screened in this study, and a prognosis risk score formula was established for the first time according to the prognosis PRGs. TME immunocyte infiltration characterization and immune response were assessed using ssGSEA, CIBERSORT algorithm, TIMER, and TISIDB database. Besides a ceRNA network was built up. This study not only laid foundations for further exploring pyroptosis and TME in improving prognosis of GBM, but also provided a new idea for more effective guidance on clinical immunotherapy to patients and developing new immunotherapeutic drugs.

Cancer Biology

Glioblastoma

Pyroptosis

Tumor microenvironment

Immunocyte infiltration

Competing endogenous RNA network

Glioblastoma (GBM) is the most common malignant tumor of the central nervous system [1]. It is characterized by invasive growth and easy recurrence [2]. Recently, despite continuous developments in treatment methods, such as surgery, radiotherapy, and chemotherapy have been reported, the median survival time of patients remains 12–15 months, and no substantial improvement has been achieved yet [3]. Some reasons are related to the potential molecular and cell heterogeneity of GBM [4]. Therefore, exploring the potential molecular mechanism of GBM and determining the key therapy targets are of important significance.

Pyroptosis is an inflammatory form of programmed cell death, and its characteristic lies in membrane perforation, cell swelling, plasma membrane rupture, and intracellular content release [5]. Pyroptosis is mainly mediated by the gasdermin family, which includes GSDMA, GSDMB, GSDMC, GSDMD, GSDME/DFNA5, and DFNB59. Except DFNB59, the other members of the gasdermin family are activated by cleaving the C-terminal and N-terminal domains. Specifically, the N-terminal segments form holes in the plasma membrane, thus resulting in the swelling and rupture of cells [6]. Pyroptosis cells release proinflammatory molecules to the extracellular environment, thus triggering inflammation and immunoreaction [7]. Pyroptosis induced by the combination of anti-tumor drugs and autophagy inhibitors can inhibit GBM growth and increase survival rate [8]. Ren et al. [9] reported that activating NF-κB/NLRP3/GSDMD pathway can trigger pyroptosis of GBM and thereby inhibit GBM growth. Pyroptosis plays a crucial role in the pathogenesis and progress of various malignant tumors, including GBM [10].

According to cumulative evidence, the development and progression of cancer occur simultaneously and accompanied by the alteration of the surrounding stroma. The microenvironment on which tumor growth and metastasis plays an important role in the occurrence and development of tumors [11]. The tumor microenvironment (TME) is a complex network formed by distinct cell types and extracellular matrix (ECM) components, including fibroblasts, vascular endothelial cells, various infiltrating immunocytes (including cytokines and chemokines released by these cells), and ECM proteins [12]. Interactions between tumor cells and immune cells can promote tumor growth and metastasis and influence tumor progression [11]. Pyroptosis can reprogram the TME to an immunostimulation state [13], and the TME plays an important role in the progression and treatment resistance of GBM [14]. Immunotherapeutic approaches targeting the TME is one of the effective means for treating GBM [15]. Therefore, immunotherapy combined with pyroptosis-related therapy offer prospects for GBM treatment.

Interactions between tumor cells and infiltrating immunocytes are generally regulated by competing endogenous RNA (ceRNA) networks, which are composed of messenger RNAs (mRNAs), long noncoding RNAs (lncRNAs), and microRNAs (miRNAs). These networks can regulate the post-transcription of oncogenes and tumor suppressor genes and the interactions between protein and genes, thus regulating the biological behavior of tumors, particularly invasion and metastasis [16].

Most studies are limited to one or two pyroptosis-related genes (PRGs) and TME cell types. The antitumor effect is a characteristic feature of highly coordinated interactions among many tumor suppressors [17]. Therefore, comprehensive understanding of TME cell infiltration characteristics related to multiple PRGs and ceRNA networks might provide a new perspective for the study of the potential mechanisms of the occurrence and development of GBM and prediction of its responses to immunotherapeutic approaches.

In the present study, potential independent prognosis pyroptosis-related genes (prognosis PRGs) of GBM were screened comprehensively, and a prognostic model was constructed. Independent PRGs and the immunoreactions of tumors are associated. Subsequently, a prognosis risk evaluation formula was established on the basis of independent prognosis PRGs, and the TME cell infiltration characterization of differential risk scoring of GBM was assessed. The results demonstrated that independent prognosis PRGs played a considerable role in the formation of the TME characterization features of GBM. These genes might be attributed to the potential immunotherapeutic targets and biomarkers of GBM. The ceRNA network of mRNA (CASP3, TP63)–microRNA (hsa-miR-519c-5p)–lncRNA (GABPB1-AS1) was constructed for the first time on the basis of a pyroptosis-related prognostic model with bioinformatics tools for the reverse prediction of target genes. The possible action mechanisms of this regulation axis in the growth, invasion, and metastasis of GBM were identified using Gene Set Enrichment Analysis (GSEA) for the interpretation of the prognosis value of PRGs in GBM and relevant molecular mechanisms.

Datasets and preprocessing

The RNA-sequencing (RNA seq) data of 155 GBM and five normal subjects and relevant clinical data were collected from The Cancer Genome Atlas (TCGA) database [18] (https://portal.gdc.cancer.gov/). The RNA seq data of the brain tissues of 115 normal subjects were obtained from the Genotype-Tissue Expression (GTEx) project [19] (https://www.gtexportal.org/). Single-nucleotide polymorphism (SNP) data were downloaded from TCGA. The unit of RNA seq data was normalized as FPKM.

Identification of differentially expressed PRGs

A total of 49 PRGs were extracted from GSEA-MSigDB (http://www.gsea-msigdb.org/gsea/login.jsp), Reactome (https://reactome.org/), Harmonizome (https://maayanlab.cloud/Harmonizome/), and Pubmed (https://pubmed.ncbi.nlm.nih.gov/). By using R software (version 4.1.1) and “ggpubr” R package, 42 PRGs were screened according to P < 0.05. The protein–protein interaction (PPI) networks of the PRGs were presented using the Search Tool for the Retrieval of Interacting Genes (STRING; version 11.5, https://string-db.org/) and Cytoscape software (version 3.8.2). The “igraph” and “reshape2” R packages were used in building a differential gene co-expression network.

Analysis of the somatic mutation characteristics of PRGs

The somatic mutation frequency of the 49 PRGs was analyzed and visualized with the “maftools” R package. A mutation abstract plot and oncoplot were then generated.

Functional enrichment analysis

We performed Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Gene and Genomes (KEGG) pathway analyses with the “clusterProfiler” R package. The biological processes, cell components, molecular functions, and KEGG pathways of PRGs were visualized using the “ggplot2” and “GOplot” R packages.

Establishment and validation of the pyroptosis-related gene prognostic model

Cox regression analysis was used in evaluating the prognosis value of the PRGs. Forest map was used to display the hazard ratio (HR), 95% confidence interval (CI), and P values of each variable. HR was used in determining the protective genes (HR < 1) and risk genes (HR > 1). Based on these PRGs related with prognosis, a Lasso (Least absolute shrinkage and selection operator) Cox regression model was constructed with the “glmnet” R package to screen optimal prognosis PRGs. The optimal value of penalty parameter (λ) was selected according to 1000 cross-validation runs. The weighted regression coefficient and expression levels of the prognosis PRGs were extracted, and the risk score was calculated using the following formula as the index for measuring the survival risk of each patient:

According to the median risk score, patients with GBM in the TCGA were divided into low-risk (risk score < median risk value) and high-risk (risk score > median risk value) groups. The data visualization algorithms of principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were performed using the “Rtsne” and “ggplot2” R packages. The survival times of the two groups were compared through Kaplan-Meier analysis. Receiver operating characteristic (ROC) analysis was carried out with the “survival,” “survminer,” and “timeROC” R packages for the evaluation of the prediction accuracy of different genes and risk scores.

Independent prognostic analysis and clinical value of the risk model

To determine whether the risk score from the gene characteristic model can be used as an independent prognosis factor for patients with GBM, we conducted a univariate and multivariate Cox regression analysis to the five prognosis PRGs. A clinical feature heatmap was built by extracting the clinical information (age and gender) of patients in the TCGA group and combining with risk score in the regression model. A nomogram can be used in predicting the prognosis of cancer. In this study, a prognostic nomogram was built by including the prognosis PRGs and used in analyzing the probability of 1-, 2-, and 3-year survival times of patients with GBM. This model was implemented using the “rms” R package. The calibration curve of the nomogram was plotted for the evaluation of its prediction accuracy.

Immunocyte infiltration analysis

The correlation between prognosis PRGs and immunocyte levels (e.g., B cell, CD8+ T cell, CD4+ T cell, macrophage, neutrophil, and dendritic cell) was detected with Tumor IMmune Estimation Resource (TIMER) database (http://timer.cistrome.org/). The correlation between the expression of prognosis PRGs and immunoregulators (including immunoinhibitor, immunostimulator and MHC molecule) was analyzed using the TISIDB database (http://cis.hku.hk/TISIDB/).

To quantify the relative proportions of immunocytes in the high-risk and low-risk groups, we calculated immunocyte infiltration with the CIBERSORT algorithm. Samples were screened according to the standard (P < 0.05), and the percentage of each immunocyte type in the samples was calculated, which were exhibited in a bar diagram. Correlation heatmap analysis was carried out with the “corrplot” R package, which disclosed the correlations among 22 types of immunocytes. Additionally, the scores of 16 immunocytes and 13 immune-related functions were calculated with the ssGSEA algorithm for the analysis of the differential expression levels in the high-risk and low-risk groups.

CeRNA network construction

To elaborate the potential functions of prognosis PRGs in GBM, a ceRNA network was built. Starbase (version 2.0, http://starbase.sysu.edu.cn/) and MirTarbase (version 9.0 beta, http://mirtarbase.cuhk.edu.cn/) were used in predicting the target miRNA of the prognosis PRGs. The top 10 highly connected miRNAs were screened by calculating the degree among various nodes with the CytoHubba plug-in of Cytoscape software (version 3.8.2), and the prognosis values of the screened miRNAs were analyzed using the Oncolnc database (http://www.oncolnc.org/). Based on above miRNA, Lncbase (version v.2, https://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=lncbasev2/) and Starbase were used in predicting the lncRNA targets that interacted with miRNAs, and the intersection of the lncRNA targets were screened with a Venn diagram (https://bioinfogp.cnb.csic.es/tools/venny/). The expression and prognosis conditions of the lncRNA targets were analyzed with GEPIA2 (http://gepia2.cancer-pku.cn/).

Gene Set Enrichment Analysis

To explore the influences of gene expression on this pathway, GSEA was carried out. The first 50% of PRG expression was set as the high-expression group, and the rest was set as the low-expression group. GSEA was carried out with GSEA software (version 4.1.0).

Defining the expression of PRGs in GBM

A total of 42 differential expressed PRGs in GBM were recognized through the comparison of the differential expressions of 49 PRGs in 120 normal subjects and 155 GBM from the TCGA and GTEx database (P < 0.05). The expression levels of 13 genes (NLRP2, NLRP7, TNF, IL1B, IL6, NLRP1, IL1A, PRKACA, GSDMB, CHMP3, CHMP4B, CHMP2B, and CHMP7) in GBM were downregulated relative to those in the brain tissues of normal subjects, whereas the expression levels of 29 genes (PLCG1, CASP9, CHMP6, SCAF11, CHMP4A, CHMP2A, HMGB1, IRF2, NLRC4, GZMB, NOD1, TP63, BAK1, CASP8, BAX, IRF1, CASP3, PYCARD, IL18, GSDMD, AIM2, GSDMA, NOD2, GZMA, CASP5, CASP6, TP53, CASP1, and CASP4) were upregulated (Figure 1A). To further explore the interactions of these differentially expressed PRGs, we carried out PPI analysis. As shown in Figure 1B, the minimum required interaction score was set at 0.4, and the degree of each node in the network was expressed and arranged according to color. The top 10 genes in the term of degree were TNF, IL1B, IL18, CASP8, IL6, PYCARD, AIM2, NLRC4, AIM2, and NLRC4. The correlation network results of the differentially expressed PRGs are shown in Figure 1C. Most differentially expressed PRGs were positively related, whereas negative correlations between GZMA and CHMP2A, between IL6 and PLCG1, and between GSDMB and CHMP4A were found.

Analysis of the somatic mutation characteristics of PRGs

Somatic mutation frequency of 49 PRGs in GBM was analyzed. The results demonstrated that missense mutation was the most common mutation classification, and SNP was the most common variant type. Single-nucleotide variation mainly occurred in the form of C > T (Figure 2A). In the GBM samples, the top 10 mutation genes related to pyroptosis were NLRP3, NLRP7, NLRP2, NOD1, CASP1, NLRP1, PLCG1, GZMB, NOD2, and CHMP4C (Figure 2B).

Functional enrichment analysis of the PRGs

GO and KEGG databases were used in analyzing the related pathways of differentially expressed PRGs and clarifying their functions. As shown in Figure 3A, the differentially expressed PRGs were mainly enriched in the positive regulation of cytokine production, regulation of interleukin-1 production, positive regulation of cysteine-type endopeptidase activity, and pyroptosis biological processes. In addition, the genes were mainly related to multiple cell components, such as inflammasome complex, ESCRT complex, multivesicular body, late endosome membrane, pore complex, and immunological synapse. Moreover, molecular function analysis verified that the PRGs were correlated to endopeptidase activity, cytokine receptor binding, protease binding, cysteine-type peptidase activity, peptidase activator activity, and CARD domain binding. KEGG pathway enrichment analysis demonstrated that the PRGs were closely related to NOD-like receptor signaling pathway, necroptosis, Salmonella infection, legionellosis, lipid and atherosclerosis, and pathogenic Escherichia coli infection (Figure 3B and C).

Construction of the prognostic PRG model

According to the survival-related genes screened through univariate Cox regression analysis, 34 differentially expressed PRGs (CHMP4A, GSDMA, CHMP3, AIM2, CHMP2B, CASP5, NLRC4, BAX, IL18, IRF1, IL1B, CASP8, NLRP1, IL6, IL1A, CASP1, CHMP6, TP63, PRKACA, NLRP2, PYCARD, GZMA, CHMP7, CHMP2A, CHMP4B, CASP6, GSDMD, CASP3, NOD2, NOD1, IRF2, CASP4, GZMB, and NLRP7) had HRs of >1, indicating that they were related to increased risk. The other genes (CASP9, BAK1, BAK1, HMGB1, GSDMB, PLCG1, SCAF11, and TNF) had HRs of <1 and were protective factors (Figure 4A).

To evaluate the prognosis value of these differentially expressed PRGs, a Lasso-Cox regression model was built with “glmnet” R package and used in further screening. According to the minimum penalty parameter (λ), five of the 42 differentially expressed PRGs were retained, namely, CASP3, NLRP2, TP63, GZMB, and CASP9 (Figure 4B and C). A prognosis risk score formula was established according to weighted regression coefficient and expression levels based on multivariate Cox regression analysis: Risk score = 0.0046 × (expression value of CASP3) + 0.0060 × (expression value of NLRP2) + (0.0227) × (expression value of TP63) + 0.1804 × (expression value of GZMB) - (0.0481) × (expression value of CASP9). According to the median score calculated in the prognosis risk score formula, 152 patients with GBM were divided into low- and high-risk groups. PCA and t-SNE analysis demonstrated that patients with different risks can be divided into two types (Figure 4D and E). The risk score distribution and survival states of the patients are shown in Figure 4F and G. As risk score increased, the risk of death in the patients increased, whereas survival time was shortened. According to the Kaplan–Meier survival curve, the low-risk group had more obvious survival advantages (P = 0.009, Figure 4H). The area under curve values of the 1-, 3-, and 5-year ROC curves were 0.675, 0.669, and 0.723, respectively (Figure 4I), indicating that the pyroptosis model can predict the 1-, 3-, and 5-year survival rates of patients with GBM.

Independent prognostic analysis and nomogram construction based on risk model

Whether risk score from the gene characteristic model can be used as an independent prognosis biomarker of GBM was determined using univariate and multivariate Cox regression analysis. Univariate and multivariate Cox regression analyses demonstrated that the risk score was an independent factor that influenced the prognosis of patients with GBM (HR = 6.917, 95% CI: 2.760–17.330 and HR = 6.168, 95% CI: 2.467–15.420; Figure 5A and B). In addition, the clinical feature heatmap in Figure 5C showed that THE gender and age of patients showed no significant differences between the low- and high-risk groups (P > 0.05).

To discuss THE importance of these five prognosis PRGs to THE prognosis prediction of patients with GBM, we constructed a nomogram to visualize the Cox regression model as the survival probability of patients. The performance of the model was evaluated using a calibration curve. Compared with the ideal model, the predicted total 3-year survival rate of patients with GBM based on prognosis PRGs was slightly different from the practical value. All the five prognosis PRGs had good prognosis prediction ability for GBM (Figure 5D and E).

The risk score was closely related to the expression levels of the five prognosis PRGs. It was an important factor for predicting the prognosis results of patients with GBM and had good prediction ability. The five prognosis PRGs might participate in the progress of GBM.

Correlation between the expression of prognostic PRGs and immunocyte levels and immunoregulators in GBM

Using the TIMER database, we detected the correlation of five prognosis PRGs with the levels of immunocyte infiltration in GBM. As shown in Figure 6A, CASP3 was significantly associated with CD8+ T cell (r = 0.109; P = 2.65e-02), macrophage (r = 0.084; P = 8.72e-02), neutrophil (r = 0.155; P = 1.50e-03), and dendritic cell (r = 0.193; P = 6.89e-05). NLRP2 was significantly associated with CD8+ T cell (r = −0.197; P = 4.96e-05). TP63 was significantly associated with CD4+ T cell (r = −0.143; P = 3.34e-03), macrophage (r = −0.167; P = 5.93e-04), and neutrophil (r = −0.196; P = 5.65e-05). GZMB was significantly associated with CD8+ T cell (r = −0.203; P = 3.00e05) and CD4+ T cell (r = −0.174; P = 350e-04). CASP9 was significantly associated with CD4+ T cell (r = 0.199; P = 4.22e-05) and macrophage (r = 0.174; P = 3.48e-04). The above results showed a significant correlation between prognosis PRGs and tumor immunocyte infiltration.

To further address the influences of prognostic PRGs on tumor immune response, the correlation between the expression of prognostic PRGs and immunoregulators was calculated according to the TISIDB database. CASP3 and CASP9 were negatively correlated with immunoinhibitor, immunostimulator and MHC molecule in GBM, whereas NLRP2, TP63, and GZMB were positively related to them (Figure 6B). Hence, the expression levels of the prognosis PRGs were closely related to the abundance of immunocyte infiltration and immunoregulators.

Immunocyte infiltration between the high-risk and low-risk groups

We used the CIBERSORT algorithm to analyze the abundance of immunocyte infiltration between the high-risk and low-risk groups in patients with GBM. The high-risk group presented higher distribution levels of monocytes, activated dendritic cells, activated mast cells, and eosinophils and lower distribution levels of resting NK cells and resting mast cells (Figure 7A). According to the correlations among 22 types of immunocytes, activated dendritic cells were positively related with naive CD4 T cells (r = 0.47), whereas M0 macrophages were negatively related with regulatory T cells (r = 0.42). The resting NK cells were negatively correlated with activated NK cells (r = −0.69), monocytes, and M2 macrophages were negatively correlated with M0 macrophages (r = −0.67; r = −0.6; Figure 7B).

The differential expression of immunocyte infiltration and immune-related function score between the high- and low-risk groups was analyzed by ssGSEA algorithm. Significant differences in B cells, CD8+ T cells, dendritic cells, macrophages, neutrophils, NK cells, T helper cells, tumor-infiltrating lymphocytes, and regulatory T cells were observed between the groups (Figure 7C). Moreover, significant differences in the scores of APC co-inhibition, APC co-stimulation, CCR, check point, cytolytic activity, HLA, inflammation-promoting process, parainflammation, T-cell co inhibition, T-cell co stimulation, type I FN response, and type II FN response were found between the groups (Figure 7D).

Significant difference in immunocyte infiltration was found between the groups. Therefore, monocytes, dendritic cells, eosinophils, NK cells, and mast cells might be potential immunocytes that play important roles and participated in the progress of GBM.

Construction of the network of mRNA–miRNA–lncRNA

To elaborate the potential molecular mechanism of the five pyroptosis-related mRNAs (CASP3, NLRP2, TP63, GZMB, and CASP9) which had prognosis values in GBM, an mRNA–miRNA–lncRNA co-expression network was built. miRNA targets that combined with these five prognostic PRGs were predicted reversely by using mirTarBase and Starbase. A total of 514 miRNA targets of CASP3, NLRP2, TP63, and CASP9 were obtained (Figure 8A). The degree values among miRNAs nodes were calculated using the CytoHubba plug-in of Cytoscape software, and the top 10 highly connected miRNAs were determined (hsa-miR-342-3p, hsa-miR-582-5p, hsa-miR-3163, hsa-miR-320a, hsa-miR-320c, hsa-miR-320d, hsa-miR-320b, hsa-miR-526a, hsa-miR-519b-5p, and hsa-miR-519c-5p; Figure 8B). The OncoLnc database determined a miRNA, hsa-miR-519c-5p, and GBM patients with high expression levels of hsa-miR-519c-5p. These patients had relatively high total survival rates (P = 0.0337; Figure 8C). Based on the above miRNA, the upstream lncRNA targets were explored for the construction of a miRNA–lncRNA axis. As shown in Figure 8D, for the intersection lncRNAs of Lncbase and Starbase, KCNQ1OT1, GABPB1-AS1, ENTPD1-AS1, XIST, LINC01018, and LINC00662 were determined as targets. According to survival analysis and differential expression of the lncRNA targets, only GABPB1-AS1 significantly decreased the survival probability of patients with GBM (P = 0.05; Figure 8E) and it was upregulated significantly in the GBM samples (|log2FC| < 1; P < 0.01; Figure 9F). These results indicated that the regulatory axis of lncRNA GABPB1-AS1/hsa-miR-519c-5p/CASP3/TP63 plays an important role in the expression of GBM.

GSEA enrichment

The GSEA results showed that CASP3 and TP63 were enriched on ECM-related signaling pathways (e.g. ECM receptor interaction, adherens junction, and TGFβ signaling pathways) and immune-related signaling pathways (e.g. Toll-like receptor and NOD-like receptor signaling pathways; Figure 9). This result showed that the two pyroptosis-related mRNAs are closely related to the growth, metastasis, and diffusion of GBM, and they play important roles in the immunoregulation of GBM TME.

Pyroptosis is closely related to malignant tumor [10], and TME plays an important role in the progress, metastasis, and treatment resistance of GBM [14]. Most studies only focused on a single PRG or TME cell type, and antitumor effect develops through interactions among tumor suppressor factors and system coordination [17]. Therefore, the effects of several PRGs on immunocyte infiltration characterization in TME were discussed by combining the ceRNA network, which is useful in elucidating the potential mechanism of initiation and development of GBM and prediction of the responses of immunotherapeutic approaches. Hence it can guide immunotherapy more effectively in clinics.

This study elaborated the differential expression of 49 PRGs in GBM samples and normal tissues and disclosed potential signaling pathways. The expression levels of NLRP2, NLRP7, TNF, IL1B, IL6, and NLRP1 in GBM were downregulated relative to those in normal tissues, whereas the expression levels of PLCG1, CASP9, CHMP6, SCAF11, and CHMP4A were upregulated. According to GSEA results, the differential genes mainly participated in pyroptosis, cytokine production, NOD-like receptor signaling pathway, necroptosis, and other signaling pathways.

Subsequently, a prediction gene model based on five prognosis PRGs (CASP3, NLRP2, TP63, GZMB, and CASP9) was constructed by using Lasso-Cox regression analysis. The results demonstrated that the five prognosis PRGs had good prognosis prediction ability for GBM. According to the TIMER and TISIDB database, the correlations between prognosis PRGs and abundance of immunocyte infiltration in GBM were evaluated. The expression levels of the five prognosis PRGs were closely related to CD8+ T cells, neutrophils, macrophages, dendritic cells, and CD4+ T cells, as well as immunoregulators. In addition, a prognosis risk score formula was established using the pyroptosis-related prognostic model. Patients with GBM were divided into low-risk and high-risk groups. The Kaplan-Meier survival curve showed that the low-risk group had obvious survival advantages. Notably, the high-risk and low-risk groups showed significant differences in immunocyte infiltration. The high-risk group presented relatively high distribution levels of monocytes, activated dendritic cells, activated mast cells, and eosinophils but relatively low distribution levels of resting NK cells and resting mast cells. Most of the non-neoplastic population of GBM was composed of infiltrating immune cells, whereas local inflammatory TME promoted tumor aggressiveness and the drug resistance of GBM [20]. For example, some study demonstrated that the frequency of M2 macrophages or microglia in the TME of recurrent GBM increased [21], and dendritic cells presented tumor cell peptides in the GBM, leading to cytotoxic T cells response and secretion of proinflammatory cytokines [20]. NK cells are activated by IL-12 and can kill GBM cells [22]. Moreover, the effectiveness of NK cells in inhibiting systemic metastasis of GBM in xenograft mouse models was reported [23]. Hence, TME immunocyte infiltration under different risk scores had different features and different tumor immune-reactions. The infiltrating immunocytes were closely related with the growth, invasion, and metastasis of GBM.

To elaborate the potential molecular mechanisms of five prognosis PRGs, a ceRNA network was built by searching the database and overlapping prediction results, which was used in recognizing the lncRNA GABPB1-AS1/hsa-miR-519c-5p/CASP3/TP63 regulatory axis. LncRNA is a type of noncoding RNAs with a length of more than 200 nucleotides [24]. Changes in lncRNA expression and its mutation can promote the occurrence and metastasis of tumors [25]. lncRNA plays a key role in TME intracellular signal transduction [26]. GABPB1-AS1 is an lncRNA in the cytoplasm. Recently, many studies have pointed out that lncRNA is a poor prognosis marker of glioma [27], cervical cancer [28], breast cancer [29], and prostate cancer [30]. According to this study, significant differences in GABPB1-AS1 expression was found between tumor and normal tissues. Similarly, Li et al. [31] found high expression of GABPB1-AS1 in the glioma tissues, and in vitro and in vivo experiments have demonstrated that GABPB1-AS1 knockdown reduced the proliferation and invasiveness of glioma cells. According to our results, Hsa-miR-519c-5p can be regulated by GABPB1-AS1 specifically, thus regulating the transcriptional levels of CASP3 and TP63. On this basis, Hsa-miR-519c-5p can control ECM-related signaling pathways (e.g., ECM receptor interaction, adherens junction, and TGFβ signaling pathways) and immune-related signaling pathways (e.g., Toll-like receptor and NOD-like receptor signaling pathways). This demonstrated that ECM and immune inflammation reactions played the key role in occurrence, invasion and metastasis of GBM. The TGF-β in the tumor tissues can inhibit immune cells, decrease local inflammation, and promote the reconstruction of ECM [32]. ECM components and high expression levels on the corresponding cell surface receptors are closely related with the progression of GBM [33]. Toll-like receptors (TLRs) can recognize non-self-molecules and activate the inflammation process. TLR expression had been observed in various samples, such as GBM cases and GBM cell lines, indicating that TLR plays an important role in tumor invasion and metastasis [34]. Therefore, we speculated that the mechanism by which GABPB1-AS1 regulates the transcriptional levels of CASP3 and TP63 and relevant signaling pathways through Hsa-miR-519c-5p might play a crucial role in the occurrence and development of GBM.

In conclusion, pyroptosis and TME-related independent prognostic markers were screened, and a prognosis risk score formula was established for the first time on the basis of prognosis PRGs. TME immunocyte infiltration characterization and immune response were assessed with ssGSEA, CIBERSORT algorithm, TIMER, and TISIDB database. A ceRNA network was then constructed. This study not only laid the foundation for the further exploration of pyroptosis and TME for the improvement of GBM prognosis but also provided novel insights for clinical immunotherapy and development of novel immunotherapeutic drugs. Deep experimental and clinical studies are still needed to elucidate the molecular mechanism of GBM and clinical applications of biomarkers and immunotherapy.

Acknowledgements

Not applicable.

Authors’ contributions

HMA and BH designed the study, coordinated technical support and funding. BH revised the manuscript. MRD performed the study and drafted the manuscript. YJQ, XP, JFC, MXZ and TZ participated the study. All authors read and approved the final manuscript.

Funding

This study was supported by Medical Innovation Research Program of Shanghai Municipality (grant number 21Y11920900), and Scientific and Technological Innovation Projects of Longhua Hospital (grant number CX202052).

Availability of data and materials

The datasets used during this study are available from the corresponding author on reasonable request.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no conflict of interest.

Ohgaki H, Kleihues P. The definition of primary and secondary glioblastoma. Clin Cancer Res. 2013;19(4):764–72.
Lucki NC, Villa GR, Vergani N, Bollong MJ, Beyer BA, Lee JW, Anglin JL, Spangenberg SH, Chin EN, Sharma A, Johnson K, Sander PN, Gordon P, Skirboll SL, Wurdak H, Schultz PG, Mischel PS, Lairson LL. A cell type-selective apoptosis-inducing small molecule for the treatment of brain cancer. Proc Natl Acad Sci USA. 2019;116(13):6435–6440.
Hu J, Xiao Q, Dong M, Guo D, Wu X, Wang B. Glioblastoma Immunotherapy Targeting the Innate Immune Checkpoint CD47-SIRPα Axis. Front Immunol. 2020;11:593219.
Bakas S, Akbari H, Pisapia J, Martinez-Lage M, Rozycki M, Rathore S, Dahmane N, O'Rourke DM, Davatzikos C. In Vivo Detection of EGFRvIII in Glioblastoma via Perfusion Magnetic Resonance Imaging Signature Consistent with Deep Peritumoral Infiltration: The φ-Index. Clin Cancer Res. 2017;23(16):4724–4734.
Fang Y, Tian S, Pan Y, Li W, Wang Q, Tang Y, Yu T, Wu X, Shi Y, Ma P, Shu Y. Pyroptosis: A new frontier in cancer. Biomed Pharmacother. 2020;121:109595.
Wang L, Qin X, Liang J, Ge P. Induction of Pyroptosis: A Promising Strategy for Cancer Treatment. Front Oncol. 2021;11:635774.
Tsuchiya K. Switching from Apoptosis to Pyroptosis: Gasdermin-Elicited Inflammation and Antitumor Immunity. Int J Mol Sci. 2021;22(1):426.
Kong Y, Feng Z, Chen A, Qi Q, Han M, Wang S, Zhang Y, Zhang X, Yang N, Wang J, Huang B, Zhang Q, Xiang G, Li W, Zhang D, Wang J, Li X. The Natural Flavonoid Galangin Elicits Apoptosis, Pyroptosis, and Autophagy in Glioblastoma. Front Oncol. 2019;9:942.
Ren LW, Li W, Zheng XJ, Liu JY, Yang YH, Li S, Zhang S, Fu WQ, Xiao B, Wang JH, Du GH. Benzimidazoles induce concurrent apoptosis and pyroptosis of human glioblastoma cells via arresting cell cycle. Acta Pharmacol Sin. 2021. <uvertical-align:super;>doi</uvertical-align:super;><uvertical-align:super;>:</uvertical-align:super;><uvertical-align:super;> </uvertical-align:super;><background-color:#CFBFB1;uvertical-align:super;>10.1038/s41401-021-00752-y</background-color:#CFBFB1;uvertical-align:super;><uvertical-align:super;>.</uvertical-align:super;>
Jiang M, Qi L, Li L, Li Y. The caspase-3/GSDME signal pathway as a switch between apoptosis and pyroptosis in cancer. Cell Death Discov. 2020;6:112.
Hinshaw DC, Shevde LA. The Tumor Microenvironment Innately Modulates Cancer Progression. Cancer Res. 2019;79(18):4557–4566.
Domingues P, González-Tablas M, Otero Á, Pascual D, Miranda D, Ruiz L, Sousa P, Ciudad J, Gonçalves JM, Lopes MC, Orfao A, Tabernero MD. Tumor infiltrating immune cells in gliomas and meningiomas. Brain Behav Immun. 2016;53:1–15.
Tsuchiya K. Switching from Apoptosis to Pyroptosis: Gasdermin-Elicited Inflammation and Antitumor Immunity. Int J Mol Sci. 2021;22(1):426.
Da Ros M, De Gregorio V, Iorio AL, Giunti L, Guidi M, de Martino M, Genitori L, Sardi I. Glioblastoma Chemoresistance: The Double Play by Microenvironment and Blood-Brain Barrier. Int J Mol Sci. 2018;19(10):2879.
Xiao K, Tan J, Yuan J, Peng G, Long W, Su J, Xiao Y, Xiao Q, Wu C, Qin C, Hu L, Liu K, Liu S, Zhou H, Ning Y, Ding X, Liu Q. Prognostic value and immune cell infiltration of hypoxic phenotype-related gene signatures in glioblastoma microenvironment. J Cell Mol Med. 2020;24(22):13235–13247.
Huang R, Meng T, Chen R, Yan P, Zhang J, Hu P, Zhu X, Yin H, Song D, Huang Z. The construction and analysis of tumor-infiltrating immune cell and ceRNA networks in recurrent soft tissue sarcoma. Aging (Albany NY). 2019;11(22):10116–10143.
Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. 2020;19(1):53.
Cancer Genome Atlas Research Network. Comprehensive genomic characterization defines human glioblastoma genes and core pathways. Nature. 2008;455(7216):1061–8.
GTEx Consortium. The Genotype-Tissue Expression (GTEx) project. Nat Genet. 2013;45(6):580–5.
DeCordova S, Shastri A, Tsolaki AG, Yasmin H, Klein L, Singh SK, Kishore U. Molecular Heterogeneity and Immunosuppressive Microenvironment in Glioblastoma. Front Immunol. 2020;11:1402.
Wang Q, Hu B, Hu X, Kim H, Squatrito M, Scarpace L, deCarvalho AC, Lyu S, Li P, Li Y, Barthel F, Cho HJ, Lin YH, Satani N, Martinez-Ledesma E, Zheng S, Chang E, Sauvé CG, Olar A, Lan ZD, Finocchiaro G, Phillips JJ, Berger MS, Gabrusiewicz KR, Wang G, Eskilsson E, Hu J, Mikkelsen T, DePinho RA, Muller F, Heimberger AB, Sulman EP, Nam DH, Verhaak RGW. Tumor Evolution of Glioma-Intrinsic Gene Expression Subtypes Associates with Immunological Changes in the Microenvironment. Cancer Cell. 2017;32(1):42-56.e6.
Castriconi R, Daga A, Dondero A, Zona G, Poliani PL, Melotti A, Griffero F, Marubbi D, Spaziante R, Bellora F, Moretta L, Moretta A, Corte G, Bottino C. NK cells recognize and kill human glioblastoma cells with stem cell-like properties. J Immunol. 2009;182(6):3530–9.
Lee SJ, Kang WY, Yoon Y, Jin JY, Song HJ, Her JH, Kang SM, Hwang YK, Kang KJ, Joo KM, Nam DH. Natural killer (NK) cells inhibit systemic metastasis of glioblastoma cells and have therapeutic effects against glioblastomas in the brain. BMC Cancer. 2015;15:1011.
Wang J, Su Z, Lu S, Fu W, Liu Z, Jiang X, Tai S. LncRNA HOXA-AS2 and its molecular mechanisms in human cancer. Clin Chim Acta. 2018;485:229–233.
Bhan A, Soleimani M, Mandal SS. Long Noncoding RNA and Cancer: A New Paradigm. Cancer Res. 2017;77(15):3965–3981.
Botti G, Scognamiglio G, Aquino G, Liguori G, Cantile M. LncRNA HOTAIR in Tumor Microenvironment: What Role? Int J Mol Sci. 2019;20(9):2279.
Luan F, Chen W, Chen M, Yan J, Chen H, Yu H, Liu T, Mo L. An autophagy-related long non-coding RNA signature for glioma. FEBS Open Bio. 2019;9(4):653–667.
Ou R, Lv M, Liu X, Lv J, Zhao J, Zhao Y, Li X, Li W, Zhao L, Li J, Ren Y, Xu Y. HPV16 E6 oncoprotein-induced upregulation of lncRNA GABPB1-AS1 facilitates cervical cancer progression by regulating miR-519e-5p/Notch2 axis. FASEB J. 2020;34(10):13211–13223.
Suvanto M, Beesley J, Blomqvist C, Chenevix-Trench G, Khan S, Nevanlinna H. SNPs in lncRNA Regions and Breast Cancer Risk. Front Genet. 2020;11:550.
Alkhateeb A, Rezaeian I, Singireddy S, Cavallo-Medved D, Porter LA, Rueda L. Transcriptomics Signature from Next-Generation Sequencing Data Reveals New Transcriptomic Biomarkers Related to Prostate Cancer. Cancer Inform. 2019;18:1176935119835522.
Li X, Wang H. Long Non-Coding RNA GABPB1-AS1 Augments Malignancy of Glioma Cells by Sequestering MicroRNA-330 and Reinforcing the ZNF367/Cell Cycle Signaling Pathway. Neuropsychiatr Dis Treat. 2021;17:2073–2087.
Bryukhovetskiy I, Shevchenko V, Arnotskaya N, Kushnir T, Pak O, Victor Z, Zaitsev S, Khotimchenko Y, Bryukhovetskiy A, Sharma A, Sharma HS. Transforming growth factor-β mimics the key proteome properties of CD133- differentiated and CD133+ cancer stem cells in glioblastoma. Int Rev Neurobiol. 2020;151:219–242.
Xiao W, Wang S, Zhang R, Sohrabi A, Yu Q, Liu S, Ehsanipour A, Liang J, Bierman RD, Nathanson DA, Seidlits SK. Bioengineered scaffolds for 3D culture demonstrate extracellular matrix-mediated mechanisms of chemotherapy resistance in glioblastoma. Matrix Biol. 2020;85-86:128–146.
Moretti IF, Franco DG, de Almeida Galatro TF, Oba-Shinjo SM, Marie SKN. Plasmatic membrane toll-like receptor expressions in human astrocytomas. PLoS One. 2018;13(6):e0199211.

Download PDF

Editor assigned by journal
20 Nov, 2021
Submission checks completed at journal
19 Nov, 2021
Editor invited by journal
19 Nov, 2021
First submitted to journal
18 Nov, 2021

You are reading this latest preprint version

Pyroptosis-related Prognosis Model, Immunocyte Infiltration Characterization, and Competing Endogenous RNA Network of Glioblastoma

Status:

Version 1

Abstract

Figures

Background

Materials And Methods

Datasets and preprocessing

Identification of differentially expressed PRGs

Analysis of the somatic mutation characteristics of PRGs

Functional enrichment analysis

Establishment and validation of the pyroptosis-related gene prognostic model

Independent prognostic analysis and clinical value of the risk model

Immunocyte infiltration analysis

CeRNA network construction

Gene Set Enrichment Analysis

Results

Defining the expression of PRGs in GBM

Analysis of the somatic mutation characteristics of PRGs

Functional enrichment analysis of the PRGs

Construction of the prognostic PRG model

Independent prognostic analysis and nomogram construction based on risk model

Correlation between the expression of prognostic PRGs and immunocyte levels and immunoregulators in GBM

Immunocyte infiltration between the high-risk and low-risk groups

Construction of the network of mRNA–miRNA–lncRNA

GSEA enrichment

Discussion

Conclusions

Declarations

Acknowledgements

Authors’ contributions

Funding

Availability of data and materials

Ethics approval and consent to participate

Consent for publication

Competing interests

References

Status:

Version 1