Evolution and heterogeneity of lethal metastatic bladder cancer subtypes

Metastatic bladder cancer (BLCA) is a common and highly lethal disease that remains a major clinical challenge owing to its marked histological and molecular heterogeneity^5,6,7,8. Approximately one in four patients with metastatic BLCA have tumours with subtype histological features such as squamous, plasmacytoid, neuroendocrine (NE) or sarcomatoid differentiation^1,2. These histological subtypes often have distinct disease trajectories compared with conventional urothelial carcinoma (UC)^1,2,3,9. However, it is not well understood how genomic, transcriptomic and microenvironmental heterogeneity evolve in these tumours and shape their distinct natural histories and responses to therapy^{10,11,12,13,14,15}. Charting the evolutionary history of metastatic tumours in patients can uncover the emergence of treatment-resistant clones and identify optimal points for intervention. However, so far, a comprehensive characterization of tumour evolution and molecular heterogeneity in histological subtypes has been limited, particularly in the metastatic setting.

The University of Washington and Fred Hutchinson Cancer Center (UW–Fred Hutch) metastatic BLCA rapid autopsy programme was designed to enrich for histological subtypes of carcinomas arising in the urothelium. This resource was developed to investigate intra-patient and inter-patient tumoural heterogeneity and to retrace the evolutionary trajectories of metastatic BLCA^16,17,18. Here using whole-genome sequencing (WGS) of primary and multiple metastatic tumours from each patient, we reconstruct tumour clonal phylogeny and metastatic seeding, which reveal insights into the nature of metastatic spread and drivers of aggressive histological subtypes and therapy resistance. We analyse bulk RNA sequencing (RNA-seq) and single-nucleus RNA sequencing (snRNA-seq) data to uncover transcriptional heterogeneity and molecular hallmarks that underlie histological subtypes, their tumour microenvironments (TMEs) and disease progression. Finally, we evaluate the potential of post-mortem cell-free DNA (cfDNA) to capture the genomic and transcriptional heterogeneity of metastases and to enable a minimally invasive approach to identify driver genomic alterations and distinguish aggressive metastatic BLCA subtypes. Our detailed analyses of the clonal dynamics of metastatic tumour spread and heterogeneity offer new insights into the temporal evolution of lethal metastatic BLCA.

The UW–Fred Hutch BLCA rapid autopsy programme was established in 2015 to collect tumours and benign tissues from patients with terminal disease¹⁶. Tissue samples (20 benign, 8 primary tumours and 80 metastases) were collected via rapid autopsy (median 3.7 h following death) from 20 patients with carcinomas of the urothelium (Fig. 1a, Extended Data Fig. 1a and Supplementary Table 1). Primary tumours were obtained from diagnostic transurethral resection of bladder tumour (TURBT; n = 5), surgery (n = 11) or autopsy (n = 8); 4 patients had both TURBT and autopsy primary tumours. The tumours represented a spectrum of histological subtypes⁴, including 50 with conventional UC, and an enrichment of subtypes, including 24 plasmacytoid UC (PUC), 19 UC with squamous differentiation (UCSD), 8 NE and 3 UC with sarcomatoid subtype (UC-Sarc). PUC and NE tumours were predominantly variant histology (mean >95%), whereas UCSD and UC-Sarc tumours had variable histological purity. Patients with PUC tumours had the shortest survival (Extended Data Fig. 1b), whereas liver metastasis presence was not associated with survival after adjusting for histology (Extended Data Fig. 1c).

WGS (mean 64.7× coverage) or whole-exome sequencing (WES; mean target coverage 43.2×) were performed on 93 tumours (Supplementary Table 1). Bulk RNA-seq (n = 83 tumours), snRNA-seq (n = 15 tumours) and cfDNA sequencing (n = 17 targeted, n = 8 WGS) were also conducted. Detailed clinical histories for each patient enabled integration of multi-omic data with survival, treatment response, imaging and pathology data (Extended Data Fig. 1a and Supplementary Tables 1 and 2).

To study the evolutionary history between primary and metastatic tumours, we reconstructed phylogenetic relationships, clonality and metastatic seeding patterns on the basis of somatic mutation analyses (Methods). For each patient, mutations were clustered by their cellular prevalence and used to infer clonal relationships between tumours. Mutation clusters were assigned to three categories: founder (highest cellular prevalence clone present in all tumours); shared (non-founder clones present in multiple tumours); and private (clones exclusive to one tumour) (Supplementary Table 3). Using these mutational phylogenetic trees, we predicted the sequence of metastatic clonal migration for 16 patients with high-quality primary tumours (Fig. 1b,c and Extended Data Fig. 1d). In most patients (9 out of 16; 56%), a single early metastatic tumour seeded all other metastases, whereas the primary tumour was the dominant seeder in 5 (31%) patients. Thirteen (81%) patients had at least one metastasis-to-metastasis seeding event. Complex seeding was observed in patient 19-004, in whom three tumours displayed seeding ability, and in patients 15-109 and 18-065, in whom a single metastatic niche was seeded by multiple tumours.

Computed tomography (CT) scans performed before autopsy were used to corroborate the ordering of inferred metastasis-to-metastasis seeding. We evaluated scans for all 20 patients and identified 5 for whom metastatic timing could be determined (Methods, Supplementary Table 2 and Supplementary Fig. 2). For patients 17-030, 19-004 and 19-037, radiographic imaging detected a single metastasis in the patient’s first staging scan that confirmed the first seeded metastasis (Fig. 1c). In patients 16-070 and 17-020, the initial scan detected two early metastatic lesions, and computational seeding analysis further resolved clonal constraints to identify the initial seeded metastasis (Extended Data Fig. 1e). Metastases of the abdominal wall (patient 17-030), liver (patient 19-004) and perinephritic site (patient 19-037) that were absent in the last scans probably emerged later in disease.

We observed polyclonal seeding in all 16 patients, including 12 (75%) with polyphyletic migration (events involving multiple distinct clonal lineages) (Supplementary Table 3). The polyclonality of seeding, defined for each patient as the average number of clones per metastatic seeding event, trended with increased mortality risk (hazard ratio (HR) = 6.8, P = 0.09; Fig. 2a and Extended Data Fig. 2a). This polyclonality metric improved survival risk stratification (log-rank P = 0.0007; Fig. 2b) beyond histology, age and sex alone (log-rank P = 0.008), whereas the purity of variant histology did not contribute to survival risk (log-rank P = 0.06; Extended Data Fig. 2b,c). Metastases seeded an increased proportion of polyclonal polyphyletic events compared with primary tumours (0.36 versus 0.65, P = 0.038; Fig. 2c). Tumours with non-UC histology also had higher proportions of polyphyletic events compared with UC tumours (P ≤ 0.04; Fig. 2d). We did not observe differences in seeding patterns between patients treated with chemotherapy or immune checkpoint inhibitors (ICIs) (Extended Data Fig. 2d). However, tumours arising after ICI treatment (identified by CT scan) were seeded via more migrating clones (P = 0.11) and polyclonal polyphyletic events (P = 0.029) than pre-ICI metastases (Fig. 2e and Extended Data Fig. 2e,f). These ICI-persistent clones were not site-specific but exhibited increased expansion in the post-ICI metastases (P = 0.022; Fig. 2f and Extended Data Fig. 2g). Together, these findings indicate that the migration of more distinct clones is linked to non-UC metastatic BLCAs and ICI-related tumour evolution.

Among 128 known BLCA-associated genes⁵, we observed the following number of alterations: 143 pathogenic mutations (single nucleotide variants (SNVs) and insertions and deletions (indels)), 2–16 per patient; 165 copy number alterations (CNAs), 3–18 per patient; and 40 transections by structural variant (SV) breakpoints, 1–7 per patient (Methods, Extended Data Fig. 3 and Supplementary Table 4). Pathogenic events were acquired earlier than non-pathogenic events across all alteration types (Fig. 2g). Pathogenic mutations had higher enrichment in founder clones (median 68.3%) than either CNA or SV events (median 33.3% and 6.2%, respectively; P < 0.040). This finding indicates that larger genomic alterations may be acquired later than mutations in BLCA, as seen in other cancers¹⁹. Moreover, positive selection (ratio between non-synonymous and synonymous variation (dN/dS) > 1) for truncating mutations across all clones was observed. This result suggests that there is continued selection for tumour-suppressor inactivation throughout tumour evolution (Supplementary Fig. 3a).

All 20 patients had late subclonal pathogenic alterations that contributed to intra-patient heterogeneity, including mutations (n = 17 patients), CNAs (n = 18) and transecting SVs (n = 10). Private mutations that are druggable were observed in FGFR3 (R248C, patient 19-022), PIK3CA (E545K, patient 15-108) and KMT2D (frameshift deletion, patient 19-007) (Supplementary Fig. 3b–d). CNAs involving PIK3CB and NECTIN4 were frequently non-founder and subclonal events, a finding consistent with previous observations of enrichment in metastases¹⁸ (Extended Data Fig. 4a–e). Conversely, CNAs affecting TP53 and NOTCH1 were primarily founder and clonal, a result that highlights the variable evolutionary timing of BLCA drivers. Founder alterations of TP53 and TERT were also enriched in patients with metastasis-to-metastasis seeding, whereas founder CDKN2A alterations were enriched in primary-to-metastasis seeding (Extended Data Fig. 4f). These results are consistent with reports that CDKN2A alterations are mutually exclusive of TP53 loss in muscle-invasive bladder cancer, enriched in metastatic (49%) versus localized (22%) BLCA and associated with high metastatic capacity^5,20,21.

For the patients with primary tumour samples collected at both TURBT and autopsy, we observed multiple distinct evolutionary patterns. In patient 17-020, the primary tumour samples were more similar to each other than to the metastases, which indicated divergent metastatic evolution (Fig. 1b and Extended Data Figs. 3 and 5a). Conversely, in patients 19-007 and 17-026, the metastases were more genomically similar to the TURBT sample than the primary tumour sample at autopsy (Extended Data Fig. 5b). For example, chromosome 6q loss of heterozygosity was exclusive to the TURBT and metastatic samples from patient 17-026. This result suggests that the primary tumour sampled at autopsy in patient 17-026 probably evolved from a minor subclone without chromosome 6q loss of heterozygosity early in disease (Extended Data Fig. 5c–f). Despite these various evolutionary trajectories, the source and timing of primary tumour collection did not lead to systematic differences in clonal composition between primary and metastatic tumours (Extended Data Fig. 5g,h).

Next, we identified 142 unique complex SVs²² (Fig. 2h and Supplementary Table 4), with only 20 (14%) classified as early founder events. However, the majority of SV-disrupted BLCA-associated genes were transected by founder complex SVs (23 out of 33; 69%) (Supplementary Fig. 4). We observed founder breakage–fusion–bridge cycle amplifications and non-cyclic SVs amplifying CCND1 and NECTIN4 (Supplementary Fig. 5). By contrast, in patient 19-007, ERBB2 was amplified via extrachromosomal DNA (ecDNA) in shared clones specific to the lymph node metastases and TURBT sample, which were confirmed to have higher HER2 protein expression than the autopsy primary tumour sample (Fig. 2i and Supplementary Fig. 6). Thus, metastatic BLCA progression is characterized by early clonal driver events, with continued evolution and heterogeneity through complex genomic structural alterations.

The burden, heterogeneity and timing of acquired genomic alterations were markedly different among histological subtypes. PUC tumours had the highest coding mutational burden (Extended Data Fig. 5i), whereas UCSD and UC tumours had moderately more CNAs and SVs, respectively (Extended Data Fig. 5j,k). Patients with PUC and, to some extent, NE tumours exhibited a significantly higher proportion of pathogenic alterations in founder clones, whereas those with UC or UCSD tumours had even distributions across early and late clones (Fig. 2j). This trend corresponded with frequent founder alterations of histology-associated drivers, including CDH1, TP53 and RB1 in PUC tumours, and TP53, PTEN and RB1 in NE tumours^23,24 (Extended Data Fig. 3). The early, stunted evolution of PUC tumours was also consistent with having smaller, less expanded subclones (Supplementary Fig. 7a–c).

The higher proportion of founder CNAs in PUC and NE tumours correlated with shorter survival (Fig. 2k and Supplementary Fig. 7d). By contrast, private alterations in UC and UCSD tumours increased with survival. This result suggests that pathogenic alterations acquired throughout tumour evolution may have contributed to progression later in disease. Collectively, these results support the conclusion that pathogenic founder alterations promote early aggressiveness in PUC and NE tumours.

We next investigated the evolution and timing of mutational processes by assessing clone-specific single-base mutational signatures (SBSs) (Methods). The most prominent mutational processes were clock-like ageing (SBS1 + SBS5), APOBEC (SBS2 + SBS13) and platinum chemotherapy (SBS31 + SBS35) (Extended Data Fig. 6a and Supplementary Table 5). Clock-like and APOBEC-mediated mutation processes were enriched in founder clones (P ≤ 0.039), whereas mutations induced by platinum chemotherapy and DNA repair deficiencies were increased in shared and private clones (P ≤ 0.042), as previously reported⁶ (Fig. 3a). Clock-like and APOBEC signatures did not differ significantly across histological subtypes (Extended Data Fig. 6b,c). However, PUC tumours displayed a significantly lower proportion and total burden of platinum-induced mutations in shared and private clones than UC or UCSD tumours (P = 0.009; Fig. 3b and Extended Data Fig. 6d), despite similar cisplatin exposure. PUC tumours also had comparable treatment-to-death intervals and an increased overall tumour mutational burden (Extended Data Fig. 6e–g).

The limited presence of the SBS31 + SBS35 mutational signature in cisplatin-treated PUC tumours was not explained by altered gene expression levels of several common resistance mechanisms²⁵ (Extended Data Fig. 6h). Although PUC tumours displayed decreased glutathione expression, this would be expected to increase cisplatin-induced mutations²⁵. However, expression of components of the Fanconi anaemia (FANC) pathway, which repair cisplatin-induced inter-strand crosslinks²⁶, trended higher in PUC than in UC or UCSD tumours (P = 0.071; Fig. 3c). Furthermore, increased FANC pathway expression was significantly associated with decreased SBS31 + SBS35 proportion (ρ = −0.7, P = 0.024; Fig. 3d). We observed similar correlations in independent cohorts of patients with metastatic BLCA treated with cisplatin⁵ (r = −0.31, P = 0.025; Fig. 3e) and patients with ovarian cancer²⁷ treated with platinum (r = −0.19, P = 0.039; Extended Data Fig. 6i and Supplementary Table 5). FANC expression also correlated with cell cycle progression scores (r = 0.93, P = 1.47 × 10⁻¹⁸), which suggests that its upregulation reflects the high proliferation of PUC tumours (Extended Data Fig. 6j).

Next, we investigated whether expression of FANC pathway genes is associated with cisplatin resistance. There are currently no models of PUC; however, publicly available data²⁸ have suggested that FANC expression is positively correlated with cisplatin resistance in UC-derived BLCA cell lines (Extended Data Fig. 6k). Therefore, we tested UC-origin BLCA cell lines with either FANC-high expression (SW780 and HT1376) or FANC-low expression (UBLC1 and BFTC905) (Extended Data Fig. 6l). First, we examined whether increased FANC gene expression was sufficient to increase pathway activity. After 24 h of 5 µM cisplatin treatment, the FANC-high lines displayed a significant increase in ubiquitinated FANCD2 (P = 0.009), a downstream indicator of FANC activity^29,30, whereas FANC-low lines did not (P = 0.699) (Fig. 3f). This difference in pathway activity translated to cisplatin resistance in FANC-high cell lines, with SW780 and HT1376 cells displaying increased cell confluence (mean half-maximum inhibitory concentration (IC₅₀) values of 13.2 versus 6.1 µM, P = 0.004) and decreased cell death (P = 0.029) after cisplatin treatment (Fig. 3g,h and Extended Data Fig. 6m,n). FANC-high cell lines also responded to mitomycin C (MMC), another crosslink-inducing chemotherapy, with increased ubiquitinated FANCD2 and resistance to cell death, although no difference was observed in cell confluence (Extended Data Fig. 6o–q). The similar effects between cisplatin and MMC support the idea that resistance stems from FANC repair of inter-strand crosslinks, whereas other resistance pathways may have contributed to variation in sensitivity.

We further validated the relationship between FANC expression and cisplatin response by knocking down or overexpressing FANCF, a key component of the FANC core complex known to modulate crosslinker sensitivity^26,31, in the highest (SW780) and lowest (UBLC1) FANC-expressing cell lines, respectively. Differential DNA repair activity in these cell lines was confirmed by assessing γH2AX foci (Extended Data Fig. 7a). FANC-low UBLC1 cells accumulated more DNA damage after cisplatin treatment, whereas SW780 cells plateaued, which indicated that there was increased DNA repair in FANC-high SW780 cells (Fig. 3i). After FANCF knockdown, SW780 cells displayed decreased ubiquitinated FANCD2 and increased cisplatin sensitivity (Fig. 3j–l and Extended Data Fig. 7b–e). Reciprocally, FANCF overexpression in FANC-low UBLC1 cells increased pathway activity and cisplatin resistance (Extended Data Fig. 7f–i). These data demonstrate that FANC pathway expression is highly linked to cisplatin responses in BLCA cells. This finding provides support for a mechanism in which increased FANC activity in PUC tumours may underlie the difficulty in treating this aggressive, often platinum-resistant subtype³².

