Comprehensive analysis of tumor immune-related gene signature for predicting prognosis, immunotherapy, and drug sensitivity in bladder urothelial carcinoma

Introduction

The incidence of bladder urothelial carcinoma (BLCA), the tenth most frequent malignancy overall, is around four times higher in males than in females (1). Tobacco use currently remains the leading risk factor of BLCA, accounting for almost 50% of cases (2). Following surgical removal and, when required, with supplementary intravesical maintenance immunotherapy or chemotherapy, most individuals present with a superficial, non-muscle-invading tumor, leading to a promising prognosis (3). Yet, bladder cancer that invades the muscle is marked by elevated morbidity and mortality rates, making it a profoundly fatal ailment. Currently, multimodal therapy that combines a cystectomy with neoadjuvant chemotherapy and immunotherapy is effective in treating metastatic bladder cancer (4).

Current treatments for advanced BLCA involve immune checkpoint inhibitors (ICIs), specifically inhibitors targeting programmed cell death protein 1 (PD-1) and programmed death ligand 1 (PD-L1) (5). ICI therapy for muscle-invasive and advanced conditions has demonstrated considerable survival advantages (6). However, immune checkpoint blockade (ICB) is often used to treat bladder cancer, with little to no therapeutic effect for the majority of patients (7). The effectiveness of ICI is influenced by various factors, including the tumor microenvironment (TME), tumor mutation burden (TMB), tumor immune dysfunction and exclusion (TIDE), and microsatellite instability (MSI) (8-11). Few accurate biomarkers can forecast how a particular patient will react to ICI treatment. The individualization of immunotherapy for BLCA patients may be made possible by the identification of suitable prognostic indicators. Regrettably, there are a paucity of established biomarkers that reliably signify a patient’s response to treatment, and we still do not fully understand the TME, TMB, and TIDE of BLCA. There is an urgent need for more accurate prognostic signatures.

In our study, we conducted an exhaustive analysis of the bladder carcinoma immune gene signature (BCIGS). Our aim was to establish a robust prognostic framework to appraise prognostic relevance, gauge immunotherapy reactions, and ascertain medication susceptibility in BLCA patients. This endeavor aids experts in initiating early diagnostic procedures and commencing treatment. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1053/rc).

Methods

Data collection and immune gene expression variability

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The Cancer Genome Atlas (TCGA) transcriptome, clinical, and mutation data can be accessed at https://portal.gdc.cancer.gov; it comprises 411 samples from patients with BLCA and 19 samples from normal bladder tissue. A list of 2,660 immune genes was compiled using information from the innateDB (https://www.innatedb.com) and ImmPort (www.immport.org/home) websites. We employed the limma R package (12) and established the criteria as |log₂ fold change (FC)| ≥1 and false discovery rate (FDR) <0.05. Subsequently, we utilized the Wilcoxon test to detect differentially expressed genes (DEGs) and immune genes (DEIGs) by comparing standard bladder samples to BLCA samples. Finally, we produced a heatmap of the DEG and DEIG using the R package pheatmap.

Weighted gene co-expression network analysis (WGCNA), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of DEIGs

We performed GO enrichment and KEGG pathway analyses for the 224 DEIGs using the “clusterProfiler” (13), “enrichplot”, “ggplot2”, and “org.Hs.eg.db” R packages, with P and Q value thresholds set at 0.05. Additionally, we also performed WGCNA for the expression data of DEIGs using the “WGCNA” (14) and “limma” R packages. Using cluster analysis, the closely related genes were separated into gene modules. Gene modules smaller than 15 were not taken into consideration. We set a cut height threshold at 0.3 to analyze the gene modules. Subsequently, we created a heatmap to visualize the Pearson correlation coefficients and their corresponding P values between module eigengenes (MEs) and tumor status. The key module linked to the BCIGS was identified for subsequent studies.

Prognosis-related DEIGs identification

We performed a univariate Cox regression analysis and conducted a Kaplan-Meier test on the majority of the significant module genes identified through WGCNA, using the “survminer” and “survival” R packages. Genes that showed significant results in the univariate Cox regression analysis, with a P value <0.05, were regarded as potential immune genes associated with prognosis. The prognostic significance of these DEIGs was assessed using the Kaplan-Meier test, with genes having a P value below 0.05 identified as potential central immune genes with prognostic relevance.

Prognostic model establishment and evaluation

For the TCGA BLCA patient cohort, we evenly partitioned them into training and testing sets. Subsequently, we employed the glmnet (15), survival, caret, survminer, and survivalROC R packages to construct a multivariate Cox proportional hazards regression model using the DEIGs associated with prognosis from the training set. This model yielded the risk score essential for patient prognosis assessment. The risk score was computed as follows: risk score = α1 × Ex1 + α2 × Ex2 + … + αi × Exi, where Ex is the expression level, and α denotes the regression coefficient. We divided the training dataset into two groups, distinguishing between low- and high-risk, using the median risk score obtained from the model’s equation. We then applied the same median risk score to classify the testing dataset and external validation datasets GSE39281 (16) and IMvigor210 (17). These datasets served as validation cohorts to confirm the predictive efficacy of our model. In both the training and testing datasets, as well as the external validation datasets GSE39281 and IMvigor210, we used the survival and survminer R packages to assess overall survival (OS) differences between risk groups via log-rank tests. To evaluate the model’s precision, we generated a risk plot and heatmap using the pheatmap package, along with a receiver operating characteristic (ROC) curve using the survivalROC package. We conducted a Cox regression analysis using the survival R package to assess the predictive capabilities of various clinical features in the BLCA patients within the training dataset. Moreover, the ComplexHeatmap R package enabled us to discern variations in numerous clinical parameters between the training dataset’s high- and low-risk categories.

