Skip to main content
Download PDF

Apoptosis related genes mediated molecular subtypes depict the hallmarks of the tumor microenvironment and guide immunotherapy in bladder cancer

BMC Medical Genomics volume 16, Article number: 88 (2023) Cite this article

Abstract

Apoptosis has been discovered as a mechanism of cell death. The purpose of this study is to identify the diagnostic signature factors related to bladder cancer (BLCA) through apoptosis related genes (ARGs). Clinicopathological parameters and transcriptomics data of 1,440 BLCA patients were obtained from 7 datasets (GSE13507, GSE31684, GSE32548, GSE32894, GSE48075, TCGA-BLCA, and IMvigor210). We first identified prognosis-related ARGs in BLCA and used them to construct two ARGs molecular subtypes by using consensus clustering algorithm. By using principal component analysis algorithms, a ARGscore was constructed to quantify the index of individualized patient. High ARGscore correlated with progressive malignancy and poor outcomes in BLCA patients. High ARGscore was associated with higher immune cell, higher estimate scores, higher stromal scores, higher immune scores, higher immune checkpoint, and lower tumor purity, which was consistent with the “immunity tidal model theory”. Preclinically, BLCA immunotherapy cohorts confirmed patients with low ARGscore demonstrated significant therapeutic advantages and clinical benefits. These findings contribute to our understanding of ARGs and immunotherapy in BLCA. The ARGscore is a potentially useful tool to predict the prognosis and immunotherapy in BLCA.

Peer Review reports

Introduction

Bladder cancer (BLCA) is a major contributor to cancer-related deaths [1]. BLCA can present as non-muscle-invasive BLCA, muscle-invasive BLCA, and metastatic disease. Currently, the efficacy of BLCA treatment strategies is not satisfactory due to high heterogeneity and drug resistance [2, 3]. Immune checkpoint inhibitors (ICIs) have shown promising therapeutic results in BLCA [4, 5]. However, due to differences in individual immunogenicity and tumor microenvironments (TME), the response rate to this therapy in patients without selective treatment is not satisfactory. Meanwhile, the biomarkers of ICIs has certain limitations [6, 7]. Thus, it is important to develop immunotherapy biomarkers from a genomic perspective to meet the needs of individualized therapy for patients with different genetic backgrounds.

Recent studies have shown that apoptosis regulates tumorigenesis and progression extensively [8]. Tumor development and immune response are all influenced by the apoptosis signaling pathway [9]. In response to stimulation by appropriate external or internal signals, cells may alter the expression of genes encoding proteins associated with the apoptotic process. Alteration of genes associated with apoptosis in bladder cancer leads to a shift from a default apoptotic state to drug-resistant cells with higher survival characteristics [10]. The TME contains tumor cells, vasculature, extracellular matrix (ECM), stromal, and immune cells [11,12,13]. In vitro and in vivo experiments have also shown that apoptotic tumor cells induce antitumor immunogenicity by cross-initiating and proliferating CD8 T cells [14, 15].

Therefore, this study focuses on the functional role of apoptosis related genes (ARGs) in BLCA. The purpose of this study is to explore the expression of ARGs in BLCA by bulk RNA sequencing, and to identify their expression. To further investigate the effect of ARGs expression in BLCA on the development of BLCA, and to explore its influence on the biological behavior of BLCA, and to provide a basis for the etiology, pathogenesis and prevention of primary BLCA. Meanwhile, we adopted the consensus clustering algorithm to construct the consensus matrix and cluster the samples’ BLCA related expression profile data. Then, a ARGscore was constructed to quantify the apoptosis index of individualized patient. Our study might provide clues for targeted therapeutics in BLCA.

Materials and methods

Data collection and pre-processing

Clinicopathological parameters and transcriptomics data were obtained from nine independent BLCA cohorts (GSE13507 [16], GSE31684 [17], GSE32548 [18], GSE32894 [19], GSE48075 [20], TCGA-BLCA, and IMvigor210 [21]). The FPKM of BLCA and normal tissues were downloaded from TCGA database website, and the raw data were extracted into matrix files by Perl software. The correspondence between human gene names and gene ids was downloaded from ensemble website, and the gene ids in the raw data were converted into gene names by Perl software. Then, the FPKM values are converted to TPM values. The TCGA data and GEO data are integrated into a meta cohort using the "sva" package. The data information is summarized in Table 1 and Additional file 1: Table S1.

Table 1 Basic information of BLCA cohorts included in this study
Full size table

Expression and prognosis of ARGs in BLCA

We obtain 91 ARGs from the previous studies (Additional file 2: Table S2) [22]. Then, we compared ARGs expression in TCGA-BLCA (|log FC|≥ 1.0 and adj. P < 0.05) by using “limma” package [23]. Univariate Cox regression analyses were performed by “survival” packages to acquire prognostic-related ARGs [24].

Consensus clustering analysis

In order to better classify the ARGs into different clusters, we adopted the ConsensusClusterPlus tool to construct the consensus matrix and cluster the ARG related expression profile data of samples. We adopted the "KM" algorithm and "1-Pearson correlation", and applied the metric distance and 500 bootstraps, with each bootstrap process including 80% of patients in the training set. We identified the cluster number as 2 to 10. The number of clusters (K) is chosen by the maximum average profile value [25].

Identification of DEGs and prognosis genes between ARGs subtypes

