Predicting immunotherapy response of advanced bladder cancer through a meta-analysis of six independent cohorts

Abstract

Advanced bladder cancer patients show very variable responses to immune checkpoint inhibitors (ICIs) and effective strategies to predict response are still lacking. Here we integrate mutation and gene expression data from 707 advanced bladder cancer patients treated with anti-PD-1/anti-PD-L1 to build highly accurate predictive models. We find that, in addition to tumor mutational burden (TMB), enrichment in the APOBEC mutational signature, and the abundance of pro-inflammatory macrophages, are major factors associated with the response. Paradoxically, patients with high immune infiltration do not show an overall better response. We show that this can be explained by the activation of immune suppressive mechanisms in a large portion of these patients. In the case of non-immune-infiltrated cancer subtypes, we uncover specific variables likely to be involved in the response. Our findings provide information for advancing precision medicine in patients with advanced bladder cancer treated with immunotherapy.

Fig. 1. Overview of the analyzed sequencing data of six metastatic bladder cancer datasets.

A Number of patients with different types of data. Sequencing data of 707 patients was downloaded. Patients with a RECIST classification of partial or complete response to anti-PD-L1/PD-1 were considered responders (R) and progressive disease were considered non-responders (NR). Stable disease and non-evaluable were not included for further analysis. B Of 466 patients with response to immunotherapy information, 348 patients have RNA-Seq data and a TMB estimate available. C Most frequently mutated BLCA genes mutated in ≥5% of the samples over all four datasets with WES data (N = 318 patients). The X axis represents the patients from the different cohorts, the corresponding TMB is indicated. Mutations are classified in different types as indicated; multi hit describes cases where the gene carries more than one mutation in the patient. R responders, NR non-responders, Mut panel DNAmutation panel.

Fig. 2. Relationship between somatic mutations and the response to ICI.

A Responders have a higher tumor mutational burden (Z score TMB) than non-responders (median Z score TMB R = 0.18, NR = −0.45). TMB is calculated as the number of non-synonymous mutations per 50 Mb. The differences were highly significant (two-sided, two-sample Wilcoxon test, p value = 2.1e-13). B Responders show a higher mean TMB than non-responders. This finding is consistent over all six datasets, and four of the six datasets display significant differences (two-sided, two-sample Wilcoxon test). C The difference between treatment response groups remains when separating the TMB into clonal and subclonal by a cancer cell fraction (CCF) cutoff of 0.9 (two-sided, two-sample Wilcoxon test). Median clonal TMB R: 2.15 mut/50 Mb, NR: 1.1 mut/50 Mb; median subclonal TMB R: 1.19 mut/50 Mb, NR: 0.8 mut/50 Mb. D Responders are significantly enriched in APOBEC-induced mutations compared to non-responders (two-sided, two-sample Wilcoxon test, p value = 9e-04; median APOBEC-enrichment score R: 3.45, NR: 2.85). E Spearman correlation between different DNA-derived variables. F Number of non-synonymous mutations by type and association of the different mutation types with response. Missense mutations, nonstop mutations, and putative neoantigens were found to be significantly associated with the response (two-sided, two-sample Wilcoxon test). Medians for number of missense mutations R: 202 muts, NR: 88 muts; nonstop mutations R: 0 muts, NR: 0 muts, frameshift insertions/deletions R: 3 muts, NR: 3 muts; in-frame insertions/deletions R: 0 muts, NR: 0 muts; putative neoantigens R: 22.5 neoantigens, NR: 12 neoantigens. Neoantigens were predicted from missense mutations applying a threshold of 500 nM IC₅₀ binding affinity in NetMHCpan 4.0. N = 234 patients. R responders (yellow), NR non-responders (turquoise), TMB tumor mutational burden, INDELs insertions and deletions. All p values are indicated in the according plots. Source data are provided as a Source Data file.

Fig. 3. Dissecting the effects of different immune-related variables.

A Normalized expression values for selected genes described in the context of Immune activity or suppression. All comparisons were significant after multiple testing adjustment, except for the genes TGFB1 and CCND1 (two-sided, two-sample Wilcoxon test). B Expression of the immune checkpoint molecules PD-L1 (gene CD274) and PD-1 (gene PDCD1) is significantly higher in responders (two-sided, two-sample Wilcoxon test). C Deconvolution analysis using CIBERSORT shows higher immune cell abundance in responders (two-sided, two-sample Wilcoxon test). D Six signatures of tumor antigen presentation and immune response are enriched among responders, while immune suppression signature scores are higher in non-responders. All comparisons were significant after multiple testing adjustment, except T cells inflamed and Stroma/EMT (two-sided, two-sample Wilcoxon test). E Gene signatures combining the expression of HLA-I and HLA-II types are higher in responders than in non-responders (two-sided, two-sample Wilcoxon test). F Responders show higher expression of a signature related to tumor-specific long non-coding RNAs (lncRNA) than non-responders (two-sided, two-sample Wilcoxon test). G Correlation matrix including TMB and RNA-Seq variables. P values obtained by two-sided, two-sample Wilcoxon test. N = 420 patients. R responders (yellow), NR non-responders (turquoise), APM antigen-presenting machinery, lncRNA long non-coding RNA. All p values are indicated in the according plots. Source data are provided as a Source Data file.