We next investigated inter-tumoural and intra-tumoural transcriptional heterogeneity of metastatic BLCA histological subtypes. First, we classified tumours into established, prognostic consensus molecular subtypes⁸ (Methods). We observed concordance between histological and molecular subtypes from bulk RNA-seq, including >90% of UCSD and UC-Sarc tumours subtyped as basal/squamous (Ba/Sq) and 71% of NE tumours as NE-like (Fig. 4a and Supplementary Table 6). UC tumours were either luminal papillary (LumP) (32%) or Ba/Sq (51%). PUC tumours were distributed into four subtypes, with mainly luminal unstable (LumU) metastases (62%), consistent with their high tumour mutational burden, but stroma-rich primary tumours (Fig. 4b). PUC tumours were the only variant to exhibit this dichotomy. Despite most (13 out of 18) patients having uniform intra-patient molecular subtyping, unbiased analyses of transcriptional variation revealed notable intra-subtype heterogeneity in UC, PUC and UCSD histological subtypes (Fig. 4c and Extended Data Fig. 8a).

To investigate transcriptional heterogeneity at single-cell resolution, we sequenced and analysed RNA from 39,888 nuclei across 15 tumours in 9 patients (Supplementary Table 1 and Supplementary Fig. 8a,b). We identified distinct epithelial, immune and stromal cell compartments (Fig. 4d and Supplementary Fig. 8c–e), which clustered by histological subtypes and individual tumours (Extended Data Fig. 8b,c). The presence of CNAs confirmed that the 24,412 epithelial cells are tumour cells (Extended Data Fig. 8d). In this epithelial compartment, PUC samples had increased expression of cell cycle and DNA repair pathways, including high FANC (Extended Data Fig. 8e–h). UCSD tumours had increased expression of epithelial-to-mesenchymal transition (EMT) signalling, whereas immune-related pathways (for example, TGFβ, WNT and IFNα) were altered in both PUC and UCSD tumours (Extended Data Fig. 8g). These results were confirmed in the bulk tumour RNA-seq samples, thereby supporting these pathways as hallmarks of these histological subtypes (Extended Data Fig. 8i).

Tumour cells clustered according to both histological and molecular subtypes, and transitional ‘bridge’ samples may indicate the presence of mixed phenotypes (Fig. 4e,f). The PUC tumours comprising these bridges (samples 15-109G5 and 17-071I1; Extended Data Fig. 8f) were classified as LumP and lacked the PUC hallmarks of ERBB2 overexpression and CDH1 loss^24,33 (Extended Data Fig. 8j,k). Furthermore, four UC tumours were classified as Ba/Sq and clustered along with the UCSD samples, which indicated a potential transcriptional-state change before squamous differentiation emerged. These mixed cell states were evident in the heterogeneous molecular subtypes observed in UC and PUC tumours, including in the mixed UC–PUC histology specimen 17-071I1 (Fig. 4g). This increased intra-patient transcriptional heterogeneity (Methods) correlated with both increased survival from first metastasis (R² = 0.60, P  = 0.041; Fig. 4h) and elevated genomic heterogeneity (R² = 0.74, P = 0.014), while being unrelated to treatment (Extended Data Fig. 8l–o). These findings highlight that although there are transcriptional features associated with histological subtypes, significant heterogeneity persists among subtypes that correlates with survival.

We compared immune and fibroblast compositions in the TME between histological subtypes using snRNA-seq. PUC tumours had the highest proportion of B cells and T cells in lymph node metastases (median 29% B cells and 15% T cells in PUC lymph nodes versus 4% B cells and 1% T cells in UC lymph nodes) (Fig. 4i, Extended Data Fig. 9a and Supplementary Table 7), which was further supported by bulk RNA-seq deconvolution (P ≤ 0.037; Extended Data Fig. 9b and Supplementary Table 6). Re-clustering of the 6,763 immune cells further indicated that these PUC lymph node metastases were immune-inflamed, with an enrichment of B cells and cytotoxic CD8 T cells, whereas UC tumours contained more exhausted CD8 T cells (Fig. 4j,k and Extended Data Fig. 9c–f). UCSD tumours were enriched in fibroblasts and myeloid cells in both snRNA-seq (P ≤ 0.029 versus UC; Fig. 4i and Extended Data Fig. 9a) and bulk RNA-seq data (P ≤ 0.036; Extended Data Fig. 9b). This expanded myeloid population in UCSD tumours was mainly composed of MRC1⁺ tumour-associated macrophages (TAMs), including SPP1⁺ TAMs, which indicate an immunosuppressive environment³⁴ (Fig. 4k).

Changes in immune cell composition can be related to intercellular signalling. Analyses of intercellular ligand–target interactions identified TGFB1, encoding a cytokine that represses T cell function³⁵, as the ligand most related to the differences between PUC and UC–UCSD T cells (Methods). Both epithelial and fibroblast cells in PUC tumours had decreased TGFB1 expression, which led to low TGFB1 target activity in T cells (Fig. 4l). Re-clustering of the 5,039 fibroblasts also revealed differences in cancer-associated fibroblasts (CAFs) across histological subtypes (Fig. 4m and Extended Data Fig. 9g–i). UCSD tumours were enriched with matrix CAFs (78%, P < 0.03) and had increased expression of fibroblast activation protein (FAP), which can support SPP1⁺ TAMs and inhibit lymphocyte infiltration^36,37 (Extended Data Fig. 9j). This result may relate to the minimal efficacy of ICI therapies (median 1.38 months before progression) for patients with UCSD tumours (Extended Data Fig. 1a). By contrast, PUC tumours were enriched for inflammatory CAFs (18%, P < 0.036) and reticular-like CAFs (24%), which are known to support T cell recruitment in lymph nodes³⁸. Notably, patient 18-120, who had PUC and an inflamed TME, experienced an exceptional response (>12 months) to anti-PD1 therapy (Extended Data Figs. 1a and 9e). These results suggest that metastatic BLCA histological subtypes vary in immune activation, composition and signalling, with potential implications for immunotherapy response in end-stage disease.

cfDNA collected contemporaneously with autopsy tissue for 17 patients provided an ideal opportunity to assess the dynamics of circulating tumour DNA (ctDNA) shedding^39,40, which can affect the detection of clinically relevant states^41,42. We observed that the cfDNA tumour fraction (TFx; median 5.9%) was increased in patients with 3+ liver metastases (P = 0.0004; Methods and Extended Data Fig. 10a). No relationship was observed between the time from death to autopsy and TFx (ρ = 0.018, P = 0.94), although total cfDNA yield decreased with increasing post-mortem interval (ρ = −0.41, P = 0.10) (Extended Data Fig. 10b,c). To assess the utility of cfDNA sequencing breadth versus depth, we performed deep WGS (mean 73× fragment coverage) for samples with TFx values of ≥8% (n = 8) and ultradeep targeted panel (TP) capture of 397 cancer genes (mean 3,469× target consensus coverage) for all samples (n = 17; Methods and Supplementary Table 1).

Overall, cfDNA captured a median of 70% of all tumour SNVs in high TFx samples (≥8%) and 34% in low TFx samples (<8%), independent of the death-to-autopsy interval (Fig. 5a, Extended Data Fig. 10d and Supplementary Table 8). SNV detection was highly dependent on clonality, with nearly all (median 98–100%) of founder, 79–90% of shared and 28–32% of private SNVs captured in samples with a high TFx (Fig. 5b). We observed similar trends for copy number gains detected by cfDNA WGS (Extended Data Fig. 10e). Notably, every tumour with a private clone was detected in cfDNA WGS at a higher rate than explained by sequencing error, a result that highlights the advantages of genome-wide analysis (Fig. 5c, Methods and Extended Data Fig. 10f). The detection of shared and private mutations varied extensively and was not explained by TFx alone (Extended Data Fig. 10g). Instead, shared clone capture correlated with the TFx-adjusted patient clone abundance (total abundance of each clone across tumours, ρ > 0.38, P < 0.018; Fig. 5d and Methods). However, the detection rates of private SNVs were not well explained by cellular prevalence, nor metastatic site, but rather by tumour-specific purity (from histology, ρ = 0.59, P = 0.0049 for WGS) and volume (from CT scans, ρ = 0.53, P = 0.02 for TP) of the individual tumour with each private clone (Fig. 5e,f and Extended Data Fig. 10h–j).

Nearly all BLCA-driver SNVs identified in tumours (median 93–100% per patient) were detected in high TFx cfDNA samples, with variable (median 64%) detection in low TFx samples (Fig. 5g,h). The ultradeep TP enabled the capture of additional potentially actionable private SNVs, including FGFR3 (R248C, patient 19-022) and NOTCH1 (C575C, patient 16-079) (Fig. 5g). BLCA-driver genes altered by CNAs were similarly well-captured in cfDNA WGS (overall median 83%) (Extended Data Fig. 10k), with 94% and 100% of founder gains and losses, respectively detected (Fig. 5i). Altogether, these results demonstrate that most early (98–100%) and driver (93%) alterations are readily captured from cfDNA, and detection of later-evolving SNVs is influenced by clonal burden, tumour volume and purity and cfDNA TFx.

Analyses of nucleosome occupancy from cfDNA WGS data can reveal differential transcriptional activity between tumour subtypes^43,44. We profiled nucleosome patterns to infer cfDNA contribution from different tissues, which enabled us to differentiate BLCA from healthy cfDNA and to identify shedding from immune cells and tissue damage (Methods and Extended Data Fig. 11a,b). Next, we analysed nucleosome accessibility at transcription-factor-binding sites (TFBSs). We observed increased transcription factor activity (reduced cfDNA coverage at TFBSs) of the known NE markers⁴⁵ ASCL1 and SOX2 specifically in patient 19-037, who had NE tumours (Fig. 5j and Extended Data Fig. 11c). Moreover, low activity of HES1 in cfDNA was exclusive to patient 15-109, a result consistent with the observed reduced expression of HES1 in PUC tumours (Extended Data Fig. 11d). The use of cfDNA may aid in the minimally invasive identification of poor-prognosis PUC tumours, as reduced HES1 activity can promote EMT and increase tumour aggressiveness⁴⁶.

We next sought to differentiate between luminal and basal subtypes in UC and UCSD tumours, as basal tumours have poorer prognosis⁸. We identified RARG, GRHL2 and PPARG among the transcription factors with increased activity in the LumP subtype (Fig. 5k,l and Supplementary Table 8). By contrast, HMGA1 activity was increased in the Ba/Sq subtype, consistent with its known function in stemness and EMT⁴⁷. Notably, analyses of cfDNA from patient 19-022 revealed partial luminal transcription factor activity (PPARG and RARG) alongside reduced ZNF217 and ZNF146 activity, which are typically suppressed in basal tumours⁴⁸, a result consistent with the co-occurrence of LumP and Ba/Sq lesions in this patient (Fig. 4b). Overall, inferred transcription factor activity was able to delineate both histological and molecular subtypes, which indicated that cfDNA can provide insights into metastatic BLCA prognosis (Extended Data Fig. 11e).

This study presents a comprehensive multi-omic analyses of a unique BLCA rapid autopsy programme featuring a broad representation of metastatic BLCA histological subtypes. Genomic evolutionary analyses revealed that metastasis-to-metastasis seeding was the dominant pattern of metastatic BLCA spread, consistent with reports from other cancers⁴⁹. CT scans corroborated metastatic seeding patterns, although the possibility of early-seeding micrometastases cannot be excluded. Complex seeding was associated with poor survival, variant histology and immunotherapy treatment. Multiregional or longitudinal sampling could refine these results, including differentiating between the effects of therapeutic pressure versus natural clonal evolution. Our findings indicate that earlier intervention for oligometastatic disease in metastatic BLCA may help suppress metastasis-to-metastasis spread, limit highly prevalent polyclonal evolution and extend survival.

Genomic phylogenetic analyses revealed enrichment of founder driver alterations in PUC and possibly NE subtypes, which may explain their aggressive disease course. However, non-founder drivers were present in every patient, which suggests that single-sample profiling may miss late events that could underlie resistance mechanisms. Analyses of additional tumours is needed to validate specific genes observed here as late driver events. We observed that the late platinum chemotherapy mutation signature recently described in BLCA⁶ was reduced in PUC tumours and associated with increased FANC pathway expression. BLCA cell line experiments revealed that high FANC expression was sufficient to reduce sensitivity to cisplatin. Although further mechanistic experimentation is required, these results highlight a possible clinical rationale for targeting FANCD2 deubiquitination, such as with USP1 inhibitors, to disrupt the FANC pathway and sensitize PUC tumours to chemotherapy^50,51.

Our results revealed heterogeneity in molecular and histological subtypes. PUC tumours are known to be luminal⁵²; however, we observed that PUC metastases were LumU, whereas primary tumours were stroma-rich. This result highlights the evolutionary progression and role of the stroma in early disease. Furthermore, patients with even minor squamous differentiation often had entirely Ba/Sq tumours, which suggests that minimal variant histology can have patient-wide significance. Despite homogeneity across bulk RNA-seq data in most patients, snRNA-seq uncovered transitional states and heterogeneous subtyping in PUC and UCSD tumours, which indicated potential transdifferentiation or mixed phenotypes^52,53. This study, to our knowledge, is the first to investigate the TME of metastatic BLCA histological subtypes at single-cell resolution. Signalling interactions between epithelial, fibroblast and immune cells cooperated to create an immune-inflamed TME in PUC lymph nodes and immunosuppressive TME in UCSD tumours. In end-stage UCSD tumours, the increased immunosuppressive myeloid populations may underlie poor ICI responses, consistent with previous reports^13,54. Despite the small sample size of our study, these findings highlight the potential for patient stratification to immunotherapy by histological subtype^9,55.

Analyses of cfDNA captured nearly all founder and pathogenic alterations, whereas detection of late-evolving mutations depended on tumour clonality, volume and purity. ctDNA and radiological imaging have previously been compared in cancer surveillance^56,57, but we highlight a direct association between cfDNA and CT volumetric data. Future study of clonal haematopoiesis of indeterminate potential mutations is required to identify ctDNA-specific alterations. Furthermore, cfDNA distinguished transcriptional patterns differentiating histological and molecular subtypes, thereby expanding on phenotypic profiling performed in other cancers^{42,44,58,59,60}. This work highlights the potential use of cfDNA in metastatic BLCA management, although studies involving larger cohorts and on-treatment collection are needed to validate its clinical utility.

This study generated a compendium of multi-omic datasets from a distinct cohort of rapid autopsy cases. Despite the modest cohort size, it is among the largest of its kind, including multiple tumours per patient and histological subtypes. Our results establish the foundational landscape of tumour evolution and heterogeneity of BLCA subtypes, uncovering molecular characteristics that shape the natural history of progression to end-stage disease.

Patient recruitment and sample collection

Specimens were obtained from 20 patients, including histopathologically normal tissue (n = 20, flash-frozen tissue), primary tumour samples (n = 24; 16 FFPE and 8 flash-frozen tissue) and 1–7 metastases per patient (n = 80, flash-frozen tissue). Details of tissue source, location, histology, among other information, for each sample are provided in Supplementary Table 1. Autopsy samples (n = 108), which included all normal, all metastatic and 8 primary tumour samples, were obtained within 11.5 h (median 3.7 h) of death as part of the rapid autopsy programme (Supplementary Table 1). For the logistical framework of the rapid autopsy programme, see Fig. 1a. A representative selection of the metastases observed at autopsy (excluding bone metastases) was chosen for sequencing, in addition to autopsy-derived primary tumour samples when available. Archival FFPE primary tumour samples (n = 16) included diagnostic (pre-treatment) TURBT (n = 5) and surgical specimens (n = 11) (Supplementary Table 1). All primary tumour and metastatic specimens selected for sequencing had more than 80% tumour cellularity based on genitourinary pathologist review (F.V.-L.). All normal tissues were confirmed to be normal by histology.

Study approval

All samples were obtained from patients with signed informed consent documents under the aegis of the Genitourinary Cancer Biorepository at the University of Washington (University of Washington IRB 2341). All 20 patients signed written informed consent for the rapid autopsy programme. Metastases and the primary tumour (if present) were identified and collected. In accordance with study protocols approved by the institutional review board, no metastatic biopsy samples collected from living patients were obtainable for use in this study.

Sectioning, H&E staining and pathologist assessment of tumour sections

All visceral metastases and matched primary tumour were embedded in Optimal Cutting Temperature compound (OCT; Tissue-Tek, Sakura Finetek) or FFPE. Haematoxylin and eosin (H&E) staining was completed as previously described^18,61. H&E-stained sections (5 µm) of the OCT-embedded and FFPE tissues were reviewed by an independent, dedicated genitourinary pathologist (F.V.-L.). The percentage of tumour, necrosis and histological subtypes were assessed and recorded for each section. Between one and ten sections per tumour or normal sample were assessed, depending on the size of the specimen available. The areas of >80% tumour were marked for macrodissection for DNA and RNA extraction purposes.