Gene set enrichment analysis and gene set variation analysis (GSVA)

We carried out a gene set enrichment analysis (GSEA) on the training dataset’s high- and low-risk groups using the limma and org packages. For this assessment, gene sets from the C2 (c2.cp.kegg.v7.4.symbols.gmt) and C5 (c5.go.v7.4.symbols) families were specifically selected. 1,000 permutations were used to compute the FDR q values; a gene set was deemed highly enriched if its normalized enrichment score (NES) was more than 1 and the FDR was less than 0.25. GSVA is a method used to measure changes in pathway activity over different sample groups (18). To calculate the relationship between KEGG pathway and model genes and risk scores, we employed the GSVA approach. We employed the GSVA package (19), along with limma, GSEABase, GSVA, reshape2, and ggplot2 for the analysis.

Immune cell infiltration, immune-related function, and survival analysis

Utilizing normalized gene expression data, the cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) analytical tool is adept at characterizing the cellular composition of multifaceted tissues (20). We employed CIBERSORT to evaluate and quantify the proportional presence of 22 immune cell types infiltrating each BLCA specimen. Additionally, we sourced the relative prevalence data of immune cell infiltrates in TCGA samples from Thorsson et al.’s documentation (21). For the Kaplan-Meier survival analysis, which involved comparing immune cell types infiltrating the high- and low-risk groups, we made use of the R packages “survminer”, “limma”, and “survival”. To further investigate the relationship between patient risk scores and the extent of immune cell infiltration, we utilized the “scales”, “ggtext”, “ggplot2”, “ggpubr”, and “tidyverse” R packages. Bindea et al.’s study (22) furnished us with 29 immunological datasets detailing specific immune cell categories, immune-focused pathways, and immune-driven functionalities. To compute the infiltration score per BLCA specimen amidst the high- and low-risk categories, using the 29 immunological datasets, we utilized the “limma”, “GSVA”, “ggpubr”, “GSEABase”, and “reshape2” R packages. With the “limma”, “survival”, and “survminer” R packages, we executed a survival analysis juxtaposing immune-associated operations between the high- and low-risk groupings.

TMB, TME, and TIDE analysis

In this study, we performed a comprehensive analysis of the TME using advanced statistical techniques and the R packages “estimate” and “limma”. The primary goal was to understand the composition of the TME and its potential association with different risk groups. To achieve this, we first assessed the presence of immune and stromal components within the TME, quantifying them using metrics such as estimation of stromal and immune cells in malignant tumor tissues using expression data score (ESTIMATEScore), StromalScore, and ImmuneScore. By doing so, we could distinguish between portions of immune-stromal, pure stromal, and solely immune elements, providing valuable insights into the TME’s heterogeneity. To discern disparities between high- and low-risk cohorts, Wilcoxon rank-sum tests were conducted on these scores, employing additional R packages such as reshape2 and ggpubr. Furthermore, external data including TIDE scores, MSI scores, T cell dysfunction scores, and T cell exclusion scores were sourced to enrich their analysis and gain a more holistic understanding of the tumors under investigation. The connection between the risk score and TMB, a crucial factor in cancer biology, was also explored. Using a combination of R packages such as limma, ggExtra, ggplot2, and ggpubr, potential differences in TMB between the risk groups were probed, shedding light on a potentially crucial aspect of cancer biology and risk assessment.

Immune subtype and immune checkpoint analysis

We examined over 10,000 TCGA cancers and identified six immune subgroups with distinct characteristics and prognoses utilizing cluster analysis of a large list of immune expression signatures reported by prior studies: C1, healing of wounds; C2, dominance of IFN-γ; C3, inflammation; and C4, lymphocyte depletion (21). In our study, we utilized the RColorBrewer package in R to examine potential differences in immune subtypes between high- and low-risk groups. Spearman correlation coefficients and their corresponding P values were calculated for the expression levels of immune checkpoint genes. To gain deeper insights into the relationships among model genes, immune checkpoint genes, and risk scores, we harnessed the power of the ggplot2, reshape2, and limma packages in R. This comprehensive analysis allowed us to unravel key associations within our dataset.

Clinical application of prognostic model