Fig. 4. The TCGA-subtypes luminal and neuronal have higher TMB, while basal-squamous and luminal-infiltrated show higher immune activity.

A Proportion of responders (yellow) and non-responders (turquoise) in each of the five TCGA subtypes. TMB was weakly correlated with the immune activation signatures. The neuronal subtype is the only one with a significant excess of responders over non-responders (two-sided Fisher’s exact test, p value = 0.014, indicated with *). B Basal-squamous and luminal-infiltrated have high values of immune cell abundance (absolute score obtained from CIBERSORT), while luminal, luminal-papillary and neuronal show low immune invasion. Immune cell abundance of CD4 T cells is similarly high in basal-squamous and luminal-infiltrated subtypes (CIBERSORT). Luminal-infiltrated shows the highest mean of CD4 memory resting T cells (CIBERSORT) (Kruskal–Wallis test). C Basal-squamous and luminal-infiltrated have the highest expression of immune biomarkers, both immune activating and suppressive markers (Kruskal–Wallis test). D. HLA expression is also highest in luminal-infiltrated and basal-squamous (Kruskal–Wallis test). E TMB is highest in the luminal and neuronal subtypes (Kruskal–Wallis test). F The tumor-infiltrated subtypes basal-squamous and luminal-infiltrated show the highest values of TGF-β gene expression (Kruskal–Wallis test). G Basal-squamous is enriched in antigen-presenting machinery (APM) compared to other subtypes (Kruskal–Wallis test). N = 420 patients. Ba/Sq: Basal-squamous (yellow), LumInf Luminal-infiltrated (orange), LumP Luminal-papillary (green), Lum Luminal (dark blue), NE Neuronal (light blue). All p values are indicated in the according plots. Source data are provided as a Source Data file.

Fig. 5. Influence of different biomarkers in the immune-infiltrated and non-immune-infiltrated subtypes.

A Overall, TMB is higher in the non-immune infiltrated group. In both groups, responders have significantly higher TMB (two-sided, two-sample Wilcoxon test). B In both groups, responders have higher enrichment scores of the APOBEC mutational signature (two-sided, two-sample Wilcoxon test). C Responders have higher IFN-γ expression values, and lower TGF-β expression values than non-responders for both subgroups (two-sided, two-sample Wilcoxon test). D CIBERSORT score is higher in responders of the infiltrated groups, but not the non-infiltrated group. Differences between responders and non-responders are significant for CD8+ T cells and macrophages M1 (two-sided, two-sample Wilcoxon test). E Only in the infiltrated group, responders have higher expression of the lncRNA signature. No difference was observed between response groups in the non-immune-infiltrative group (two-sided, two-sample Wilcoxon test). F Relationship between CD8 T cell and TGF-β gene expression. Immune-infiltrated samples tend to have high CD8+ T cell abundance, and in many cases also high levels of the TGF-β signature. The dashed line marks the optimal cutoff for each gene signature obtained by ROC (TGF-β: 0.0163; CD8 t effector cells: 0.246), responder are marked with circles, non-responder with filled circles. N = 420 patients. R: responders, NR: non-responders. All p values are indicated in the according plots. Source data are provided as a Source Data file.

Fig. 6. Predictive models of the response to immunotherapy.

A ROC curves for the complete model, TMB and TMB + RNA (TMB and RNA-seq derived variables), with AUC and number of samples (N). FPR: false positive rate; TPR: true positive rate. The average of 1000 runs is shown. B Feature positively or negatively associated with response. The length of the bar represents feature importance from random forest. Color reflects association with response taken from the previous manuscript sections (green: positively associated, red: negatively associated). C ROC curve for testing the TMB + RNA model in the validation cohort JAVELIN Bladder 100 trial. AUC and N values are indicated. D ROC curve of the complete model and TMB + RNA for the subset immune-infiltrated tumors, with AUC and N values. E ROC curve of the complete model and TMB + RNA for the subset of non-immune-infiltrated tumors, with AUC and N values. F ROC curves for TMB + RNA model when removing one of the cohorts at a time.