A Newly Defined Pyroptosis-Related Gene Signature for the Prognosis of Bladder Cancer

Introduction

Bladder cancer (BC), as the most Common urinary system tumor type and the main cause of tumor-related death, causes 165 thousand deaths every year and leads to high morbidity and mortality in China.¹ Kinds of intrinsic and extrinsic factors could affect the occurrence and progression of BC (ethnic group, geographical area, treatment modalities, cell type of the bladder cancer, etc). The most common risk factors contain effects of environmental carcinogens, occupational exposure and smoking.² Although a large percentage of BC is non-muscle invasive bladder cancer (NMIBC) at the first diagnosis, the percent of reoccurrence is almost 30%, needing transurethral resection of bladder tumor for a second time.³ A big part of NMIBC patients will finally change into invasive disease.⁴ Neoadjuvant chemotherapy with following radical cystectomy is a standard treatment for muscle invasive bladder cancer (MIBC). However, the OS rate of 5 years is unsatisfactory (38.6%) among cT3-4aN0 patients after above treatments.⁵ In view of the limitations of BC therapy, new therapeutic targets are needed to improve the survival outcome of BC; therefore, there is an urgent need for reliable new prognostic models to make targeted therapy more flexible.

Pyroptosis, as a new form of programmed cell death, is known as cellular inflammatory necrosis.^6,7 The gasdermin family is the main executor of pyroptosis.^8–12 Pyroptosis was firstly found to be an important mechanism for combating infection. Afterwards, more and more researches found that it was related to the occurrence and development of tumors. It has been published that proinflammatory cytokines, gasdermin proteins, and inflammatory vesicles, which are significant components of pyroptosis, are related with invasion, oncogenesis, and metastasis.¹³ Dupaul-Chicoine et al¹³ found that the transgenic mice knocked out inflammatory vesicle-related genes (CASP1 and NLRP3) were more probably to develop colon cancers. Besides, the process of pyroptosis always be accompanied by inflammatory response and immune system activation, which is unlike apoptosis.¹⁴ Some studies have found that the potent pro-inflammatory effect of pyroptosis plays an important role in the regulation of tumor immune microenvironment. In addition, some studies reported that pyroptosis also plays an important role in the activity of CD8+T lymphocytes¹⁵ and anti-tumor function of NK cells.¹⁶

Taking into account the existing research results, we found that pyroptosis plays a vital role in tumor and anti-tumor process. However, its specific functions in BC have been less studied. Therefore, we conducted this study in order to clarify the expression levels of pyroptosis-related genes between normal and bladder tumor tissues, explore the potential prognostic value of these genes, and try to discover the correlation between them and the tumor immune microenvironment.

Materials and Methods

Datasets