We used the rms R package to create a nomogram based on our prognostic model, allowing us to predict the probability of OS. To assess the potential efficacy of immunotherapy, we referenced The Cancer Immunome Atlas (TCIA) database (https://tcia.at/) for BLCA’s cytotoxic T-lymphocyte antigen 4 (CTLA-4) and PD-1 blockade treatment details. To identify potential differences in treatments targeting CTLA-4 and PD-1 across BLCA patient groups delineated by risk scores (high/low), we used the ggpubr R package. The “pRRophetic” package was used for drug response prediction. This package utilizes half maximum inhibitory concentration (IC₅₀) values specific to each BLCA patient, derived from the Genomics of Drug Sensitivity in Cancer (GDSC) database (23). Then, using freely accessible drug sensitivity databases, we identified potential medications that target the immune gene profile associated with bladder cancer.

Statistical analysis

Statistical evaluations were executed with R version 4.0.2 and relevant packages. For categorical data, we employed the Chi-squared test, whereas the independent t-test was utilized for contrasting continuous variables across two groups. The log-rank test facilitated the Kaplan-Meier survival analysis. For regression studies, hazard ratios (HRs) and their associated 95% confidence intervals (CIs) were determined. The significance was denoted as: *, P<0.05; **, P<0.01; and ***, P<0.001.

Results

DEIGs identification and GO, KEGG analyses

The detailed study design is shown in Figure 1. Relevant R programs were utilized for data processing and the selection of DEGs and DEIGs. From this analysis, 3,914 DEGs and 224 DEIGs satisfied the criteria set for the study. Figure 2A,2B show the heatmaps for DEGs and DEIGs, respectively. Using R and correlation tools, we delved deep into the GO and KEGG enrichment studies of DEIGs. The GO analysis underscored a notable enrichment in DEIGs for functionalities such as inflammation modulation, chemokine-mediated signaling, cell chemotaxis, cellular calcium ion homeostasis, external stimuli response amplification, and responses to chemokines, as illustrated in Figure 2C. Meanwhile, the KEGG exploration pinpointed significant DEIG engagement in pathways such as cytokine-cytokine receptor interactions, calcium signaling, neuroactive ligand-receptor interplays, interactions with viral proteins and cytokine receptors, melanoma, and the pathways of PI3K-Akt, Ras, and MAPK, all shown in Figure 2D.

Figure 1 Flowchart of this study. TCGA, The Cancer Genome Atlas; BLCA, bladder urothelial carcinoma; DEGs, differentially expressed genes; DEIGs, differentially expressed immune genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; WGCNA, weighted gene co-expression network analysis; ROC, receiver operating characteristic; GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; TIDE, tumor immune dysfunction and exclusion; TME, tumor microenvironment; TMB, tumor mutation burden.

Figure 2 Analysis of differentially expressed immune genes and WGCNA. (A) The heatmap of differentially expressed genes. (B) The heatmap of DEIGs. (C) Circular plot of GO enrichment analysis of DEIGs. (D) Circular plot of KEGG pathway analysis of DEIGs. (E) The merged graphical result shows the final clustering of samples under different network modules. (F) Gene modules related to BLCA obtained by WGCNA. WGCNA, weighted gene co-expression network analysis; GO, Gene Ontology; FC, fold change; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEIGs, differentially expressed immune genes; BLCA, bladder urothelial carcinoma.

WGCNA of DEIGs

We carefully evaluated the soft-thresholding power to maintain scale independence and mean linkage in co-expression networks. The network topology was characterized by a soft-threshold power of 0.9 for the scale-free index and 3 for mean connectivity, resulting in a nearly scale-free state (Figure S1). Five integrated co-expression modules were subsequently built (Figure 2E). The relationship between normal, malignant, and MEs was subsequently established. The results showed a direct association between the tumor and the blue and brown modules, as depicted in Figure 2F. In addition, we found that the tumor was adversely linked with the turquoise and yellow modules (Figure 2F). The WGCNA found 155 DEIGs in the turquoise module.

Prognosis-related DEIGs selection

Results from the univariate Cox regression analysis (shown in Figure 3A) yielded 32 potential hub DEIGs linked to prognosis. These were then further analyzed for their prognostic relevance. Cox regression analysis involving multiple variables was conducted on the 32 DEIGs, leading to the identification of 11 hub DEIGs that were deemed independent predictors of BLCA, as illustrated in Figure 3B and outlined in Table S1. The Kaplan-Meier survival analysis (represented in Figure 3C) indicated 32 hub DEIGs related to prognosis.

Figure 3 Survival analysis of DEIGs. (A) Univariate and (B) multivariate analyses of the gene expression level and OS. (C) Kaplan-Meier survival analysis for OS between high and low gene expression groups. CI, confidence interval; DEIGs, differentially expressed immune genes; OS, overall survival.

Prognosis-related model construction and validation

Out of the 404 BLCA patients, 204 were randomly assigned to the training dataset, with the rest allocated to the testing dataset. The predictive model was constructed using the 11 hub DEIGs for prognosis with the training data. To determine the risk score for every patient, we used the subsequent formula: risk score = (−0.4595 × ExpCXCL12) + (1.5529 × ExpLCN1) + (−2.6484 × ExpFGF10) + (0.8514 × ExpSEMA6D) + (−0.8561 × ExpCORT) + (−1.4912 × ExpNTF3) + (0.7531 × ExpPGR) + (0.5880 × ExpPTGER3) + (0.5979 × ExpUSP2) + (0.4599 × ExpJAM3) + (0.4272 × ExpC7).

Next, our aim was to assess the prediction potential. Figure 4A,4B depict patient survival status, risk scores, and an expression heatmap based on the 11 DEIGs for both low- and high-risk groups in both the training and testing datasets. Results from the training and testing set showed that the high-risk group had a reduced OS relative to the low-risk group, as depicted in Figure 4C,4D. The 11 hub DEIGs’ predictive usefulness was verified by ROC analysis. In the training and testing dataset, the model achieved an area under the ROC curve (AUC) of 0.712 and 0.631, respectively, indicating improved diagnostic accuracy, as demonstrated in Figure 4E,4F.

Figure 4 Validation of the BCIGS predicting model and performance analysis in BLCA. (A,B) Risk score distribution, survival status, and the expression of 11 DEIGs for patients in low- and high-risk groups from TCGA-BLCA training dataset (A) and TCGA-BLCA testing dataset (B). (C,D) Kaplan-Meier survival curves showing the comparison of OS between the low- and high-risk groups from TCGA-BLCA training dataset (C) and TCGA-BLCA testing dataset (D). (E,F) ROC curves of TCGA-BLCA training dataset (E) and TCGA-BLCA testing dataset (F). (G) Univariate Cox analysis of clinical factors and the risk score. (H) Multivariable Cox analysis of clinical factors and the risk score. (I) Clinical characteristics of different risk subgroups. *, P<0.05. ROC, receiver operating characteristic; AUC, area under the receiver operating characteristic curve; CI, confidence interval; BCIGS, bladder carcinoma immune gene signature; BLCA, bladder urothelial carcinoma; DEIGs, differentially expressed immune genes; TCGA, The Cancer Genome Atlas; OS, overall survival; T, tumor; N, node; M, metastasis.

In this study, Cox regression analysis played a pivotal role in shedding light on the prognosis of BLCA patients. Through the use of this statistical method, we meticulously examined various clinical attributes and their impact on OS. Our findings, as demonstrated in Figure 4G, indicated that risk score, age, and stage were significantly associated with OS when assessed individually through univariate Cox regression. Moreover, comprehensive multivariate Cox analysis, as depicted in Figure 4H, confirmed that these factors independently influenced the OS of BLCA patients in our training dataset. Furthermore, Figure 4I illuminated substantial variations in T stage, underscoring the crucial role of clinical traits in distinguishing high- and low-risk groups among BLCA patients. This robust statistical approach provides valuable insights into the complex interplay of clinical attributes and their implications for patient prognosis.

To validate the predictive capability of our risk score model, the same approach was applied to the testing dataset and external validation datasets GSE39281 and IMvigor210. In the external validation datasets GSE39281 and IMvigor210, the OS of the high-risk group was significantly lower than that of the low-risk group (Figure 5A,5B), with AUCs of 0.609 and 0.563 (Figure 5C,5D), respectively. The patient survival status, risk scores, and expression heatmap based on the 11 DEIGs for both the high- and low-risk groups are shown in Figure 5E,5F. These suggest potential avenues for enhancing the model’s prognostic accuracy and reliability.

Figure 5 External validation of the BCIGS predicting model in BLCA. (A,B) Kaplan-Meier survival curves showing the comparison of OS between the low- and high-risk groups from GSE39281 dataset (A) and IMvigor210 dataset (B). (C,D) ROC curves of GSE39281 dataset (C) and IMvigor210 dataset (D). (E,F) Risk score distribution, survival status, and the expression of 11 DEIGs for patients in low- and high-risk groups from GSE39281 dataset (E) and IMvigor210 dataset (F). ROC, receiver operating characteristic; AUC, area under the receiver operating characteristic curve; BCIGS, bladder carcinoma immune gene signature; BLCA, bladder urothelial carcinoma; OS, overall survival; DEIGs, differentially expressed immune genes.

The results of GSEA and GSVA

The KEGG GSEA highlighted key pathways in our study. GRAFT_VERSUS_HOST_DISEASE was strongly linked to the high-risk group, whereas LINOLEIC_ACID_METABOLISM was prominent in the low-risk group (Figure 6A,6B). REGULATION_OF_HUMORAL_IMMUNE_RESPONSE was significant in the high-risk group, and AROMATASE_ACTIVITY stood out in the low-risk group, providing insights into their distinct molecular profiles (Figure 6C,6D). We used GSVA to compare KEGG pathways, model gene expression levels, and risk scores to investigate the related cancer hallmark pathways in relation to BCIGS. A total of 21 KEGG pathways showed good correlations with model risk ratings above the predetermined limit (Figure 6E). These findings supported further investigation by confirming that BCIGS was connected to the majority of tumor-related pathways.

Figure 6 GSEA and GSVA. (A) GSEA shows top enriched pathways in high-risk group. (B) GSEA shows top enriched pathways in low-risk group. (C) GSEA shows top enriched GO terms in the high-risk group. (D) GSEA shows top GO terms in the low-risk group. (E) Heatmap shows the correlation between KEGG pathway and model genes and risk scores. *, P<0.05; **, P<0.01; ***, P<0.001. KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; GO, Gene Ontology.

BCIGS was associated with BLCA immune status

In order to gain a deeper understanding of the immunological landscape among patients at varying levels of risk, the researchers in this study employed a comprehensive analysis approach. They utilized the single-sample gene set enrichment analysis (ssGSEA) method to assess the infiltration scores of multiple immune cell subtypes and to explore immune-related functions and pathways. The findings, as depicted in Figure 7A,7B, provide valuable insights into the intricate relationship between the patients’ risk profiles and their immune systems. One of the key observations highlighted in the study was the significant disparity in immune cell activity between high- and low-risk patients. Various immune cell subtypes, such as CD8⁺ T cells, macrophages, and regulatory T cells, exhibited notably elevated scores in high-risk individuals, as illustrated in Figure 7C. This heightened immune activity in high-risk patients suggests a potential link between immune system activation and the underlying risk factors. However, the most striking aspect of the study emerged when considering patient outcomes. The Kaplan-Meier survival analysis revealed that patients with increased levels of certain immune cell types, such as macrophages M0 and M2, experienced shorter survival spans. Conversely, patients with elevated expressions of other immune cell types, such as dendritic cells activated and T cells follicular helper, exhibited longer survival durations (Figure 7D). Similarly, immune-related functions and pathways also played a crucial role in determining patient survival, with some showing a clear association with reduced survival yet others being linked to prolonged survival (Figure 7E). In summary, this research sheds light on the intricate interplay between the immune system and patients’ risk profiles, providing valuable insights that may have implications for prognosis and potential therapeutic interventions.

Figure 7 The difference in immune cell infiltration and immune-related function/pathways between low- and high-risk groups. (A) The relative percent of 22 immune cell types in the low- and high-risk groups of BLCA patients. (B) The correlation analysis of different immune cell types and risk score of patients by various statistical algorithms. (C) Box plot showing the scores of immune-related function/pathways between low- and high-risk groups in BLCA patients. (D) Kaplan-Meier survival analysis of the different immune cells between high- and low-expression groups. (E) Kaplan-Meier survival analysis of the different immune-related function/pathways between high- and low-expression groups. *, P<0.05; **, P<0.01; ***, P<0.001. BLCA, bladder urothelial carcinoma.

The results of TME, TIDE, TMB, and immune subtype

A comprehensive study was conducted comparing high- and low-risk groups concerning the development of metastasis in cancer patients, which focused on various scoring systems to evaluate the potential for metastasis development. Notably, the high-risk group exhibited significantly higher scores in key areas, including StromalScore, ImmuneScore, and ESTIMATEScore, compared to the low-risk group (Figure 8A). These findings imply that individuals with elevated scores in these categories have an increased likelihood of developing metastasis. Moreover, the study delved into several other factors that could influence metastasis and the effectiveness of immunotherapy. Interestingly, although factors such as T cell dysfunction, TIDE score, and MSI score did not exhibit significant differences between the high- and low-risk groups (Figure 8B-8D), the T cell exclusion score was higher in the high-risk group (Figure 8E). This suggests that individuals in the high-risk category, as determined by the BCIGS, might be particularly suitable candidates for immunotherapy, especially ICIs. Furthermore, the research revealed distinct immunological clusters within the two risk groups, with the C2 subtype dominating the high-risk group and C1 being more prevalent in the low-risk category (Figure 8F). This information could have important implications for tailoring treatment approaches based on the patient’s risk profile. Additionally, the study explored the relationships between immune checkpoint genes and other genes within the model. Positive correlations were observed between several immune checkpoint genes and specific genes such as CXCL12, LCN1, FGF10, SEMA6D, PGR, PTGER3, JAM3, and C7, whereas CORT showed an inverse relationship with most immune checkpoint genes (Figure 8G). Lastly, the study investigated the TMB and its association with the risk score, revealing a negative correlation between the two in bladder cancer (Figure 8H,8I). This suggests that as the risk score increases, the TMB tends to decrease, providing valuable insights into the underlying mechanisms of metastasis in this context.

Figure 8 Characterization of prognostic model. (A) TME scores, (B) TIDE, (C) dysfunction, (D) MSI scores, (E) exclusion, and (F) immune subtype were analyzed in the low- and high-risk groups of BLCA patients. (G) The correlation analysis between 11 signature genes and immune checkpoint genes. (H) Comparison of TMB between low- and high-risk groups. (I) The correlation analysis of risk score and TMB. ns, P≥0.05; *, P<0.05; **, P<0.01; ***, P<0.001. TME, tumor microenvironment; TIDE, tumor immune dysfunction and exclusion; MSI, microsatellite instability; TCGA, The Cancer Genome Atlas; TMB, tumor mutation burden; BLCA, bladder urothelial carcinoma.

Clinical application

Using the aforementioned findings, a nomogram for predicting OS (Figure 9A) was constructed by incorporating the expression of CXCL12, LCN1, FGF10, SEMA6D, CORT, NTF3, PGR, PTGER3, USP2, JAM3, and C7. We also assessed the model’s propensity for prediction. According to Figure 9B, the AUC for the 1-, 2-, and 3-year OS stood at 0.727, 0.772, and 0.765 respectively. Figure 9C,9D demonstrate that the predictive capability of our model surpasses that of other models, underscoring its superior diagnostic and predictive accuracy. The results of predicting the response to anti-CTLA4 and anti-PD-1 immunotherapy (Figure 9E-9G) indicated that anti-CTLA4 treatment had a more significant impact on the low-risk patients compared to the high-risk patients. This suggests that low-risk patients may potentially exhibit a more favorable response to anti-CTAL4 immunotherapy. In our study of patients with BLCA, we employed the pRRophetic technique to forecast the response to chemotherapy by utilizing the IC₅₀ data from the GDSC database. Our analysis identified a total of 20 chemical compounds with significantly diverse responses between the low- and high-risk groups, as depicted in Figure 10. The most significant distinction was observed in the IC₅₀ value of the TW.37 drug between the low- and high-risk patient groups. These findings suggest that these little molecular compounds might serve as BLCA therapy agents; however, more research is still required. Our findings provide promising molecular chemotherapeutic drugs for BLCA patients.

Figure 9 The prognostic prediction of nomogram and the evaluation of immunotherapy response. (A) The nomogram for the prediction of OS. (B) ROC curves were constructed for 1-, 2-, and 3-year OS of the prognostic model. (C) Comparison of ROC curves among prognostic model, TIDE model, and TIS model. (D) Comparison of ROC curves between prognostic model and other clinical features. (E-G) The difference of predicting anti-CTLA4 and anti-PD-1 immunotherapy response between high-risk and low-risk group. AUC, area under the receiver operating characteristic curve; TIDE, tumor immune dysfunction and exclusion; TIS, tumor inflammation signature; OS, overall survival; ROC, receiver operating characteristic; CTLA4, cytotoxic T-lymphocyte antigen 4; PD-1, programmed cell death protein 1.

Figure 10 Analysis of GDSC database identifies novel candidate drugs targeting the bladder carcinoma immune genes signature. IC₅₀, half maximum inhibitory concentration; GDSC, Genomics of Drug Sensitivity in Cancer.

Discussion

The most prevalent cancer in the genitourinary system is bladder cancer, sometimes referred to as urothelial carcinoma (24). Bladder cancer is the ninth most common cancer worldwide and ranks sixth in prevalence in the United States (25). The treatment strategy for bladder cancer is determined by its pathological severity during transurethral resection of the bladder tumor (TURBT) and the ensuing staging, which considers tumor lymph node metastases (26). For patients with pure carcinoma in situ of the bladder, those diagnosed at an older age have an increased risk of recurrence and progression, and elderly patients may not respond to Bacillus Calmette-Guérin (BCG) therapy (27). For high-grade and high-risk non-muscle invasive bladder cancer (NMIBC), intravesical immunotherapy with BCG is considered the gold standard treatment to reduce the risk of recurrence and progression following the initial TURBT (28). For patients with high-grade NMIBC undergoing intravesical BCG therapy, the study by Ferro et al. indicated that the predictive factors for response to BCG in high-grade NMIBC patients include multifocality, lymphovascular invasion, and high-grade on repeat TURBT (re-TURBT) (29). The current standard treatment for bladder cancer patients with stage IV or metastatic disease is cisplatin-based chemotherapy, and cisplatin-based combination chemotherapies have grown to be crucial adjuvant treatments (30). Nonetheless, cisplatin-based therapy often leads to the development of chemoresistance, leading to treatment failures and limiting its effectiveness in treating bladder cancer (31). Shi et al. reported that the risk factors influencing the prognosis of patients with metastatic bladder cancer include age, chemotherapy, histological type, bone metastasis, lung metastasis, and liver metastasis (32). Immunotherapy can only be used as a first-line treatment for patients who are not eligible for cisplatin-based therapies (33). ICIs typically have a low rate of success in bladder cancer management. Therefore, it is crucial to identify patients who are most likely to benefit from ICI treatment to avoid unfavorable treatment outcomes. It is now understood that several immune environment subtypes affect tumor responses to ICIs (34). The importance of establishing the prognostic and immunotherapy-associated gene profiles of ICI in BLCA is demonstrated by this.

In this research, we developed an immune-related gene signature specific to bladder cancer, consisting of 11 DEIGs: CXCL12, LCN1, FGF10, SEMA6D, CORT, NTF3, PGR, PTGER3, USP2, JAM3, and C7, all of which have prognostic implications. It was determined by the drawing of risk curve, risk heatmap, ROC curve, and survival curve, and that the risk model did, in fact, have a good predictive effect a finding that was in line with that discovered when it was evaluated on the testing cohort, the external validation datasets GSE39281 and IMvigor210. The results of GSVA showed a good correlation between model risk scores and 21 KEGG pathways. According to a report by the C1 subtype, suppressive TME was implied by intermediate immune infiltration with M2 macrophages and active transforming growth factor (TGF)-β signaling pathway. The C2 subtype, meanwhile showed the strongest CD8⁺ signal and the widest variety of T cell receptors, suggesting a favorable immune-activated phenotype. In our approach, C1 and C2 subtypes account for a combined 87% of instances. The high-risk category had a greater prevalence of C2, whereas the low-risk group exhibited a higher percentage of C1. According to these findings, high-risk individuals had an active immune system whereas low-risk patients had an immunosuppressive immune system. Furthermore, patients in both the low- and high-risk groups exhibited statistically significant differences in immune function and immune-related survival. According to the correlation study results, most immune checkpoint genes displayed positive associations with CXCL12, LCN1, FGF10, SEMA6D, PGR, PTGER3, JAM3, and C7.

A nomogram was used to illustrate the AUC of the risk score in predicting 1-, 2-, and 3-year OS. The AUC values were 0.727, 0.772, and 0.765, respectively, indicating the robust predictive capability of the signature. According to the findings, the model is more precise than other models, including those based on age, gender, grade, the TNM (tumor-node-metastasis) staging system, tumor inflammation signature (TIS), and TIDE. The high-risk group displayed elevated T cell exclusion scores, which are typically indicative of a poor prognosis as they suggest that tumor cells can evade immune surveillance. This further validates the predictive ability of our model for patient responses to immunotherapy (35). In contrast, the low-risk group showed higher expression of CTLA4, potentially enhancing their susceptibility to ICIs. Consequently, patients in the low-risk group were more likely to benefit from immunotherapy compared to those in the high-risk group. Consistent with Wu et al.’s study, observations indicated that the high-risk prognostic group had heightened CTLA4 expression and demonstrated enhanced responsiveness to anti-PD-1 therapy (36). It is well-established that immune and stromal cells associated with tumors can have either tumor-inhibiting or tumor-promoting properties (37). The TME comprises immune cells that play a pivotal role in cancer development. In this study, the high-risk group exhibited elevated levels of immune score, stromal score, and ESTIMATE score, underscoring the BCIGS’s ability to effectively distinguish between different immune states. In terms of the prediction of drug sensitivity, 20 small molecular compounds were discovered, and the difference in the IC₅₀ value of the drug TW.37 between low- and high-risk patient groups was particularly pronounced. These substances may one day be used to make medications to treat patients with BLCA, but further research is required to confirm this. This study offers solid justification for the therapeutic application of BCIGS. In the future, we can analyze samples of BLCA patients for adjuvant therapy to give patients more individualized care.

The study is subject to a variety of limitations. As the model validation data were only gathered from public sources, multicenter clinical studies are still needed to corroborate them. There is a lack of immunotherapy data especially for BLCA patients. Therefore, it is imperative to validate the predictive effectiveness of BCIGS for immunotherapy by utilizing data from other cancer types. The accuracy and validity of our conclusions about the individual or combined effects of the 11 genes that make up BCIGS must be confirmed by further research.

Conclusions

We analyzed the role of immune-related genes in BLCA and developed BCIGS to predict prognostic significance, immunotherapy response, and drug sensitivity. BCIGS can effectively distinguish between the survival rate, immune subtype, and immune function of BLCA patients. Patients in the BCIGS high group exhibited high TME scores and T cell exclusion scores. In contrast, patients in the BCIGS low group displayed significantly higher CTLA4 expression levels. These findings hold significant potential for guiding the development of personalized immunotherapy and chemotherapy strategies, ensuring that BLCA patients receive more effective and precisely targeted treatments, ultimately improving their clinical outcomes.

Acknowledgments

Funding: This study was supported by National Natural Science Foundation of China (No. 82361148724) and Inner Mongolia Natural Science Foundation (No. 2023LHMS08072).

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1053/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1053/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1053/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Lobo N, Afferi L, Moschini M, et al. Epidemiology, Screening, and Prevention of Bladder Cancer. Eur Urol Oncol 2022;5:628-39. [Crossref] [PubMed]
2. Jubber I, Ong S, Bukavina L, et al. Epidemiology of Bladder Cancer in 2023: A Systematic Review of Risk Factors. Eur Urol 2023;84:176-90. [Crossref] [PubMed]
3. Teoh JY, Kamat AM, Black PC, et al. Recurrence mechanisms of non-muscle-invasive bladder cancer - a clinical perspective. Nat Rev Urol 2022;19:280-94. [Crossref] [PubMed]
4. Ulamec M, Murgić J, Novosel L, et al. New Insights into the Diagnosis, Molecular Taxonomy, and Treatment of Bladder Cancer. Acta Med Acad 2021;50:143-56. [Crossref] [PubMed]
5. Stühler V, Maas JM, Bochem J, et al. Molecular predictors of response to PD-1/PD-L1 inhibition in urothelial cancer. World J Urol 2019;37:1773-84. [Crossref] [PubMed]
6. Crispen PL, Kusmartsev S. Mechanisms of immune evasion in bladder cancer. Cancer Immunol Immunother 2020;69:3-14. [Crossref] [PubMed]
7. Chen X, Xu R, He D, et al. CD8(+) T effector and immune checkpoint signatures predict prognosis and responsiveness to immunotherapy in bladder cancer. Oncogene 2021;40:6223-34. [Crossref] [PubMed]
8. Genova C, Dellepiane C, Carrega P, et al. Therapeutic Implications of Tumor Microenvironment in Lung Cancer: Focus on Immune Checkpoint Blockade. Front Immunol 2022;12:799455. [Crossref] [PubMed]
9. Samstein RM, Lee CH, Shoushtari AN, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet 2019;51:202-6. [Crossref] [PubMed]
10. Lin A, Zhang J, Luo P. Crosstalk Between the MSI Status and Tumor Microenvironment in Colorectal Cancer. Front Immunol 2020;11:2039. [Crossref] [PubMed]
11. Yang C, Zhang Z, Tang X, et al. Pan-cancer analysis reveals homologous recombination deficiency score as a predictive marker for immunotherapy responders. Hum Cell 2022;35:199-213. [Crossref] [PubMed]
12. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
13. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
14. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008;9:559. [Crossref] [PubMed]
15. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw 2010;33:1-22. [Crossref] [PubMed]
16. Riester M, Werner L, Bellmunt J, et al. Integrative analysis of 1q23.3 copy-number gain in metastatic urothelial carcinoma. Clin Cancer Res 2014;20:1873-83. [Crossref] [PubMed]
17. Mariathasan S, Turley SJ, Nickles D, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 2018;554:544-8. [Crossref] [PubMed]
18. Amirouchene-Angelozzi N, Nemati F, Gentien D, et al. Establishment of novel cell lines recapitulating the genetic landscape of uveal melanoma and preclinical validation of mTOR as a therapeutic target. Mol Oncol 2014;8:1508-20. [Crossref] [PubMed]
19. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
20. Gentles AJ, Newman AM, Liu CL, et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat Med 2015;21:938-45. [Crossref] [PubMed]
21. Thorsson V, Gibbs DL, Brown SD, et al. The Immune Landscape of Cancer. Immunity 2019;51:411-2. [Crossref] [PubMed]
22. Bindea G, Mlecnik B, Tosolini M, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 2013;39:782-95. [Crossref] [PubMed]
23. Geeleher P, Cox NJ, Huang RS. Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol 2014;15:R47. [Crossref] [PubMed]
24. Xu Y, Luo C, Wang J, et al. Application of nanotechnology in the diagnosis and treatment of bladder cancer. J Nanobiotechnology 2021;19:393. [Crossref] [PubMed]
25. Compérat E, Amin MB, Cathomas R, et al. Current best practice for bladder cancer: a narrative review of diagnostics and treatments. Lancet 2022;400:1712-21. [Crossref] [PubMed]
26. Ouzaid I, Panthier F, Hermieu JF, et al. Contemporary surgical and technical aspects of transurethral resection of bladder tumor. Transl Androl Urol 2019;8:21-4. [Crossref] [PubMed]
27. Ferro M, Chiujdea S, Musi G, et al. Impact of Age on Outcomes of Patients With Pure Carcinoma In Situ of the Bladder: Multi-Institutional Cohort Analysis. Clin Genitourin Cancer 2022;20:e166-72. [Crossref] [PubMed]
28. Lopez-Beltran A, Cookson MS, Guercio BJ, et al. Advances in diagnosis and treatment of bladder cancer. BMJ 2024;384:e076743. [Crossref] [PubMed]
29. Ferro M, Barone B, Crocetto F, et al. Predictive clinico-pathological factors to identify BCG, unresponsive patients, after re-resection for T1 high grade non-muscle invasive bladder cancer. Urol Oncol 2022;40:490.e13-20. [Crossref] [PubMed]
30. Teply BA, Kim JJ. Systemic therapy for bladder cancer - a medical oncologist's perspective. J Solid Tumors 2014;4:25-35. [Crossref] [PubMed]
31. Li F, Zheng Z, Chen W, et al. Regulation of cisplatin resistance in bladder cancer by epigenetic mechanisms. Drug Resist Updat 2023;68:100938. [Crossref] [PubMed]
32. Shi S, Peng G, Luo L, et al. Predictive nomograms for risk and prognostic factors in metastatic bladder cancer: a population-based study. Transl Cancer Res 2023;12:3284-302. [Crossref] [PubMed]
33. Konala VM, Adapa S, Aronow WS. Immunotherapy in Bladder Cancer. Am J Ther 2022;29:e334-7. [Crossref] [PubMed]
34. Teng MW, Ngiow SF, Ribas A, et al. Classifying Cancers Based on T-cell Infiltration and PD-L1. Cancer Res 2015;75:2139-45. [Crossref] [PubMed]
35. Wang Z, Moresco P, Yan R, et al. Carcinomas assemble a filamentous CXCL12-keratin-19 coating that suppresses T cell-mediated immune attack. Proc Natl Acad Sci U S A 2022;119:e2119463119. [Crossref] [PubMed]
36. Wu X, Lv D, Cai C, et al. A TP53-Associated Immune Prognostic Signature for the Prediction of Overall Survival and Therapeutic Responses in Muscle-Invasive Bladder Cancer. Front Immunol 2020;11:590618. [Crossref] [PubMed]
37. Lei X, Lei Y, Li JK, et al. Immune cells within the tumor microenvironment: Biological functions and roles in cancer immunotherapy. Cancer Lett 2020;470:126-33. [Crossref] [PubMed]