Different expression of ARGs were analyzed in 2 ARGs subtypes using the "limma" package, and the screening criteria were | log2Fold Change (FC)|> 1, FDR < 0.05. Univariate Cox regression analyses were performed by “survival” packages to acquire prognostic-related DEGs.

GO and KEGG enrichment analysis

To explore the biological functions mainly performed by DEGs and prognosis genes between ARGs subtypes, then the up-regulated and down-regulated genes were analyzed for GO and KEGG enrichment in metascape [26,27,28]. The screening criteria were P < 0.05 and FDR < 0.05[29].

Evaluation of immune infiltrating cells in BLCA

The ESTIMATE algorithm can calculate the immune and stromal fractions of the tumor microenvironment of tumor samples, and a larger fraction indicates a larger proportion of immune cells or stromal cells in the tumor microenvironment. ssGSEA and TIMER were utilized to reveal the immune cell infiltration patterns, determine the proportions of circulating immune cells and analyze the correlation between these immune cells [30].

Tracking tumor immunophenotype analysis

Tracking Tumor Immunophenotype (TIP) analysis was done according to website protocol. Transcriptomic expression profile was used as the input for analysis [31].

Construction of ARGscore

Based on the DEGs and prognosis genes between ARGs subtypes, we established a ARGscore by using PCA method [32,33,34,35]. The final ARGscore was calculated as follows:

$$\mathrm{ARGscore }=\Sigma \left(P\mathrm{C}{1}_{i}+P\mathrm{C}{2}_{i}\right),$$

where i is the expression value of each gene in the signature. A high ARGscore group and a low ARGscore group were divided base on the median ARGscore.

GSEA analysis, clinical relevance, and immune correlation analysis of ARGscore

The differentially expressed genes in the ARGscore high and low groups in liver cancer samples were analyzed by GSEA to explore the impact of the difference in ARGscore levels on cell signaling pathways. Enrichment analysis was performed according to the default reference settings, and q < 0.05 was selected as the set of significantly enriched genes [36].

Immunohistochemistry (IHC) staining

IHC was performed on formalin fixed, paraffin-embedded tissue sections using a two-step protocol. Then using traditional method to evaluate under the optical microscope (N = 10). The paired BLCA tissues and non-tumor tissues were collected from the patients at Wuming Hospital, Guangxi Medical University. This study was evaluated and reviewed by the Ethics Committee of Wuming Hospital, Guangxi Medical University (WM-2022-15). Written consent was obtained from all patients prior to the study.

Statistical methods

Data analysis and graphing were carried out using R 4.2.2 software. The Wilcoxon assay was used to assess the difference in ARGs between normal and tumor tissues. The association between ARGs regulators expression and prognosis was evaluated using univariate Cox regression analysis. Wilcoxon analysis was performed to assess the correlation between ARGscore and immune cell infiltration levels and immune checkpoint inhibitors expression. The predictive efficiency of the ARGscore model was tested by the ROC curve. The difference was considered statistically significant at P < 0. 05 (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001).

Results

Differential expression and prognosis of ARGs

The workflow of our analysis was shown in Fig. 1. Firstly, GSEA showed that apoptosis pathway was significant enriched in the BLCA (Fig. 2A). Then, we compared the expression of 91 ARGs, of which the expression level of 10 ARGs were significantly higher in BLCA than in the normal tissues (Fig. 2B, C, Additional file 3: Table S3). In total, 219 nonsmokers and 59 smokers had adequate knowledge of age, sex, grade, clinical staging, survival status, and follow-up time. The mean age of all participants at diagnosis was 68.97 ± 10.57 years, with 71.22% of male participants, almost 2.5 times more than females. Next, univariate Cox regression analysis was performed based on TCGA data, 24 ARGs coorelated with outcome in BLCA (Fig. 2D, Additional file 4: Table S4). Based on these results, 4 ARGs (CAPN14, EGFR, YWHAQ, MAP2K1) were considered as risk genes. The 4 ARGs had positive correlation with each other (Fig. 2E). Moreover, CAPN14, EGFR, YWHAQ, MAP2K1 expression were higher in stage III & stage IV than stage I & stage II, higher in radiation therapy group than in no radiation therapy group, higher in PD&SD group than in PR&CR, higher in dead group than in alive, higher in non-papillary group than in papillary, higher in high Grade group than in low Grade, higher in smoker group than in no smoker (Fig. 2F–L). Meanwhile, these four genes were significantly associated with immune infiltration (Additional file 8: Fig. S1A–D). The above results suggest that these 4 ARGs may be associated with the metastasis of BLCA.

Fig. 1
figure 1

Flow chart of the study

Full size image
Fig. 2
figure 2

Identification of differential expression and prognostic of ARGs in TCGA-BLCA. A GSEA showed that the apoptosis pathway are differentially enriched in BLCA patients. B Venn diagram to identify differentially expressed genes between tumor and normal tissue that were correlated with overall survival. C Volcano plot of differentially expressed ARGs between tumor and normal tissue. Dot stands for gene. Red dots represent up-regulated genes. D Forest plots showing the results of the univariate Cox regression between 8 ARGs and overall survival in BLCA. E Correlation heat map of 4 ARGs. F The illustration shows the expression distribution of 4 ARGs between stage 3&4 (red) and stage 1&2 (blue). G The illustration shows the expression distribution of 4 ARGs between PR&CR group (red) and PD&SD group (blue). H The illustration shows the expression distribution of 4 ARGs between radiation therapy group (red) and no radiation therapy group (blue). I The illustration shows the expression distribution of 4 ARGs between dead (red) and alive (blue). J The illustration shows the expression distribution of 4 ARGs between papillary group (red) and non-papillary group (blue). K The illustration shows the expression distribution of 4 ARGs between high grade group (red) and low grade (blue). L The illustration shows the expression distribution of 4 ARGs between smoker group (red) and no smoker group (blue). (*P < 0.05; **P < 0.01; ***P < 0.001)