Histology was defined both at the patient level and the tumour (sample) level. All classifications were based on pathology assessment of FFPE H&E-stained specimens by pathologist (F.V.-L.) review based on the WHO classifications of ‘Urinary and Male Genital Tumours’ (5th edition)⁴. NE histology was further verified by synaptophysin positivity. PUC, UC-Sarc and UCSD tumours are UC with plasmacytoid, sarcomatoid or squamous differentiation, respectively. Patients and tumours were designated as non-UC if there was any evidence of variant histology. For patient-level histology, the histology was defined as a variant if there was evidence of subtype histology in any tumour in that patient (for example, if a patient had any tumour exhibiting any squamous differentiation, they were classified as UCSD). Patients 16-070, 17-030, 18-101 and 19-022 were designated UCSD because they had a mix of UC and UCSD tumours, whereas all tumours in patients 16-097 and 17-026 exhibited UCSD histology. Sample-level histology was defined for each tumour individually (for example, the primary tumour for patient 16-070 is UCSD, whereas their metastases are UC at the sample level). Further details on the percentage of histological subtypes in each tumour are given in Supplementary Table 1 (first sheet).

HER2 immunohistochemistry

FFPE sections (5 μm) were deparaffinized and rehydrated in sequential xylene and graded ethanol series¹⁷. HER2 protein expression was evaluated using a PATHWAY anti-HER2/neu (clone 4B5) kit on a Ventana BenchMark ULTRA platform (Ventana/Roche Tissue Diagnostics) following the manufacturer’s standardized protocols. Staining intensity (scored 0 to 3+) and the percentage of positive tumour cells were assessed by two pathologists (M.C.H. and E.S.).

Radiological assessment of CT scans

Clinical staging CT scans were requested for all 20 patients. A total of 59 scans (range 0–19 per patient) were available and acquired (Supplementary Table 2). These were scrubbed of all identifying metadata and the de-identified DICOM data were then uploaded to MIM (MIM Software, v.7.3.7, build O606-01). All CT scans were centrally reviewed by an independent radiologist as part of standard of care, and the de-identified radiology reports and scans were further evaluated by a radiation oncologist (O.Y.M.) to confirm the presence of the metastatic sites collected at autopsy. Primary and metastatic lesions corresponding to pathologically confirmed sites of disease at autopsy were circumscribed on all axial slices. Anterior–posterior plane projection digitally reconstructed radiographs were generated in MIM (Supplementary Fig. 2), with superimposed two-dimensional colour projections of all circumscribed radiographically identifiable malignant lesions (Fig. 1c and Extended Data Fig. 1e).

Corroboration of metastatic seeding

For each patient with metastasis-to-metastasis seeding (n = 9), the first available clinical staging CT scan after diagnosis was assessed for visible tumours to determine whether an initiating metastasis could be identified. For five patients in whom the early CT scan identified a first metastasis, subsequent CT scans were assessed to follow metastatic evolution. In patient 19-037, the sequenced peri-pancreatic LN corresponded to a portion of the RP LN visible in the CT scan. The remaining four patients were unable to be further assessed because of a lack of any staging CT scans (patient 15-108), all sequenced metastases were visible at the first CT scan (patients 16-097, 19-001) or there was a lack of clear correspondence between lesions visible on the first scan and metastatic sites collected at autopsy (patient 17-071) (Supplementary Fig. 2).

Determination of pre-ICI versus post-ICI metastasis

For each patient who underwent at least one full cycle of ICI treatment (n = 9; Supplementary Table 2), the CT scan that immediately preceded the start of ICI treatment was assessed for visible tumours. The time between the CT scan preceding ICI and the start of ICI was variable (median 16 days, range 2–260 days). For patient 19-022, the most recent CT scan available before ICI was 260 days before treatment, but clinical records confirmed that only the lung metastasis was observed before starting ICI treatment. The metastases already in place before ICI were designated as ‘pre-ICI metastases’, whereas metastases collected and sequenced at autopsy but not visible by CT before ICI were designated ‘post-ICI metastases’. Of the nine ICI-treated patients, there were three for whom all metastases were visible by CT before the start of ICI therapy. Six patients had both metastases visible pre-ICI and metastases that only appeared post-ICI. In total, 17 metastases were designated pre-ICI and 15 were post-ICI (Supplementary Table 3). Clonal composition (number and cellular prevalence of clones), metastatic seeding and metastatic site were compared between metastases seeded pre-ICI versus post-ICI.

We explored whether we could perform a similar retrospective seeding analysis for patients who received chemotherapy as we did for ICI treatment. Unfortunately, only three patients in our cohort (15-109, 17-020 and 17-026) received chemotherapy alone. Moreover, these three patients received chemotherapy early in the disease course; therefore, no metastases could be confirmed by CT scans before treatment.

Tumour volume calculations

Tumour volumes (in cubic centimetres) were calculated in MIM from patients’ last CT scan before death and exported in site-specific tabular format for subsequent correlative analyses.

Blood collection and fractionation from autopsy and healthy donors

Cardiac or venous blood collection was performed within 2.8–11.5 h of death, at the start of the rapid autopsy, using a 16–18 gauge needle and 20–50 cc syringes into Tiger Top (serum) or Purple Top (plasma) blood collection tubes and spun at 1,600g for 15 min. Supernatant plasma or serum were aliquoted using a Serum Bank pipette into 560–600 μl aliquots and stored at –80 °C until cfDNA extraction. Healthy donor plasma was purchased from Research Blood Components (no. 1-009, Fresh Plasma Single Donor 10 ml). Blood was collected in Cell-free DNA BCT Streck tubes and plasma was separated by centrifugation of whole blood at 1,600g for 10 min at room temperature. The upper plasma layer was removed, leaving the buffy coat undisturbed. Plasma was shipped on dry ice and stored at –80 °C until cfDNA extraction.

cfDNA extraction from plasma and serum

cfDNA was extracted using a QIAamp Circulating Nucleic Acid kit (Qiagen). After thawing 0.5–1.5 ml plasma or serum on ice, samples were spun at 15,000g for 10 min at room temperature. Supernatant was removed, avoiding disrupting any remaining pellet, and processed according to the manufacturer’s instructions, eluting cfDNA with 55 μl Buffer AVE. cfDNA concentration was quantified using a Qubit 1× dsDNA High Sensitivity Assay kit (Thermo Fisher Scientific, Q33230) on a Qubit 4 Fluorometer. cfDNA quality and fragment-size distribution were assessed from 3 μl (range 100–4,000 pg μl^–1) cfDNA run with an Agilent Cell-free DNA ScreenTape assay (Agilent Technologies) on an Agilent 4200 TapeStation system following the manufacturer’s instructions.

cfDNA sequencing

Seventeen patients had post-mortem cfDNA samples available from plasma (n = 7) or serum (n = 10), with high cfDNA yields (median 440 ng ml^–1) and detectable ctDNA of ≥3% TFx (median 5.9%). Ultralow-pass (ULP) WGS was performed for all samples. Deeper WGS and/or TP was subsequently performed on selected samples based on the estimated TFx.

ULP library preparation

ULP library preparation and sequencing were performed by Broad Clinical Labs for both plasma and serum samples. Initial cfDNA input was normalized to be within 25–52 ng in 50 μl TE buffer (10 mM Tris-HCl 12 mM EDTA, pH 8.0) according to PicoGreen quantification, with no shearing of cfDNA before library construction. Library preparation was performed using a commercially available kit (KAPA HyperPrep kit with Library Amplification product KK8504, KAPA Biosystems) and xGen UDI-UMI duplex adapters (Integrated DNA Technologies (IDT)). Unique 8-bp dual index sequences embedded in the p5 and p7 primers (IDT) were added during PCR. Enzymatic clean-up was performed using Beckman Coulter AMPure XP SPRI beads (Beckman Coulter, A63880), with elution volumes reduced to 30 μl to maximize library concentration. Library quantification was performed using an Invitrogen Quant-It broad range dsDNA quantification assay kit (Thermo Fisher Scientific, Q33130) with 1:200 PicoGreen dilution. Following quantification, each library was normalized to a concentration of 35 ng µl^–1 using Tris-HCl, 10 mM, pH 8.0.

ULP library pooling and sequencing

In preparation for ULP libraries, approximately 4 µl of the normalized library was transferred into a new receptacle and further normalized to a concentration of 2 ng µl^–1 using Tris-HCl, 10 mM, pH 8.0. Following normalization, up to 95 ULP WGS samples were pooled together using equivolume pooling. The pool was quantified by qPCR and normalized to the appropriate concentration to proceed to sequencing. Cluster amplification of library pools was performed according to the manufacturer’s protocol (Illumina) using Exclusion Amplification cluster chemistry and NovaSeq SP flow cells. Flow cells were sequenced on a NovaSeq 6000 Sequencer using the XP workflow and a v.1.5 300 cycle NovaSeq SP kit. Each pool of ULP whole-genome libraries was run on one lane using paired 151 bp runs.

TFx estimation using ichorCNA