The RNA sequencing (RNA-seq) data and clinical information of 414 BC patients and 19 normal human bladder samples were downloaded from the TCGA database on 29 July 2021 (https://portal.gdc.cancer.gov/repository). We obtained the RNA-seq data and the clinical information from TCGA database and GEO database (https://www.ncbi.nlm.nih.gov/geo/, ID: GSE31684).

Distinguish the Differentially Expressed Genes Related to Pyroptosis

Fifty-two pyroptosis-related genes from prior reviews were extracted,^17–20 and the list of full names is presented in the supplementary file (Table S1). We used “limma” package to distinguish DEGs with a P value <0.05. Then, we constructed a DEGs-related PPI network by Search Tool (STRING, version 11.5, https://string-db.org/) for the Retrieval of Interacting Genes.

Development and Validation of the Prognostic Model Related to Pyroptosis Gene

We further took Cox regression analysis to excavate the relationships between each gene and survival status in the TCGA group. The cut-off P-value was 0.2, and 9 survival-related genes were selected. The LASSO Cox regression model was used by R package “glmnet” to improve the prognostic model. Finally, the eight genes and their coefficients were gotten. The risk score was calculated by “scale” function in R in the TCGA data. Then, The BC patients were divided into high- and low-risk groups according to the median risk score, and the Kaplan–Meier analysis of OS time, 1,3,5-year ROC curve, PCA and t-SNE was performed. Distribution and survival status of patients was showed by the “bioRiskPlot” function based on the risk score. For the validation studies, the same procedure operates in GEO data (GSE31684).

Independent Prognostic Analysis of the Risk Score

We got the clinical information (age, gender, stage, grade, TNM) of patients in the TCGA group. Then, we employed Univariate and multivariable Cox regression models for the analysis. Due to the lack of above clinical information in GEO data, Univariate and multivariable Cox regression analysis was not continued.

Functional Enrichment Analysis of the DEGs

BC patients were divided into two cohorts by the median risk score. The DEGs between the high- and low-risk groups were filtered (FDR < 0.05 and |log2FC| ≥1). Based on these DEGs, GO and KEGG analyses were performed by applying the “clusterProfiler” package. The “gsva” package was utilized to conduct the ssGSEA to calculate the scores of infiltrating immune cells and to evaluate the activity of immune-related pathways.

Statistical Analysis

Single-factor analysis of variance was used for comparing the gene expression levels, while the Pearson chi-square test was applied for categorical variables. The Kaplan–Meier method with a two-sided Log rank test, univariate and multivariate Cox regression models was employed for comparing the OS of patients between subgroups and assessing the independent prognostic value of the risk model. The Mann–Whitney test was used when comparing the immune cell infiltration and pathway activation between the two cohorts. All statistical analyses were accomplished with R software (v4.1.0).

Results

Identification of DEGs Between Normal and Tumor Tissues

The 53 gene expression levels related to pyroptosis were compared in TCGA data from 414 tumor and 19 normal tissues, and finally we found 29 DEGs (P < 0.05). Among these genes, 6 genes (CHMP7, NLRP1, CHMP3, NLRP3, IL6 and ELANE) were downregulated, and 23 genes (CASP5, AIM2, NLRP2, IL1A, NLRP7, CHMP4C, CASP6, GSDMB, BAX, TP63, PYCARD, BAK1, CASP3, CASP8, GSDMD, TP53, PLCG1, GPX4, CHMP2A, HMGB1, CHMP4B, CHMP4A and CYCS) were enriched in the tumor cohort. The heatmaps of RNA levels and the PPI analysis are shown in Figures 1A and B. The minimum required interaction score was 0.4, and we determined that AIM2, BAK1, BAX, CASP3, CASP5, CASP6, CASP8, CHMP2A, CHMP3, CYCS, ELANE, GSDMD, HMGB1, IL1A and IL6 were hub genes. The correlation network including all pyroptosis-related genes is showed in Figure 1C.

Figure 1 Expressions of the 29 pyroptosis-related genes and the interactions among them. (A) Heatmap of the pyroptosis-related genes between the normal and the tumor tissues. (B) PPI network showing the interactions of the pyroptosis-related genes (interaction score = 0.4). (C) The correlation network of the pyroptosis-related genes. P values were showed as: **P < 0.01; ***P < 0.001.

Tumor Classification Based on the DEGs

We used a consensus clustering analysis with 414 BC patients in the TCGA cohort for exploring the connections and we found that the intragroup correlations were the highest and the inter-group correlations were low when k = 2 (Figure 2A). The gene expression and the clinical information including age (≤65 or >65 years), gender (female or male), the degree of tumor differentiation (high grade, low grade or unknown), T stage (T0,1,2,3,4 or unknown), N stage (N0, 1, 2, 3 or unknown) and M stage (M0,1 or unknown) were showed in a heatmap, and we found there was little differences in clinical features, except gender (Figure 2B). The OS was also compared between the two cohorts, but no significant differences were discovered (P = 0.293, Figure 2C).

Figure 2 Tumor classification based on the pyroptosis-related DEGs. (A) 414BC patients were divided into two clusters by the consensus clustering matrix (k = 2). (B) Heatmap and the clinicopathologic characters of the two clusters classified by these DEGs. (C) Kaplan–Meier curves for the OS of these two clusters.

Development of a Prognostic Gene Model in the TCGA Cohort

We used univariate Cox regression analysis for survival-related genes. The 9 genes (BAK1, CASP1, GZMB, IRF1, TP63, AIM2, NOD 2, TNF and GZMA) were selected for further analysis (P < 0.2) (Figure 3A). After LASSO Cox regression analysis, we constructed an 8-gene signature (Figure 3B and C). The result was showed as follows: risk score = (−0.010*BAK1 exp.) + (−0.012*GZMB exp.) + (−0.073*IRF1 exp.) + (−0.063*TP63 exp.) + (−0.033*AIM2 exp.) + (0.233*NOD2 exp.) + (−0.114*TNF exp.) + (−0.032*GZMA exp.). 414 patients were equally divided into high-risk and low-risk subgroups, and a significant difference of OS was detected between these two groups (P < 0.001, Figure 3D). Time-dependent ROC analysis was used to evaluate the specificity and sensitivity of the prognosis model, and the AUC was 0.631, 0.626, and 0.631 for 1-, 3- and 5-year survival (Figure 3E). High-risk group had more patient deaths and shorter survival time than low-risk group (Figure 3F and G). The PCA and t-SNE (Figure 3H and I) showed that patients were well divided into two clusters.

Figure 3 Construction of risk signature in the TCGA cohort. (A) Univariate cox regression analysis of OS for selected pyroptosis-related gene with P < 0.2. (B) LASSO regression of the 9 OS-related genes. (C) Cross-validation for tuning the parameter selection in the LASSO regression. (D) Kaplan–Meier OS curves of patients in the high- and low-risk groups. (E) ROC curves demonstrated the predictive efficiency of the risk score. (F) Distribution of patients based on the risk score. (G) The survival status for each patient. (H) PCA plot for BCs based on the risk score. (I) t-SNE plot for BCs based on the risk score.

External Validation of the Risk Signature

Ninety-three BC patients from GEO cohort (GSE31684) were used as the validation cohort. Forty-nine patients in the GEO cohort were separated into the low-risk group, and 44 patients were separated into the high-risk group by the median risk score in the TCGA cohort. A significant difference in OS was found between the high- and low-risk groups (P = 0.012, Figure 4A). Time-dependent ROC analysis was used to assess the sensitivity and specificity of the prognosis model, and we detected that the area under the AUC was 0.562, 0.578 and 0.653 for 1-, 3- and 5-year survival (Figure 4B).

Figure 4 Validation of the risk model in the GEO cohort. (A) Kaplan–Meier curves for comparison of the OS between low- and high-risk groups. (B) Time-dependent ROC curves for BCs. (C) Distribution of patients in the GEO cohort based on the median risk score in the TCGA cohort. (D) The survival status for each patient (low-risk population: on the left side of the dotted line; high-risk population: on the right side of the dotted line). (E) PCA plot for BCs. (F) t-SNE plot for BCs based on the risk score.

High-risk group had more patient deaths and shorter survival time than low-risk group (Figure 4C and D). The PCA and t-SNE (Figure 4E and F) showed that patients were well divided into two clusters.

Independent Prognostic Value of the Risk Model

Univariate and multivariable Cox regression analyses were used to find the independent prognostic factor, and finally indicated that the risk score was an independent factor in both univariate and multivariate Cox regression analysis (HR = 5.214, 95% CI: 1.993–13.644, Figure 5A; HR = 3.629, 95% CI: 1.334–9.873, Figure 5B). Besides, we created a heatmap of clinical features for the TCGA cohort and detected that the M stage was significantly different between the two subgroups (P < 0.01, Figure 5E).

Figure 5 Univariate and multivariate Cox regression analyses for the risk score. (A) Univariate analysis for the TCGA cohort. (B) Multivariate analysis for the TCGA cohort. (C) Heatmap for the connections between clinicopathologic features and the risk groups. P values were showed as: **P < 0.01.

Functional Analyses Based on the Risk Model

One hundred and fourteen DEGs between the high and low groups in the TCGA cohort were identified after the use of “limma” R package. Among these genes, 20 genes were up-regulated in the high-risk group, and 94 genes were down-regulated (showed in the supplementary file (Table S2)). The results of GO enrichment analysis and Kyoto KEGG pathway analysis indicated that the DEGs play an important role in the immune response, response to virus and antigen processing and presentation (Figure 6A and B).

Figure 6 Functional analysis based on the DEGs between the two-risk groups in the TCGA cohort. (A) Bubble graph for GO enrichment (the bigger bubble means the more genes enriched, and the increasing depth of red means the differences were more obvious; q-value: the adjusted p-value). (B) Barplot graph for KEGG pathways (the longer bar means the more genes enriched, and the increasing depth of red means the differences were more obvious).

Comparison of the Immune Activity Between Subgroups

We employed the ssGSEA by comparing the enrichment scores of the activity of 13 immune-related pathways and 16 types of immune cells between the different groups in TCGA and GEO groups. Except for mast cells, the other immune cells showed lower levels of infiltration in the high-risk group in the TCGA cohort and GEO cohorts (Figures 7A and C). In the TCGA and GEO cohorts (Figures 7B and D), the high-risk subgroup generally had lower activity of all immune pathways, except the type-2 IFN response pathway.

Figure 7 Comparison of the ssGSEA scores for immune cells and immune pathways. (A–D) Comparison of the tumor immunity between low- and high-risk group in the TCGA and GEO cohort. *P < 0.05; **P< 0.01; ***P < 0.001.

Discussion

In the study, we researched the mRNA levels of 53 recently known genes related to pyroptosis in normal bladder and BC tissues and discovered that more than half of these genes were differentially expressed. Then, we used the consensus clustering analysis to divide BC patients into two clusters based on the DEGs, but there was no significant difference in the survival. In order to further explore the prognostic role of these genes, we created an 8-gene risk signature. The following functional analyses suggested that the DEGs between the high- and low-risk groups were bound up with immune response. We compared the immune cells and immune pathways in the low- and high-risk groups and discovered that the high-risk group had less infiltrating immune cells and less active immune-related pathways.

Pyroptosis is a new type of programmed death, which has now been found to have an important double-sided role in the development of tumors. On the one hand, it can stimulate normal cells to transform into tumor cells by releasing a large number of inflammatory factors. On the other hand, the process of pyroptosis could be a potential therapeutic target. In patients with bladder cancer, the correlation between these pyroptosis-genes and the relationship between them and the survival is still unclear. The study constructed a signature featuring 8 pyroptosis-related genes (BAK1, GZMB, IRF1, TP63, AIM2, NOD 2, TNF and GZMA) and found it could forecast survival in BC patients.

BAK1 is initially established as an ultimate death effector of apoptosis triggering cues. One study demonstrated that BAK1 could be directly targeted by miR-125b to mediate breast cancer progression.²¹ It showed that BAK1 could significantly enhance Taxol-induced cytotoxicity and apoptosis. Similarly, another study indicated that BAK1 suppression caused by specific miRNA could lead to androgen-independent growth of prostate cancer cells.²² However, BAK1 is regarded as an unfavorable prognostic factor in some studies for Hepatocellular Carcinoma and Head and Neck Squamous Cell Carcinoma.^23,24 In the study, BAK1 seems to be a cancer-promoting gene, because its expression in tumor tissues is significantly higher than that in normal tissues; however, it helps to prolong patient survival for that it is enriched in the low-risk group. Given that the data provided by BC is limited and the results of different tumors are often contradictory, our results on BAK1 provide some insights for further research.

Granzyme B (GZMB) is normally produced as a toxic granule-secreted enzyme with killing activity in the tumor microenvironment by natural killer cells and cytotoxic T lymphocytes.^25,26 But some reports have reported that low GZMB expression may be related to early CRC metastasis.²⁷ Besides, GZMB can cut downstream caspase-3 and bid to activate them with hydrolytic activity as a serine protease, resulting in apoptosis of target-cells.^28,29 Moreover, a study detected that, in Colorectal cancer cells, LINC02474 was associated with the transcriptional repression of GZMB determined by the dual-luciferase reporter assay and then interference of GZMB could increase the metastatic abilities of CRC cells by reducing apoptosis. In our study, we found that the GZMB gene is enriched in the low-risk group, which is an indicator of good prognosis. However, more studies are demanded to determine how this gene is involved in tumor suppression in the pyroptosis program.

Unlike GZMB, the actual function and mechanism of action of GZMA (another member of the granzyme family) remain unclear, likely because of the lack of sufficiently relevant experimental systems. A recent study discovered that GZMA from cytotoxic lymphocytes cleaves and activates GSDMB to induce target cell pyroptosis.³⁰ This immune effector mechanism promotes CTL-mediated tumor clearance in mice. We found that the high expression of GZMA is associated with better survival outcomes. Therefore, whether its mechanism in BC is the same as the above research is worthy of our follow-up study.

AIM2 could use ASC-mediated junctional proteins to activate CASP-1 and to promote pyroptosis.³¹ AIM2 was considered to be a tumor suppressor for that its mutation or inactivation was discovered in kinds of tumors, including gastric, endometrial and colon cancers, but it was different in non-small-cell lung, nasopharyngeal, and oral cancer.³² Therefore, AIM2 may play a different role in different kinds of cancers. In this study, AIM 2 is likely a cancer-promoting gene and may be associated with a good prognosis. We hope that this discovery can pave the way for subsequent research on bladder cancer.

Interferon regulatory factor 1 (IRF1), as a transcription factor induced in response to interferons,³³ plays crucial roles in various physiological and pathological processes, including inflammation, microbial infection, tumor immune surveillance, and innate and adaptive immune responses.^34,35 IRF1 controls the gene expression of guanylate-binding proteins, the inducible nitric oxide synthase, and caspase-1 (CASP1),^34,36 which have been implicated in various inflammatory diseases. Besides, IRF1 was reported that it could prevent tumorigenesis in a spontaneous mouse model of colorectal cancer due to defective pyroptosis, apoptosis, and necroptosis.³⁷ More studies are needed focusing on how/whether IRF1 takes part in tumor suppression and pyroptosis.

TP63 belongs to the TP53 family of tumor suppressor genes,³⁸ and these genes induce cell cycle arrest and apoptosis.^39,41 In particular, TP63, which encodes 2 subtypes of TAP63 and carcinogenic deltaNp63 subtypes,⁴⁰ 14 has been involved in many cancers such as bladder, uterine, and breast cancer, as well as squamous cell carcinoma of the head and neck.^42–45 Somerville et al⁴⁶ found that the deltaNp63 subtype-driven enhancer reprogramming enhances tumor growth and the invasive ability of cells, promoting epithelial differentiation and metastasis of cancer cells. Moreover, they showed that the continuous expression of TP63 is important for the growth of cancer epithelial cells. However, the correlation between TP63 and cell pyroptosis or tumor prognosis is less researched, and it is worthy of our continued exploration.

The TNF family, as one of the best-studied protein families, includes the 29 TNF receptor superfamily and 19 TNF ligands superfamily.⁴⁷ Some studies have found that TNF can through TNFR1 and TNFR2 and then activate immune cells, as an inflammatory regulator.⁴⁸ Moreover, some studies have detected that many cancers are closely associated with the TNF family.⁴⁹ In our study, we found that TNF is associated with the prognosis of patients with bladder cancer, and its high expression means a better prognosis. However, further experimentally studies are needed for the specific mechanism.

NOD2 belongs to the NOD like receptor family and can interact with many proteins and regulate immune responses. A study found that activated NOD2 can recruit receptor-interacting protein 2, and conducts signal through nuclear factor (NF)-κB and MAPK pathway.^50,51 Moreover, NOD2 can promote myocardial apoptosis⁵² and promote renal injury in diabetic nephropathy.⁵³ A recent study found that LPS can induce Extracellular histone H3 and then cause pyroptosis via VSIG4/NLRP3 and NOD2 pathway during sepsis. However, the mechanism of tumor-related pyroptosis induced by NOD2 has rarely been studied. Our study found that high expression of NOD2 is associated with a worse prognosis. This may remind us that NOD2 may affect the prognosis of cancer patients by mediating pyroptosis.

In summary, all eight of the above-mentioned pyroptosis-related genes have some positive or negative effects on the prognosis of bladder cancer patients in our study. The interrelationship of these genes in the pyroptosis process is worthy of our further study.

Nowadays, the research on pyroptosis is not thorough enough. Pyroptosis, as one of the processes of programmed death, may co-occur with apoptosis and other mechanisms that promote cell death and interlock with each other.⁵⁴ For example, BAK1 and GZMB are the key regulator in the apoptosis mechanism. Different from apoptosis, pyroptosis could lead to the release of cell contents and then activates a strong inflammatory response, which can also be regarded as an important natural immune response of the body. Through the analysis of patients in different risk groups, we found that DEGs mainly covers cellular immunity and inflammation. Through the following GO and KEGG analysis, we can speculate that pyroptosis is related to the cellular immune microenvironment.

In addition, we compared the levels of key antitumor infiltrating immune cells between high- and low-risk cohorts from TCGA and GEO datasets, and finally found an overall suppression of immune functions in the high-risk group. Previous studies reported that treg cells could suppress antitumor immunity and are associated with poor clinical outcomes. In the study, we found that treg cells are enriched in a low-risk group. The reason for different conclusions may be associated with GZMB. GZMB has been detected to be expressed by human and mouse Treg cells and suppress the immune function and cancer immune-evasion by inhibiting NK and Tc cell responses to control viral-induced lung inflammation or in tumor models.^55–58 Besides, it may be related to different Treg cell subtypes.⁵⁹ Through these discoveries, the poor survival outcome of high-risk BCs may be generated by decreased levels of antitumor immunity, and the specific mechanism needs our future research.

There are also some limitations with our study. This is a preliminary research for the prognostic value of pyroptosis-related genes, and we hope to provide some theoretical basis for future research. Because of the lack of relevant research, we cannot confirm whether the above regulators play similar roles in pyroptosis pathways in bladder cancer patients, so in-depth researches are needed to confirm this question. In addition, there were potentially important confounders, such as treatment measures, for which data were not available in our database. Therefore, we hope that further prospective studies could be conducted to confirm our findings. Last but not least, we believe that the real-life research by our institution would contribute and validate these results in the future.

In conclusions, we discovered that pyroptosis is related to BC closely, because a big part of the pyroptosis-related genes were differently expressed between normal and BC tissues. Besides, the score based on 8 pyroptosis-related genes generated from our risk signature was an independent risk factor for predicting OS. The DEGs between the high and low-risk groups were connected to tumor immunity. Our study offers a new gene signature for predicting the prognosis of BC patients and provides an important foundation for future researches of the relations between pyroptosis-related genes and immunity in BC.