Full size image

Identification of ARGs subtypes in BLCA

Then, we adopted the ConsensusClusterPlus tool to construct the consensus matrix and cluster the 4 ARGs related expression profile data of samples. The unsupervised method NMF analysis results in 2 subtypes ARGs.cluster.A (n = 548) and ARGs.cluster.B (n = 544) (Additional file 8: Fig. S2A–F). The overall survival (OS) time in ARGs.cluster.A group was significantly longer than that in ARGs.cluster.B group (HR = 1.48, P < 0.001, Fig. 3A). The PCA demonstrated there is significant DEGs between the two clusters, and 2859 prognostic ARGs-related DEGs were identified (Fig. 3B, Additional file 5: Table S5). Then, we applied GO and KEGG, which showed that ARGs.cluster.B was significantly enriched in immune-related pathways (T cell activation, T cell migration, T cell apoptotic process, MHC protein binding, and cytokine-cytokine receptor interaction) and cancer-related pathways (DNA replication, cell growth, apoptosis, and PI3K-Akt signaling pathway) (Fig. 3C, D), suggesting that ARGs.cluster.B may play an important role in tumor development and immune regulation.

Fig. 3
figure 3

The molecular subtypes categorization of BLCA base on 4 ARGs. A Kaplan–Meier curve showed a significant difference between the two ARGs clusters. B UMAP analysis for the transcriptome profiles of ARGs.cluster.A and ARGs.cluster.B, showing a remarkable difference on transcriptome between different group. C GO enrichment analysis, D KEGG enrichment analysis for the different expression that were correlated with OS genes. E The correlation of tumor microenvironment condition and ARGs.cluster. F The abundance of each TME infiltrating cell in two ARGs.clusters. The upper and lower ends of the boxes represented the interquartile range of values. The lines in the boxes represented median value, and black dots showed outliers. G–J The ESTIMATE score, stromal score, immune score, and tumor immunity levels in the ARGs.cluster.A and ARGs.cluster.B groups by using ESTIMATE algorithm. K–M The PD-1, PD-L1 and CATA4 levels in the ARGs.cluster.A and ARGs.cluster.B groups. (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001)

Full size image

In order to clarify the differences in TME between the 2 established clusters, we investigated the gene expression and infiltration level of immune cells by ssGSEA. The results suggested that the infiltration level of most immune cells was significantly lower in ARGs.cluster.A group (Fig. 3E,F). The ESTIMATE results showed that the stromal score, immune score and ESTIMATE score in ARGs.cluster.B group were much higher than those in ARGs.cluster.A group (Fig. 3G–J). Immune checkpoint blockade-mediated immunotherapy has shown promising effects against tumors, and has proved as effective in a significant proportion of refractory patients undergoing standard therapy. In this analysis, we evaluated some representative targets, and found that immune checkpoint were highly expressed in ARGs.cluster.B group (Fig. 3K–M, Additional file 8: Fig. S3). The above results suggest that ARGs.cluster.B was consistent with the “immunity tidal model theory”, thus leading to a worse prognosis for ARGs.cluster.B.

Construction of ARGscore

A ARGscore was constructed based on 2859 prognostic ARGs-related DEGs. Alluvial diagram confirmed the relationship between the ARGs.cluster and ARGscore (Fig. 4A). ARGs.cluster.B group had a higher proportion of high ARGscore (Fig. 4B). Next, we analyzed the distribution of ARGscore in the ARGs.cluster, and ARGscore was lower in the ARGs.cluster.A group (Fig. 4C). The results also showed that patients with high ARGscore had a worse overall survival (HR = 1.80, P < 0.001, Fig. 4D). To test the prediction performance of ARGscore, we plotted ROC curves with AUC values of 0.841, suggesting ARGscore had a good predictive performance (Fig. 4E). In addition, the PD-1 gene expression was highly expressed in high ARGscore group (Fig. 4F). The above results suggest that high ARGscore may account for the worse prognosis of patients in the ARGs.cluster.B group.

Fig. 4
figure 4

Construction of ARGscore. A Alluvial diagram showing the changes of ARGs.cluster and ARGscore. B The proportion of patients with high and low ARGscore in ARGs.cluster.A and ARGs.cluster.B groups. C Differences in ARGscore among two ARGs.cluster. D Kaplan–Meier curves for high and low ARGscore groups. E The predictive value of ARGscore. F The PD-1 levels in the high and low ARGscore groups. G Differences in various steps of Cancer-Immunity Cycle. (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001)

Full size image

