Molecular correlates of cisplatin-based chemotherapy response in muscle invasive bladder cancer by integrated multi-omics analysis

Abstract

Overtreatment with cisplatin-based chemotherapy is a major issue in the management of muscle-invasive bladder cancer (MIBC), and currently none of the reported biomarkers for predicting response have been implemented in the clinic. Here we perform a comprehensive multi-omics analysis (genomics, transcriptomics, epigenomics and proteomics) of 300 MIBC patients treated with chemotherapy (neoadjuvant or first-line) to identify molecular changes associated with treatment response. DNA-based associations with response converge on genomic instability driven by a high number of chromosomal alterations, indels, signature 5 mutations and/or BRCA2 mutations. Expression data identifies the basal/squamous gene expression subtype to be associated with poor response. Immune cell infiltration and high PD-1 protein expression are associated with treatment response. Through integration of genomic and transcriptomic data, we demonstrate patient stratification to groups of low and high likelihood of cisplatin-based response. This could pave the way for future patient selection following validation in prospective clinical trials.

Conflict of interest statement

L.D. has sponsored research agreements with C2i genomics, AstraZeneca, Natera and Ferring, and has an advisory/consulting role at Ferring. J.B.J. is proctor for Intuitive Surgery, member of advisory board for Olympus Europe, Cephaid and Ferring, and has sponsored research agreements with Medac, Photocure ASA, Cephaid and Ferring. The following authors declare no competing interests: A.T., E.C., P.L., I.N., F.F.P., S.V.L., K.B., T.L.H.O., M.K., J.S.P., M.A., and T.S.

Figures

Fig. 1. Overview of the genomic alterations correlated to chemotherapy response.

a Oncoplot showing the significantly mutated genes in 165 tumors annotated by exome coverage, mutation load stratified by impact (as defined by SnpEff) and mutational signature deconvolution (top panels) and by clinical response, neoantigen load, number of damaging mutations in DDR genes, percentage of genome in allelic imbalance, expression subtypes, regulon cluster, RNA immune score, hypermethylation cluster and immune phenotype (bottom panel). b–e Violin plots showing the total number of SNVs, indels, neoantigens or the percentage of the genome in allelic imbalance compared to chemotherapy response. f Presence of damaging DDR mutations compared to chemotherapy response. g–i Violin plots showing the total number of SNVs, indels or the percentage of the genome in allelic imbalance compared to the two main mutational subtypes. j Presence of DDR mutations compared to the two main mutational subtypes. k Mutational subtypes compared to chemotherapy response. All p-values were calculated using a Wilcoxon rank-sum test for continuous variables and a Fisher’s exact test for categorical variables. Source data are provided as a Source data file.

Fig. 2. BRCA2 mutated tumors are associated with mutations in an SBS5 context and response to chemotherapy.

a Number of mutations in an SBS5 context in relation to response to chemotherapy. b Mutation status of ERCC2 compared to chemotherapy response. c, d Volcano plots showing the difference between the median number of mutations for mutated tumors and the median number of mutations for wild-type tumors for all genes mutated in more than 5% of the patient cohort (only mutations with moderate- or high protein impact are considered). The red dashed lines indicate significance levels at p = 0.05. c Represents the number of mutations in an SBS5 context and d represents the number of mutations in an SBS2+13 (APOBEC) context. P-values were calculated using a permutation test (n = 100,000) that controls for mutation burden per sample and gene. e Mutation status of BRCA2 compared to chemotherapy response. f Kaplan–Meier survival analysis depicting the probability of overall survival stratified by BRCA2 mutation status. g–j Based on TCGA data. g Number of mutations in an SBS5 context for ERCC2 and BRCA2 wild-type and mutant samples in TCGA data. h Number of mutations in an SBS2+13 (APOBEC) context for ERCC2 and BRCA2 wild-type and mutant samples in TCGA data. i, j Mutation status of ERCC2 and BRCA2 compared to response to chemotherapy in TCGA data. P-values were calculated using a Fisher’s exact test for categorical variables, a Wilcoxon rank-sum test for continuous variables and a log-rank test for comparing survival curves. For all boxplots, the center line represents the median, box hinges represent first and third quartiles, whiskers represent ±1.5 × interquartile range (IQR) and points represent outliers. Source data are provided as a Source data file.

Fig. 3. Delineation of metastatic evolution before chemotherapy.

All six patients were treated with cisplatin-based chemotherapy upon detection of metastatic disease. a Location of metastasis and the time to recurrence (mo = months). P = primary tumor, M = metastasis. The images were created using Biorender.com. b Clonal relationships between primary tumor samples and metastatic lesions depicted by phylogenetic trees. Trunk/branch lengths correspond to the number of SNVs. Mutations in genes involved in DDR, frequently mutated in TCGA or identified as drivers in BC (IntOGen) are indicated. Green = trunk, yellow = branch, blue = primary tumors, pink = metastatic lesions. c Variant allele frequencies for mutations identified in either of the available samples per patient. Identified mutations were subjected to read-counting in processed bam files to enable identification of mutations called in one sample, and present, but not called in another sample. The required read depth for identifying a given mutation was calculated for every position based on the lowest observed allele frequency. Only positions with sufficient read depth in all investigated samples were included. d Box plots depicting the observed allele frequencies for trunk and branch mutations. Asterisks indicate p-values below 2.2e−16. e Relative contribution of mutational signatures in the trunks (left circles) and branches (right circles). P-values were calculated using a Wilcoxon rank-sum test. For all boxplots, the center line represents the median, box hinges represent first and third quartiles, whiskers represent ±1.5 × interquartile range (IQR) and points represent outliers. Source data are provided as a Source data file.

