1 Introduction

Bladder cancer (BLCA) is the tenth most common cancer globally, divided into non-muscle-invasive (NMIBC) and muscle-invasive (MIBC) types¹. NMIBC patients often have high recurrence after surgery, so postoperative treatment like bladder instillation with BCG or chemotherapy is recommended¹. MIBC is treated with radical cystectomy, often combined with platinum-based neoadjuvant chemotherapy for better results². New immunotherapy using PD-1/PD-L1 inhibitors is showing promise for BLCA. In 2017, the FDA approved atezolizumab and pembrolizumab for advanced cases intolerant to platinum-containing chemotherapy³. Response rates of 11.68% and 24.05% were seen with atezolizumab and pembrolizumab^4,5. Chemotherapy can affect the tumor immune environment⁶, and ongoing trials are exploring combining immune checkpoint inhibitors (ICIs) with chemotherapy (ClinicalTrials.gov Identifier: NCT04383743, NCT04630730)⁷. Research in BLCA immunotherapy is in early stages, highlighting the need for innovative approaches combining chemotherapy and immunotherapy.

Immunogenic cell death (ICD) is a regulated cell death that, in an active immune system, triggers an adaptive immune response by exposing antigens from dying cells⁸. Certain chemotherapy drugs, like cisplatin⁹ and gemcitabine¹⁰, can induce ICD. While studied in preclinical BLCA models, the potential benefits of ICD-based therapies in BLCA lack conclusive evidence¹¹. It is crucial to explore the impact of ICD in clinical settings for BLCA treatment, be it through chemotherapy or immunotherapy. In this study, we discovered a connection between ICD-related genes and the prognosis and immune infiltration in BLCA patients using TCGA-BLCA, GEO datasets, and tissue microarray staining. Our validated risk-scoring model effectively evaluates immune infiltration, prognosis, immunotherapy response, and drug sensitivity in BLCA, providing guidance for personalized treatment and future research.

2 Materials and Methods

2.1 Identification of Differentially Expressed ICD-related genes

ICDs were sourced from a previous article¹². We found differentially expressed genes (DEGs) in TCGA-BLCA using the “Limma” R package. This involved comparing gene expression in normal and tumor samples. DEGs interacting with ICDs were identified. Enrichment scores of ICDs were assessed using ssGSEA analysis with the “GSVA” R package. Patients were then grouped as “high” or “low” based on these scores. Finally, DEGs were visualized through volcano and heatmaps using R packages “pheatmap” and “ggplot2”.

2.2 Functional Enrichment Analysis

We compared ICD-high and ICD-low groups by analyzing Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) using the “clusterProfiler” R package. Gene Set Enrichment Analysis (GSEA) assessed differences in expressed gene sets between the groups using the MSigDB collection.

2.3 Analysis of Somatic Mutation

BLCA patient mutation data from TCGA were visualized using “Maftools” in waterfall plots. Boxplots summarized tumor mutational burden (TMB), and the 20 most frequently mutated genes were shown in waterfall plots.

2.4 Analysis of Immune Infiltration in the TME

Used the “ESTIMATE” R package to predict the tumor microenvironment (TME) scores for BLCA samples and assess immunocyte proportions. XCELL, TIMER, QUANTISEQ, MCPcounter, EPIC, CIBERSORT and CIBERSORT-ABS were used to explore immune cell infiltration, visualizing differences with a violin diagram and analyzing with the Wilcoxon signed-rank test. Evaluated immunotherapeutic response using the IPS algorithm and compared immune checkpoint and HLA-related gene expression in different groups. Also, accessed IPS data from The Cancer Immunome Atlas (TCIA) to compare responses to ICIs across BLCA patient groups.

2.5 Construction and validation of risk-scoring model