To further explore the relationship between ARGscore and tumor immunology, we performed TIP analysis. To further explore the relationship between ARGscore and tumor immunology, we performed TIP analysis. The results suggest that steps 1, 2, 3 and 4 were more abundant in high ARGscore group, whilesteps 6 and 7 were more abundant in low ARGscore group. This result indicated that although high ARGscore group was helpful for initiation and processing phases of immune response, effective antitumor immunity was still suppressed (Fig. 4G).

The results also showed that cell cycle-related pathways (cell cycle checkpoint signaling, cell cycle phase transition, chromosome segregation, cell cycle) and tumor-related pathways (bladder cancer, DNA replication, P53 signaling pathway) were enriched in the high ARGscore group (Fig. 5A–L). These results suggested that high ARGscore may have a significant impact in tumor development.

Fig. 5
figure 5

GSEA on the meta cohort to explore mechanisms underlying ARGscore. A–F GSEA GO identified high and low ARGscore groups related signaling pathways in BLCA. G–L GSEA KEGG identified high and low ARGscore related signaling pathways in BLCA

Full size image

ARGscore is a robust prognosis factor in BLCA

To further validate the robustness of ARGscore, the prognostic implication of ARGscore was examined in multiple independent datasets. In both patient sets, patients in the high ARGscore group had a worse OS than in the low ARGscore group (GSE13507: HR = 2.35, 95% CI = 1.45–3.81, P = 0.001, GSE31684: HR = 2.02, 95% CI = 1.04–3.93, P = 0.038, GSE32548: HR = 2.30, 95% CI = 1.00–5.25, P = 0.049, GSE32894: HR = 2.17, 95% CI = 1.05–4.47, P = 0.037, GSE48075: HR = 1.93, 95% CI = 1.02–3.66, P = 0.044, TCGA-BLCA: HR = 1.37, 95% CI = 1.01–1.85, P = 0.043) (Fig. 6A, C, E, G, I, K). In addition, the ROC curve analysis results demonstrated that the AUC of ARGscore in the GSE13507, GSE31684, GSE32548, GSE32894, GSE48075, TCGA-BLCA set were 0.799, 0.860, 0.873, 0.847, 0.802, and 0.807, respectively (Fig. 6B, D, F, H, J, L).

Fig. 6
figure 6

Internal and external validation of ARGscore. A Kaplan–Meier curves for high and low ARGscore groups in GSE13507. B The predictive value of ARGscore in GSE13507. C Kaplan–Meier curves for high and low ARGscore groups in GSE31684. D The predictive value of ARGscore in GSE31684. E Kaplan–Meier curves for high and low ARGscore groups in GSE32548. F The predictive value of ARGscore in GSE32548. G Kaplan–Meier curves for high and low ARGscore groups in GSE32894. H The predictive value of ARGscore in GSE32894. I Kaplan–Meier curves for high and low ARGscore patient groups in GSE48075. J The predictive value of ARGscore in GSE48075. K Kaplan–Meier curves for high and low ARGscore patient groups in TCGA-BLCA. L The predictive value of ARGscore in TCGA-BLCA

Full size image

High ARGscore suppress anti-PD1/L1 immunotherapy

To further validate the robust of ARGscore in immunotherapy, 2 cluster in IMvigor210 cohort was identified by consensus clustering algorithm (Additional file 8: Fig. S4A–F). Then, a ARGscore was constructed based on 13 prognostic ARGs-related DEGs (Additional file 6: Table S6). Patients with high ARGscore had a worse OS (HR = 2.14, 95% CI = 1.47–3.12, P < 0.001, Fig. 7A). The AUC values 1, 3, and 5-survival rates in the cohort were 0.807, 0.767, 0.770, suggesting ARGscore had a good predictive performance (Fig. 7B). The ARGscore was lower in CR/PR group than in PD/SD groups (Fig. 7C). Moreover, the number of CR/PR patients was markedly higher in low ARGscore group than that in high ARGscore group (38% vs. 19%) (Fig. 7D, F, G), suggesting low ARGscore was sensitive to immunotherapy than high ARGscore group. The low ARGscore group identified twice as many patients as the high ARGscore group who responded to immunotherapy. In addition, ARGscore was lower in CR and PR group than in PD and SD groups (Fig. 7E). Multivariate Cox indicated also indicated that ARGscore could be acted as independent risk factor of BLCA patients in IMvigor210 cohort (HR = 1.84, CI = 1.37–2.33, P < 0.001) (Fig. 7H, Additional file 7: Table S7). The results showed that ARGscore was lower in desert group, TCGA I and II subtype than that other groups (Fig. 7I, J).

Fig. 7
figure 7

ARGscore in the role of anti-PD-1/L1 immunotherapy. A Survival analyses for low and high ARGscore groups in IMvigor210 cohort. B The predictive value of ARGscore in IMvigor210 cohort. C Differences in ARGscore among distinct anti-PD-1 clinical response groups. D The proportion of patients with response to PD-L1 blockade immunotherapy in low or high ARGscore groups. E Distribution of ARGscore in distinct anti-PD-L1 clinical response groups. F, G The number of CR, PR, PD, and SD patients in high and low ARGscore groups. H Multivariate Cox regression analysis for ARGscore in IMvigor210 cohort shown by the forest plot. I, J Differences in ARGscore between immune subtypes and TCGA subtypes

Full size image

Validation 4 ARGs in BLCA tissues