Fig. 4. Gene expression characteristics and relation to chemotherapy response.

a Visualisation of four identified consensus gene expression subtypes by selected gene sets. b Gene expression subtypes compared to chemotherapy response. c Kaplan–Meier survival analysis showing the probability of overall survival for patients with and without Ba/Sq gene expression subtype. d Number of mutations for gene expression subtypes. e Heatmap showing relative expression values for identified regulons and deconvoluted immune cells. f Regulon clusters compared to gene expression subtypes. g Regulon clusters in relation to response to chemotherapy. h Immune score across the identified gene expression subtypes. i Immune score compared to response to chemotherapy. j Summarized expression of the antigen-presenting machinery for immune hot (above median immune score) samples stratified by RECIST 1.1 response values. Only one tumor was classified as NE-like and was therefore omitted from this figure. Missing data is depicted in gray. P-values were calculated using a Fisher’s exact test for categorical variables, a Wilcoxon rank-sum test for continuous variables and a log-rank test for comparing survival curves. For all boxplots, the center line represents the median, box hinges represent first and third quartiles, whiskers represent ±1.5 × interquartile range (IQR) and points represent outliers. Source data are provided as a Source data file.

Fig. 5. DNA methylation subtypes based on hypermethylated cancer-specific CpG sites.

a Clustering of samples based on hypermethylation events (n = 5000). Heatmap shows beta values and the right panel shows normal bladder and leukocyte beta values for comparison. b Methylation clusters compared to gene expression subtypes. c Gene set scores calculated using xCell and stratified by methylation clusters (HMC2: n = 17; HMC1: n = 31; HMC3: n = 16). d RECIST 1.1 response measurements stratified by methylation clusters. MEscore = Microenvironment score. P-values were calculated using a Wilcoxon rank-sum test. For all boxplots, the center line represents the median, box hinges represent first and third quartiles, whiskers represent ±1.5 × interquartile range (IQR) and points represent outliers. Source data are provided as a Source data file.

Fig. 6. Immune tumor microenvironment analysis by spatial proteomics.

Immunostaining performed on bladder cancer tissue microarrays from 184 patients. All protein measurements were performed once for each distinct sample. a Staining results shown for multiplex immunofluorescence (mIF) panel 1 and PD-1 with corresponding image analysis application (APP) of four tumors representing each immune subtype and high PD-1 expression, respectively. Red dashed lines divide tissue into intratumoral and peritumoral regions of interest. Scale bar: 20 µm. b Spatial organisation of immune cells and immune evasion mechanisms stratified by immune subtypes. Heatmap shows z-scores and the barplot the mean + SD immune scores (IS), cell percentages, or H-score [1 × (% cells low intensity) + 2 × (% cells moderate intensity) + 3 × (% cells high intensity)] for MHC-expression on carcinoma cells. Asterisks denote the barplot representing the H-score. Points represent the corresponding data points. c Immune subtypes compared to chemotherapy response. d Kaplan–Meier survival curves showing overall survival (OS) stratified by immune subtypes. e–g Intratumoral and peritumoral fractions of immune evasion mechanisms (PD-1, PD-L1, MHC) compared to chemotherapy response. h Intratumoral combined PD-1/PD-L1 expression stratified by chemotherapy response. PD-1 high/PD-L1 low was compared to PD-1 high/PD-L1 high, PD-1 low/ PD-L1 low and PD-1 low/PD-L1 high combined. i–k Relationship between immune evasion mechanisms and immune subtypes. Statistical significance was assessed using a chi-square test for categorical variables, a Wilcoxon rank-sum test for continuous variables and a log-rank test for comparing survival curves. For all boxplots, the center line represents the median, box hinges represent first and third quartiles, whiskers represent ±1.5 × interquartile range (IQR) and points represent outliers. Source data are provided as a Source data file.

Fig. 7. Integrative analysis.

a Chemotherapy response stratified by patient groups of high genomic instability (HGI) or low genomic instability (LGI). b Integration of genomic, transcriptomic and proteomic data for all patients with genomic data available and a gene expression subtype assigned (n = 121), stratified by high and low genomic instability. c Chemotherapy response for patients with Ba/Sq gene expression subtype and patients with non-Ba/Sq gene expression subtype (i.e., LumP, LumU, NE-like or Stroma-rich). d Overview of odds ratios (OR) calculated for molecular and clinical features for all patients (black), NAC (green) and first-line treated patients (yellow). Continuous variables were dichotomized based on the median and high vs. low is presented. Dots indicate odds ratios and horizontal lines show 95% confidence intervals (CI). Open lines indicate that the full range of the 95% CI is not shown. Asterisks denote p-values below 0.05. Currency signs denote p-values below 0.05 in a logistic regression model using continuous data. e Likelihood of cisplatin-based chemotherapy response stratified by genomic instability and gene expression subtypes. P-values were calculated using a chi-square test. Source data are provided as a Source data file.