To explore the prognostic value of ICD-related genes in BLCA. Using TCGA-BLCA data, we conducted Cox regression analysis and LASSO Cox regression to build a prognostic model. Risk-score=CALR*(0.5378) + IFN1B*(-0.6349) + IFNG*(-0.2028) + IL1R1*(0.0919). Risk-scores were calculated, and patients were divided into low-risk and high-risk groups. We validated the model using the GSE13507 data and BLCA tissue arrays (HBlaU079Su01, Shanghai Outdo Biotech Co., Ltd). Kaplan-Meier curves for overall survival (OS) in low-risk and high-risk groups were created with the R packages “Survminer” and “Survival”. Nomogram was created using “rms” package to integrate risk-scores and other clinicopathological information to assess survival probability. ROC analysis was used to evaluate the nomogram’s ability to predict survival. The R package “pRRophetic” was used to evaluate the drug sensibility. The Remmele and Stegner’s semiquantitative immunoreactive score (IRS) scale was employed to assess the expression levels of each marker, as detailed in earlier study¹³. The R umap package was utilized to conduct the UMAP analyses of the individual cells. A two-sided P <0.05 was considered valuable. The reagents used are as follows: anti-IFNB1 (cat #: bs-0787R, BIOSS), anti-IL1R1 (cat #: bs-20697R, BIOSS), anti-IFNG (cat #: MAB48116, Bioswamp), and anti-CALR (cat #: A1066, ABclonal).

3 Results

3.1 Stratification and pathway enrichment analysis of TCGA-BLCA based on ICD-related genes

We conducted unsupervised clustering analysis on 34 ICD-related genes to identify two ICD-associated subtypes. Subsequently, we defined ICD-high and ICD-low groups based on the expression levels of these genes and compared their ICD gene expression profiles (Figure 1A). A clinical heatmap demonstrated differences between the two subtypes in terms of Grade (Figure 1B). Furthermore, survival analysis (Figure 1C) indicated that the ICD-high group exhibited significantly longer survival times compared to the ICD-low group (p=0.039). Subsequently, we found that patients in the ICD-high group exhibited a significantly higher TMB (p < 0.05) (Figure S1A). We also found that TP53, TTN, KMT2D, MUC16, and ARID1A are the most frequently mutated genes in both groups (Figure S1B, C).

Figure 1.

Next, we examined the differential gene expression between the two subgroups, identifying 4,321 DEGs, consisting of 2,177 downregulated and 2,144 upregulated genes. The heatmap displayed the top 50 upregulated and 50 downregulated genes (Figure 1D, E). To further identify pathways associated with immune activation specific to ICD-high and ICD-low groups, GO, KEGG, and GSEA were conducted. The results indicate a close association between elevated ICD expression and the activation of the immune microenvironment (Figure1F-L).

3.2 The immune characteristics within the TME in ICD-high and ICD-low group

We employed the ESTIMATE algorithm for analysis. The ICD-high group exhibited higher ESTIMATE scores, immune scores, and stromal scores compared to the ICD-low group, while the ICD-high group demonstrated lower tumor cell purity (Figure 2A-D). Further analysis of immune cell infiltration between the two subtypes revealed significantly higher levels of T cells CD8+, T cells CD4+ memory activated, and T cells follicular helper in the ICD-high group, while Tregs displayed the opposite trend (Figure 2E). Furthermore, we compared the expression of immune function-related pathways, immune checkpoint genes, and HLA-related genes between the two subtypes. The results indicate that the ICD-high group exhibited elevated expression in nearly all 13 immune pathways (Figure 2F). In the 8 immune checkpoint genes, PDCD1LG2, CTLA4, PDCD1, CD274, HAVCR2, LAG3, TIGIT had higher expression levels, while SIGLEC15 had lower expression (Figure 2G). Among the 24 HLA-related genes, the ICD-high group had higher expression levels for almost all genes (Figure 2H). These results suggest that ICD-high patients have a significantly higher quantity of tumor-infiltrating immune cells, but they appear to be in a state of functional exhaustion. Thus, we evaluated the effects of ICIs PD1 and CTLA4 on the ICD-high and ICD-low groups. Treated with PD1 and CTLA4 or PD1 alone, ICD-high group had better outcomes. With CTLA4 alone, both subtypes had similar results. Without PD1 or CTLA4 treatment, ICD-high group had significantly worse outcomes than ICD-low group (Figure 2I). Demonstrating that targeting PD-1 is crucial to enhance immune cell functionality and improve the prognosis of ICD-high patients.

Figure 2.

3.3 Construction and validation of risk-scoring model

In Cox univariate analysis, we identified 4 ICD-related genes significantly associated with OS of patients: CALR (p=0.003), IFNB1 (p=0.037), IFNG (p=0.022), IF1R1 (p=0.047) (Figure 3A). Subsequently, these 4 ICD-related genes were tested and selected in the LASSO regression analysis for constructing the predictive model (Figure 3B, C). Subsequently, patients are categorized into high-risk and low-risk groups according to the median (the boundary value is 4.372) of the risk-scores. Survival analysis indicated that, in both the TCGA-BLCA dataset and the GSE13507 dataset, the high-risk group exhibits significantly shorter survival times than the low-risk group (P<0.05) (Figure 3D, E). Next, we analyzed somatic mutations and generated waterfall plots for the top 20 mutated genes (Figure S2A, B) and calculated the TMB for both groups. Among the ICD-high-risk group and ICD-low-risk group, TP53, TTN, KMT2D, MUC16, and ARID1A were the most frequently mutated genes, but the TMB does not differ between the ICD high-risk group and the ICD low-risk group (Figure S2C). Survival analysis also indicated longer survival times in the high TMB group (P<0.001), with the best prognosis observed in the low-risk + high-TMB group and the worst prognosis in the high-risk + low-TMB group (Figure S2D, E).

Figure 3.

We further conducted a comprehensive analysis using the TCGA-BLCA dataset (Figure 3F-H) as the training set and the GSE13507 dataset (Figure 3I-K) as the validation set. Based on the median of the risk-scores, we categorized the datasets from TCGA and GEO into high-risk and low-risk groups. We analyzed the correlation between risk-scores and survival status, showing that higher risk-scores increased the risk of patient mortality. Additionally, we examined the expression of the four genes in the risk signature, finding higher expression of CALR and IL1R1 in the high-risk group, and higher expression of IFNB1 and IFNG in the low-risk group (Figure 3F-K). Overall, these results suggest a less favorable prognosis for high-risk patients.

3.4 The association between risk-scores and clinical features

We analyzed the association between the four genes composing the risk-score and patient prognosis (Figure 4A-D). Patients with low expression of CALR, low expression of IF1R1, high expression of IFNB1, and high expression of IFNG had better prognosis, which is consistent with CALR and IF1R1 being lowly expressed and IFNB1 and IFNG being highly expressed in the low-risk group (Figure 3H, K). The clinical heatmap demonstrates differences between the two subtypes in terms of Stage and Grade (Figure 4E). In the high-risk group, there is a higher proportion of Stage III and Stage IV cases, while Stage I and Stage II cases are less common (Figure 4F). Additionally, the high-risk group has a higher proportion of high-grade cases and a lower proportion of low-grade cases (Figure 4G).

Figure 4.

To aid clinicians in making more intuitive predictions of individual patient survival probabilities at 1, 3, and 5 years, we created a nomogram that incorporated gender, grade, age, and risk-score (Figure 4H). The AUC values for the 1-year, 3-year, and 5-year ROC curves are 0.630, 0.635, and 0.651, respectively (Figure 4I). The AUC values for risk-score, age, gender, and stage are 0.630, 0.669, 0.488, and 0.637, respectively (Figure 4J). Additionally, we found that the risk-score can serve as an independent prognostic factor (p < 0.001) (Figure 4K, L). This suggests that the risk-score is an effective indicator for predicting patient prognosis.

We checked risk-scores and survival in different subgroups. Higher risk-scores were seen in advanced stages (high-grade, stage III, stage IV, and M1 cases), suggesting an association with advanced disease (Figure S3A-G). Analyzing the relationship with patient prognosis, significant survival differences were noted between high-risk and low-risk groups in most subgroups (Figure S3H-S). This implies that high-risk classification is linked to shorter survival, especially in advanced stages compared to early stages.To explore the link between risk-scores and the immune microenvironment, we used the ESTIMATE algorithm. The high-risk group showed significantly lower immune scores and ESTIMATE scores than the low-risk group, while there was no notable difference in stromal scores between the two groups (Figure 4M). In the context of 13 immune-related pathways, the high-risk group demonstrated lower pathway activity (Figure S4A). Further analysis of immune cell infiltration between the two subtypes revealed that the high-risk group had significantly lower levels of T cells CD8, T cells CD4 memory activated, M1, and NK cells resting (Figure 4N). Correlation analysis between risk-scores and immune cells indicated a negative correlation with anti-tumor immune cells. Conversely, it showed a positive correlation with immune inhibitory cells (Figure S5A-P, Figure S4B). In addition, a significant correlation was observed between the 4 individual genes (CALR, IFN1B, IFNG and IL1R1) comprising the risk-score and immune cells, immune checkpoints or immune regulatory pathways (Figure S6A-P, Figure S7A-C). These analyses suggest that higher risk-scores are linked to weaker anti-tumor immune function and more tumor-promoting elements, partly explaining the shorter survival in the high-risk group.

3.5 The association between the risk-scores and the sensitivity to chemotherapy and immunotherapy

One of the primary factors inducing ICD is chemotherapy drugs, and targeted therapy plays a crucial role in BLCA treatment. Therefore, we further explored the correlation between the risk-scores and the sensitivity to chemotherapy drugs and targeted therapies. The high-risk group exhibited higher sensitivity (lower IC50 values) to sorafenib, epothilone B, docetaxel, elesclomol, MG132, lapatinib, and others. Conversely, the ICD-high risk group showed increased resistance (higher IC50 values) to lenalidomide, sunitinib, gefitinib, methotrexate, and camptothecin (Figure 5A-P). Furthermore, we investigated the prognostic characteristics of both groups in response to immunotherapy. The outcomes indicate that without anti-CTLA4 and anti-PD1 treatment, there was no significant difference in survival time between the groups. However, with at least one form of immunotherapy, the low-risk group had a significantly higher response rate than the high-risk group (Figure 5Q-T). These findings highlight the clinical potential of the risk-score in predicting sensitivity to various therapies and guiding treatment decisions for better outcomes.

Figure 5.

3.6 Validate the risk-scoring model through BLCA tissue arrays and analyze the gene expression with single-cell sequencing

To validate the reliability of the aforementioned results, we conducted an analysis of the expression levels of CALR, IFN1B, IFNG and IL1R1 and their associations with prognosis and immune function using BLCA tissue arrays. Immunofluorescence staining revealed that CALR expression in tumor tissues was higher than in normal tissues, while the other three genes, IL1R1, IFNB1, and IFNG, showed no significant differences (Figure 6A, B, Figure S8A). This observation is consistent with the results from the TCGA-BLCA dataset exported using the GEPIA 2 tool (Figure S8B-E). Patients were divided into high and low expression groups based on the median expression levels (IRS scores) of CALR, IFN1B, IFNG, and IL1R1. It was found that patients in the CALR high-expression group and IL1R1 high-expression group had shorter survival times, while those in the IFNB1 high-expression group and IFNB high-expression group had longer survival times (Figure 6C), consistent with the results of TCGA-BLCA dataset (Figure 4A-D). Patients’ risk-scores were similarly calculated using the risk-scoring model: Risk-score = CALR-IRS*(0.5378) + IFN1B-IRS *(-0.6349) + IFNG-IRS *(-0.2028) + IL1R1-IRS *(0.0919). The high-risk group (50% vs. 50%) had significantly shorter survival times than the low-risk group (Figure 6D), corroborating the previous results (Figure 3D, E). More importantly, we found a significant positive correlation between the risk-score and CD39, a molecule we previously confirmed to be associated with immune suppression¹⁴. Additionally, the risk score showed a significant positive correlation with immune exhaustion (CD8+LAG3+), while having a negative correlation with CD8+ T cell infiltration (the IRS scores are based on our previous research¹⁴). However, there was no correlation with tumor volume and T-stage (Figure 6E-I). This suggests a strong link between higher risk scores and immune suppression.

Figure 6.

Finally, using our previous single-cell sequencing data¹⁴ and the public dataset GSE135337, we analyzed the cellular expression patterns of CALR, IFN1B, IFNG, and IL1R1. CALR is expressed in nearly all cells in the TME, IL1R1 is mainly expressed in tumor-associated fibroblasts and endothelial cells, and the expression levels of IFNG and IFNB1 are minimal, primarily by T lymphocytes. This suggests a crucial role of cancer-associated fibroblasts in remodeling the immune microenvironment of BLCA (Figure 6J, K, Figure S8F, G).

4 Disscussion

ICD, a regulated cell death type, can be triggered by internal or external antigens, activating an adaptive immune response⁸. Some chemotherapies, such as cisplatin and gemcitabine, can induce ICD, thereby activating the immune microenvironment to achieve antitumor effects^9,10. ICIs are a key focus in cancer research, blocking tumors by reactivating immune cells through targeting checkpoints like PD1/PD-L1 and CTLA4¹⁵. Combining chemotherapy with ICI treatment is expected to enhance immune function and inhibit tumors, given ICD’s role in activating the immune microenvironment, and some studies targeting key genes in the ICD process or combining chemotherapy with immunotherapy have also achieved encouraging results^6,16,17. However, the results of clinical trials indicate that the response rate of BLCA to ICIs is significantly lower than expected. The tumor immunosuppressive microenvironment is undoubtedly an important influencing factor. Considering the complexity of the TME and the impact of ICD, we investigate and confirm the clinical utility of ICD-related genes using BLCA public databases, BLCA tissue microarrays, and single-cell sequencing data. This offers theoretical guidance for applying ICD in clinical settings.

Through unsupervised clustering analysis, we categorized TCGA-BLCA patients into two groups: ICD-high and ICD-low. Pathway enrichment analysis revealed that the ICD-high group predominantly clustered in pathways related to the immune system, consistent with previous literature reports^18,19. In the analysis of immune infiltration, we found that the ICD-high group showed higher immune cell infiltration, coupled with lower tumor purity, indicating a lower invasive ability of the tumor and an activated immune state²⁰. We checked if ICD boosts ICI effectiveness by studying the connection between ICD-related gene expression and immunotherapy response rates. Turns out, the ICD-high group showed a much better response compared to the ICD-low group. This suggests that ICI treatment could work better in people with high ICD. Thus, we created a risk-scoring model using four ICD-related genes (CALR, IL1R1, IFNB1, and IFNG) identified through Cox regression and LASSO regression analyses. This model helps classify BLCA and predict immune infiltration, prognosis, and response to immunotherapy. Of them, IFNB1 and IFNG, both belonging to the interferon family, play crucial roles in anti-tumor immunity^21,22. IL1R1 is a receptor for IL1α and IL1β, a member of the interleukin family, and plays an immunosuppressive role in immune regulation²³. In ICD, the CALR calcium-binding protein moves to the cell membrane to aid antigen presentation, a widely accepted process²⁴. Yet, the general increase in CALR levels does not necessarily have a positive regulatory effect. Studies indicate that elevated CALR expression in breast cancer can enhance tumor cell metastasis and resistance to chemotherapy²⁵. Similarly, our study found that high CALR expression in BLCA has an inhibitory effect. Survival analysis showed that high CALR and IL1R1 expression was detrimental to patient survival, while high IFNB1 and IFNG expression was beneficial.

To verify our risk-scoring model, we used univariate and multivariate Cox regression analyses. The results indicate that the ICD scoring system can independently predict the prognosis of BLCA patients. Subgroup analysis based on clinical characteristics showed no difference between high-risk and low-risk groups in the early stages of the disease. However, in later stages, the high-risk group had shorter survival times. This may be due to changes in the TME at different disease stages, with a decrease in immune-killing molecules and an increase in inhibitory molecules as the tumor advances²⁶. Additional analysis of the correlation between the risk-scores and immune cell infiltration confirmed this: as the risk-score increased, the infiltration of cells with anti-tumor effects decreased, while the infiltration of M2-type macrophages and cancer-associated fibroblasts increased. These findings confirm the efficacy of our risk-scoring model in assessing patients’ immune infiltration, prognosis, and tumor progression.

In treating BLCA, chemotherapy, targeted therapy, and immunotherapy are effective but have limitations in sensitivity and specificity²⁷. Our risk signature predicts drug sensitivity. High-risk group responds better to chemotherapy, and low-risk group shows higher response in immunotherapy (PD1 or CTLA4). This matches the higher immune infiltration in the low-risk patient group, emphasizing that immune infiltration level is crucial for successful immunotherapy. However, we find here that the quantity of immune cells in the tumor does not directly decide the patient’s prognosis. This is because many infiltrating immune cells are functionally exhausted, and only immunotherapy can activate the immune response within the tumor, improving prognosis. This is confirmed in BLCA tissue microarray stains. Specifically, the risk-score shows a negative correlation with CD8+ T cell infiltration but a positive correlation with T cell exhaustion (CD8+LAG3+ or CD39). This shows our risk scoring model can categorize patients based on immune infiltration, specifically identifying a high-risk group (lower immune infiltration) and a low-risk group (higher immune infiltration). Personalized treatment is crucial for patients in this context. Additionally, in BLCA single-cell sequencing, we found IL1R1 is abundantly expressed in cancer-associated fibroblasts, linked to tumor promotion and immunosuppression²⁸. This suggests a crucial role for these fibroblasts in shaping the immune microenvironment of BLCA, guiding our next research.

5 Conclusions

ICD, a type of cell death triggered in an immunogenic environment, influences immune response. Our study links ICD-related genes to BLCA prognosis and tumor immune infiltration. We created a validated risk-scoring model with four key genes (CALR, IL1R1, IFNB1, and IFNG) using TCGA-BLCA, GEO datasets, and tissue microarray staining. This model effectively assesses immune infiltration, prognosis, immunotherapy response, and drug sensitivity in BLCA, guiding personalized treatment and future research.

Acknowledgements

None.

CRediT authorship contribution statement

Lilong Liu: Funding acquisition;Data curation; Formal analysis;Writing - original draft; and Writing - review & editing.Zhenghao Liu, Lei Fan and Zhipeng Yao: Data curation; Formal analysis; Writing - original draft; and Writing - review & editing. Junyi Hu, Yaxin Hou,Yang Li, Yingchun Kuang, and Yuhong Ding: Data curation and Writing - review & editing. Yi Hao played crucial roles in Funding acquisition,conceptualization and Writing - review & editing. Zheng Liu, and Ke Chen played crucial roles in conceptualization and Writing - review & editing. All authors approved the final manuscript for publication.

Ethics approval and consent to participate

Patients provided informed consent for this work. All experimental procedures were approved by the Institutional Review Board of Tongji Medical College, Huazhong University of Science and Technology (IRB: S116).

Consent for publication

Not applicable.

Declaration of interest statement

None.

Funding

This research was supported by the National Natural Science Foundation of China (No.82303623), the National Key Research and Development Project of China (No.2019YFA0905600), the China Postdoctoral Science Foundation (No.2023M731199), and the Postdoctoral Innovation Research Post of Hubei Provincial Department of Human Resources and Social Security (No.331048). The authors acknowledge that the views expressed in this work are personal and may not necessarily represent the perspectives of official institutions.

Data Availability Statement

We obtained TCGA data from https://portal.gdc.cancer.gov/ and GEO data from https://www.ncbi.nlm.nih.gov/geo/. Access our single-cell sequencing data¹⁴ at https://ngdc.cncb.ac.cn/gsa-human/browse/HRA000212 or https://www.ncbi.nlm.nih.gov/bioproject/PRJNA662018/.