We used IHC and qPCR to determine the protein level and mRNA expression of 4 ARGs with a tissue array of 10 BLCA patient samples. The IHC and qPCR result showed the expression of EGFR, YWHAQ, and MAP2K1 were significantly higher in BLCA than in their corresponding adjacent nontumor tissues, CAPN14 was significantly lower in BLCA than in their corresponding adjacent nontumor tissues (Fig. 8A, B).

Fig. 8
figure 8

Representative immunohistochemistry images and qPCR of 4 ARGs expression in normal tissues and BLCA tissue. (A) Immunohistochemistry of 4 ARGs expression in normal tissues and BLCA tissue. (B) qPCR of 4 ARGs expression in normal tissues and BLCA tissue

Full size image

Discussion

Currently, chemotherapy remains one of the main treatment modalities for the vast majority of patients with advanced cancer, but tumor drug resistance has become a bottleneck limiting the efficacy of chemotherapy. Cell death can overcome the resistance of tumor cells to conventional chemotherapeutic drugs and promote the clearance of defective cells. Therefore, the induction of apoptosis in tumor cells may be a new approach for tumor treatment. Chen et al. revealed that guggulsterone inhibits BLCA growth and metastasis through caspase-dependent apoptosis [37]. Yin et al. showed that epigallocatechin-3-gallate can inhibited bladder cancer cells via inducing autophagy-related apoptosis [38]. We constructed a ARGscore prognostic model based on 1092 BLCA patients. Subsequently, the role of TME in BLCA was explored. This study revealed potential biomarkers and therapeutic targets in apoptosis-related signaling pathways.

In this study, we first downloaded transcriptomic information and clinical data of BLCA patients from the TCGA database. By comparing the expression of 91 ARGs in the tumor group and normal group, we showed that CAPN14, EGFR, YWHAQ, MAP2K1 were prognostic predictors of BLCA by univariate Cox regression analysis, and high expression of EGFR, YWHAQ, MAP2K1 were associated with poor prognosis of BLCA, low expression of CAPN14 was associated with poor prognosis of BLCA. Then, we divided BLCA into two molecular subtypes based on the expression of 4 ARGs. Cluster.B group had a worse prognosis compared to cluster.A, and patients in this group had a higher immune infiltration and higher ICI expression. This type is known as immunity tidal model type. Next, ARGscore model was constructed based on prognostic ARGs-related DEGs. ROC curve analysis confirmed that the model showed good predictive power for patients with BLCA.

Tregs and macrophage M2 have been shown to exert pro-tumor functions, while activated NK cells have the potential to kill tumor cells in different ways and are important mediators of cancer immune detection. In addition, Th1 and Th2 balance is associated with anti-tumor immunity, and in breast cancer patients, Th1/Th2 imbalance and elevated cytokine release from Th2 have been observed, and Caihu saponin can inhibit breast cancer growth by shifting Th1/Th2 balance to Th1 [39,40,41]. In the present study, immunosuppressive cells such as Tregs, Th2 and macrophages were more abundant in the high ARGscore group with poorer prognosis. The results suggest that the model is to some extent related to the immune landscape of the BLCA microenvironment, and the ARGscore model we constructed may play a key role in immune cell infiltration in BLCA.

M2-like macrophages, a major part of TAMs, have a strong influence on the formation of the immunosuppressive microenvironment [42, 43]. TAMs secrete a variety of growth factors, pro-inflammatory factors and enzymes. IL-6 and TNF-α can activate signaling pathways such as signal transduction molecules Smad and Wnt, which in turn activate transcription factors Snail family and Twist, prompting tumor cells to shift from epithelial to mesenchymal cell type, leading to a decrease in intercellular adhesion, shedding of tumor cells into the interrogative stroma and gaining the ability to easily metastasize to distant sites [44]. In vitro experiments showed that culturing BLCA cells with medium culturing M2-type macrophages resulted in activation of the TLR4/STAT3 signaling pathway in BLCA cells, increased expression levels of EMT markers, and enhanced metastatic ability of tumor cells [45].

In recent years, ICIs have significantly transformed oncology treatment with the rise of immune checkpoint blockade therapies. Advances in checkpoint blockade therapies have deepened our understanding of the interactions between the immune system, cancer cells, and the tumor microenvironment. Tumor mutational burden and immune checkpoint-associated gene expression do not accurately predict response to ICIs therapy. Therefore, identifying biomarkers that can accurately predict the clinical efficacy of ICIs therapies is important to further advance precision immunotherapy. Among the many immunotherapeutic strategies, immune checkpoint inhibitors have shown great benefits in the treatment of a range of cancer types, such as the use of CTLA-4 and PD-1 or PD-L1 to enhance anti-tumor immunity. At present, the biomarkers of ICIs therapy that show the most promise are PD-L1 expression and tumor mutation burden (TMB), and microsatellite instability (MSI) [46, 47]. However, each of these biomarkers has certain limitations. The expression of PD-L1 is regulated by many mechanisms, and as a result, it has low reliability since there are often great differences in its expression at different times and conditions. Meanwhile, the positivity rate of PD-L1 in BLCA was less than 10% [48]. Study also showed a significant correlation between positive expression of PD-L1 and immune efficacy of immunotherapy. TMB is greatly influenced by the equipment of different laboratories, and the level of TMB in BLCA was not significant [49]. In contrast, MSI or MMR occurs in only 2–3% of BLCA patients [50]. Through our research, we have confirmed the influence of ARGscore on the prognosis and immunotherapy of BLCA patients who received ICIs treatment. Taking the expression of 4 ARGs as the significant criteria, this paper expounds the influence of ARGscore on the immune microenvironment and immunogenicity, so as to create accurate guidelines for the use of ICIs treatment in BLCA patients. The low ARGscore group identified twice as many patients as the high ARGscore group who responded to immunotherapy.