ULP WGS data were analysed using ichorCNA⁶² (https://github.com/GavinHaLab/ichorCNA; commit ID: d31ed52) to estimate the TFx. Settings for hg38 build and 1 Mbp bin size were used. The final configuration file is available from GitHub (https://github.com/GavinHaLab/BLCA-subtype-evolution-paper). The TFx estimate was used to inform follow-on deeper WGS and TP sequencing.

Deep WGS from ULP libraries for high TFx ctDNA samples

Following initial ULP sequencing, selected ULP libraries with a TFx of ≥8% underwent further sequencing. Libraries were initially normalized to 2 ng µl^–1, pooled, then quantified by qPCR using a KAPA Biosystems kit that uses probes specific to the ends of the adapters. On the basis of the qPCR results, the libraries were adjusted to 2.2 nM before proceeding to the next sequencing stage using an automated Agilent Bravo liquid-handling platform. Pools were denatured with sodium hydroxide, diluted using an Illumina-provided pre-load buffer and transferred to a uniquely-barcoded 8-lane strip tube with a Hamilton Starlet liquid handler. Strip tubes were loaded into a 300 cycle NovaSeq X 25B kit and the run was initiated with a 151-bp end, dual-indexed read structure. On the basis of the pool size, the number of lanes was calculated to ensure samples reached the desired mean coverage.

PanCancer TP sequencing

After ULP library construction, in-solution hybridization and TP capture were performed using the relevant components of a XGen hybridization and wash kit (IDT) following the manufacturer’s suggested protocol, but with several exceptions. A set of 12-plex pre-hybridization pools were created. These pre-hybridization pools were created by equivolume pooling of the normalized libraries, human COT-1 and IDT XGen blocking oligonucleotides. The pre-hybridization pools underwent lyophilization using Biotage SPE-DRY. After lyophilization, custom PanCancer bait (Twist Biosciences) along with hybridization master mix were added to the lyophilized pool before resuspension. Samples were incubated overnight. Library normalization and hybridization setup were performed on a Hamilton Starlet liquid-handling platform, whereas target capture was performed on an Agilent Bravo automated platform. After capture, PCR was performed to amplify the capture material. After post-capture enrichment, library pools were quantified by qPCR (automated assay on the Agilent Bravo) using a kit purchased from KAPA Biosystems with probes specific to the ends of the adapters. On the basis of qPCR quantification, pools were normalized using a Hamilton Starlet to 2 nM and sequenced using Illumina sequencing technology.

TP cluster amplification and sequencing

Cluster amplification of library pools was performed according to the manufacturer’s protocol (Illumina) using Exclusion Amplification cluster chemistry and HiSeqX flow cells. Flow cells were sequenced on v.2 Sequencing-by-Synthesis chemistry for HiSeqX flow cells. The flow cells were then analysed using RTA (v.2.7.3 or later). Each pool of libraries was run on paired 151 bp runs, reading the dual-indexed sequences to identify molecular indices and sequenced across the number of lanes needed to meet coverage for all libraries in the pool.

Tissue WES and WGS

WES was performed as previously described¹⁶. WGS of flash-frozen normal control tissue (30× WGS) and FFPE or flash-frozen tumour samples (60× WGS) were performed by Broad Clinical Labs. gDNA derived from FFPE samples for a subset of primary tumours was sequenced using Human Whole Genome Sequencing PCR Plus (v.1.1–v.1.3), whereas a PCR-free method was used for gDNA from flash-frozen tissue. An aliquot of gDNA (100 ng in 50 µl for FFPE samples, 350 ng for flash-frozen samples) was used as input into DNA fragmentation by acoustic shearing using a Covaris focused-ultrasonicator, targeting 385 bp and 350 bp fragments, respectively. Following fragmentation, additional size selection was performed using SPRI cleanup. Library preparation was performed using a commercially available kit (KAPA Hyper Prep with Library Amplification Primer Mix, product KK8504, KAPA Biosystems) and with palindromic forked adapters using unique 8-bp index sequences embedded in the adapter (purchased from Roche). For FFPE samples, the libraries were then amplified using 10 cycles of PCR. Following sample preparation, libraries were quantified by qPCR (kit purchased from KAPA Biosystems) with probes specific to the ends of the adapters on an automated Bravo liquid-handling platform (Agilent). On the basis of qPCR quantification, libraries were normalized to 2.2 nM and pooled into 24-plexes. Sample pools were combined with NovaSeq Cluster Amp Reagents DPX1, DPX2 and DPX3 and loaded into single lanes of a NovaSeq 6000 S4 flow cell using a Hamilton Starlet liquid-handling system. Cluster amplification and sequencing were performed on NovaSeq 6000 instruments using sequencing-by-synthesis kits to produce 151 bp paired-end reads. Output from Illumina software was processed using the Picard data-processing pipeline to generate BAM files containing demultiplexed, aggregated aligned reads.

TURBT WGS library preparation and sequencing

Sample processing for FFPE TURBT specimens was performed in similar manner to the FFPE primary samples. Tumour regions were macrodissected to achieve >80% tumour content, and dual DNA and RNA extraction was carried out at the Fred Hutchinson Cancer Center core facility. For two patients (17-020 and 17-026), normal tissues that were originally sequenced by WES were re-processed for WGS to serve as controls for the WGS TURBT samples, following the same protocol. FFPE tumour DNA and flash-frozen normal DNA were quantified using an Invitrogen Qubit 2.0 Fluorometer (Thermo Fisher Scientific). For library preparation, 100 ng DNA per sample was fragmented to a target size of 400 bp using a Covaris LE220-plus focused ultrasonicator. Sequencing libraries were then prepared from the fragmented DNA using a xGen cfDNA & FFPE DNA Library Prep v2 MC kit (IDT) in combination with xGen UDI indexing primers (IDT). Library quantification was performed using an Invitrogen Qubit 2.0 Fluorometer, and fragment size distribution was assessed using an Agilent 4200 TapeStation (Agilent Technologies). Individual libraries were pooled in 8-plex at weighted molar concentrations based on DNA type. Sequencing was conducted on an Illumina NovaSeq X Plus system (Illumina) across two lanes of a 25B-300 flow cell using a paired-end 150 bp read configuration, with an average sequencing output per library of 837 million read pairs (range 711–1,107 million).

Bulk tumour RNA-seq

For flash-frozen autopsy samples (n = 74), RNA was isolated from OCT-embedded specimens using a combination of RNA STAT-60 reagent (Tel-Test) and a RNeasy Mini kit (Qiagen), with an in-solution DNase treatment step included before purification. For FFPE surgical samples (n = 8), RNA was isolated using an AllPrep DNA/RNA FFPE kit (Qiagen). RNA quality was assessed by measuring the RNA integrity number using an Agilent Bioanalyzer (Agilent Technologies). Sample 18-016_G2 was excluded from bulk RNA-seq analysis owing to a low RNA integrity number score of extracted RNA. For RNA-seq, libraries were generated from 300 ng total RNA using an Illumina TruSeq RNA Exome Sample Prep kit, following the standard protocol (Illumina). The resulting barcoded libraries were pooled and sequenced on an Illumina HiSeq 2500 platform to produce 50-bp paired-end reads.

snRNA-seq

snRNA-seq was performed by Singulomics as follows. Nuclei were isolated from OCT-embedded frozen human BLCA tissue samples, and 3′ single-cell gene expression libraries (Next GEM v.3.1) were constructed using the 10x Genomics Chromium system. Libraries were sequenced with around 200 million 150-bp paired-end reads per sample on an Illumina NovaSeq X Plus. The sequencing reads were analysed using human reference genome GRCh38 with Cell Ranger (v.7.1.0). Introns were included in the analyses.

Cell culture experiments

Cell lines used in this study were SW780, HT1376, BFTC905 and UBLC1 BLCA cells. All cell lines were sourced from the American Type Culture Collection, authenticated using STR profiling and tested mycoplasma negative. Each cell line was cultured in DMEM high-glucose and l-glutamine medium (Gibco 11965-092) supplemented with 10% FBS, 1% penicillin–streptomycin, 1% l-glutamine and 1% MEM non-essential amino acids. Cell lines were cultured in an incubator at 37 °C with 5% CO₂.

Western blotting

Cells were treated in 10 cm dishes with 5 µM cisplatin or 1 µM MMC (Sigma, M5353) for 24 h before collection for western blotting. Cell pellets were lysed using RIPA buffer (Fisher Scientific) supplemented with protease and phosphatase inhibitors (Sigma). Equivalent amounts of protein per sample (5–40 μg depending on the target) were run on a polyacrylamide gel and transferred to a PVDF membrane. Proteins were detected using respective primary and HRP-conjugated secondary antibodies, including rabbit anti-FANCD2 (Abcam 108928 1:1,000), mouse anti-FANCF (Santa Cruz sc-271952 1:500), mouse anti-tubulin (Sigma-Aldrich T8203 1:1,000), rabbit anti-vinculin (Cell Signaling Technology 13901S 1:1,000), goat anti-rabbit (Invitrogen 31460 1:5,000) and goat anti-mouse (Invitrogen 31430 1:5,000). Quantification of images was performed in ImageJ.

siRNA-mediated FANCF knockdown

FANCF knockdown was performed using siRNA transfection in the FANC-high BLCA cell line SW780. Cells were seeded in parallel in 96-well plates (3,000 cells per well) for live-cell imaging assays and in 6-well plates (2.5 × 10⁵ cells per well) for RNA and protein analyses. Cells were transfected the following day using Lipofectamine RNAiMAX (Thermo Fisher Scientific) according to the manufacturer’s instructions (forward transfection), with either non-targeting control (Dharmacon D-001810-10-05) or FANCF-targeting (Dharmacon L-014206-00-0005) siRNA SMARTpools. At 24 h after transfection, cells plated in 6-well plates were collected for RNA extraction to assess knockdown efficiency by qPCR. Protein lysates were collected at 48 h after transfection for immunoblot validation.

Lentiviral overexpression of FANCF

FANCF overexpression was achieved by lentiviral transduction in the FANC-low BLCA cell line UBLC1. Plasmids encoding empty vector control (pCR1265) or FANCF (pCR1446) were obtained from Addgene⁶³. Lentivirus was produced in HEK293T cells via calcium phosphate (CaCl₂–HBSS) transfection using third-generation packaging plasmids. Viral supernatant was collected 72 h after transfection, filtered through a 0.45 μm membrane, aliquoted and stored at −80 °C before use. UBLC1 cells were transduced with lentiviral supernatant overnight without transduction enhancers. Following transduction, cells were selected with geneticin (G418; 1 mg ml^–1) for 7 days to establish stable expression. Stable expression of FANCF was confirmed by both qPCR and immunoblotting.

qPCR

Total RNA was extracted from cell pellets using a RNeasy Plus Mini kit (Qiagen). iScript Reverse Transcription Supermix (Bio-Rad) was used to reverse-transcribe cDNA from equal volumes of RNA across samples. qPCR was performed using SsoAdvanced Universal SYBR Green Supermix (Bio-Rad) for FANCF (5′-GGTGGCGGCTAGTCACTAAA-3′, 5′-GCTAGTCCACTGGCTTCTGG-3′) and GAPDH (5′-ACCCACTCCACCTTTGAC-3′, 5′-ATGAGGTCCACCACCCTGTTG-3′).

Immunofluorescence staining for γH2AX

Cells (5 × 10⁴) were seeded onto poly-l-lysine-coated glass coverslips in 24-well plates and allowed to adhere overnight. The following day, cells were treated with cisplatin at doses approximating the IC₅₀ for each cell line. FANC-high cells (SW780) were treated with 10 μM cisplatin, whereas FANC-low cells (UBLC1) were treated with 5 μM cisplatin. Cells were exposed to cisplatin for 3 h, after which the drug-containing medium was removed and replaced with fresh medium. Cells were fixed at 0, 4, 24 and 48 h following drug removal. For immunofluorescence staining, cells were rinsed briefly with PBS and fixed with 4% paraformaldehyde in PBS for 10 min at room temperature. Fixation was quenched using 500 mM Tris and 125 mM glycine for 5 min, followed by washing with PBS. Cells were permeabilized with 0.5% Triton X-100 in TBS for 5 min and washed once with PBS. Samples were blocked using MaxBLOCK (Active Motif) for 1 h at 37 °C, followed by 2 washes with TBS containing 0.1% Tween-20 (TBST). Cells were incubated with a primary antibody against γH2AX (mouse; Sigma, 05-636; 1:4,000) diluted in 1% BSA in TBST for 1 h at room temperature. Coverslips were washed twice with TBST and incubated with Alexa Fluor 633 goat anti-mouse secondary antibody (ThermoFisher, A-21052; 1:1,000) for 1 h at room temperature in the dark. Following secondary incubation, cells were stained with DAPI (2 μg ml^–1 in TBST) for 10 min, washed twice with TBST and rinsed once with water to remove residual salts. Coverslips were mounted onto glass slides using ProLong Diamond Antifade mountant (Thermo Fisher Scientific) and allowed to cure at room temperature before imaging. Slides were stored at 4 °C. Images were acquired with a Leica Stellaris 8 scanning confocal microscope equipped with a 40×/1.3 NA oil HC PL APO CS2 objective. A total of 100 z stacks were acquired for each slide. DAPI and Alexa Fluor 594 were imaged using 420–504 nm and 600–750 nm detection windows, respectively.

Images were analysed using Imaris (v.11) by first segmenting nuclei based on DAPI staining using the machine-learning surface tool with an estimated diameter of 8 μm in UBLC1 cells and 7 μm in SW780 cells. Nuclei with a median γH2AX intensity of more than 15 were excluded from further analyses, as these were probably dying cells. The median number of nuclei per image after filtering was 1,130 (minimum = 537, maximum = 1,395). γH2AX foci were identified in nuclei using the spots tool with an estimated diameter of 0.4 μm, minimum mean intensity of 15 and minimum quality of 3.5. The number of foci per nucleus was calculated using the Imaris cells tool.

Cell growth and death assays

Cell lines were plated at between 5,000 and 15,000 cells per well in 96-well plates to achieve starting densities of approximately 30%. The next day, the cell medium was changed to include 0–20 µM cisplatin (0, 1, 2.5, 5, 10 and 20 µM doses) or 0–10 µM MMC (0, 0.05, 0.10, 0.5, 1, 2.5, 5 and 10 µM doses) and, for cell death experiments, 250 nM of IncuCyte Cytotox dye (Sartorius). Each cell line was plated in four or five technical replicate wells per condition. Cells were moved immediately after treatment to an IncuCyte S3 or IncuCyte SX5 platform for continuous imaging and confluence analyses. Cells were imaged over the course of 72 h, and IncuCyte software (v.2023A Rev1, Sartorius) was used to measure changes in confluence and to count dead and dying cells (fluorescent puncti from Cytotox dye) in each well at each time point.

For cell growth analyses, the cell confluence was normalized to the starting confluence for each well. To calculate IC₅₀ values, dose–response curves were fit using a four-parameter log-logistic (LL.4) function and applied to these normalized growth values as a function of cisplatin dose. Curve fitting was performed using the drm function of the R package drc (v.3.0-1). The auc function from the R package MESS (v.0.6.0) was used to calculate AUC values from the normalized confluence values.

For cell death analyses, the count of dead and dying cells (Cytotox fluorescent puncti) in each well was normalized to the cell confluence of that well to obtain a representation of the proportion of dying cells. The auc function of the R package MESS (v.0.6.0) was used to calculate AUC values from the normalized dead cell values.

Sequence alignment

DNA sequencing reads were aligned to the GRCh38 human genome using BWA-MEM (v.0.7.17). Alignments were then sorted and indexed using samtools (v.1.10), and duplicates were marked using picardtools MarkDuplicates (v.2.18.29). Finally, alignments were subjected to base quality score recalibration using BaseRecalibrator and ApplyBQSR from GATK (v.4.1.8.1). Read counts and the per cent of properly mapped reads were collected using CollectAlignmentSummaryMetrics, and CollectWgsMetrics was used to gather the mean coverage (GATK v.4.1.8.1). These steps generated recalibrated BAM files that were used for further analyses. For the WES samples, metrics were collected using the tool HsMetrics and an exome capture bedfile. The healthy donor plasma samples were processed using the same approach as described above to generate the final BAM files for analyses. The sequencing metrics were computed using CollectWgsMetrics.

BLCA-driver gene list creation

We compiled a total of 128 known potential metastatic BLCA-driver genes. First, we used the IntOGen database⁶⁴ to identify a list of 95 genes from 8 cohorts comprising 867 samples. Next, to enhance this list, we manually incorporated 24 additional genes identified by the Hartwig Medical Foundation. Finally, we included nine genes reported in BLCA literature as commonly affected in metastatic BLCA^{6,65,66,67,68}. This curated list of 128 driver genes was used for all downstream analyses.

Next, we annotated the functional roles of these 128 genes in BLCA by using the Cancer Gene Census⁶⁹ and OncoKB⁷⁰ databases. For each gene, if both databases listed the gene as a tumour suppressor or as an oncogene, then this annotation was used. If the databases differed with each other, did not list an annotation or listed multiple conflicting annotations, then the gene was designated ‘Other’. The driver genes and annotations that were used for all analyses are provided in Supplementary Table 4.

Somatic mutation consensus calling

To identify high-confidence somatic mutations from tumours, we implemented a consensus-based variant calling strategy using four somatic callers: Mutect2 (ref. ⁷¹), Strelka2 (ref. ⁷²), VarScan2 (ref. ⁷³) and MuSE⁷⁴. All tools were run in paired tumour–normal mode, leveraging matched normal samples to exclude germline variants. For Mutect2, we also constructed a panel of normals to remove recurrent technical artefacts. Different panels of normals were generated for WES (n = 7) and WGS (n = 13) cohorts using the ‘Create a Panel of Normals’ workflow in GATK⁷⁵.

Genomic variant annotation was conducted using ANNOVAR (v.2020-06-07)⁷⁶ with the table_annovar.pl script to functionally interpret variants identified in the study. The analysis was performed using the human genome build GRCh38 (Broad version). The input VCF file was annotated against a comprehensive set of databases, including RefSeq genes (refGene) for gene-based annotation and cytogenetic band information (cytoBand) for region-based annotation. The variant frequency in the population and functional scores were annotated using clinical databases, including the NHLBI-ESP 6500 dataset (esp6500siv2_all), dbSNP (v.144 and v.150; avsnp144 and avsnp150), the 1000 Genomes Project (ALL.sites.2015_08), gnomAD genome and exome datasets (gnomad_genome and gnomad_exome), ExAC (v.0.3; exac03), ClinVar (release 20190305), InterVar (release 20180118) and dbNSFP (v.3.3a) for functional prediction scores, and gnomAD (v.3.1.2) genome data (gnomad312_genome) and COSMIC (v.70) for somatic mutation data. The annotation process used the --vcfinput flag to accept VCF format input and the -polish option to refine annotations. Variants with missing annotations were marked with a dot using the -nastring option, and intermediate files were removed after processing using the -remove flag.

Criteria for high-confidence consensus mutation filtering

A SNV was considered high-confidence if it was identified by at least two out of the four callers and met the following somatic call filtering criteria: (1) the matched normal sample contained ≥8 reference reads; (2) the tumour sample had ≥14 total reads and ≥5 mutant reads; (3) the tumour variant allele frequency (VAF) was ≥1%; and (4) the normal VAF was ≤5%. The frequency of SNVs greater than 10% in gnomAD or ExAC databases were filtered out to generate the final high-confidence somatic SNV set for downstream analyses.

For calling indels, we used, Mutect2, Strelka2 and SvABA⁷⁷, also in paired mode with the corresponding matched normal. Indels were retained only if they were detected by all three callers and the Mutect2 tumour VAF was >1%. Mutations with a frequency greater than 10% in gnomAD or ExAC databases (to exclude common population variants) were filtered out to generate the final consensus indel call set for downstream analyses. Last, for all analyses involving somatic mutations, four primary tumour samples, 17-047pD10, 18-101M3, 18-120pA14 and 19-044pB11, were excluded because they had a tumour purity of <20% or had poor-quality FFPE DNA sequencing data.

Analysis of positive and negative selection

We used the dNdScv R package (v.0.010, https://github.com/im3sanger/dndscv) to estimate dN/dS ratios across the cohort, which enabled detection of positive or negative selection in a set of mutations⁷⁸. This method applies maximum-likelihood modelling of synonymous and non-synonymous mutations, accounting for gene-specific background mutation rates, sequence context and trinucleotide mutational signatures. All somatic mutations were included, stratified by founder, shared and private status. Analyses were performed separately for BLCA-driver genes and all other genes (‘All genes’) and only the global dN/dS estimates are reported (Supplementary Fig. 3).

CNA analysis

Copy number calls

To identify CNAs, we used the TITAN⁷⁹ (v1.15.0) pipeline (https://github.com/GavinHaLab/TitanCNA_SV_WGS; commitID: bedbd76). TITAN-derived results were manually curated to confirm the optimal solution for each sample and are shown in Supplementary Table 1. The ploidy and purity of each tumour were derived from these TITAN optimal solutions. Default settings were used, except the bin size for WES data, which was set to 50 kb windows, whereas for WGS data, 10 kb windows were used.

Gene overlap of tumour CNAs

To construct a gene-level CNA matrix, a copy number status was assigned to each gene on the basis of its overlap with TITAN-derived copy number segments and gene coordinates from Ensembl BioMart (GRCh38.p12). To enable consistent comparisons across tumour samples, gene-level copy numbers were further normalized to account for sample-specific ploidy. First, the approximate ploidy was computed as the median copy number across all genes in the sample. Then, for each gene, its copy number was divided by this approximate ploidy to produce a normalized copy ratio (c_gene). CNAs were classified as gains if c_gene > 1.3, losses if c_gene < 0.66 and neutral if 0.66 ≤ c_gene ≤ 1.33. The resulting gene-level, ploidy-adjusted copy number matrix was used for downstream analyses. The same ploidy-based normalization and classification approach was applied to the validation cohorts (Hartwig Medical Foundation and The Cancer Genome Atlas (TCGA)) for copy number comparisons.

Founder, shared and private copy number clonality status

To classify the clonality of CNAs as being in founder, shared or private clones, we analysed gene-level CNAs across multiple tumours in each patient. Each gene was assigned either loss, neutral or gain status (as described above), and a founder, shared or private status based on whether the event was present with the same alteration status in all samples of the patient (founder), in multiple but not all samples (shared) or unique to a single sample (private).

Patient-level copy-number status and fraction of genome altered estimation

For each tumour, we divided the genome into 1 Mb bins and, for each bin, calculated the median copy number across all overlapping CNA calls. This size of binning enabled joint analyses of CNA results for WGS and WES, which were originally analysed using different bin sizes. Next, the copy number for each bin was divided by the median corrected copy number from TITAN across all 1 Mb bins in the sample, which produced a normalized copy ratio (c_bin). CNAs were classified as gains if c_bin > 1.3, losses if c_bin < 0.66 and neutral if 0.66 ≤ c_gene ≤ 1.33. The fraction of genome altered was calculated for each patient by counting the number of 1 Mb bins showing a copy number gain or loss and dividing the total number of all bins in that patient.

For all CNA-based analyses, six samples were excluded, including 17-030pM9 and 19-022pJ13-B, as well as the four samples already excluded for the mutation analyses, owing to high variability in copy number profiles.

Mutation clonality analysis

To assess cellular prevalence, defined as the proportion of tumour cells with a given somatic mutation, we used PyClone-VI⁸⁰ (v.0.1.1; https://github.com/Roth-Lab/pyclone-vi;commit ID: 6607ea1), a Bayesian clustering framework used to infer clonal population structure by grouping mutations with similar cellular prevalence. This analysis incorporated consensus SNV, indels, CNAs and matched tumour purity and ploidy estimates from TITAN. For each patient, we began by merging SNVs and indels from all available tumour samples in a patient into a unified mutation list. In cases when the primary tumour samples were FFPE, we applied an additional filtering step whereby we retained only the mutations that were also present in at least one metastatic tumour from the same patient. This step was done to mitigate the impact of sequencing artefacts and false positives in FFPE-derived data. As a result, mutations found exclusively in FFPE primary tumours were not considered in the mutation clonality analyses.

Following the generation of the unified mutation list, we extracted allelic read counts for each SNV using the following prioritization: (1) if an SNV was present in the tumour sample of the consensus somatic SNV set, then read counts were obtained from the corresponding variant caller in the following order of priority: Mutect2, Strelka2, MuSE and VarScan; or (2) if a SNV was not present in the tumour sample of the consensus set, allelic read counts were obtained using CollectAllelicCounts from GATK. For indels, read counts were extracted using bam-readcount (https://github.com/genome/bam-readcount, commit ID: c7c76e6). A final file was created to contain read counts for both SNVs and indels, as well as major and minor allele copy numbers across all tumours for a patient.

This final file was used as input for the fit function in PyClone-VI. The beta-binomial model was used with default parameters, except we specified a maximum of 7 clusters, 10,000 iterations and a burn-in of 1,000.

Phylogenetic tree reconstruction using LICHeE

Phylogenetic reconstruction was performed using LICHeE⁸¹ (v.1.0) (https://github.com/viq854/lichee; commitID:26c2a70) with the mutation clusters derived directly from PyClone-VI. For each patient, we generated two input files: a mutation file with SNV and indel presence or absent, and a cluster file with the cellular prevalence (C.P. from PyClone-VI) profiles across tumour samples. We ran LICHeE in cellular prevalence mode (-cp) with an error margin of 0.4 to retain all clusters. We specified the following parameters: -maxVAFAbsent 0.0, -minVAFPresent 0.5, -maxClusterDist 0.2, -minRobustNodeSupport 2, -minClusterSize 20, -minPrivateClusterSize 1 and -e 0.4. Before tree construction, we excluded clusters containing fewer than 20 mutations or with a cellular prevalence below 0.03 to omit small clonal clusters with less mutation support. For each patient, we selected the top-scoring phylogenetic tree based on the tree with the minimum sum of squared deviations.

Founder, shared and private clonal cluster definitions

The resulting trees represented subclonal evolutionary relationships^82,83. Clusters were classified on the basis of their average cellular prevalence and distribution across samples (Supplementary Table 3). A cluster was designated as founder if it had the highest average cellular prevalence (median = 0.99) in the patient and was consistently detected across all samples (minimum average cellular prevalence = 0.61). Shared clusters were present in more than one sample but had a lower average cellular prevalence than the founder cluster in that patient. Private clusters were restricted to a single sample and exhibited a non-zero cellular prevalence in that sample.

Visualization of clones

We visualized the trees using the LICHeE Java viewer and exported them as annotated DOT images and pdfs for downstream analyses. To visualize subclonal architecture in the context of phylogenetic structure and cellular prevalence, we used the R package cloneMap⁸⁴ to make octagon-shaped visualizations of each tumour (https://github.com/amf71/cloneMap; commitID: 7628d70). Input data included mutation clusters and lineage relationships inferred by LICHeE, along with cellular prevalence estimates from PyClone-VI. These were formatted according to cloneMap specifications, and the tool was executed using default parameters. Clone sizes were scaled according to sample-specific cellular prevalence values, whereas spatial positioning of clones reflected their inferred phylogenetic relationships. Final clone maps were exported as high-resolution vector graphics for integration into downstream figure panels, including those illustrating metastatic seeding patterns.

Metastatic seeding analysis using MACHINA

To infer metastatic seeding patterns and clonal migration histories, we applied MACHINA⁸⁵ (https://github.com/raphael-group/machina; commit ID: ac71282; v.1.2), a computational tool designed to reconstruct tumour seeding across primary and metastatic sites. MACHINA metastatic seeding analysis was performed on 16 out of of the 20 patients; 4 patients (17-047, 18-101, 18-120 and 19-044) were excluded from this analysis as their primary tumours either had estimated tumour purity <20% or poor-quality FFPE DNA sequencing data. The input data included mutation clusters and clonal assignments derived from LICHeE-generated trees, formatted according to the input requirements of MACHINA.

MACHINA was run using both the parsimonious migration history (pmh) and pmh with tumour reseeding (pmh_tr) models to account for single and multiple seeding events between tumour sites. Migration inference was performed under an unrestricted migration model that allowed the following seeding types: primary seeding (seeding from the primary tumour); metastatic seeding (seeding to other sites); metastasis-to-metastasis seeding; and reseeding.

For trees that included polytomies (that is, LICHeE clusters with three or more descendant clones), we favoured solutions from the pmh_tr model. To select the most parsimonious solution, we first identified those with the lowest total number of migration and comigration events. If multiple solutions had the same total, we prioritized the one with fewer migration events, assuming that migration is less frequent. If a tie remained, we selected the solution in which the most recent common ancestral node (in the LICHeE tree) had the highest cellular prevalence, thereby reflecting a more dominant seeding clone. Candidate clone trees were evaluated for each patient and, when available, CT scans were used to help confirm the selected solutions.

Migration trees and clonal relationships were visualized using Graphviz. MACHINA was run on a high-performance computing cluster using Gurobi (v.11.0) with a WLS licence.

Multivariate Cox regression models that included known prognostic factors of age, sex and histology were used to assess the relationship between polyclonality of seeding (the mean number of migrating clones across all metastatic seeding events in each patient) and patient survival from first metastasis to death. The time that the patient’s first metastasis was observed by CT scan was set to t₀ to reduce confounding from the variable interval of disease-free survival from diagnosis to first metastasis. For one patient who presented with metastases at diagnosis, we use diagnosis as the time of first metastasis, although there is uncertainty as to when these tumours may have first arisen in the patient. To compute the mean seeding polyclonality metric, we used the following equation. Let c_m,p be the number of clones required to seed metastasis m in patient p. Then, the mean number of clones per seeding event across M total metastases in the patient p was defined as \(\overline{{C}_{p}}=\left({\sum }_{m}^{M}{c}_{m,p}\right)/M\). Patient 16-070 was excluded from this analysis owing to being a statistical outlier in survival from first metastasis to death (z score = 3.5).

SV analysis

SV consensus calling

To establish a set of high-confidence consensus SVs, we used three tools, SvABA⁷⁷, Manta⁸⁶ and GRIDSS⁸⁷, to generate a consensus list. A SV event was retained for downstream analyses only if it was detected by at least two out of the three tools. A SV event identified by different callers were considered the same event if each of their two breakpoints overlapped in a 1 kbp window. This created a set of consensus SVs that were merged across callsets made by each tool. Next, SVs were classified as insertions, deletions, duplications, translocations or inversions on the basis of their breakpoint orientations. We further refined this list by retaining only SVs greater than 1 kbp in length, thereby focusing on larger, potentially more impactful genomic alterations and not accounting for small deletions.

SVs with breakpoints located in the ENCODE⁸⁸ blacklist (https://github.com/Boyle-Lab/Blacklist) were excluded. A gene was considered transected (broken) by an SV if at least one of its breakpoints overlapped the gene based on the coordinates from GENCODE release 44 (GRCh38.p14; basic gene annotation, CHR regions).

For all SV-based analyses, six samples were excluded, including 17-030pM9 and 19-022pJ13-B, as well as the four samples already excluded for the mutation analysis, owing to high variability in copy number profiles.

Complex SV analysis

Complex SVs were detected using two tools: Junction Balance Analysis (JaBbA)²² and Amplicon Architect (AA)⁸⁹. For JaBbA, the tool was installed from GitHub (v.1.1, https://github.com/mskilab-org/JaBbA; commitID: e3dc184). JaBbA was run on a high-performance computing cluster using a Gurobi WLS licence and Gurobi (v.11.0). The consensus SV data in BEDPE format were used as junction inputs. For genome-wide coverage input, we used the tumour-to-normal coverage ratio generated by ichorCNA as part of the TITAN pipeline (bin size of 10 kbp) and used the default algorithm CBS for tumour–normal copy number segmentation input. The purity and ploidy estimates were provided from TITAN-curated solutions (Supplementary Table 1). All other parameters were kept at their default settings.

Complex SVs such as breakage–fusion–bridge cycles, double minutes, rigma, pyrgo and templated insertion chains were called on junction-balanced genome graphs and rds files from JaBbA using the companion package gGnome (v.0.1, https://github.com/mskilab-org/gGnome; commitID: ee19bea). Any simple SVs such as inversions, translocations and duplications detected from JaBbA were not used for analyses. gGnome was used to visualize certain complex SV events and to save the rds files generated by JaBbA as.txt files for further downstream analyses. We further manually curated the presence of these events by looking into the chromosome-level plots (Fig. 2i and Supplementary Fig. 4) and confirmed the presence or absence of complex SVs.

AA was run using the end-to-end wrapper AmpliconSuite-pipeline (v.1.2.2, https://github.com/AmpliconSuite/AmpliconSuite-pipeline; commitID: 7ae7f32), which enabled the detection and classification of focal copy number amplifications such as ecDNA and breakage–fusion–bridge events from WGS data. All relevant libraries and repositories for the GRCh38 reference genome were used.

AA identified focal amplifications by first defining seed intervals as genomic regions larger than 50 kb with a copy number greater than 4.5 and using a downsampling factor of 20. For FFPE samples, a downsampling factor of 1 was used to mitigate FFPE-associated sequencing artefacts. The resulting breakpoint graph was decomposed into simple and complex cycles to detect potential circular DNA structures, including ecDNA. We classified tumours by the presence or absence of ecDNA or breakage–fusion–bridge or complex SVs. Any linear amplifications (chromosomal) detected were not considered for analyses. If both chromosomal and ecDNA were present, patients were included in the ecDNA group.

Annotation and merging JaBbA and AA results

SV footprints identified by JaBbA and amplicons detected by AA were merged and unique events were retained. Merging was performed on the basis of a reciprocal overlap criterion: SVs were considered concordant if they shared at least 75% overlap in breakpoint coordinates or if the overlapping region encompassed at least 50% of the genes reported by either tool. To consolidate SVs across multiple tumour samples from the same patient, we applied the same merging criteria—75% breakpoint overlap or 50% gene overlap—thereby generating patient-level SV event counts. Each SV event was annotated for overlap with known BLCA driver genes: an event was classified as a driver complex SV if any breakpoint intersected a driver gene.

SV clonality analysis

Clonality of simple SVs such as insertions, deletions, duplications, inversions and translocations was assessed using SVclone⁹⁰ (v.1.1.3, https://github.com/mcmero/SVclone; commit ID: 93dafb2), a tool designed to estimate the cancer cell fraction (CCF) of somatic SVs from WGS data. The pipeline was run using default parameters unless otherwise specified in a single sample mode. Input data included high-confidence consensus somatic SV calls in VCF format and matched tumour–normal BAM files.

SVclone first annotates and filters SVs, then calculates VAFs by counting supporting and non-supporting reads at each breakpoint. These counts were adjusted for copy number and tumour purity (from TITAN) to estimate the CCF of each SV. SVs were then clustered on the basis of their CCFs to infer subclonal architecture.

Only SVs passing all quality filters and with sufficient read support (a minimum of 1 split and 1 discordant read) were included in downstream SV clonality analyses. Clusters with CCFs > 0.90 were considered clonal, whereas those with lower CCFs were classified as subclonal. To further obtain a SV CCF status per patient (SVclone is for single sample runs) we used our custom script, which merges the different SVs in tumours of a patient to obtain a final list of SVs with CCF for each patient. For each patient, we merged overlapping SV calls across that patient’s tumour samples using a 500 bp window in BEDPE format to produce a unified per-patient SV list. Each SV was then classified as founder if it was clonal and detected in every tumour, shared if it occurred in two or more but not all tumours, or private if it was unique to a single tumour.

Visualization of somatic alterations (comut)

Somatic alterations for curated BLCA-driver genes were visualized using the Python package comut⁹¹ (v.0.0.3, https://github.com/vanallenlab/comut, commitID: 0ee1408). Mutations were shown on the basis of founder, shared or private clonal status. For each SNV and indel, mutations that did not qualify as being present were defined as one of the following: (1) below threshold (when a mutation is clustered by Pyclone-VI as being in a shared or founder cluster present in the tumour, but the sample itself does not have sufficient evidence for the mutation (it was filtered out based on the criteria listed in criteria for high-confidence consensus mutation filtering)); (2) no mutation (for any sample in a patient for whom the mutation was not found in the consensus SNV call set, despite having a read depth providing sufficient statistical power (≥80%)); or (3) no power (the depth at the mutation site did not provide sufficient power). The power was estimated for SNVs based on a previously used approach¹⁶ that incorporates the total read depth at each SNV locus, tumour purity and ploidy. SNVs with a power estimate below 80% were considered to have insufficient power.

Classifying BLCA pathogenic alterations

To determine whether a mutation is pathogenic, we used a consensus-based approach that incorporated multiple predictive algorithms. Specifically, if any 3 or more of the following 11 predictors classify a mutation as deleterious, it was labelled as pathogenic: SIFT, PolyPhen-2 HDIV, PolyPhen-2 HVAR, MutationTaster, MutationAssessor, FATHMM, PROVEAN, MetaLR, M-CAP, ClinVar and FATHMM-MKL. In addition to these computational predictions, we also considered other forms of supporting evidence: the mutation is listed in the COSMIC (v.70) database and occurs in one of the128 BLCA driver genes; it is a splicing mutation affecting BLCA-driver gene; or it is a frameshift indel in a BLCA-driver gene.

This approach improved the identification of potentially pathogenic mutations by combining multiple lines of evidence, including predictive algorithms and biological context. CNAs were considered pathogenic only if the event was a gain in an oncogene or a loss in a tumour suppressor from the list of BLCA-driver genes. Any CNA in a BLCA gene that could not be classified as clearly oncogene or tumour suppressor was not considered as pathogenic. We classified any SV intersecting a BLCA-driver gene as pathogenic, even if only a single breakpoint fell in the driver gene. For analyses of clinical associations between pathogenic alterations and survival (time from diagnosis to death or metastasis), patient 17-047 was excluded as an outlier owing to markedly longer survival from diagnosis to death (9.5 years, z score = 2.8). Patient 19-004 was also excluded from these analyses as the sole case with a UC-Sarc histological subtype.

Mutational signature analysis of patient phylogenetic trees

Phylogenetic trees were constructed for each patient using LICHeE (v.1.0, see above). SNVs assigned to each node in the resulting phylogeny were then analysed using SigProfiler suite of tools to estimate the SBS mutational signature composition of each cluster or clone in the tree.

Mutational signature analysis was performed using SigProfilerAssignment⁹² (v.0.1.1, https://github.com/AlexandrovLab/SigProfilerAssignment; commit ID: 327fb93) SigProfilerAssignment was run on individual mutation clusters (clones) identified by PyClone-VI, which resulted in a distinct set of SBS signatures for each clone in a patient.

To achieve this, SNVs from each cluster were formatted into VCF files using a five-column format: chromosome, genomic coordinate, sample ID, reference allele and alternative allele. All sample VCFs were then processed collectively using SigProfilerAssignment⁹² (v.1.1.23) to extract mutational signatures and to generate a matrix of mutation counts across the 96 SBS patterns. Signature extraction was conducted across all samples and matched to known COSMIC (v.3.2) SBS reference signatures. A final signature refitting step was performed using a curated subset of COSMIC SBS signatures. Non-relevant chemotherapy-associated signatures (for example, SBS11 (temozolomide), SBS86 (unknown chemotherapy), SBS87 (thiopurine chemotherapy) and SBS42 (haloalkane exposure)) and UV-related signatures (SBS7a–d and SBS38) were excluded, which resulted in a set of 51 COSMIC (v.3.2) SBS signatures. One of the seven patients who did not receive platinum-based chemotherapy displayed a minor proportion of platinum-associated signatures (SBS31 + SBS35), which is probably an artefact owing to the low mutation count in this mutation cluster. Therefore, for these seven patients, SBS31 and SBS35 were also excluded, and refitting was performed separately to avoid such artefacts, which generated a final set of 49 signatures used for downstream analyses in this group (Supplementary Table 5).

Mutation signatures were grouped into aetiologies for analyses as follows: APOBEC (SBS2 + SBS13); Clock-like (SBS1 + SBS5); platinum chemotherapy (SBS31 + SBS35); and DNA repair deficiencies (SBS3, SBS6, SBS14, SBS15, SBS20, SBS21, SBS26, SBS30, SBS36 and SBS44).

Bulk tumour RNA-seq analysis

Raw sequencing reads were assessed for quality using FastQC (v.0.12.1). Using STAR2 (v.2.7.3a), RNA-seq reads were aligned to the UCSC hg38 genome assembly and quantified for gene-level expression using HTSeq (v.0.11.1) against the gencode (v.22) gene annotation database. Variance stabilizing transformation (vst) expression values from DESeq2 (v.1.48.1) were used for all RNA-seq analyses. Dimension reduction was performed with the R package umap (v.0.2.10.0), using the top 1,000 most variably expressed genes.

FANC pathway analysis

The following genes were used to define FANC and cisplatin-related pathways⁹³: FANC: FANCD1, BRCA2, FANCD2, BRCA1, FANCI, RAD51, RAD51C, FANCA, FANCC, FANCE, FANCF and FANCG; cisplatin import: SLC31A1 and SLC31A2; cisplatin export: ATP7A, ATP7B, ABCC1, ABCC2, ABCC3 and ABCC5; GSH metabolism: GCLC, GCLM, GSTM4, GSTP1 and GSTT1; and nucleotide-excision repair: ERCC1, ERCC4, ERCC2, ERCC3, ERCC5, ERCC6 and ERCC8.

Pathway scores for each of these gene sets were calculated by first z-score-normalizing vst expression values for each gene across all high-quality flash-frozen bulk RNA-seq samples. These normalized expression values were then summed across all genes in the pathway of interest to generate a composite activity score for each sample. A LMM with the patient ID as a random effect (to account for intra-patient variability) was fit to assess the association between FANC pathway activity scores and histology using the R package lme4 (v.1.1-36). Cell line FANC scores were calculated in the same way as for bulk tumour RNA-seq samples using RNA-seq data downloaded from DepMap (release 23Q4)^94,95 (https://depmap.org/portal, ‘OmicsExpressionProteinCodingGenesTPMLogp1.csv’). snRNA-seq FANC scores were calculated per cell using the AddModuleScore function from Seurat (v.5.2.1).

Gene set variation analysis pathway scores

The gene set variation analysis (GSVA; v.1.52.3) package in R was used to create RNA expression scores for different pathways of interest. For the pathways identified as differentially expressed in snRNA-seq analyses, the top 25 ‘leading edge’ genes of each pathway from the fGSEA analyses were used to generate GSVA scores in each bulk RNA-seq tumour. The cell cycle proliferation pathway has been previously defined with the following genes⁹⁶: FOXM1, ASPM, TK1, PRC1, CDC20, BUB1B, PBK, DTL, CDKN3, RRM2, ASF1B, CEP55, CDK1, DLGAP5, SKA1, RAD51, KIF11, BIRC5, RAD54L, CENPM, PCLAF, KIF20A, PTTG1, CDCA8, NUSAP1, PLK1, CDCA3, ORC6, CENPF, TOP2A and MCM10.

Consensus molecular subtype classification

The consensusMIBC (v.1.1.0) package in R was used on each bulk RNA-seq tumour sample to determine its consensus molecular subtype as previously described⁸. For each tumour, the top scoring molecular subtype was used, except if the separationLevel was ≤0.2, in which case the molecular subtype was deemed ‘indeterminate’.

Deconvolution of cell types in bulk tumour RNA-seq data

The online CIBERSORTx⁹⁷ tool was used to deconvolute cell types from bulk RNA-seq data (log₂ TPM values). For deconvolution of immune cell types, an ‘Impute Cell Fractions’ job was run using the website’s standard ‘LM22.update-gene-symbols.txt’ immune gene signatures. For deconvolution of broad fibroblast, epithelial, immune and endothelial cell types, an Impute Cell Fractions job was run using previously described TR4 gene signatures (supplementary table 2L in ref. ⁹⁷). Both jobs were run without batch correction or quantile normalization in relative run mode with 500 and 1,000 permutations, respectively.

Hartwig Medical Foundation data analysis

Hartwig Medical Foundation data were accessed following approval of Data Access Request (DR-250). Data from 97 patients with metastatic BLCA⁵ and 133 patients with metastatic ovarian cancer²⁷ were included, comprising clinical metadata, somatic CNAs, mutation calls (VCF/TXT) and RNA-seq (TPM-normalized; TXT) data, all aligned to GRCh38. A subset of 52 patients with BLCA that were treated with cisplatin and had matched RNA-seq and WGS mutation data were used as a validation cohort for mutational signature and FANC analyses. For ovarian cancer, 115 patients had matched RNA-seq and tumour–normal WGS data, of whom 106 received carboplatin and 9 received cisplatin; these data were used as an additional validation cohort. Somatic mutation calls (VCF) were generated by Hartwig Medical Foundation using the SAGE caller, and we filtered the calls to retain PASS variants, excluding panel-of-normals, sequencing artefacts and common germline variants (gnomAD), followed by additional read-level filtering (see the section ‘Somatic mutation consensus calling’ above). Specifically, SNVs were required to have a tumour read depth of ≥14, ≥5 alternative reads and a tumour VAF ≥ 1%, with a matched normal read depth of ≥8 and VAF ≤ 0.05. SBS mutational signatures were assigned at the sample level using SigProfilerAssignment, applying the same reference signatures and versions as in the UW–Fred Hutch Rapid Autopsy programme (see the section ‘Mutational signature analysis of patient phylogenetic trees’ above). FANC pathway activity scores were computed using the same z score normalization approach as in the UW–Fred Hutch cohort. BLCA samples from the cohort described above were used to assess gene-level CNAs for genes of interest and compared with the TCGA cohort. Copy number data were normalized and thresholded using a previously described approach (see the section ‘Gene overlap of tumour CNAs’ above).

TCGA data analysis

Copy number data for BLCA were obtained from TCGA at the gene level via the Genomic Data Commons. Associated clinical and sample metadata were curated from previously published sources^20,65. The analysis was restricted to primary solid tumour samples (n = 344). Copy number profiles were processed in the GRCh38 (hg38) reference genome coordinate system and analysed using a ploidy-aware framework consistent with previously described approaches (see the section ‘Gene overlap of tumour CNAs’ above) for gene-level CNA normalization and thresholding.

cfDNA analysis

Overview and rationale for cfDNA sequencing

WGS (around 70× fragment coverage) and TP (about 3,000× fragment coverage) were performed on post-mortem plasma cfDNA samples. Each sequencing data type provides unique molecular features for assessment of cfDNA. WGS of cfDNA enabled whole-genome assessment of copy number. WGS also enabled the assessment of tumour phenotype and detection of tissue damage by using cfDNA nucleosome profiling analysis, which we and others have previously used^43,44,98,99. The breadth of WGS has been proposed as advantageous for ctDNA detection owing to the higher number of variants that can be observed genome-wide, thereby increasing the power to detect the presence of tumours¹⁰⁰. However, the depth of coverage for WGS data was limiting for samples with low ctDNA fraction (that is, TFx < 8%). Ultradeep TP sequencing provided the statistical power to detect mutations in samples with a low TFx (<8%), although it is limited to interrogating mutations and indels only in the targeted regions. For samples with a high TFx (≥8%), we had the opportunity to compare both WGS and TP to evaluate their ability to capture mutation clonality from the primary and metastatic tumours.

Duplex consensus calling of cfDNA TP sequencing

TP sequencing data were received as BAM files aligned to GRCh37. Duplex consensus calling was performed on these sequences using the fgbio tool suite (v.2.1.1, with Java v.11.0.2), and data were processed as follows: first, reads were grouped by the unique molecular identifier (stored in the RX tag of each read) using the GroupReadsByUmi function in fgbio, then consensus sequences were constructed using CallDuplexConsensusReads in fgbio, with a minimum read count of 1. Consensus reads were then filtered using FilterConsensusReads, with cutoffs of a minimum base quality of 20 and minimum number of reads of 2 being applied. Next, the unaligned consensus sequences (stored in a ubam) were aligned to GRCh38 using BWA MEM (v.0.7.17), after which alignments were sorted using samtools sort (v.1.10). Insert size metrics were collected using CollectInsertSizeMetrics in Picard (v.2.25.1). This above consensus called and re-aligned BAM files were used for all downstream analyses of the TP.

Detection of tumour SNVs in cfDNA

To determine the per cent capture of tumour SNVs in cfDNA sequencing, the number of reference and variant reads for each tumour SNV locus was assessed in the cfDNA using CollectAllelicCounts in GATK. If the cfDNA sequencing contained three or more variant reads matching the tumour mutation at the SNV locus, this SNV was considered to be captured in the cfDNA. For TP cfDNA sequencing, only tumour SNVs within 150 bp of the TP capture region were considered in calculations of the proportion of tumour SNVs captured.

For the comparison of observed to expected private mutations captured in cfDNA, the expected false positives were calculated on the basis of the probability of observing three or more variant reads matching each private tumour SNV based on sequencing error alone. The false positive probability for each tumour SNV was computed using a binomial error model as follows:

where n is the total number of cfDNA reads for the given locus and P = 0.001 is the sequencing error rate as published for the NovaSeq 6000 and NovaSeq X sequencers^101,102. The sum of such probabilities for all the private SNVs in the tumour was used as the total expected false positive value for that tumour.

Detection of tumour CNAs in cfDNA

CNAs were identified using TITAN for the eight cfDNA WGS data as described above for tumour WGS data, with the following details or changes. For four cfDNA samples with matched normal (tissue) WGS data, the standard tumour–normal paired pipeline was used. For the other four cfDNA samples with normal (tissue) WES data, the tumour-only pipeline was used because the normal WES sample was not suitable as a proper normal in the paired analysis. The full pipeline configurations are provided in GitHub (https://github.com/GavinHaLab/BLCA-subtype-evolution-paper).

To determine the per cent capture of tumour CNAs in cfDNA sequencing, we used the gene-level CNA matrix of tumours as a reference (see the section ‘Gene overlap of tumour CNAs’ above). A corresponding matrix was generated for cfDNA samples by intersecting gene coordinates with TITAN-derived segmented copy number regions. Genes with a copy number >2 were annotated as gains, and those with <2 as losses. As with the tumour analysis, cfDNA events were classified as founder, shared or private on the basis of their presence and directionality (gain or loss) across samples from the same patient. The per cent capture in cfDNA was then calculated separately for founder, shared and private status. For each status, it represents the proportion of tumour events (gain or loss) that were also detected in cfDNA with same directionality.

TFx estimation

The final TFx values used for cfDNA samples in this study were derived on a per-patient basis. As all sequencing samples (ULP WGS, deep WGS or TP) were derived from a single blood sample per patient, one TFx estimate was used for all samples of each patient. For patients with a high TFx (≥8%) and had both ULP and deep WGS, TITAN TFx estimates from the deep WGS were used. For patients with a low TFx and only ULP WGS, ichorCNA TFx estimates were used.

Griffin analysis and nucleosome profiling

Griffin⁹⁹ (v0.2.0) is a computational tool for profiling nucleosome protection and chromatin accessibility at genomic sites. For this study, Griffin was installed from GitHub (https://github.com/GavinHaLab/Griffin; commitID: b624c7a) and analysis was performed following the guidelines provided in the Griffin documentation.

After GC-bias correction, nucleosome profiling was performed on a per-sample basis. Only nucleosome-sized fragments (100–200 bp) were retained for downstream analyses. The human genome reference used for all analyses was GRCh38. Two features were derived from normalized nucleosome coverage profiles by aggregating sites (n = 10,000) for each entity (for example, transcription factor binding sites or tissue-specific accessible chromatin): (1) mean central coverage (−30 to +30 bp) and (2) mean window coverage (−990 to +990 bp) relative to the site centre.

All the figures for nucleosome profiling were based on mean central coverage. The TFBSs were taken from the gene transcription regulation database, and additional details on the source of the ChIP–seq (chromatin immunoprecipitation followed by sequencing) sites and their filtering were adopted from previously published studies^44,99. cfDNA from healthy donors (n = 6) was processed using the same Griffin pipeline as that applied to WGS cfDNA from eight patients.

Tissue of origin analysis

Griffin analysis was performed to determine tissues contributing to each cfDNA sample (see Griffin details above). Tissue-specific sites were obtained from the ENCODE regulatory index and CATtlas (https://catlas.org/humanenhancer/) single-cell ATLAS of ATAC-seq (assay for transposase-accessible chromatin with sequencing) data¹⁰³. The top 10,000 peaks were selected by peak scores. These sites were used to generate the Griffin composite site profile to extract mean central coverage and mean window coverage.

BLCA-specific ATAC-seq sites

Griffin analysis was performed to determine BLCA active sites contributing to each cfDNA sample (see Griffin details above). BLCA-specific chromatin accessibility sites were obtained from a published ATAC-seq dataset¹⁰⁴. To focus on BLCA-relevant regulatory regions, we applied a filtering strategy consistent with previous work^44,105, removing haematopoietic and ubiquitously accessible sites. From the remaining set, the top 10,000 peaks were selected on the basis of peak scores. These curated sites were then used to generate Griffin composite site profiles, from which mean central coverage and mean window coverage were extracted.

snRNA-seq analysis

Filtering and clustering

CellRanger (v.7.1.0 from 10× Genomics) was used to align, quantify and provide basic quality-control metrics for the snRNA-seq data. Cells with fewer than 1,000 detected reads, greater than 50,000 detected reads or greater than 20% mitochondrial reads were filtered from subsequent analyses. Using the standard Seurat (v.5.2.1) pipeline, the snRNA-seq data were normalized (NormalizeData) and clustered (FindVariableFeatures, FindNeighbors with dims = 1:30, and FindClusters with resolution = 2). CCAIntegration was used to batch correct between the two sequencing batches of samples, and preliminary cell types were determined by applying the Human Primary Cell Atlas Data reference (celldex::HumanPrimaryCellAtlasData()) to the dataset using SingleR (v.2.6.0). Doublets were removed separately for each sequencing batch (re-clustered each using dims = 1:20 and resolution = 1) using DoubletFinder (v.2.0.4). DoubletFinder parameters for batch 1 samples were as follows: PCs = 1:23, pN = 0.25, pK = 0.01, nExp = 1131. Parameters for batch 2 samples were as follows: PCs = 1:23, pN = 0.25, pK = 0.005, nExp = 723. After doublets were removed, the fully filtered dataset was re-normalized, re-clustered and re-batch-corrected as above with dims = 1:30 and resolution = 1.75. Seurat’s UMAP dimension reduction technique (dims = 1:30) was used for visualization.

Cell-type annotation

The main seven cell types were annotated using a combination of the Human Primary Cell Atlas¹⁰⁶ re-applied to the finalized, filtered dataset as above and manual validation with cell-type marker genes. The HPCA dataset (using ‘label.main’) was subsetted to relevant cell types (B cells, dendritic cells, endothelial cells, epithelial cells, erythroblasts, fibroblasts, hepatocytes, HSC-G-CSF cells, keratinocytes, macrophages, monocytes, mesenchymal stem cells, neuroepithelial cells, neutrophils, natural killer cells, platelets and T cells) for this application, and the most commonly called cell type in each cell cluster was applied to all the cells in that cluster. Because marker genes were not able to distinguish between clusters called monocyte versus macrophage in this manner, both were re-labelled as ‘myeloid’.

The specific fibroblast cell types were annotated using supplementary table 3 in ref. ¹⁰⁷ as a reference dataset. The top 50 differentially expressed genes (DEGs) for each fibroblast type in this dataset were used to create a cell-specific score for each fibroblast type (AddModuleScore). After preliminary calling of fibroblast types using the max score for each cell, we removed all calls for fibroblast types with fewer than 100 total cells called (apCAF, dCAF, hsp tCAF, tCAF and ifnCAF) and finalized CAF annotations as the top scoring fibroblast type across iCAF, vCAF, rCAF, mCAF and pericyte.

Fine-grain immune cell types were annotated primarily using manual marker gene assessment (Extended Data Fig. 9d), with annotations from HPCA, the Tumour Immune Cell Atlas (TICA)¹⁰⁸ and top DEGs for each immune cell cluster (FindAllMarkers) as a guide. HPCA predictions were re-done as described above using the ‘label.fine’ annotations. TICA predictions were performed using the Annotation function of the web app (https://singlecellgenomics-cnag-crg.shinyapps.io/TICA/) with input of the top 100 DEGs (FindAllMarkers) from each cluster of the subsetted immune dataset. One immune cell type was chosen for each immune cell cluster as the final annotation.

Epithelial, immune and fibroblast subsets

Epithelial, immune and fibroblast cells were subset from the filtered dataset using the seven main cell-type annotations and re-processed independently. Each subset underwent the same normalization, batch correction, clustering and UMAP workflow described for the full dataset above, with the following parameters: epithelial cells (‘Epithelial cells’) used dims = 1:20, resolution = 1.25; immune cells (‘B cell’, ‘T cells’ and ‘myeloid’) used dims = 1:28, resolution = 1.5; fibroblasts (‘fibroblasts’) used dims = 1:25, resolution = 1.5.

Differential gene expression and GSEA of histological subtypes

DEGs for PUC and UCSD tumours were identified using FindMarkers on the epithelial-cell-specific dataset. PUC DEGs were called between the group of samples 19_007K3, 19_007D5 and 18_120L1 versus all other samples. UCSD DEGs were called between the group of samples 19-022_H1 and 17-026_K2, 17-026_J3 versus all other samples. Gene set enrichment analysis (GSEA) was performed on these DEG lists using fgsea (v1.30.0), with the MSigDB Hallmark gene sets and log₂[fold change] × –log₁₀[adjusted P] × |pct.1 – pct.2| as the ranking metric, where pct.1 and pct.2 denote the proportions of cells with detectable expression of each gene in the first and second groups of the differential-expression comparison, respectively.

Calling consensus molecular subtypes in snRNA-seq

Per-cell and per-cluster consensus molecular subtypes were called using the R package consensusMIBC (v.1.1.0). For cell-specific subtyping, getConsensusClass was applied to the log-transformed normalized expression matrix of the epithelial-cell-specific dataset. For cell-cluster subtyping, each cell cluster in the epithelial-cell-specific dataset was first pseudobulked into normalized expression profiles for the cluster, and getConsensusClass was applied to each cluster profile.

Shannon entropy intra-patient heterogeneity index

Intra-patient heterogeneity scores were computed by first grouping cells from the epithelial-specific dataset by patient and cell cluster identity to obtain the proportion of each cluster in each patient. The Shannon entropy was then calculated per patient using these proportions with the entropy() function of the entropy package (v.1.3.1) in R. Intra-tumoural heterogeneity scores were computed in the same fashion by grouping instead by individual sample and cell-cluster identity.

NicheNet cell–cell signalling analysis

Cell–cell signalling analysis was performed using nichenetr¹⁰⁹ (v.2.1.7) in R using the sender-focused approach for human ligand–receptor–target networks as described in the GitHub seurat_steps.md vignette (https://github.com/saeyslab/nichenetr/blob/master/vignettes/seurat_steps.md). Epithelial cells and fibroblasts were used as the sender cell types, and T cells were set as the receiver cell type, using the full, filtered snRNA-seq Seurat object. The DEG set of interest was defined from the comparison between UC and PUC T cells using FindMarkers. Predicting ligand activities in this manner nominated TGFB1 as the ligand with the highest area under the precision–recall curve. The receptors and targets of TGFB1 were then identified using NicheNet’s interaction network, and the expression of these in UC versus PUC T cells was assessed using dot plots.

Numbat CNA analysis

The R package Numbat was used to identify CNAs from snRNA-seq data on both clone and cell-specific levels¹¹⁰. This analysis was performed according to the user guide on the author’s GitHub page (https://kharchenkolab.github.io/numbat/articles/numbat.html). First, the data were prepared using the given pileup_and_phase.R script for each patient (for patients with multiple tumours, all BAM and barcode files were run in one merged script).

CNA identification was performed for each patient. For patients with multiple tumours in snRNA-seq, the multiple tumours were combined, except for 19-022. Because the tumours from patient 19-022 were so dissimilar, they were run separately through the CNA identification step of Numbat. For each patient, the expression matrix was created from the finalized Seurat object, subsetted to the patient of interest, using GetAssayData for the counts layer. The reference expression dataset for each patient was created from that patient’s T cells, B cells, myeloid cells, endothelial cells and hepatocytes.

Bulk DNA sequencing (WES or WGS) CNA profiles were input into the Numbat run to inform CNA identification. The TITAN segment files (titan.ichor.seg.txt) were first reformatted for Numbat input, with cnv_state defined by the Corrected_MajorCN and Corrected_MinorCN (neutral: major=1 and minor=1; bdel: major=0 and minor=0; del: major+minor=1; loh: major≥1 and minor=0; amp: major+minor>2; bamp: major>1 and minor>1). For patients with only one tumour in snRNA-seq, the TITAN calls for the matching tumour in bulk DNA sequencing were used. For patients with multiple tumours in snRNA-seq, a combined CNA profile of these tumours was created. All segment breakpoints across all tumours were used, which created a union set of segments. For each of these segments, if all samples agreed on the copy number state, that state was used. In cases of disagreement, if one sample indicated an amplification (amp or bamp) and the other a deletion (del or bdel), the segment was assigned neutral to reflect this conflict. Otherwise, the merged state was determined by a priority hierarchy: amp > loh > del > bamp > bdel > neu, whereby the highest-priority state among the conflicting calls was used. Patient 15-109G5 did not have a matching bulk DNA sequencing copy number profile; therefore, a combination of CNAs from other tumours (15-109I10, 15-109M2 and 15-109N1) of this patient was used as a reference profile. Moreover, the ncut = 6 parameter was used to limit the number of clones identified for patient 15-109G5.

For determination of tumour versus normal from Numbat results, the clone_post_2.tsv output compartment_opt column was used. For identifying clonality and cellular prevalence of driver gene CNAs, copy number states for each segment were pulled from the bulk_clones_final.tsv.gz output file. For each driver gene of interest, the segment overlapping the gene was identified and the CNA calls for each tumour clone were pulled from the cnv_state_post column. CNAs were clonal if all tumour clones in a patient had the same CNA state. Cellular prevalence of subclonal CNAs was calculated as the proportion of tumour cells in a patient having a certain CNA state out of all the tumour cells of that patient.

Graphing and statistics

A combination of R (v.4.4.3) and Python (v.3.19, pandas v.1.40) was used for scripting, graphs and statistics. In R, ggplot2 (v.3.5.1) was used for graphing, rstatix (v.0.7.2) for t-tests and Mann–Whitney U-tests, and lme4 (v.1.1-36) for LMMs. For the snakemake pipelines, Python (v.3.7.4) was used. All statistical tests were two-sided unless otherwise specified, with the specific statistical tests used for each analysis detailed in the respective figure legends. All box plots display the median (horizontal line) and IQR (box bounds are 25th and 75th percentiles). Whiskers extend to the most extreme data point within 1.5× the IQR of the box; individual data points are overlaid.

LMMs were applied for repeated measures in patients (for example, multiple tumours in a patient), incorporating patient ID as a random effect. For intra-patient comparisons (for example, proportions of founder, shared and private alterations), we used paired Wilcoxon signed-rank tests. Paired t-tests were used for paired cell line experiments (for example, siNTC versus siFANCF). All linear fits shown in the figures were obtained using simple linear least-squares regression.

Survival analyses were performed using the Cox proportional hazards model implemented in the R package survival (v.0.5.0), with survival time defined as the interval from metastasis to death. HR values and corresponding P values were derived from the fitted models, and results are presented as forest plots generated with the forestmodel package. For risk stratification, patient-specific risk scores were calculated as the linear predictor from the Cox model and dichotomized at the median into high-risk and low-risk groups. Kaplan–Meier survival curves were generated using the survminer package.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

WES and WGS data are accessible at dbGAP under accession phs001797.v2.p1. Bulk RNA-seq and snRNA-seq data are accessible from the Gene Expression Omnibus under accession GSE302273. The Hartwig Medical Foundation WGS data and mutation call set were obtained via controlled access, available only with permission from the Hartwig Medical Foundation data access request (https://www.hartwigmedicalfoundation.nl/en/data/data-access-request/). Publicly available BLCA datasets from TCGA (https://portal.gdc.cancer.gov/projects/TCGA-BLCA) were accessed through the Genomic Data Commons Data Portal and analysed in this study. Publicly available cisplatin sensitivity and transcriptomic data for BLCA cell lines analysed as part of this study were obtained from the Cancer Dependency Map (DepMap release 23Q4; https://doi.org/10.25452/figshare.plus.24667905.v2)⁹⁵, including the Genomics of Drug Sensitivity in Cancer (GDSC2) dataset (https://depmap.sanger.ac.uk/documentation/datasets/drug-sensitivity/). All other processed data necessary to reproduce the results are found in the Supplementary Tables and at GitHub (https://github.com/GavinHaLab/BLCA-subtype-evolution-paper). All analyses were performed using the human reference genome GRCh38 (Broad version; https://console.cloud.google.com/storage/browser/gcp-public-data–broad-references/hg38/v0). Source data are provided with this paper.

Software and pipeline configurations and custom code written for this study are accessible at GitHub (https://github.com/GavinHaLab/BLCA-subtype-evolution-paper).

1. Black, A. J. & Black, P. C. Variant histology in bladder cancer: diagnostic and clinical implications. Transl. Cancer Res. 9, 6565–6575 (2020).

   Article  PubMed  PubMed Central  Google Scholar 

2. Chalasani, V., Chin, J. L. & Izawa, J. I. Histologic variants of urothelial bladder cancer and nonurothelial histology in bladder cancer. Can. Urol. Assoc. J. 3, S193 (2009).

   Article  PubMed  PubMed Central  Google Scholar 

3. Alderson, M., Grivas, P., Milowsky, M. I. & Wobker, S. E. Histologic variants of urothelial carcinoma: morphology, molecular features and clinical implications. Bladder Cancer 6, 107–122 (2020).

   Article  Google Scholar 

4. WHO Classification of Tumours Editorial Board. Urinary and Male Genital Tumours 5th edn, Vol. 8 (IARC Press, 2022).

5. Nakauma-González, J. A. et al. Comprehensive molecular characterization reveals genomic and transcriptomic subtypes of metastatic urothelial carcinoma. Eur. Urol. 81, 331–336 (2022).

   Article  PubMed  Google Scholar 

6. Nguyen, D. D. et al. The interplay of mutagenesis and ecDNA shapes urothelial cancer evolution. Nature 635, 219–228 (2024).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

7. Warrick, J. I. et al. Intratumoral heterogeneity of bladder cancer by molecular subtypes and histologic variants. Eur. Urol. 75, 18–22 (2019).

   Article  CAS  PubMed  Google Scholar 

8. Kamoun, A. et al. A consensus molecular classification of muscle-invasive bladder cancer. Eur. Urol. 77, 420–433 (2020).

   Article  PubMed  Google Scholar 

9. Teo, M. Y. et al. Natural history, response to systemic therapy, and genomic landscape of plasmacytoid urothelial carcinoma. Br. J. Cancer 124, 1214–1221 (2021).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

10. Taber, A. et al. Molecular correlates of cisplatin-based chemotherapy response in muscle invasive bladder cancer by integrated multi-omics analysis. Nat. Commun. 11, 4858 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

11. Marusyk, A., Janiszewska, M. & Polyak, K. Intratumor heterogeneity: the Rosetta stone of therapy resistance. Cancer Cell 37, 471–484 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

12. Beckabir, W. et al. Immune features are associated with response to neoadjuvant chemo-immunotherapy for muscle-invasive bladder cancer. Nat. Commun. 15, 4448 (2024).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

13. Bakaloudi, D. R. et al. Response and survival with immune checkpoint inhibitor in patients with advanced urothelial carcinoma and histology subtypes. Clin. Genitourin. Cancer 23, 102356 (2025).

    Article  PubMed  Google Scholar 

14. Koehne, E. L. et al. Adjuvant chemotherapy and survival after radical cystectomy in histologic subtype bladder cancer. Clin. Genitourin. Cancer 22, 102100 (2024).

    Article  PubMed  Google Scholar 

15. Yang, H. et al. Bladder cancer variants share aggressive features including a CA125⁺ cell state and targetable TM4SF1 expression. Nat. Commun. 16, 5312 (2025).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

16. Winters, B. R. et al. Genomic distinctions between metastatic lower and upper tract urothelial carcinoma revealed through rapid autopsy. JCI Insight 4, e128728 (2019).

    Article  Google Scholar 

17. Lee, H. J. et al. Patterns of HER2 expression in metastatic prostate and urothelial cancers: implications for HER2-targeted therapies. Cancer Res. Commun 5, 1419–1428 (2025).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

18. Ghali, F. et al. Metastatic bladder cancer expression and subcellular localization of Nectin-4 and Trop-2 in variant histology: a rapid autopsy study. Clin. Genitourin. Cancer 21, 669–678 (2023).

    Article  PubMed  PubMed Central  Google Scholar 

19. Zhao, K. et al. Longitudinal and multisite sampling reveals mutational and copy number evolution in tumors during metastatic dissemination. Nat. Genet. 57, 1504–1511 (2025).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

20. Robertson, A. G. et al. Comprehensive molecular characterization of muscle-invasive bladder cancer. Cell 171, 540–556 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

21. Rebouissou, S. et al. CDKN2A homozygous deletion is associated with muscle invasion in FGFR3-mutated urothelial bladder carcinoma. J. Pathol. 227, 315–324 (2012).

    Article  CAS  PubMed  Google Scholar 

22. Hadi, K. et al. Distinct classes of complex structural variation uncovered across thousands of cancer genome graphs. Cell 183, 197–210 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

23. Wang, L. et al. A genetically defined disease model reveals that urothelial cells can initiate divergent bladder cancer phenotypes. Proc. Natl Acad. Sci. USA 117, 563–572 (2020).

    Article  ADS  CAS  PubMed  Google Scholar 

24. Al-Ahmadie, H. A. et al. Frequent somatic CDH1 loss-of-function mutations in plasmacytoid-variant bladder cancer. Nat. Genet. 48, 356–358 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

25. Galluzzi, L. et al. Molecular mechanisms of cisplatin resistance. Oncogene 31, 1869–1883 (2012).

    Article  CAS  PubMed  Google Scholar 

26. Ceccaldi, R., Sarangi, P. & D’Andrea, A. D. The Fanconi anaemia pathway: new players and new functions. Nat. Rev. Mol. Cell Biol. 17, 337–349 (2016).

    Article  CAS  PubMed  Google Scholar 

27. de Witte, C. J. et al. Distinct genomic profiles are associated with treatment response and survival in ovarian cancer. Cancers 14, 1511 (2022).

    Article  PubMed  PubMed Central  Google Scholar 

28. Yang, W. et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 41, D955–D961 (2013).

    Article  CAS  PubMed  Google Scholar 

29. van Twest, S. et al. Mechanism of ubiquitination and deubiquitination in the Fanconi anemia pathway. Mol. Cell 65, 247–259 (2017).

    Article  PubMed  Google Scholar 

30. Taniguchi, T. et al. Disruption of the Fanconi anemia–BRCA pathway in cisplatin-sensitive ovarian tumors. Nat. Med. 9, 568–574 (2003).

    Article  CAS  PubMed  Google Scholar 

31. Chen, Q., Van der Sluis, P. C., Boulware, D., Hazlehurst, L. A. & Dalton, W. S. The FA/BRCA pathway is involved in melphalan-induced DNA interstrand cross-link repair and accounts for melphalan resistance in multiple myeloma cells. Blood 106, 698–705 (2005).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

32. Diamantopoulos, L. N. et al. Plasmacytoid urothelial carcinoma: response to chemotherapy and oncologic outcomes. Bladder Cancer 6, 71–81 (2020).

    Article  PubMed  PubMed Central  Google Scholar 

33. Kim, B. et al. HER2 protein overexpression and gene amplification in plasmacytoid urothelial carcinoma of the urinary bladder. Dis. Markers 2016, 8463731 (2016).

    Article  PubMed  PubMed Central  Google Scholar 

34. Lyu, A. et al. Evolution of myeloid-mediated immunotherapy resistance in prostate cancer. Nature 637, 1207–1217 (2025).

    Article  ADS  CAS  PubMed  Google Scholar 

35. Thomas, D. A. & Massagué, J. TGF-β directly targets cytotoxic T cell functions during tumor evasion of immune surveillance. Cancer Cell 8, 369–380 (2005).

    Article  CAS  PubMed  Google Scholar 

36. Forsthuber, A. et al. Cancer-associated fibroblast subtypes modulate the tumor-immune microenvironment and are associated with skin cancer malignancy. Nat. Commun. 15, 9678 (2024).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

37. Qi, J. et al. Single-cell and spatial analysis reveal interaction of FAP⁺ fibroblasts and SPP1⁺ macrophages in colorectal cancer. Nat. Commun. 13, 1742 (2022).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

38. Link, A. et al. Fibroblastic reticular cells in lymph nodes regulate the homeostasis of naive T cells. Nat. Immunol. 8, 1255–1265 (2007).

    Article  CAS  PubMed  Google Scholar 

39. Pereira, B. et al. Cell-free DNA captures tumor heterogeneity and driver alterations in rapid autopsies with pre-treated metastatic cancer. Nat. Commun. 12, 3199 (2021).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

40. Takai, E. et al. Post-mortem plasma cell-free DNA sequencing: proof-of-concept study for the “liquid autopsy”. Sci. Rep. 10, 2120 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

41. Abbosh, C. et al. Phylogenetic ctDNA analysis depicts early stage lung cancer evolution. Nature 545, 446–451 (2017).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

42. Herberts, C. et al. Deep whole-genome ctDNA chronology of treatment-resistant prostate cancer. Nature 608, 199–208 (2022).

    Article  ADS  CAS  PubMed  Google Scholar 

43. Ulz, P. et al. Inference of transcription factor binding from cell-free DNA enables tumor subtype prediction and early detection. Nat. Commun. 10, 4666 (2019).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

44. De Sarkar, N. et al. Nucleosome patterns in circulating tumor DNA reveal transcriptional regulation of advanced prostate cancer phenotypes. Cancer Discov. 13, 632–653 (2023).

    Article  PubMed  PubMed Central  Google Scholar 

45. Nouruzi, S. et al. ASCL1 activates neuronal stem cell-like lineage programming through remodeling of the chromatin landscape in prostate cancer. Nat. Commun. 13, 2282 (2022).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

46. Maraver, A. et al. NOTCH pathway inactivation promotes bladder cancer progression. J. Clin. Invest. 125, 824–830 (2015).

    Article  PubMed  PubMed Central  Google Scholar 

47. Pegoraro, S. et al. HMGA1 promotes metastatic processes in basal-like breast cancer regulating EMT and stemness. Oncotarget 4, 1293–1308 (2013).

    Article  PubMed  PubMed Central  Google Scholar 

48. Reddy, J. et al. Predicting master transcription factors from pan-cancer expression data. Sci. Adv. 7, eabf6123 (2021).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

49. Koyyalagunta, D., Ganesh, K. & Morris, Q. Inferring cancer type-specific patterns of metastatic spread using Metient. Nat. Methods 23, 574–584 (2026).

    Article  CAS  PubMed  Google Scholar 

50. Nijman, S. M. B. et al. The deubiquitinating enzyme USP1 regulates the Fanconi anemia pathway. Mol. Cell 17, 331–339 (2005).

    Article  CAS  PubMed  Google Scholar 

51. Liang, Q. et al. A selective USP1–UAF1 inhibitor links deubiquitination to DNA damage responses. Nat. Chem. Biol. 10, 298–304 (2014).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

52. Kossaï, M. et al. Plasmacytoid urothelial carcinoma (UC) are luminal tumors with similar CD8⁺ T cell density and PD-L1 protein expression on immune cells as compared to conventional UC. Urol. Oncol. 40, 12.e1–12.e11 (2022).

    Article  PubMed  Google Scholar 

53. Sfakianos, J. P. et al. Epithelial plasticity can generate multi-lineage phenotypes in human and murine bladder cancers. Nat. Commun. 11, 2540 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

54. Jindal, T. et al. Impact of squamous histology on clinical outcomes and molecular profiling in metastatic urothelial carcinoma patients treated with immune checkpoint inhibitors or enfortumab vedotin. Clin. Genitourin. Cancer 21, e394–e404 (2023).

    Article  PubMed  Google Scholar 

55. Cristescu, R. et al. Transcriptomic determinants of response to pembrolizumab monotherapy across solid tumor types. Clin. Cancer Res. 28, 1680–1689 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

56. Eraky, A. et al. Bladder cancer with undetectable circulating tumor DNA after radical cystectomy may be amenable to a less intense imaging surveillance protocol: a diagnostic accuracy study. Eur. Urol. Oncol. 8, 1505–1512 (2025).

    Article  PubMed  Google Scholar 

57. Tran, H. T. et al. Circulating tumor DNA and radiological tumor volume identify patients at risk for relapse with resected, early-stage non-small-cell lung cancer. Ann. Oncol. 35, 183–189 (2024).

    Article  CAS  PubMed  Google Scholar 

58. Franceschini, G. M. et al. Noninvasive detection of neuroendocrine prostate cancer through targeted cell-free DNA methylation. Cancer Discov. 14, 424–445 https://doi.org/10.1158/2159-8290.CD-23-0754 (2024).

59. Hiatt, J. B. et al. Molecular phenotyping of small cell lung cancer using targeted cfDNA profiling of transcriptional regulatory regions. Sci. Adv. 10, eadk2082 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

60. Peneder, P. et al. Multimodal analysis of cell-free DNA whole-genome sequencing for pediatric cancers with low mutational burden. Nat. Commun. 12, 3230 (2021).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

61. Jana, S. et al. Transcriptional–translational conflict is a barrier to cellular transformation and cancer progression. Cancer Cell 41, 853–870 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

62. Adalsteinsson, V. A. et al. Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors. Nat. Commun. 8, 1324 (2017).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

63. Richardson, C. D. et al. CRISPR–Cas9 genome editing in human cells occurs via the Fanconi anemia pathway. Nat. Genet. 50, 1132–1139 (2018).

    Article  CAS  PubMed  Google Scholar 

64. Gonzalez-Perez, A. et al. IntOGen-mutations identifies cancer drivers across tumor types. Nat. Methods 10, 1081–1082 (2013).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

65. Weinstein, J. N. et al. Comprehensive molecular characterization of urothelial bladder carcinoma. Nature 507, 315–322 (2014).

    Article  ADS  CAS  Google Scholar 

66. Choi, W. et al. Identification of distinct basal and luminal subtypes of muscle-invasive bladder cancer with different sensitivities to frontline chemotherapy. Cancer Cell 25, 152–165 (2014).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

67. Majumdar, S. et al. Loss of Sh3gl2/Endophilin A1 is a common event in urothelial carcinoma that promotes malignant behavior. Neoplasia 15, 749–760 (2013).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

68. Yun, S. J. & Kim, W.-J. Role of the epithelial–mesenchymal transition in bladder cancer: from prognosis to therapeutic target. Korean J. Urol. 54, 645–650 (2013).

    Article  PubMed  PubMed Central  Google Scholar 

69. Sondka, Z. et al. The COSMIC Cancer Gene Census: describing genetic dysfunction across all human cancers. Nat. Rev. Cancer 18, 696–705 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

70. Chakravarty, D. et al. OncoKB: a precision oncology knowledge base. JCO Precis. Oncol. 1, PO.17.00011 (2017).

    Google Scholar 

71. Benjamin, D. et al. Calling somatic SNVs and Indels with Mutect2. Preprint at bioRxiv https://doi.org/10.1101/861054 (2019).

72. Kim, S. et al. Strelka2: fast and accurate calling of germline and somatic variants. Nat. Methods 15, 591–594 (2018).

    Article  CAS  PubMed  Google Scholar 

73. Koboldt, D. C. et al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 22, 568–576 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

74. Ji, S., Zhu, T., Sethia, A. & Wang, W. Accelerated somatic mutation calling for whole-genome and whole-exome sequencing data from heterogenous tumor samples. Genome Res. 34, 633–641 (2024).

    CAS  PubMed  PubMed Central  Google Scholar 

75. Van der Auwera, G. A. et al. From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline. Curr. Protoc. Bioinformatics 11, 11.10.1–11.10.33 (2013).

    Google Scholar 

76. Wang, K., Li, M. & Hakonarson, H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res 38, e164 (2010).

    Article  PubMed  PubMed Central  Google Scholar 

77. Wala, J. A. et al. SvABA: genome-wide detection of structural variants and indels by local assembly. Genome Res. 28, 581–591 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

78. Martincorena, I. et al. Universal patterns of selection in cancer and somatic tissues. Cell 171, 1029–1041 (2017).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

79. Ha, G. et al. TITAN: inference of copy number architectures in clonal cell populations from tumor whole-genome sequence data. Genome Res. 24, 1881–1893 (2014).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

80. Gillis, S. & Roth, A. PyClone-VI: scalable inference of clonal population structures using whole genome data. BMC Bioinformatics 21, 571 (2020).

    Article  PubMed  PubMed Central  Google Scholar 

81. Popic, V. et al. Fast and scalable inference of multi-sample cancer lineages. Genome Biol. 16, 91 (2015).

    Article  PubMed  PubMed Central  Google Scholar 

82. Al Bakir, M. et al. The evolution of non-small cell lung cancer metastases in TRACERx. Nature 616, 534–542 (2023).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

83. Birkbak, N. J. & McGranahan, N. Cancer genome evolutionary trajectories in metastasis. Cancer Cell 37, 8–19 (2020).

    Article  CAS  PubMed  Google Scholar 

84. Frankell, A. M., Colliver, E., Mcgranahan, N. & Swanton, C. cloneMap: a R package to visualise clonal heterogeneity. Preprint at bioRxiv https://doi.org/10.1101/2022.07.26.501523 (2022).

85. El-Kebir, M., Satas, G. & Raphael, B. J. Inferring parsimonious migration histories for metastatic cancers. Nat. Genet. 50, 718–726 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

86. Chen, X. et al. Manta: rapid detection of structural variants and indels for germline and cancer sequencing applications. Bioinformatics 32, 1220–1222 (2016).

    Article  CAS  PubMed  Google Scholar 

87. Cameron, D. L. et al. GRIDSS: sensitive and specific genomic rearrangement detection using positional de Bruijn graph assembly. Genome Res. 27, 2050–2060 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

88. Amemiya, H. M., Kundaje, A. & Boyle, A. P. The ENCODE Blacklist: identification of problematic regions of the genome. Sci. Rep. 9, 9354 (2019).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

89. Deshpande, V. et al. Exploring the landscape of focal amplifications in cancer using AmpliconArchitect. Nat. Commun. 10, 392 (2019).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

90. Cmero, M. et al. Inferring structural variant cancer cell fraction. Nat. Commun. 11, 730 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

91. Crowdis, J., He, M. X., Reardon, B. & Van Allen, E. M. CoMut: visualizing integrated molecular information with comutation plots. Bioinformatics 36, 4348–4349 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

92. Díaz-Gay, M. et al. Assigning mutational signatures to individual samples and individual somatic mutations with SigProfilerAssignment. Bioinformatics 39, btad756 (2023).

    Article  PubMed  PubMed Central  Google Scholar 

93. Huang, D. et al. A highly annotated database of genes associated with platinum resistance in cancer. Oncogene 40, 6395–6405 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

94. Arafeh, R., Shibue, T., Dempster, J. M., Hahn, W. C. & Vazquez, F. The present and future of the Cancer Dependency Map. Nat. Rev. Cancer 25, 59–73 (2025).

    Article  CAS  PubMed  Google Scholar 

95. DepMap & Broad. DepMap Public 23Q4. Figshare https://doi.org/10.25452/figshare.plus.24667905.v2 (2023).

96. Cuzick, J. et al. Prognostic value of an RNA expression signature derived from cell cycle proliferation genes in patients with prostate cancer: a retrospective study. Lancet Oncol. 12, 245–255 (2011).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

97. Newman, A. M. et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat. Biotechnol. 37, 773–782 (2019).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

98. Snyder, M. W., Kircher, M., Hill, A. J., Daza, R. M. & Shendure, J. Cell-free DNA comprises an in vivo nucleosome footprint that informs its tissues-of-origin. Cell 164, 57–68 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

99. Doebley, A.-L. et al. A framework for clinical cancer subtyping from nucleosome profiling of cell-free DNA. Nat. Commun. 13, 7475 (2022).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

100. Zviran, A. et al. Genome-wide cell-free DNA mutational integration enables ultra-sensitive cancer monitoring. Nat. Med. 26, 1114–1124 (2020).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

101. Illumina. High-accuracy next-generation sequencing with the NovaSeq X Series. Data concordance application note (Illumina, 2023).

102. Stoler, N. & Nekrutenko, A. Sequencing error profiles of Illumina sequencing instruments. NAR Genom. Bioinform. 3, lqab019 (2021).

     Article  PubMed  PubMed Central  Google Scholar 

103. Zhang, K. et al. A single-cell atlas of chromatin accessibility in the human genome. Cell 184, 5985–6001 (2021).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

104. Corces, M. R. et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018).

     Article  ADS  PubMed  PubMed Central  Google Scholar 

105. Adil, M. et al. Preeclampsia risk prediction from prenatal cell-free DNA screening. Nat. Med. 31, 1312–1318 (2025).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

106. Mabbott, N. A., Baillie, J. K., Brown, H., Freeman, T. C. & Hume, D. A. An expression atlas of human primary cells: inference of gene function from coexpression networks. BMC Genom.14, 632 (2013).

     Article  CAS  Google Scholar 

107. Cords, L. et al. Cancer-associated fibroblast classification in single-cell and spatial proteomics data. Nat. Commun. 14, 4294 (2023).

     Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

108. Nieto, P. et al. A single-cell tumor immune atlas for precision oncology. Genome Res. 31, 1913–1926 (2021).

     Article  CAS  PubMed  PubMed Central  Google Scholar 

109. Browaeys, R., Saelens, W. & Saeys, Y. NicheNet: modeling intercellular communication by linking ligands to target genes. Nat. Methods 17, 159–162 (2020).

     Article  CAS  PubMed  Google Scholar 

110. Gao, T. et al. Haplotype-aware analysis of somatic copy number variations from single-cell transcriptomes. Nat. Biotechnol. 41, 417–426 (2023).

     Article  CAS  PubMed  Google Scholar 

Download references

Author notes

1. These authors contributed equally: Pushpa Itagi, Samantha L. Schuster

2. These authors jointly supervised this work: Omar Y. Mian, Hung-Ming Lam, Andrew C. Hsieh, Gavin Ha

Authors and Affiliations

1. Public Health Sciences Division, Fred Hutchinson Cancer Center, Seattle, WA, USA

   Pushpa Itagi, Samantha L. Schuster, Thomas W. Persse, Michael Yang, Alan Min, Pooja Chandra, Mohamed Adil, Patricia C. Galipeau, Yixin Lin, Minjeong Ko, Robert D. Patton, A. Patrick McDeed, Manasvita Vashisth, Daniel W. Lin, Peter S. Nelson & Gavin Ha

2. Human Biology Division, Fred Hutchinson Cancer Center, Seattle, WA, USA

   Pushpa Itagi, Samantha L. Schuster, Sonali Arora, Jennifer A. Waters, Mohamed Adil, Patricia C. Galipeau, Dmytro Rudoy, Erolcan Sayar, Robert D. Patton, Cynthia L. Wladyka, A. Patrick McDeed, Claire B. Mills, Manasvita Vashisth, Jin Yeong Kim, Peter S. Nelson, John K. Lee, Michael C. Haffner, Omar Y. Mian, Andrew C. Hsieh & Gavin Ha

3. Department of Genome Sciences, University of Washington, Seattle, WA, USA

   Allie S. Kreitman, Peter S. Nelson, Andrew C. Hsieh & Gavin Ha

4. Department of Urology, University of Washington, Seattle, WA, USA

   Lori Kollath, Abby Meis, Samuel Lindergren, Nathan Ji, Khursheed Ali, Sarah P. Psutka, John L. Gore, Daniel W. Lin, Peter S. Nelson, Eva Corey, Colm Morrissey, Jonathan L. Wright & Hung-Ming Lam

5. Immunotherapy Integrated Research Center, Fred Hutchinson Cancer Center, Seattle, WA, USA

   Hrishi Venkatesh & Lawrence Fong

6. Department of Oral Health Sciences, University of Washington, Seattle, WA, USA

   Claire B. Mills

7. Division of Hematology and Oncology, Department of Medicine, University of Washington, Seattle, WA, USA

   Rosa Nadal, Jessica E. Hawley, Todd A. Yezefski, Heather H. Cheng, Michael T. Schweizer, Evan Y. Yu, Petros Grivas, Robert B. Montgomery & Andrew C. Hsieh

8. Clinical Research Division, Fred Hutchinson Cancer Center, Seattle, WA, USA

   Rosa Nadal, Jessica E. Hawley, Peter S. Nelson, Heather H. Cheng, Michael T. Schweizer, Evan Y. Yu, Petros Grivas, Michael C. Haffner & Andrew C. Hsieh

9. Division of Hematology/Oncology, Department of Medicine and Institute of Urologic Oncology, Department of Urology, David Geffen School of Medicine at UCLA, Los Angeles, CA, USA

   John K. Lee

10. Department of Laboratory Medicine and Pathology, University of Washington, Seattle, WA, USA

    Michael C. Haffner & Funda Vakar-Lopez

11. Radiation Oncology Division, Fred Hutchinson Cancer Center, Seattle, WA, USA

    Omar Y. Mian

12. Department of Radiation Oncology, University of Washington, Seattle, WA, USA

    Omar Y. Mian

Contributions

P.I., S.L.S., O.Y.M., H.-M.L., A.C.H. and G.H. conceived and designed the study. L.K., E.C., N.J., J.E.H., C.M., A. Meis, S.L., K.A., R.N., T.A.Y., S.P.P., J.L.G., P.S.N., H.H.C., M.T.S., L.F., E.Y.Y., P.G., J.K.L., D.W.L., R.B.M., J.L.W., H.-M.L. and A.C.H. recruited patients, obtained patient consent, coordinated acquisition of clinical specimens and contributed clinical data and interpretations. P.C.G., A. Meis, K.A., C.L.W. and S.L. performed sample processing, facilitated sequencing studies and obtained CT scans. P.I., S.L.S., S.A., O.Y.M., H.-M.L., A.C.H. and G.H. designed and performed the computational analyses and performed data interpretation. T.W.P., M.Y., A. Min, P.C., P.C.G., M.A., A.S.K., Y.L., M.K., R.D.P., H.V., C.B.M., M.V., J.Y.K. and A.P.M. assisted with computational analyses or provided interpretation. S.L.S., J.A.W., D.R. and A.C.H. designed the experiments. S.L.S. performed experiments in Fig. 3f–h. J.A.W. performed experiments in Fig. 3i–l and Extended Data Fig. 7. D.R. performed experiments in Extended Data Fig. 6o–q. E.S., M.C.H., F.V.-L. and O.Y.M. provided pathology assessment, imaging data and interpretation. S.L.S., P.I., S.A., O.Y.M. and G.H. performed data visualization. S.L.S., P.I., O.Y.M., H.-M.L., A.C.H. and G.H. wrote the manuscript. P.I., S.L.S., S.A., M.A., P.C.G., R.D.P., A.P.M., C.B.M., C.M., P.S.N., M.T.S., M.C.H., O.Y.M., H.-M.L., A.C.H. and G.H. reviewed and edited the manuscript. O.Y.M., H.-M.L., A.C.H. and G.H. supervised the study. R.B.M., J.L.W., J.K.L., H.-M.L., G.H. and A.C.H. acquired funding for the study.

Corresponding authors

Correspondence to Omar Y. Mian, Hung-Ming Lam, Andrew C. Hsieh or Gavin Ha.

Competing interests

G.H. receives research support from Pfizer and consults with Quest Diagnostics; all activities unrelated to this work. A.C.H. is a consultant to the SAB of Interdict Bio; all activities unrelated to this work. C.B.M. has received consulting fees and travel support from Illumina unrelated to the current work. C.M. has received research funding from Genentech, Astra Zeneca, Janssen and Novartis, unrelated to the current work. E.C. served as a paid consultant to DotQuant and received Institutional sponsored research funding unrelated to this work from Astra Zeneca, AbbVie, Gilead, Sanofi, Zenith Epigenetics, Bayer Pharmaceuticals, Forma Therapeutics, Genentech, GSK, Janssen Research, Kronos Bio, Foghorn Therapeutics, K36 Therapeutics, BiondlessBio, and MacroGenics. E.Y.Y. has received consulting from Astellas, Johnson & Johnson, AZ, Tolmar, Merck, Bayer, Oncternal, Lantheus, Novartis, Samsung; and has received research funding to the institution provided by Dendreon, Merck, SeaGen, Blue Earth, Bayer, Lantheus, and Tyra. J.L.W. has received research support from Merck, Seagen, Pacific Edge, and Veracyte. J.K.L. serves as a research consultant for Lyell Immunopharma and Xilio Therapeutics; is a CMO and serves on the scientific advisory board for PromiCell Therapeutics. L.F. received research support from Abbvie, Merck, and Roche/Genentech; and has served on scientific advisory boards for Abbvie, Actym, Bioatla, Boehringer Ingelheim, Bristol Myer Squibb, Daiichi Sankyo, Immunogenesis, Innovent, Merck, Nutcracker, RAPT, Senti, Sutro, and Roche/Genentech. M.C.H. consults or has received honoraria from Pfizer and Astra Zeneca and has received research funding from Merck, Novartis, Genentech, Promicell and Bristol Myers Squibb. M.T.S. consults or has received honoraria from JNJ, Pfizer, Daiichi Sankyo, K36 Therapeutics and Bayer; and has received research support to his institution from Novartis, Zenith Epigenetics, BMS, merck, Epigenetix, Xencor, Ambrx, Oric Pharmaceuticals, AZ, Lightspeed and Kyntra. O.Y.M. consults with Veracyte, RiboX Therapeutics, and Electa; all activities unrelated to this work. P.S.N. served as a paid consultant to AstraZeneca, Genentech, and Pfizer and received research support from Janssen for work unrelated to the current study. P.G., for the past 2 years, has consulted with MSD, Bristol Myers Squibb, AstraZeneca, EMD Serono, Pfizer, Janssen, Roche, Astellas Pharma, Gilead Sciences, AbbVie, Bicycle Therapeutics, Replimune, Daiichi Sankyo, Foundation Medicine, Eli Lilly, Urogen, Tyra Biosciences, Natera, Ottimo Pharma, Bayer; and received research funding from MSD, EMD Serono, Gilead Sciences, Acrivon Therapeutics, ALX Oncology, Genentech (paid to institution). S.P.P. consults with J&J, CG Oncology, and Merck. T.A.Y. consults with Fennec, and Dendreon.

Peer review information

Nature thanks Joost L. Boormans, Jeffrey Damrauer, Kent Mouw, J. Alberto Nakauma-González and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, 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 you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. 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-nc-nd/4.0/.

Reprints and permissions