The present study also has some limitations; our study is based on publicly mined databases, no clinical experimental studies have been conducted, and their role and function in apoptosis processes are yet to be validated using more advanced methods and techniques. Second, this study was retrospective, so the risk score model needs to be further validated in prospective studies and multicenter clinical trials.

Conclusion

In summary, we identified that SCAPN14, EGFR, YWHAQ, MAP2K1 could be used as prognostic factors for BLCA by comprehensive analysis of ARGs. Meanwhile, 2 molecular subtypes were constructed based on 4 ARGs, and ARGscore was constructed based on the ARGs-related DEGs between subtypes. ARGscore can be used as BLCA independent predictor with good predictive ability. Also, high ARGscore group was sensitive to immunotherapy. In conclusion our study provides a promising avenue for individualized survival prediction and clinical outcome prediction of ICIs antitumor immunotherapy in patients with BLCA.

Availability of data and materials

All data used in the study can be downloaded from the TCGA data repository (https://portal.gdc.cancer.gov/) and the GEO database (https://www.ncbi.nlm.nih.gov/geo/).

References

  1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2021. CA Cancer J Clin. 2021;71(1):7–33.

    Article  PubMed  Google Scholar 

  2. Lenis AT, Lec PM, Chamie K, Mshs MD. Bladder cancer: a review. JAMA. 2020;324(19):1980–91.

    Article  CAS  PubMed  Google Scholar 

  3. Alifrangis C, McGovern U, Freeman A, Powles T, Linch M. Molecular and histopathology directed therapy for advanced bladder cancer. Nat Rev Urol. 2019;16(8):465–83.

    Article  CAS  PubMed  Google Scholar 

  4. Patel VG, Oh WK, Galsky MD. Treatment of muscle-invasive and advanced bladder cancer in 2020. CA Cancer J Clin. 2020;70(5):404–23.

    Article  PubMed  Google Scholar 

  5. Witjes JA, Bruins HM, Cathomas R, Compérat EM, Cowan NC, Gakis G, Hernández V, Linares Espinós E, Lorch A, Neuzillet Y, et al. European association of urology guidelines on muscle-invasive and metastatic bladder cancer: summary of the 2020 guidelines. Eur Urol. 2021;79(1):82–104.

    Article  CAS  PubMed  Google Scholar 

  6. Pfail JL, Katims AB, Alerasool P, Sfakianos JP. Immunotherapy in non-muscle-invasive bladder cancer: current status and future directions. World J Urol. 2021;39(5):1319–29.

    Article  CAS  PubMed  Google Scholar 

  7. Wilkins A, Ost P, Sundahl N. Is there a benefit of combining immunotherapy and radiotherapy in bladder cancer? Clin Oncol. 2021;33(6):407–14.

    Article  CAS  Google Scholar 

  8. Wong RS. Apoptosis in cancer: from pathogenesis to treatment. J Exp Clin Cancer Res. 2011;30(1):87.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Vaghari-Tabari M, Ferns GA, Qujeq D, Andevari AN, Sabahi Z, Moein S. Signaling, metabolism, and cancer: an important relationship for therapeutic intervention. J Cell Physiol. 2021;236(8):5512–32.

    Article  CAS  PubMed  Google Scholar 

  10. McKnight JJ, Gray SB, O’Kane HF, Johnston SR, Williamson KE. Apoptosis and chemotherapy for bladder cancer. J Urol. 2005;173(3):683–90.

    Article  CAS  PubMed  Google Scholar 

  11. Hatogai K, Sweis RF. The tumor microenvironment of bladder cancer. Adv Exp Med Biol. 2020;1296:275–90.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Ghasemi H, Mousavibahar SH, Hashemnia M, Karimi J, Khodadadi I, Tavilani H. Transitional cell carcinoma matrix stiffness regulates the osteopontin and YAP expression in recurrent patients. Mol Biol Rep. 2021;48(5):4253–62.

    Article  CAS  PubMed  Google Scholar 

  13. Chen X, Chen H, He D, Cheng Y, Zhu Y, Xiao M, Lan H, Wang Z, Cao K. Analysis of tumor microenvironment characteristics in bladder cancer: implications for immune checkpoint inhibitor therapy. Front Immunol. 2021;12: 672158.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Yan X, Wu C, Chen T, Santos MM, Liu CL, Yang C, Zhang L, Ren J, Liao S, Guo H, et al. Cathepsin S inhibition changes regulatory T-cell activity in regulating bladder cancer and immune cell proliferation and apoptosis. Mol Immunol. 2017;82:66–74.

    Article  CAS  PubMed  Google Scholar 

  15. Yang M, Wang B, Hou W, Yu H, Zhou B, Zhong W, Liu Z, Li J, Zeng H, Liu C, et al. Negative effects of stromal neutrophils on T cells reduce survival in resectable urothelial carcinoma of the bladder. Front Immunol. 2022;13: 827457.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Lee JS, Leem SH, Lee SY, Kim SC, Park ES, Kim SB, Kim SK, Kim YJ, Kim WJ, Chu IS. Expression signature of E2F1 and its associated genes predict superficial to invasive progression of bladder tumors. J Clin Oncol. 2010;28(16):2660–7.

    Article  CAS  PubMed  Google Scholar 

  17. Riester M, Werner L, Bellmunt J, Selvarajah S, Guancial EA, Weir BA, Stack EC, Park RS, O’Brien R, Schutz FA, et al. Integrative analysis of 1q233 copy-number gain in metastatic urothelial carcinoma. Clin Cancer Res. 2014;20(7):1873–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Lindgren D, Sjödahl G, Lauss M, Staaf J, Chebil G, Lövgren K, Gudjonsson S, Liedberg F, Patschan O, Månsson W, et al. Integrated genomic and gene expression profiling identifies two major genomic circuits in urothelial carcinoma. PLoS ONE. 2012;7(6): e38863.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Sjödahl G, Lauss M, Lövgren K, Chebil G, Gudjonsson S, Veerla S, Patschan O, Aine M, Fernö M, Ringnér M, et al. A molecular taxonomy for urothelial carcinoma. Clin Cancer Res. 2012;18(12):3377–86.

    Article  PubMed  Google Scholar 

  20. Choi W, Porten S, Kim S, Willis D, Plimack ER, Hoffman-Censits J, Roth B, Cheng T, Tran M, Lee IL, et al. Identification of distinct basal and luminal subtypes of muscle-invasive bladder cancer with different sensitivities to frontline chemotherapy. Cancer Cell. 2014;25(2):152–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE III, Koeppen H, Astarita JL, Cubas R, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554(7693):544–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Yuan X, Zhou J, Zhou L, Huang Z, Wang W, Qiu J, Yang Q, Zhang C, Ma M. Apoptosis-related gene-mediated cell death pattern induces immunosuppression and immunotherapy resistance in gastric cancer. Front Genet. 2022;13: 921163.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7): e47.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Emura T, Matsui S, Chen HY. compound.Cox: univariate feature selection and compound covariate for predicting survival. Comput Methods Programs Biomed. 2019;168:21–37.

    Article  PubMed  Google Scholar 

  25. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics (Oxford, England). 2010;26(12):1572–3.

    CAS  PubMed  Google Scholar 

  26. Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Prot Sci Publ Prot Soc. 2019;28(11):1947–51.

    Article  CAS  Google Scholar 

  28. Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51(D1):D587-d592.

    Article  PubMed  Google Scholar 

  29. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation. 2021;2(3):100141.

    CAS  PubMed  PubMed Central  Google Scholar 

  30. Finotello F, Trajanoski Z. Quantifying tumor-infiltrating immune cells from transcriptomics data. Cancer Immunol Immunother CII. 2018;67(7):1031–40.

    Article  CAS  PubMed  Google Scholar 

  31. Xu L, Deng C, Pang B, Zhang X, Liu W, Liao G, Yuan H, Cheng P, Li F, Long Z, et al. TIP: a web server for resolving tumor immunophenotype profiling. Can Res. 2018;78(23):6575–80.

    Article  CAS  Google Scholar 

  32. David CC, Jacobs DJ. Principal component analysis: a method for determining the essential dynamics of proteins. Methods Mol Biol (Clifton, NJ). 2014;1084:193–226.

    Article  CAS  Google Scholar 

  33. Liu Y, Wang J, Li L, Qin H, Wei Y, Zhang X, Ren X, Ding W, Shen X, Li G, et al. AC010973.2 promotes cell proliferation and is one of six stemness-related genes that predict overall survival of renal clear cell carcinoma. Sci Rep. 2022;12(1):4272.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Yu L, Shen H, Ren X, Wang A, Zhu S, Zheng Y, Wang X. Multi-omics analysis reveals the interaction between the complement system and the coagulation cascade in the development of endometriosis. Sci Rep. 2021;11(1):11926.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Zhang D, Xu X, Wei Y, Chen X, Li G, Lu Z, Zhang X, Ren X, Wang S, Qin C. Prognostic role of DNA damage response genes mutations and their association with the sensitivity of olaparib in prostate cancer patients. Cancer Control J Moffitt Cancer Center. 2022;29:10732748221129452.

    Google Scholar 

  36. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005;102(43):15545–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Chen Y, Wang HH, Chang HH, Huang YH, Wang JR, Changchien CY, Wu ST. Guggulsterone induces apoptosis and inhibits lysosomal-dependent migration in human bladder cancer cells. Phytomed Int J Phytother Phytopharmacol. 2021;87: 153587.

    CAS  Google Scholar 

  38. Yin Z, Li J, Kang L, Liu X, Luo J, Zhang L, Li Y, Cai J. Epigallocatechin-3-gallate induces autophagy-related apoptosis associated with LC3B II and Beclin expression of bladder cancer cells. J Food Biochem. 2021;45(6): e13758.

    Article  CAS  PubMed  Google Scholar 

  39. Wang W, Green M, Choi JE, Gijón M, Kennedy PD, Johnson JK, Liao P, Lang X, Kryczek I, Sell A, et al. CD8(+) T cells regulate tumour ferroptosis during cancer immunotherapy. Nature. 2019;569(7755):270–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Ohue Y, Nishikawa H. Regulatory T (Treg) cells in cancer: can Treg cells be a new therapeutic target? Cancer Sci. 2019;110(7):2080–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Sakaguchi S. Regulatory T cells: key controllers of immunologic self-tolerance. Cell. 2000;101(5):455–8.

    Article  CAS  PubMed  Google Scholar 

  42. Pan Y, Yu Y, Wang X, Zhang T. Tumor-associated macrophages in tumor immunity. Front Immunol. 2020;11: 583084.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Cassetta L, Pollard JW. Tumor-associated macrophages. Curr Biol. 2020;30(6):R246-r248.

    Article  CAS  PubMed  Google Scholar 

  44. Chen D, Zhang X, Li Z, Zhu B. Metabolic regulatory crosstalk between tumor microenvironment and tumor-associated macrophages. Theranostics. 2021;11(3):1016–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Leblond MM, Zdimerova H, Desponds E, Verdeil G. Tumor-associated macrophages in bladder cancer: biological role, impact on therapeutic response and perspectives for immunotherapy. Cancers. 2021;13(18):4712.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Chan TA, Yarchoan M, Jaffee E, Swanton C, Quezada SA, Stenzinger A, Peters S. Development of tumor mutation burden as an immunotherapy biomarker: utility for the oncology clinic. Ann Oncol. 2019;30(1):44–56.

    Article  CAS  PubMed  Google Scholar 

  47. Lopez-Beltran A, Cimadamore A, Blanca A, Massari F, Vau N, Scarpelli M, Cheng L, Montironi R. Immune checkpoint inhibitors for the treatment of bladder cancer. Cancers. 2021;13(1):131.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Mancini M, Righetto M, Noessner E. Checkpoint inhibition in bladder cancer: clinical expectations, current evidence, and proposal of future strategies based on a tumor-specific immunobiological approach. Cancers. 2021;13(23):6016.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Zhu J, Armstrong AJ, Friedlander TW, Kim W, Pal SK, George DJ, Zhang T. Biomarkers of immunotherapy in urothelial and renal cell carcinoma: PD-L1, tumor mutational burden, and beyond. J Immunother Cancer. 2018;6(1):4.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Moon C, Gordon M, Moon D, Reynolds T. Microsatellite instability analysis (MSA) for bladder cancer: past history and future directions. Int J Mol Sci. 2021;22(23):12864.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was received Guangxi Natural Science Foundation Project, China (2018GXNSFAA138061).

Author information

Author notes

  1. Liquan Zhou and Guanglong Xu have contributed equally to this work and share the first authorship

  2. Liquan Zhou and Yong Tang have contributed equally to this work and share the corresponding authorship

Authors and Affiliations

  1. Department of Urology, The Second Affiliated Hospital, Guangxi Medical University, Nanning, 530006, Guangxi, China

    Liquan Zhou, Guanglong Xu, Fu Huang, Wenyuan Chen & Jiange Zhang

  2. Department of Urology, Wuming Hospital, Guangxi Medical University, Nanning, 530199, Guangxi, China

    Yong Tang

Authors
  1. Liquan Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Guanglong Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Fu Huang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Wenyuan Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Jiange Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Yong Tang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

LZ, GX and YT wrote the main manuscript text and FH, WC, and JZ prepared Figs. 1, 2, 3, 4, 5, 6, 7 and 8. All authors reviewed the manuscript.

Corresponding author

Correspondence to Yong Tang.

Ethics declarations

Ethics approval and consent to participate

The study was approved by the Ethical Committee of the Wuming Hospital, Guangxi Medical University (WM-2022-15) and followed the guidelines of the Declaration of Helsinki. Written informed consent was obtained from all patients after explanation of the study and possible consequences of the study.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1. Supplementary table 1.

Clinical characteristics of the BLCA patients used in this study.

Additional file 2. Supplementary table 2.

List of 91 apoptosis related genes used in this study.

Additional file 3. Supplementary table 3.

Differentially expressed apoptosis related genes between tumor and normal tissue inTCGA-BLCA dataset.

Additional file 4. Supplementary table 4.

The results of the univariate Cox regression analysis between gene expression and OS inTCGA-BLCA dataset.

Additional file 5. Supplementary table 5.

List of 2859 ARGs phenotype-related prognostic genes used in this study.

Additional file 6. Supplementary table 6.

List of 13 ARGs phenotype-related prognostic genes used in this study in IMvigor210cohort.

Additional file 7. Supplementary table 7.

Summary of detailed clinical information and apoptosisscore in IMvigor210 (mUC) cohort.

Additional file 8. Supplementary figure 1.

(A-D) The relationship between four genes and immune infiltration. Supplementary figure 2. (A-F) Unsupervised consensus clustering based on 4 ARGs prognostic genes in a meta cohort. Supplementary figure 3. Differences in checkpoint expression between ARGs.cluster.A and ARGs.cluster.B groups. Supplementary figure 4. (A–F) Unsupervised consensus clustering based on 4 ARGs in IMvigor210 cohort.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhou, L., Xu, G., Huang, F. et al. Apoptosis related genes mediated molecular subtypes depict the hallmarks of the tumor microenvironment and guide immunotherapy in bladder cancer. BMC Med Genomics 16, 88 (2023). https://doi.org/10.1186/s12920-023-01525-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12920-023-01525-8

Keywords

BMC Medical Genomics

ISSN: 1755-8794

Contact us