Article The Genomic and Immune Landscapes of Lethal Metastatic Breast Cancer Graphical Abstract Authors Leticia De Mattos-Arruda, Stephen-John Sammut, Edith M. Ross, ..., Florian Markowetz, Joan Seoane, Carlos Caldas Correspondence jseoane@vhio.net (J.S.), carlos.caldas@cruk.cam.ac.uk (C.C.) In Brief De Mattos-Arruda et al. profiled multiple metastases from autopsies of patients with therapy-resistant breast cancer, showing that multi-clonal spreading occurs in a small number of founder events. The analysis characterizes predicted neo-antigen landscapes, tumor microenvironments, and accumulation of HLA LOH. T cell immune responses appear to co-evolve with metastatic cancer genomes. Highlights d Genomic and transcriptomic landscapes for 10 lethal breast cancers d Within a patient, metastases group in limited clades with shared genomic ancestry d Tumor immune microenvironments across metastases are not uniform d Phylogenetic trees are correlated with TIL-TCR trees across metastases De Mattos-Arruda et al., 2019, Cell Reports 27, 2690–2708 May 28, 2019 ª 2019 The Authors. https://doi.org/10.1016/j.celrep.2019.04.098 Cell Reports Article The Genomic and Immune Landscapes of Lethal Metastatic Breast Cancer Leticia De Mattos-Arruda,1,2,20 Stephen-John Sammut,1,20 Edith M. Ross,1 Rachael Bashford-Rogers,3 Erez Greenstein,4 Havell Markus,5,6 Sandro Morganella,7 Yvonne Teng,8 Yosef Maruvka,9,10 Bernard Pereira,1 Oscar M. Rueda,1 Suet-Feung Chin,1 Tania Contente-Cuomo,5,6 Regina Mayor,2,11 Alexandra Arias,2,11 H. Raza Ali,1 Wei Cope,1 Daniel Tiezzi,1 Aliakbar Dariush,12 Tauanne Dias Amarante,7 Dan Reshef,4 Nikaoly Ciriaco,13 Elena Martinez-Saez,11,13 Vicente Peg,11,13,14 Santiago Ramon y Cajal,11,13,14 Javier Cortes,2,15 George Vassiliou,8,16,17 Gad Getz,9,10 Serena Nik-Zainal,7 Muhammed Murtaza,5,6 Nir Friedman,4 Florian Markowetz,1 Joan Seoane,2,11,18,21,* and Carlos Caldas1,19,21,22,* 1Department of Oncology and Cancer Research UK Cambridge Institute, Li Ka Shing Centre, University of Cambridge, Cambridge CB2 0RE, UK 2Vall d’Hebron Institute of Oncology (VHIO), Vall d’Hebron University Hospital, Barcelona 08035, Spain 3Department of Medicine, University of Cambridge, Cambridge CB2 0QQ, UK 4Department of Immunology, Weizmann Institute of Science, Rehovot 76100, Israel 5Center for Noninvasive Diagnostics, Translational Genomics Research Institute, Phoenix, AZ 85004, USA 6Mayo Clinic Center for Individualized Medicine, Scottsdale, AZ, USA 7Department of Medical Genetics, The Clinical School, University of Cambridge, Cambridge CB2 0QQ, UK 8Cancer Molecular Diagnosis Laboratory, NIHR Cambridge Biomedical Research Centre, Cambridge, UK 9The Broad Institute, Cambridge, MA 02142, USA 10Massachusetts General Hospital Cancer Center and Department of Pathology, Charlestown, MA 02129, USA 11Spanish Biomedical Research Network Centre in Oncology (CIBERONC), Madrid, Spain 12Institute of Astronomy, University of Cambridge, Cambridge CB3 0HA, UK 13Department of Pathology, Vall d’Hebron University Hospital, 08035 Barcelona, Spain 14Translational Molecular Pathology, Vall d’Hebron Research Institute (VHIR), Universitat Autònoma de Barcelona, 08035 Barcelona, Spain 15Ramon y Cajal Hospital, 28034 Madrid, Spain 16Wellcome Trust Sanger Institute, Wellcome Trust Genome Campus, Hinxton, Cambridge, UK 17Wellcome Trust/MRC Cambridge Stem Cell Institute, Cambridge, UK 18Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona 08010, Spain 19Breast Cancer Programme, Cancer Research UK Cambridge Cancer Centre, Cambridge CB2 2QQ, UK 20These authors contributed equally 21Senior author 22Lead Contact *Correspondence: jseoane@vhio.net (J.S.), carlos.caldas@cruk.cam.ac.uk (C.C.) https://doi.org/10.1016/j.celrep.2019.04.098 SUMMARY INTRODUCTION The detailed molecular characterization of lethal The genomic characterization of large numbers of primary breast cancers is a prerequisite to understanding resis- tumors has revealed significant inter-tumor heterogeneity and tance to therapy and escape from cancer immunoe- unraveled an increasingly refined molecular taxonomy of early diting. We performed extensive multi-platform breast cancer with profound implications for prognostication profiling of multi-regional metastases in autopsies and therapeutic stratification (Cancer Genome Atlas, 2012; Curtis et al., 2012; Dvinge et al., 2013; Nik-Zainal et al., 2012a, 2012b, from 10 patients with therapy-resistant breast can- 2016). Intra-tumor genomic heterogeneity is also seen in early cer. The integrated genomic and immune land- breast cancers, highlighting that complex clonal architectures scapes show that metastases propagate and evolve are already present in primary tumors (Nik-Zainal et al., 2016; as communities of clones, reveal their predicted Pereira et al., 2016; Shah et al., 2012; Yates et al., 2015). The neo-antigen landscapes, and show that they can tumor microenvironment (TME) in primary tumors is also different accumulate HLA loss of heterozygosity (LOH). The and distinctive across breast cancer subtypes, in particular with data further identify variable tumor microenviron- regards to adaptive immunity (Ali et al., 2016a; Rooney et al., ments and reveal, through analyses of T cell recep- 2015). The nature of the adaptive immune response, the status tor repertoires, that adaptive immune responses of immunoediting, and the diversity of the T cell receptor (TCR) appear to co-evolve with the metastatic genomes. repertoire have been analyzed in some early breast cancers These findings reveal in fine detail the landscapes (Munson et al., 2016; Park et al., 2016; Wang et al., 2017), but such information for metastatic lesions is lacking. of lethal metastatic breast cancer. 2690 Cell Reports 27, 2690–2708, May 28, 2019 ª 2019 The Authors. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Large-scale studies reporting the genomic and transcriptomic sequencing (sWGS) to obtain CNA profiles in 168 samples characterization of breast cancer metastasis (Robinson et al., from the 10 cases: 122 tumor biopsies (109 metastasis at 2013) and whole-genome sequencing of matched primary tu- autopsy, 4 metastatic biopsies during treatment, 8 primary mors and metastases (Yates et al., 2017) have identified targets breast surgical or diagnostic biopsy specimens from 6 cases, that are enriched in metastases compared with primary tumors. and a nasopharyngeal olfactory neuroblastoma) and 46 body Despite their size, these studies were generally limited to single fluid samples (24 plasma samples, 5 ascites samples, 9 cerebro- metastatic samples. spinal fluid [CSF] samples, 7 pleural fluid samples, and 1 pericar- Genomic evolution is seen in breast cancer metastases dial fluid sample) (Table S1). For 64 of these metastatic samples compared with their matched primary tumors. This was first re- (from 9 of the 10 cases), we also performed RNA sequencing ported in single cases (Shah et al., 2009; Ding et al., 2010). A (RNA-seq). The combined sWGS and RNA-seq data were used couple of warm autopsy case reports have also revealed hetero- to classify individual metastasis into one of the 10 genome geneity of genomic resistance mechanisms to targeted therapies driver-based subtypes, called integrative clusters (IntClust) (Ali across metastases (Juric et al., 2015; Murtaza et al., 2015). More et al., 2014; Curtis et al., 2012). recently two small autopsy studies with multiple-metastases The tumor CNA profiles were remarkably similar across metas- profiling have confirmed significant inter-metastasis heteroge- tases in 9 of the 10 cases, and, as expected, all metastases were neity (Hoadley et al., 2016; Savas et al., 2016). classified into the same IntClust (Figure 1B; SI1 in https://doi.org/ However, a comprehensive molecular analysis of lethal breast 10.17632/6cv77bry6m.1). An exception was case 288, an ER- cancers, interrogating both the malignant and TME compartments positive lobular breast cancer, where, besides a 1q gain and and TCR repertoires, across multiple metastases in several cases 16q loss seen in all metastases, there were additional and mutually is still lacking. Here we report extensive multi-platform molecular exclusive CNAs: amplification of 11q13/14, including CCND1 and profiling (DNA sequencing, RNA sequencing, TCR sequencing, PAK1, in lymph nodes (288-005 and 288-006), the ascites fluid digital pathology of H&E sections, and immunohistochemistry) of cell pellet, and ovaries (classified as IntClust2) and 8q and 10q multiple individual metastases from 10 warm autopsies of patients amplifications in brain and lung and pleura (classified as with lethal multi-therapy-resistant breast cancers. This collection IntClust10). These data suggest that all metastases shared a com- allowed us to characterize the mutational and copy number aber- mon ancestor with 1q gain and 16q loss, and early sub-clonal evo- ration (CNA) landscapes across the individual metastasis to infer lution with remarkable genomic divergence then occurred. The the clonal ancestries of metastases, assess the TME in each indi- PAM50 intrinsic subtypes (Parker et al., 2009) were less consistent vidual metastasis, characterize the predicted neo-antigens, and across metastases and failed to capture the clade segregation in assess the TCR repertoires across metastases, providing a case 288 (SI1 in https://doi.org/10.17632/6cv77bry6m.1). detailed molecular characterization of lethal breast cancers that We used whole-exome sequencing (WES) at a median of had been subjected to multiple lines of systemic therapies. 1323 depth, including, in each case, DNA extracted from the buffy coat as the matched germline reference to characterize RESULTS the somatic mutational landscape across 79 metastases and 7 body fluid samples, with a range of 2 to 19 metastatic samples Multi-site Genomic and Transcriptomic Landscapes of per patient sequenced. We analyzed the WES data with rigorous Lethal Metastatic Breast Cancers filters (STAR Methods) based on a recently described pipeline We performed warm autopsies in 10 patients with metastatic (Callari et al., 2017). To further validate the WES mutation calls breast cancer that had become resistant to multiple lines of ther- (Table S2), we generated ultra-deep targeted sequencing (TS) apy (Figure 1A; Figure S1). The cohort is fairly representative of (mean depth of 7,57029,8913; mean coverage > 1,0003 for the major subtypes: 8 were diagnosed with estrogen receptor 71%–100% of samples) for amplicons across 464 mutations (ER)-positive disease, and three of these were HER2-positive; (average, 46.4/case; range, 16–127) (Table S3). Matched WES one was ER-negative and HER2-positive; and one was triple- and TS data were available from 79 samples, and they validated negative. In total, 181 samples from multiple metastatic sites in the robustness of the WES data analysis pipeline we used: sensi- each patient (mean, 18.5/patient; range, 5–37) were collected tivity, 0.85; specificity, 0.99; accuracy, 0.91; precision, 0.99 and either fresh-frozen (for DNA and RNA extraction) or (STAR Methods). formalin-fixed and paraffin-embedded (FFPE). FFPE samples The WES data identified 15,430 somatic mutations across the from the original breast surgery or diagnostic biopsy were avail- 86 samples: 7,809 missense, 1258 truncating, 10 nonstop, 234 able from 6 of the patients, and metastatic biopsies during treat- in-frame deletions or insertions, and 6,119 other (Table S2). ment were also collected from 3 patients. The FFPE samples The mutation burden varied between cases, with a median of were used for histological and immunohistochemistry analysis 507 mutations per case (range, 113–997), and across metasta- and for DNA extraction. Plasma samples and a selection of ses within each case, with a median mutation load of 146 per body fluids, collected during the patient’s life or at autopsy, metastasis (interquartile range, 86.25). These numbers are were available from all 10 patients (mean, 4.7/patient; range, significantly greater (p % 2.2e16, Wilcoxon rank-sum test) 1–9/patient) and used for cell-free DNA (cfDNA) extraction. than mutation burdens reported previously for primary breast Comprehensive clinical information for the patients and analyzed tumors within The Cancer Genome Atlas (TCGA), with a median samples can be found in Table S1. mutation load of 63.5 (interquartile range, 68). The genomic landscapes of breast cancers are dominated by We classified somatic mutations found in WES as ‘‘metastatic CNAs (Ciriello et al., 2013). We used shallow whole-genome stem’’ when present in all metastases from the same case, Cell Reports 27, 2690–2708, May 28, 2019 2691 Figure 1. Molecular Profiling of 10 Lethal Metastatic Breast Cancers (A) Silhouettes representing the 10 patients with number and type of samples profiled using each platform. Patients are grouped as indicated above the silhouettes according to ER and HER2 status. Samples profiled are labeled according to the color key panel on the right. WES, whole-exome sequencing; sWGS, shallow whole-genome sequencing; TS, targeted sequencing; RNA-seq, RNA sequencing; RNA-TCRseq, targeted TCR sequencing in RNA; IHC, immunohis- tochemistry; TNBC, triple-negative breast cancer. (B) IGV plot showing the copy number aberration (CNA) landscapes across 168 metastases, with samples grouped by patient. The IntClust bar shows individual sample assignment to one of the 10 integrative clusters (Curtis et al., 2012). Copy number gains and amplifications are indicated in shades of red; copy number losses are indicated in shades of blue (see bar for corresponding Log ratio). IntClust as per color bar. NC- not classified; ON- olfactory neuroblastoma. 2692 Cell Reports 27, 2690–2708, May 28, 2019 Figure 2. Mutational Landscape of 10 Lethal Metastatic Breast Cancers (A) Mutational burden barplots across 86 metastatic samples using WES. Colors indicate mutations classified as metastatic stem, metastatic clade, and metastatic private. (B) Oncoprint plot showing the mutations in breast cancer driver genes identified by WES across 84 metastases for the 10 patients. (legend continued on next page) Cell Reports 27, 2690–2708, May 28, 2019 2693 ‘‘metastatic clade’’ when present in at least two but not all metas- mutations, whereas metastatic private driver mutations were tases, and ‘‘metastatic private’’ when present in a single metas- uncommon. tasis (Figure 2A). This revealed that the vast majority of the TS data were also obtained from 40 additional samples for mutations in individual metastasis were either stem or clade. which only FFPE blocks were available, bringing the total number We focused our analysis on mutation drivers (Table S4). We of metastatic samples, primary tumors, and liquid biopsies with considered a gene a mutation driver (e.g., associated with a TS data to 159 (average per patient, 16 samples; range, 4–25). fitness advantage when somatically mutated) using the widely The TS validated and extended the WES findings, and this was accepted framework that is based on analyses of large cancer particularly informative in case 288, showing that the bilateral mutation datasets (Lawrence et al., 2013; Vogelstein et al., ovarian metastases shared the driver mutations with the lymph 2013). For our analyses, we defined a list of 109 breast cancer nodes and ascites (Figure 2C; SI2 in https://doi.org/10.17632/ driver genes identified from reviewing the data in three large 6cv77bry6m.1). In 6 cases, FFPE blocks from the original primary cohorts (Lefebvre et al., 2016; Nik-Zainal et al., 2016; Pereira breast tumor were available, and TS data confirmed that all of et al., 2016) and a list of 527 non-breast cancer drivers (non-over- these contained the clonal ancestors of the metastases, but a lapping) from the Cancer Gene Census (https://cancer.sanger.ac. percentage of stem mutations and an even larger fraction of clade uk/census). Mutations were significantly (p = 0.0002, chi-square mutations were not detected (Figure 2D; SI2 in https://doi.org/10. test) more common in breast cancer driver genes (58 of 109, 17632/6cv77bry6m.1). This included some metastatic stem 53.21%) than in non-breast cancer driver genes (142 of 527, driver mutations (case 290, BUB1B absent in two FFPE blocks; 26.94%). The total number of driver mutations per metastasis case 308, ESR1 and PTEN absent in the two FFPE blocks; averaged 11.44 (range, 2–30) (Figure 2B; Figure S2), which is DET52, ATM not detected in ductal carcinoma in situ (DCIS) or higher than the estimated number per primary tumor (Martincor- metastatic axillary lymph nodes) and most metastatic clade driver ena et al., 2017). Metastatic stem driver mutations were identified mutations (SI2 in https://doi.org/10.17632/6cv77bry6m.1). in all 10 cases: 2 in case 288 (CDH1; non-breast driver: CALR), 9 in We next asked whether the expression of the mutant allele case 290 (BUB1B, MAP2K4, MAP3K1, NCOR1, TBX3, and TP53; was similar across mutations. A combined analysis of WES non-breast drivers: ELL, MET, and FLT3), 1 in case 291 (ERBB3), 3 and RNA-seq data were possible in 8 cases (case 291 with a sin- in case 308 (ESR1 and PTEN; non-breast driver: NKX2-1), 4 in gle metastasis with combined data and case DET52 without case 315 (ATM and GATA3, one missense and one truncating; RNA-seq data were excluded) and revealed that the normalized non-breast driver: CDC73), 5 in case 323 (CDH1 and PTEN; expression of the mutant allele was highest in stem, lower in non-breast drivers: MET, MYCL, and SDC4), 1 in case DET52 clade, and lowest in private mutations (Figure 2E). (ATM), 8 in case 328 (FGFR2, NOTCH1, and TP53; non-breast In summary, metastases keep accumulating mutations, drivers: AXIN1, CRTC1, LRIG3, SMARCE1, and WHSC1L1), and including mutations in known cancer driver genes, but an 5 in case 330 (TP53; non-breast drivers: CARD11, CNTRL, apparent hierarchy of expression (stem-clade-private) of mutant FBXO11, and PTPRK). A Li-Fraumeni syndrome patient, case alleles suggests that, as more mutations accumulate in metasta- 298, had two known metastasized malignancies: HER2-positive ses, these are increasingly passengers (e.g., not expressed). breast cancer and olfactory neuroblastoma. Previously, we A fraction of mutations (including drivers) shared across metas- showed that the brain metastasis originated from the breast can- tases (stem and clade) were not detectable in the available cer (De Mattos-Arruda et al., 2015). The HER2-amplified brain primary tumor tissue blocks, suggesting either their origin from metastases had 5 mutation drivers (non-breast drivers: a minor clone in the primary tumor or their acquisition in metasta- CACNA1D, DCTN1, FAT1, RAD21, and WHSC1) in addition to tic cells that had already left the breast. the germline TP53 mutation with associated somatic loss of heterozygosity. The WES data also revealed that the olfactory Metastases Are Initiated and Maintained as neuroblastoma had an RB1 stem driver mutation. This RB1 muta- Communities of Clones tion was detected in one of the two breast cancer brain metasta- A monoclonal origin of metastases should result in uniformly high ses (298-009), likely because of contamination by CSF cfDNA. variant allelic fractions (VAFs) of stem and clade mutations Indeed, mutations arising from both leptomeningeal neuroblas- across all metastases in a case. Plots of allelic fractions of these toma and from HER2-positive brain metastases had been de- mutations across individual metastases revealed a scatter of tected in CSF cfDNA (De Mattos-Arruda et al., 2015). In case allelic fractions using WES data (Figure 3A), and this was vali- 290, an ovarian tumor sample originally presumed to be metasta- dated using deep TS data. This observation was confirmed by tic lacked all 6 breast cancer stem mutations, including a TP53 calculating the cancer cell fraction (CCF), which is the VAF of frameshift mutation. The sample had a different TP53 p.Y220D each somatic mutation corrected for copy number and purity missense mutation, and this prompted a pathology review, which estimates, across all of the individual metastases. The probability confirmed that it was an independent primary ovarian adenocar- distributions of the CCFs were then used to classify each somatic cinoma. Most cases had also several metastatic clade driver mutation as clonal or sub-clonal (McGranahan et al., 2015; (C) Oncoprint plot showing driver mutations validated by TS (allelic fraction [AF] R 0.1%) for case 288. (D) Boxplot showing the percentage of stem and clade mutations identified as present by TS (AF R 3 SD from AF in matched normal). DNA was extracted from FFPE blocks from primary surgery samples, except for case DET52, where P1 and P3 were diagnostic biopsies (breast and axillary lymph node, respectively). (E) Boxplots of Z score-normalized mutant allele expression from RNA-seq data in metastatic stem, metastatic clade, and metastatic private mutations. TPM, transcripts per million. Bars indicate significance of difference (p values < 0.05 are considered statistically significant). 2694 Cell Reports 27, 2690–2708, May 28, 2019 Figure 3. Breast Cancer Metastases Are Communities of Clones (A) Pairwise comparisons of raw VAFs from WES data across 10 pairs of metastases from case 290. Metastatic stem and metastatic clade mutations are colored as indicated. (B) Cancer cell fraction (from WES data) of metastatic stem and metastatic clade mutations across metastases in case 290. Each symbol represents a somatic mutation in an individual metastasis. Probability distributions of the CCFs were used to classify each mutation as either clonal (blue) or subclonal (red). Error bars represent the 95% confidence interval. Plots for all remainder cases are shown in SI3 in https://doi.org/10.17632/6cv77bry6m.1. (C) Mean cellular prevalence of mutation clusters identified by PyClone from WES data across metastases in case 290. Metastatic stem (clusters 6, 9, 10, and 15) and metastatic clade (clusters 2, 5, 7, 8, 11, 13, and 14) mutation clusters are shown. (D) Boxplots showing the distribution of mutation AFs in TS data in case 290. Amplicons representative of PyClone exome-derived mutation clusters were analyzed. Plots for all remainder cases are shown in SI3 in https://doi.org/10.17632/6cv77bry6m.1. Cell Reports 27, 2690–2708, May 28, 2019 2695 (legend on next page) 2696 Cell Reports 27, 2690–2708, May 28, 2019 Figures 3B and S3; https://doi.org/10.17632/6cv77bry6m.1), and resulting from high-quality WES being available from a single the results were consistent with a fraction of the metastatic stem metastasis (291-015). In case 298, there were two breast cancer and clade mutations being sub-clonal. metastases analyzed, with the remainder of metastases We also analyzed the WES data with PyClone (Roth et al., sequenced from neuroblastoma, and metastases segregated 2014), a Bayesian clustering method for grouping sets of somatic by tumor of origin (Figure S4). Case 288 was an ER-positive mutations and estimating their cellular prevalence (Figure 3C; Fig- lobular cancer with a classical somatic truncating CDH1 muta- ure S3; SI3 in https://doi.org/10.17632/6cv77bry6m.1). In each tion (Figure 4A). This patient had bone metastases 10 years after case, a cluster constituted by metastatic stem mutations had diagnosis, followed by contralateral axillary lymph node metas- the highest predicted cellular prevalence (mean, 0.94; range, tases and, later, lung and CNS metastases. The metastases 0.81–0.99) across metastases, as would be expected (288, clus- detected clinically were sampled at autopsy and, in addition, ter 2; 290, cluster 6; 308, cluster 1; 315, cluster 3; 323, cluster 2; metastases found at autopsy in both ovaries and the uterus, 328, cluster 4; 330, cluster 5; Det52, cluster 1). However, there and ascites was also collected. The genomic phylogeny was were also stem mutation clusters with lower predicted cellular clear, with two separate clades: one consisting of lymph nodes prevalence, indicating that these were probably sub-clonal (290: (288-005 and -006), ascites (288-022), and bilateral ovarian Kru- clusters 9, 10, and 15; 308: clusters 5 and 9; 328: cluster 5; 330: kenberg metastases (13A-0267A10, -A26, and -A27) that were cluster 6; Det52: cluster 6). Amplicons representative of the clonal all ER-positive and had a truncating KMT2C R3868Kfs*13 muta- clusters identified from WES were validated using TS, and tion and 11q13 amplification and a second clade with brain (288- this confirmed that a fraction of stem and clade mutations ap- 008), leptomeningeal (CSF pellet, 288-023), and lung metastases peared to be sub-clonal (Figure 3D; SI3 in https://doi.org/10. (288-016, -017, -020, and -021) that were ER-negative, shared 17632/6cv77bry6m.1). TP53, MAP3K1, and ARID5B mutations, and lacked KMT2C- In summary, these analyses are not compatible with all meta- R3868Kfs*13 (Figure 4A). Case 290 (Figure 4B), an ER-positive/ static stem and clade mutations being clonal and support the HER2 ductal cancer, relapsed with bone metastases 2 years hypothesis that metastases are initiated and maintained as after diagnosis, followed by liver metastases 13 years later and groups of cellular clones. death shortly after development of CNS metastases 20 years after the original diagnosis. At autopsy, several liver metastases Ancestries of Multiregional Metastases Defined by were sampled and carefully mapped, in addition to brain and Phylogenetic Analyses stomach metastases being collected. The ovarian sample We aimed to reconstruct the metastatic seeding patterns with a collected at autopsy was proven genomically and upon histolog- series of phylogenetic methods. These included OncoNEM ical review to be a separate primary adenocarcinoma (see (Ross and Markowetz, 2016) and Treeomics (Reiter et al., above). The genomic phylogeny of the breast cancer metastases 2017), which employ nucleotide substitutions and short inser- revealed 3 clades. One clade, defined by the presence of an tions and deletions (from WES or TS data) as input data; Super- ESR1 Y537S mutation, in ER-positive metastases (290-014, Freq (Savas et al., 2016), which employs both single-nucleotide -015, -016A, -016B, and -021) mapped to the right inferior lobe variants (SNVs) and CNAs from WES as input data; and MEDICC of the liver. A second clade of ER-negative metastases in the (Schwarz et al., 2014), which employs CNAs from sWGS and brain (290-005, -007, and -008) and left (290-018 and -019) WES data as input (STAR Methods). Overall, the results were and upper right (290-017) lobes of the liver was defined by the consistent across these different methods, but for clarity, we presence of a KMT2A mutation. A third clade was defined by present below the results from OncoNEM (Figure 4; Figure S4A) the absence of both ESR1 and KMT2A mutations in a stomach using either WES and/or TS data (all other results are provided in metastasis (290-024). Case 308 (Figure 4C), an ER-positive/ SI4 in https://doi.org/10.17632/6cv77bry6m.1). The OncoNEM HER2 ductal cancer, relapsed with bone metastases 3 years phylogenetic trees of metastases had branched structures in after surgery, followed 3 years later by lung metastases and nine cases (Figure 4; Figure S4A; SI4 in https://doi.org/10. then, in quick succession, CNS and liver metastases and a medi- 17632/6cv77bry6m.1). These trees had a limited number of astinal mass (formed by pericardial, lung, and rib metastases) main branches forming separate clades of distinct but geneti- developing in the year prior to death. All 19 metastases analyzed cally related, by common ancestry, metastatic samples. The by WES shared activating ESR1 and truncating PTEN mutations. exception was case 291, where all metastases appeared linearly A genomic clade was defined by a PIK3CA mutation, shared by related using TS data (Figure S4A; SI4 in https://doi.org/10. pericardial (308-015, -016, and -017) and bone metastases (308- 17632/6cv77bry6m.1), but this appearance could be an artifact 018, -019, -020, -021, and -022) forming the mediastinal mass. Figure 4. Phylogenetic Ancestries of Breast Cancer Metastases (A–D) Phylogenetic trees from the OncoNEM algorithm. Shown are cases 288 (A), 290 (B), 308 (C), and 315 (D). Metastatic stem driver mutations and selected metastatic clade mutations are shown. Boxes identify clades. Tree branches are proportional to the number of mutations. (E) Phylogenetic tree from the LICHeE algorithm for case 290. Circles represent the mutation clonal clusters and digits within each circle the number of mutations for each cluster. Squares represent each individual metastasis, with colored rectangles representing the cellular prevalence of the respective clonal cluster. Cross-seeding from the KMT2A clade to 3 metastases (014, 016-A, and 016-B-WT/muc/IDC) in the ESR1 clade can be seen. Similar plots for all remainder cases are shown in SI4 in https://doi.org/10.17632/6cv77bry6m.1. (F) Mutation barplots colored according to mutational signatures for cases 288 and 298. Case 288: all mutations (left panel) and clade mutations (right panel). Case 298: all mutations across samples (left panel) and private mutations of Her2+ breast cancer metastases (right panel). Cell Reports 27, 2690–2708, May 28, 2019 2697 Most of the remainder of the metastases (except 3 of the menin- cellular prevalence. This revealed that cross-seeding between geal metastases: 308-001, -002, and -003) formed a separate clades occurs in a particular pattern: a common seeder for one clade defined by a BRCA2 missense mutation. Case 315 (Fig- clade can cross-seed metastases belonging to a separate clade ure 4D), an ER-positive/HER2+ ductal cancer treated with neo- (Figure 4E; SI4 in https://doi.org/10.17632/6cv77bry6m.1). adjuvant chemotherapy and anti-HER2 therapy followed by Cross-seeding within a clade was rare. surgery, relapsed 3 years after diagnosis with bone and liver In summary, the genomic phylogenies of metastases are com- metastases, received multiple lines of therapy mostly targeting plex but consistently show that, within a patient, individual HER2, and within 1 year of death (8 years after diagnosis) had metastasis can be grouped in phylogenic clades that share a progressive liver and CNS metastases. At autopsy, all 12 metas- common genomic ancestry, and this ancestry is mutually exclu- tases analyzed by WES shared ATM and two different GATA3 sive with the genomic ancestries of other clades. Each clade mutations. Two genomic clades were identified: one defined by group of metastases is therefore likely seeded by a common a truncating mutation in PPP2R1A, comprised of metastases in ancestor, but cross-seeding between metastases can happen, the liver (315-001, and -003), peri-pancreatic lymph node (315- although it appears that this occurs to a rather limited degree 015), para-tracheal lymph nodes (315-007 and -014), and and mostly between clades. meninges (315-017 and -018) and a second clade defined by an ESRRB mutation in liver (315-002, -004, and -005) and lung Mutational Signatures across Metastases (315-009 and -012) metastases. Case 323 (Figure S4A; SI4 in Mutational signatures, generated by different mutational pro- https://doi.org/10.17632/6cv77bry6m.1), an ER-positive/HER2 cesses, are engraved in the genomes of breast cancers (Nik-Zainal lobular cancer, relapsed 2 years after diagnosis with metastases et al., 2016). To extract these signatures accurately, direct applica- in bone, pleurae, and lymph nodes. The patient died a year later tion of non-negative matrix factorization (NNMF) on our 86 WES with CNS involvement. Genomic phylogeny showed that the two samples would lack power; hence, data from 240 additional lymph node metastases formed a single clade with CDH1 and WES single metastatic breast samples (Lefebvre et al., 2016) PTEN stem mutations. Case 328 (Figure S4A), an ER-negative/ were included in the analysis. To identify canonical signatures, HER2+ ductal carcinoma, presented with breast primary and the rank and number that could be extracted by NNMF were metastatic disease in the liver, and the patient died 19 months allowed to vary, and four were robustly seen: the APOBEC-related later. At autopsy, three anatomically distinct brain metastases signatures 2 and 13 (cosine similarity of 0.98 and 0.96, respec- were sampled; all shared TP53 and FGFR2 mutations and formed tively); signature 17, comprising mostly T > G mutations (cosine a single clade. Case 330 (Figure S4A), a triple-negative ductal similarity of 0.95); and signature 1, associated with demethylation cancer, was originally treated with neo-adjuvant chemotherapy, of cytosines (cosine similarity of 0.91) (Figure 4F; Figure S4B). with pathological complete response at surgery. The patient We explored whether signatures stratified across stem, clade, relapsed after 12 months with CNS metastases, followed months and private mutations. In case 288, considering all mutations re- later by breast metastases. Phylogenetic analysis revealed two vealed signatures 1, 2, and 13 in all metastases and signature clades: breast metastases (330-001 and -002) with GPS2, 17 only in ER-negative metastases, whereas considering clade MSH2, and FSIP2 mutations and brain (330-005) and meningeal mutations revealed that signature 13 was exclusive of ER-positive (330-003 and -004) metastases with AR and MLLT4 mutations. metastases (Figure 4F, top panels). In case 298, signatures 2 and Case DET52 (Figure S4A), an ER-positive/HER2+ ductal carci- 13 were seen only in Her2+ breast cancer brain metastases, with noma, presented with widespread metastatic disease and had enrichment of signature 13 in private mutations (Figure 4F, bottom a tree rooted in a single brain metastasis resected surgically panels). In case 328, signature 17 was seen across all 3 brain me- 18 months prior to death (Det52, mt4), with a clade formed tastases but restricted to private mutations in only two (SI4 in by ovary (Det52, mt7), liver (Det52, mt5), and lung (Det52, mt6) https://doi.org/10.17632/6cv77bry6m.1), suggesting enrichment metastases. The OncoNEM tree constructed from TS data (SI4 with later tumor evolution. In case DET52, a single lung metastasis in https://doi.org/10.17632/6cv77bry6m.1) showed one clade (DET52-mt6) with signature 17 (Figure S4B) carried in that context formed by diagnostic DCIS and axillary lymph node (Det52, the ERBB4 mutation, believed to be the driver of resistance to mt2) biopsies and invasive breast cancer sampled at autopsy lapatinib, which ultimately killed the patient. (Det52, mt3), and another clade formed by the distal metastases The assignment of mutations to the four canonical signatures collected at autopsy (bone, mt8; liver, mt5; ovary, mt7) and the revealed a considerable number of ‘‘residual’’ mutations (SI4 in surgically resected brain metastasis (mt4). https://doi.org/10.17632/6cv77bry6m.1). Mathematically small Because metastases were grouped in clades, we asked fluctuations in either direction of these residual mutations may whether cross-seeding occurred both within and between reflect the lack of power in WES data. However, a consistent clades. TS data revealed individual instances of cross-seeding excess of positive residuals indicates that many mutations may (Table S3; SI2 in https://doi.org/10.17632/6cv77bry6m.1), and be due to mutational processes previously unaccounted for. therefore we used a systematic approach to quantify these. To exclude the possibility that these residual mutations arise Cross-seeding is the result of a clone (or a group of clones) simply because of fitting to fewer canonical signatures than truly from one metastasis recirculating and seeding another metas- present, the dataset was fitted to the 12 breast cancer-associ- tasis at a different site. We used the PyClone mutation clusters ated signatures and all 30 canonical signatures in the Catalogue (see above) as a surrogate of metastatic cellular clones and of Somatic Mutations in Cancer (COSMIC). This revealed entered these into the LICHeE phylogenetic tool (Popic et al., improvement in overall fitting (as expected) at the expense of 2015), which orders clones across samples by comparing their increasing negative residuals (p = 4e–16, Wilcoxon signed-rank 2698 Cell Reports 27, 2690–2708, May 28, 2019 Figure 5. Neo-antigen Landscape across Breast Cancer Metastases (A) Bar plots of the neo-antigen landscape across cases (top panel) and LOH at the HLA allelic locus across metastases (bottom panel). (B) Violin plots of observed/expected neoantigen ratios across individual metastases. For each metastasis, 100 replicate expected mutation simulations were used, and each violin plot shows the distribution of the log2-transformed ratio. The ratio represents the relative deviation of the neo-epitope rate from expectation. test), suggesting mis-assignment to signatures that are unlikely diversity within metastatic samples compared with primary sam- to be present. Next, to demonstrate that the excess of residuals ples (Wilcoxon rank-sum test, p = 1.6e16). Moreover, the evo- was not simply due to using WES data, the same analysis was lution of diversity through the phylogenetic trees revealed a done in 640 WES primary cancers (SI4 in https://doi.org/10. greater Shannon entropy index when all mutations were consid- 17632/6cv77bry6m.1). Comparison of residuals (root-mean- ered (Wilcoxon signed-rank test, p = 1e14) versus stem muta- square error [RMSE]) between primary and metastatic samples tions (SI4 in https://doi.org/10.17632/6cv77bry6m.1). revealed that the fitting was much worse for metastatic cancers In summary, although metastases carry the same mutational (p < 2.2e16, Wilcoxon rank-sum test). Hence, we hypothesize signatures described in primary breast cancers, the increased that metastases with a longer history of exposure to mutational residuals observed suggest that they probably carry additional processes may carry additional signatures that are detected as mutation patterns. Remarkably mutation signatures are either the excess residuals. Indeed, cosine similarities between muta- shared across all metastases or across metastases within a tional profiles of metastatic samples and primary samples (with clade, suggesting that they are scars of mutational processes additional bootstrapping performed 10,0003; median p = operative in the metastatic founder clones. 3e4) revealed greater interpatient similarity in primary versus metastatic cancers. To further support this finding, we used The Predicted Neo-antigen Landscape across the Shannon entropy index to quantify the diversity within the Metastases normalized mutational profile (considering the 96 permutations Neo-antigens encoded by tumor-mutated genes result in neo- of triplet mutation context) of each sample and found greater peptides that can be presented by the major histocompatibility Cell Reports 27, 2690–2708, May 28, 2019 2699 class I complex. These neo-peptides have the potential for binding Heterogeneity of the TME across Metastases TCRs and eliciting anti-tumor adaptive immune responses (Brown We characterized the TME across individual metastases using a et al., 2014). We integrated WES and RNA-seq data to predict in combination of computational pathology of digitally scanned silico putative neo-antigens (STAR Methods). Across the metasta- H&E slides (64 frozen sections of the tissue used for RNA extrac- ses, around 16% of expressed non-silent metastatic mutations tion and 102 FFPE tumor biopsies), manually scored immunohis- yielded 1 or more predicted neo-epitopes (with IC50 < 500 nM), tochemistry (IHC) of a set of immune markers (n = 102), and gene but only a small fraction (3%) of predicted neo-antigens originated expression data (RNA-seq, n = 64). from cancer drivers. Most predicted neo-antigens arise from We previously reported the use of digitized H&E slides and metastatic stem (56%) and metastatic clade mutations (36%) machine learning methods to classify cells within a tumor as can- (Figure 5A, top panel). Recently it was reported that loss of hetero- cer, stromal, or lymphocytes (Ali et al., 2016b). Using this zygosity (LOH) at the HLA locus occurs as a mechanism of immune approach, we analyzed 166 frozen and paraffin-embedded met- escape (McGranahan et al., 2017). Using the same method, we astatic and primary tissue sections (STAR Methods), and the identified clonal LOH (present in every single metastasis) in cases data revealed significant heterogeneity of cell numbers, frac- 330 (in all three HLA class I alleles) and 315 (in HLA-C) and subclo- tions, and densities across individual metastases, showing the nal LOH (present in a fraction of the metastases) in cases 288, 290, variable spatial architecture of the TME (Figures 6A and S6A; 298, and 308 (Figure 5A, bottom panel). In case 330, the neo-anti- SI6 in https://doi.org/10.17632/6cv77bry6m.1). In parallel IHC gens were significantly (p = 0.01, Wilcoxon rank-sum test) more (102 FFPE sections), semiquantitative analysis by expert pathol- commonly predicted to be presented by the lost HLA alleles. In ogists revealed variable numbers of CD4 and CD8 T cells per the remainder of cases, this difference was not significant for any surface area (Figure 6A, bottom panel). of the lost HLA alleles. Nevertheless, on average, 55.4% of pre- The patterns of immune infiltration can also be inferred using dicted neo-antigens associated with the lost HLA allele and, deconvolution of bulk gene expression (Hackl et al., 2016). In pri- hence, could not be presented directly by tumor cells. mary tumors, these patterns are variable across subtypes and We next asked, using a previously reported approach (Rooney associated with response to therapy and survival (Ali et al., et al., 2015), whether there was evidence for selected elimination 2016a). We performed these analyses across metastases from of immunogenic sub-clones across individual metastases. The 9 cases with available RNA-seq data. The immune cytolytic activ- method relies on determining the ratio of observed-to-expected ity score (Rooney et al., 2015) was highly variable across metas- neo-epitopes, and to estimate the distribution of the number of tases (Figure 6B). We further characterized TME expression expected neo-epitopes, we used simulated mutations that signatures using a recently reported deconvolution methodology mimic the observed mutations (STAR Methods). We generated (Charoentong et al., 2017). This tool provides normalized 100 datasets of simulated mutations for each sample and Z scores for a list of cancer immunity parameters, including 20 calculated the corresponding observed-to-expected ratios single factors (major histocompatibility complex [MHC] mole- both for individual metastases and after combining all metasta- cules, immunoinhibitors, and immunostimulators) and six cell ses from the same case or from the same organ across cases. types (STAR Methods). These Z scores were visualized as The results (Figure 5B) showed that only one single metastasis immunophenograms for each individual metastasis or used to (328-003, brain metastasis) had all 100 ratios below one (i.e., generate clustered heatmaps across metastases for each case. empirical p z 0.01), suggesting immunoediting. We also Inspection of the immunophenograms revealed variability of the estimated a null distribution for the mean observed-to-expected TME in metastases both between and within each of the cases ratio using 20 of the simulated mutation datasets, generating (Figure 6C; SI6 in https://doi.org/10.17632/6cv77bry6m.1). This for each one of them 100 simulated datasets and calcu- was mirrored in IHC analysis of a total of 14 TME markers, which lating the mean of the ratios (SI5 in https://doi.org/10.17632/ also revealed heterogeneous TME across metastases (SI6 in 6cv77bry6m.1). When grouping all metastases per case, none https://doi.org/10.17632/6cv77bry6m.1). Unsupervised hierar- had a mean ratio lower than expected (compared with the 20 chical clustering based on the Euclidean distance matrix of the replicates) (Figure S5A), and when lumping metastases by target Z scores across metastases showed that the immune parame- organ across cases, no organ site showed a trend suggesting ters tended to cluster naturally into major functional groups: immunoediting (Figure S5B). immunogenic or immune-suppressive (Figure 6D; Figure S6B). In summary, in disseminated lethal breast cancer, most of the The clustering of the individual metastasis in each case also predicted neo-antigens originate from mutations shared across had interesting features. In cases 288, 290, and 330, it appeared metastases, with only a small number being private to individual as if the TME clustering of the metastases mirrored their genomic metastasis. LOH of HLA alleles was observed in many metasta- clades (Figure 6D; Figure S6B). The TME clustering seen could ses (clonal in two cases), and non-synonymous mutations pre- simply reflect the metastatic organ site, but in other cases, dicted to be neo-antigenic were frequently associated with the metastases to the same organ had very distinct TMEs: meningeal lost HLA allele, suggesting tumor cell immune escape. Overall, and bone metastases in case 308, liver metastases in case 315, the number of predicted neo-antigens in each metastasis was and brain metastases in case 328 (Figure S6). We also examined only exceptionally lower than expected, suggesting that, in late the expression of immunomodulators (Tang et al., 2018; Thors- metastatic breast cancer, tumor cells are already in the escape son et al., 2018), and this also revealed metastases largely segre- phase of the immunoediting process, where cancer cells have gated by target organ and by genomic clade (Figure 6E). acquired the ability to circumvent immune recognition or In summary, multiple orthogonal methods congruently destruction (Schreiber et al., 2011). demonstrated that the immune TME is not uniform across 2700 Cell Reports 27, 2690–2708, May 28, 2019 A D B C E Figure 6. The Tumor Microenvironment Is Heterogeneous across Metastases (A) Median lymphocyte density (computational pathology of digitally scanned H&E slides) (top panel) and CD4 and CD8 T number per square millimeter (IHC staining) (bottom panel). (B) Cytolytic activity score across metastases based on transcript levels of granzyme A (GZMA) and perforin (PRF1). (C) Immunophenograms across metastases of case 288. Each immunophenogram is color-coded in the outer part of the wheel (red, positive Z score; blue, negative Z score) for each of the parameters and gray-scaled in the inner part of the wheel, with a weighted averaged Z score within the respective category. Z scales are shown in the bars. MHC, antigen processing; CP, checkpoints/immunomodulators; EC, effector cells; SC, suppressor cells. (D) Heatmaps depicting two-way unsupervised hierarchical clustering of immune parameters and metastases for patients 288 and 330. (E) Gene expression of immunomodulators from RNA-seq gene expression (76 genes from Thorsson et al., 2018). Z-scored transformed TPMs are plotted across all 64 RNA-seq metastases from 9 patients. Cell Reports 27, 2690–2708, May 28, 2019 2701 Figure 7. Analysis of the TCR Repertoire across Metastases (A) Boxplots of proportions of TCR reads classified as metastatic stem, metastatic clade, and metastatic private. Bars indicate significance of differences (not significant [NS], p > 0.05; ***p < 0.0005). (B) Boxplots of overlap coefficients between metastatic sites of TCR b chain nucleotide sequence repertoires. Data for case 308 is shown. (C) Boxplots showing the TCR clone sizes according to their stem, clade, or private status. *p < 0.05. (D) Clustering of TCR b chain CDR3 amino acid sequences using Jaccard distance across metastases. (E) Jaccard tree for the TCR b chain CDR3 amino acid sequence (top panel) and the WES phylogenetic tree from OncoNEM (bottom panel) for case 308. 2702 Cell Reports 27, 2690–2708, May 28, 2019 metastatic sites within a patient, although it can be both relatively the total repertoire (STAR Methods). The T cell clonality in total homogeneous in metastatic clades in some cases and different mononuclear cells, CD3+ T cells, naive T cells, and central across metastases to a particular organ in other cases. memory T cells from a healthy individual’s blood sample, pro- filed using the same method, were very distinct from breast The Repertoire of TCRs across Metastatic Sites cancer metastases T cells; although the former exhibited only The landscape of neo-antigens is thought to determine the low levels of expanded T cell clones with mean largest cluster immunogenicity of cancers and, in particular, the anti-tumor re- sizes of 3.93% (range, 1.03%–7.01%), most metastases sponses mediated by T cells. We therefore characterized the showed higher levels of specific T cell clonal expansion (mean repertoires of TCRs in tumor-infiltrating lymphocytes (TILs) largest cluster size, 14.97) (Figure S7C). Interestingly stem across metastases and integrated this with the genomic data and clade TCR clones were significantly more expanded than (Figure 7; Figure S7). TCR clones from a single metastasis (Figure 7C; SI7 in https:// We sequenced the antigen binding regions of TCRs of TILs us- doi.org/10.17632/6cv77bry6m.1), and the shared TCR clones ing direct amplification of the a and b TCR chains in RNA isolated showed significant differences in clonal representation be- from 70 metastases from 8 patients (STAR Methods). Each RNA tween metastatic sites, contrasting with relatively uniform sample was amplified and high-throughput-sequenced sepa- clones in normal samples (Figure S7D; SI7 in https://doi.org/ rately for a and b TCR chains. These sequences represented a 10.17632/6cv77bry6m.1). The fact that T cell clones that are diverse set of TCRs with a mean of 5,204 unique TCR sequences shared are also significantly enlarged compared with site-spe- per sample (a, 5,551; b, 4,857). The a and b chain V-J gene usages cific clones is suggestive of immune surveillance between were highly correlated in each metastatic sample and across all metastases. metastases, attesting to the quality of the data (Figure S7A). Clustering of TCR repertoires of TILs by the level of sharing of Public TCRs are defined as TCRs that are shared between un- CDR3 amino acid sequences (using the Jaccard index) revealed related individuals, and these are thought to bind shared pep- that TCR repertoires were distinct between each metastatic tides; for example, of viral origin (Li et al., 2012). Therefore, we breast cancer patient (Figure 7D). There was only a small degree focused our analyses on non-public TCRs by first removing of sharing of TCR CDR3 sequences between unrelated patients, both the TCR b chain sequences derived from unrelated individ- which may occur by chance at low frequencies, and high levels uals in the Adaptive database (Dean et al., 2015) and the TCR a of TCR CDR3 sequence sharing between metastases within and b chain sequences shared between any of the 8 cases pro- each patient. Indeed, within each patient, the unsupervised filed. This filtering step enriched for TCRs that bind patient- CDR3 TCR clustering nearly accurately segregated the metasta- specific antigens, including tumor neo-antigens (Figure S7B). ses by organ, and clustering was consistent with both a and b The non-public TIL TCRs were classified in each case as stem TCR chains (SI7 in https://doi.org/10.17632/6cv77bry6m.1). when shared between all metastatic sites in a patient, clade Using two methods for evaluating hierarchical tree structure, when shared between some sites, or private when found in a cophenetic correlation and Robinson-Foulds metric (STAR single metastasis. Interestingly, a significant fraction of the Methods), we showed that the TCR repertoires clustered signif- TCR repertoire was comprised of stem (mean a/b, 21.60%/ icantly by organ (p < 2.2e16) rather than by chance or by differ- 24.53%) and clade (mean a/b, 58.13%/49.67%) clonotypes (Fig- ences in TCR repertoire sampling depth. ure 7A). These data indicate that a significant proportion of T cell We noted that, in 4 of the autopsy cases (288, 290, 308, and clones in TILs from an individual metastasis either recirculate 315) where both WES and RNA TCR sequencing data were to other metastatic sites or that T cell clones recognizing the available across several metastases, the tumor OncoNEM same neo-antigen are recruited to several metastatic sites phylogenetic trees and the TCR Jaccard phylogenetic trees independently. Given the evidence for T cell sharing across indi- had remarkably similar structures (Figure 7E; SI7 in https://doi. vidual metastases, we quantified the degree of clonal sharing org/10.17632/6cv77bry6m.1). To formally test this, we used (STAR Methods) of CDR3 a and b clonotypes across metastases the cophenetic statistic to assess the correlation. This analysis and showed significant variation between sites and patients. confirmed that, in 3 of the 4 cases (288, 315, and 308), the This is exemplified by case 308, which revealed significant differ- TCR a chain trees significantly correlated with the genomic trees ences in TCR clonal sharing between and within metastatic (cophenetic correlation of 0.460, 0.235 and 0.518; p % 0.05), and target organs (Figure 7B) (see SI7 in https://doi.org/10.17632/ in 2 of the 4 cases (288 and 308), the TCR b chain trees also 6cv77bry6m.1 for the remainder of cases). Indeed, the TCR correlated with the genomic trees (cophenetic correlation of clonal overlap coefficient was significantly higher between me- 0.38 and 0.598, p % 0.08 [borderline] and p < 0.01, respectively) tastases within the same target organ than in metastases to (Figure 7E; SI7 in https://doi.org/10.17632/6cv77bry6m.1). different target organs. In addition, there was a high TCR clonal To further corroborate the robustness of the findings, the cor- overlap coefficient between bone and pericardium metastases, relation between the RNA-seq and the TCR repertoire datasets which formed a mediastinal mass in case 308. (STAR Methods; Figure S7E) was tested. The sum of the log10 We next assessed the clonal architecture of tumor-infiltrating transcripts per million (TPM) values for the four CD3 complex T cells at each site using sequence diversity measures (Bash- genes (a unique marker for all T cells) and the number of TCR ford-Rogers et al., 2013; Madi et al., 2017). This measures the reads (the sum of a and b chains) were highly statistically corre- unevenness of TCR clone sizes (clonality) within each site, lated (Pearson correlation, 0.717; p < 1.684e10). where each TCR clone is defined by a unique TCR sequence, In summary, these data provide a detailed view of the adaptive and its size is defined by the frequency of that sequence within immune response in metastatic cancer and reveal that the TCR Cell Reports 27, 2690–2708, May 28, 2019 2703 repertoire of TILs is highly diverse between cases; within each metastasis within a clade, the mutations were either shared case, a significant fraction of TIL TCRs are shared between me- with other clade members or private, and this suggests that tastases, suggesting immune surveillance between metastatic seeding occurs most likely in a single spreading event. Further- sites. The clonal prevalence of shared TIL TCRs in each individ- more, the genomic segregation of the metastases was nearly ual metastasis can be very different; the TCR repertoire clusters complete, and only a limited amount of cross-seeding between metastases by target organ within each case, and tumor phylo- clades was observed. Cross-seeding within a clade was rare. genetic trees appear to be correlated with TIL-TCR trees across These data suggest that metastatic spread occurs as a result metastases within a case. This correlation suggests co-evolution of a limited number of seeding events. between tumor diversity and T cell response across metastases. The classical view of the metastatic cascade has focused on seeding from single cells (Lambert et al., 2017). However, circu- DISCUSSION lating tumor cell (CTC) clusters occur in the blood of patients with metastatic breast cancer, and mouse models show that CTC The comprehensive molecular analysis of multi-regional metas- clusters are oligoclonal and, although rare compared with single tases collected from 10 autopsies of breast cancer patients CTCs, have 23- to 50-fold increased metastatic potential (Aceto subjected to multiple lines of therapy described here details et al., 2014). Our analysis using PyClone and CCF is consistent the heterogeneous landscape of genomic aberrations, TME fea- with metastases often being composed of communities of tures, and T cell adaptive immune responses in lethal cancer. genomic clones, as indicated by metastatic stem mutations The only comparably sized study (10 autopsy cases) was and metastatic clade mutations frequently being sub-clonal. recently published, but it limited its analysis to DNA sequencing Our previously published data regarding primary tumors (Shah of the available 41 paraffin-embedded tumor samples (Brown et al., 2012) and patient-derived tumor xenografts (Bruna et al., et al., 2017). 2016; Eirew et al., 2015) also appear to show that cancers are The genomic landscapes revealed that metastatic private communities of genomic clones. The polyclonal origin of breast driver mutations are relatively uncommon and that nearly all cancer metastases has also been reported by others (Hoadley driver CNAs and SNVs are shared across all (stem) or a subset et al., 2016) and has profound implications for the study of met- (clade) of the metastases. The normalized expression of the astatic biology and for devising therapeutic strategies. mutant allele is progressively lower in metastatic clade and met- The mutation signatures across the metastases are a reflec- astatic private mutations compared with metastatic stem muta- tion of the mutational processes operative during the life history tions, suggesting that, as metastases evolve from common of the cancer. We identified previously reported mutational sig- ancestors, they accumulate higher proportions of mutant alleles natures in all cases and evidence of residual mutations not that are passengers and, therefore, have lower or suppressed explained by any of the canonical mutation signatures described expression. to date (which have been almost exclusively derived from pri- In 6 cases, targeted deep sequencing of available primary mary cancers). Larger autopsy series and WGS data will be samples confirmed that a fraction of tested mutations were de- required to definitively establish whether metastases can accu- tected at high CCFs, confirming that the metastases originated mulate novel mutation signatures reflecting both their longer from the surgically resected tumor. However, in five of these, natural history and the combined scars of therapies and the deep sequencing failed to identify some metastatic stem muta- effects of immunoediting. tions and most of the metastatic clade mutations tested, The 10 cases analyzed were patients subject to multiple lines including cancer driver mutations. Given the high depth and of targeted therapy (hormone therapy and anti-Her2) and quality of the TS data, it is exceedingly unlikely that a trivial tech- chemotherapy, to which each patient had developed resistance. nical issue (allele dropout) could explain this result. The most The small cohort and the diversity of cancer treatments the pa- likely scenarios are either that metastases originate from a minor tients received do not allow confident identification of the mech- clone not sampled in the primary tumor because of spatial het- anisms of resistance. Nevertheless, case 308 had a canonical erogeneity or that metastatic stem and clade mutations could activating mutation of ESR1 (Robinson et al., 2013; Toy et al., be acquired by metastatic cancer cells after these leave the 2013) across all metastases, and this mutation, a likely mecha- breast. Indeed, in case Det52, we reported previously that nism of resistance to hormone therapy, would have been identi- metastatic stem mutations could be identified as a minor clone fied with a single metastatic biopsy. This contrasts with cases in the axillary lymph node, consistent with an original metastatic where two different mechanisms for hormone therapy resistance ancestor being present at that site (Murtaza et al., 2015). were identified: losing ER expression in some metastases and Whether that is the rule or whether other sites (e.g., the bone 11q13 amplification in ER-positive metastases in case 288 and marrow) (Harper et al., 2016; Hosseini et al., 2016) could have losing ER expression in some metastases and an activating a similar metastases seeder role will require studying larger ESR1 mutation in ER-positive metastases in case 290. These numbers of cases. distinct forms of resistance in both cases correspond to their The genomic phylogenetic analysis showed clearly that the metastatic phylogenetic clades and imply that for both to be multiple metastases in each case (with one exception) grouped identified would require at least two metastatic biopsies. into a small number of clades (up to three). These clades were The combined analysis of somatic cancer aberrations, TME populated by a common seeder, itself a descendant of the deconvolution, predicted tumor neo-antigens, and TCR reper- original metastatic ancestor, and clades were anatomically toire from lethal metastatic breast cancer autopsies afforded distributed to one or more target organs. In each individual us a unique opportunity to analyze the interactions between 2704 Cell Reports 27, 2690–2708, May 28, 2019 the malignant and TME compartments across multi-regional me- implications for T cell-based immunotherapies. For example, tastases. Most of the predicted neo-antigens arose from meta- peptides that are derived from neo-antigens and are used for static stem and clade mutations, and this result concurs with vaccination should be tailored differently for different metastatic what has been reported in lung cancer (McGranahan et al., sites. Similarly, adoptive T cell therapies based on patient T cells 2016). Immune selection, as evidenced by depletion of neo-epi- or engineered receptors should also take into account differ- topes, was scarcely seen across metastases, suggesting that in ences in TCR composition and reactivity between metastatic advanced breast cancer, most of the metastases are in the sites. On the other hand, the presence of TCRs shared across escape phase of immunoediting (Schreiber et al., 2011). metastases and the fact that most neo-antigens are also shared TME composition and its spatial architecture were also het- across all or subsets of metastases could be translated into erogeneous across metastases from individual cases. This either T cell-based therapies or modulation of immune check- TME diversity had no direct correlation with evidence for differ- points that would elicit an effective anti-tumor response across ential immunoediting and has also profound implications for all metastases. immune checkpoint inhibitor therapy because, for example, The extensive profiling of multi-regional metastases in lethal PD1 and PDL1 expression was different across metastases. breast cancers resistant to several lines of therapy analyzed The occurrence of HLA LOH suggests that tumor cells also here has provided a unique glimpse into how metastases prop- evolve to avoid presenting neo-antigens, and this might also agate and evolve, how drug resistance mechanisms vary, how contribute to immune escape. their predicted neo-antigen landscapes look, how they shape A major determinant of the ability of the adaptive immune sys- and/or are shaped by the TME, and how the T cell adaptive tem to eliminate tumors is the diversity of the TCR repertoire in immune response appears to co-evolve with the metastatic ge- TILs. Although the level of sharing of TCRs between patients nomes. These data should motivate the research community to was minimal, an important observation was the large fraction of consider launching a lethal cancer genome project. This project, TCRs either shared across all the metastases (stem TCRs) or in across common human cancers in a sufficiently large number of at least 2 metastases (clade TCRs) within a case. This finding sug- autopsy cases from patients with detailed therapy exposure gests that specific TCRs reacting to tumor neo-antigens, which histories, will produce the detailed integrated maps required are mostly stem and clade, are present across metastatic sites. for understanding resistance to therapy and escape from cancer Shared TCRs could be clonally dominant in an individual metas- immunoediting. tasis and minor clones at other metastatic sites. We could not identify any correlation between clonal inequality (measured by STAR+METHODS either Gini or Shannon index) or intra-tumor heterogeneity, muta- tion, or neo-antigen burden (data not shown); however, other Detailed methods are provided in the online version of this paper factors, such as chemokines, may influence the migration and and include the following: proliferation of T cell clones between sites. TCR similarity was higher in metastases within a given meta- d KEY RESOURCES TABLE static organ. This observation led us to cluster the metastases d CONTACT FOR REAGENT AND RESOURCE SHARING across all patients and organs based on TCR diversity and/or d EXPERIMENTAL MODEL AND SUBJECT DETAILS similarity, revealing each patient clustered separately from all B Patients others and, within an individual patient, nearly perfect metastatic d METHOD DETAILS organ segregation. Very distinct TCR repertoires between pa- B Sample Nomenclature tients were an expected result. However, within an individual B Nucleic acid extraction patient, TCR repertoires in metastases to the same organ were B WES AND sWGS ANALYSES more similar. This observation was robust and highly statistically B RNA-sequencing analysis significant, suggesting that TCR repertoires in TILs are tuned to B Neo-antigen simulations target organs where metastases tend to share a common B T cell receptor (TCR) analysis genomic ancestor. Indeed, genomic phylogenetic trees and d QUANTIFICATION AND STATISTICAL ANALYSIS TCR repertoire Jaccard trees showed remarkably similar archi- d DATA AND SOFTWARE AVAILABILITY tectures in the 4 cases where we had parallel WES and TCR B Software sequencing data from more than 4 metastases. We tested this B Data Resources more formally using tools developed by ecologists, and the correlation was statistically significant in 3 of the cases. This is SUPPLEMENTAL INFORMATION suggestive of co-evolution of cancer genomes and the TCR rep- ertoires of the same metastases, providing unique patient-based Supplemental Information can be found online at https://doi.org/10.1016/j. evidence for the cancer immunoediting hypothesis. An alterna- celrep.2019.04.098. tive explanation that is not contradictory is that tissue-resident T cells, which have tissue-specific TCR repertoires, infiltrate ACKNOWLEDGMENTS metastases, giving rise to the observed TCR similarities. Such This work was supported by Cancer Research UK. L.D.M.A. was partly funded differences in the composition of TCR repertoires between tis- by Asociación Española contra el Cáncer. The authors acknowledge Asocia- sues can be part of the forces that influence clonal evolution of ción Española contra el Cáncer, Cellex foundation, and the clinical team at metastatic cancer cells. These findings can have important the Breast Cancer Unit of Vall d’Hebron University Hospital/Institute of Cell Reports 27, 2690–2708, May 28, 2019 2705 Oncology. We are very grateful to the 10 women and their families who Callari, M., Sammut, S.J., De Mattos-Arruda, L., Bruna, A., Rueda, O.M., Chin, donated the autopsy and tumor samples with linked clinical data. S.F., and Caldas, C. (2017). Intersect-then-combine approach: improving the performance of somatic variant calling in whole exome sequencing data using AUTHOR CONTRIBUTIONS multiple aligners and callers. Genome Med. 9, 35. Cancer Genome Atlas, N.; Cancer Genome Atlas Network (2012). Comprehen- L.D.M.-A., S.-J.S., J.S., and C.C. conceived the experiments and wrote the sive molecular portraits of human breast tumours. Nature 490, 61–70. manuscript. L.D.M.-A. and S.-J.S. performed the experiments and analyses Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, with assistance from E.M.R. and F.M. (tumor phylogenies); R.B.-R., E.G., Y.T., D., Hackl, H., and Trajanoski, Z. (2017). Pan-cancer Immunogenomic Analyses G.V., D.R., and N.F. (TCR analysis); H.M., B.P., O.M.R., S.-F.C., T.C.-C., D.T., Reveal Genotype-Immunophenotype Relationships and Predictors of and M.M. (shallow WGS, whole-exome sequencing, RNA sequencing, and tar- Response to Checkpoint Blockade. Cell Rep. 18, 248–262. geted deep sequencing data generation and analysis); S.M., T.D.A., and S.N.-Z. (mutation signatures); Y.M. and G.G. (neo-antigen observed over expected Cibulskis, K., Lawrence, M.S., Carter, S.L., Sivachenko, A., Jaffe, D., Sougnez, analysis); and R.M., A.A., H.R.A., W.C., A.D., N.C., E.M.-S., V.P., S.R. y C., C., Gabriel, S., Meyerson, M., Lander, E.S., and Getz, G. (2013). Sensitive and J.C. (clinical and pathology data). detection of somatic point mutations in impure and heterogeneous cancer samples. Nat. Biotechnol. 31, 213–219. DECLARATION OF INTERESTS Ciriello, G., Miller, M.L., Aksoy, B.A., Senbabaoglu, Y., Schultz, N., and Sander, C. (2013). Emerging landscape of oncogenic signatures across human C.C. is a member of the AstraZeneca External Science Panel and has research cancers. Nat. Genet. 45, 1127–1133. grants from Roche, Genentech, AstraZeneca, and Servier that are adminis- Curtis, C., Shah, S.P., Chin, S.F., Turashvili, G., Rueda, O.M., Dunning, M.J., tered by the University of Cambridge. Speed, D., Lynch, A.G., Samarajiwa, S., Yuan, Y., et al.; METABRIC Group (2012). The genomic and transcriptomic architecture of 2,000 breast tumours Received: January 22, 2018 reveals novel subgroups. Nature 486, 346–352. Revised: December 13, 2018 De Mattos-Arruda, L., Mayor, R., Ng, C.K., Weigelt, B., Martı́nez-Ricarte, F., Accepted: April 22, 2019 Torrejon, D., Oliveira, M., Arias, A., Raventos, C., Tang, J., et al. (2015). Cere- Published: May 28, 2019 brospinal fluid-derived circulating tumour DNA better represents the genomic alterations of brain tumours than plasma. Nat. Commun. 6, 8839. REFERENCES Dean, J., Emerson, R.O., Vignali, M., Sherwood, A.M., Rieder, M.J., Carlson, C.S., and Robins, H.S. (2015). Annotation of pseudogenic gene segments by Aceto, N., Bardia, A., Miyamoto, D.T., Donaldson, M.C., Wittner, B.S., massively parallel sequencing of rearranged lymphocyte receptor loci. Spencer, J.A., Yu, M., Pely, A., Engstrom, A., Zhu, H., et al. (2014). Circulating Genome Med. 7, 123. tumor cell clusters are oligoclonal precursors of breast cancer metastasis. Cell 158, 1110–1122. DePristo, M.A., Banks, E., Poplin, R., Garimella, K.V., Maguire, J.R., Hartl, C., Philippakis, A.A., del Angel, G., Rivas, M.A., Hanna, M., et al. (2011). A frame- Alexandrov, L.B., Nik-Zainal, S., Wedge, D.C., Campbell, P.J., and Stratton, work for variation discovery and genotyping using next-generation DNA M.R. (2013). Deciphering signatures of mutational processes operative in sequencing data. Nat. Genet. 43, 491–498. human cancer. Cell Rep. 3, 246–259. Ali, H.R., Rueda, O.M., Chin, S.F., Curtis, C., Dunning, M.J., Aparicio, S.A., and Ding, L., Ellis, M.J., Li, S., Larson, D.E., Chen, K., Wallis, J.W., Harris, C.C., Caldas, C. (2014). Genome-driven integrated classification of breast cancer McLellan, M.D., Fulton, R.S., Fulton, L.L., et al. (2010). Genome remodelling validated in over 7,500 samples. Genome Biol. 15, 431. in a basal-like breast cancer metastasis and xenograft. Nature 464, 999–1005. Ali, H.R., Chlon, L., Pharoah, P.D., Markowetz, F., and Caldas, C. (2016a). Pat- Dvinge, H., Git, A., Gräf, S., Salmon-Divon, M., Curtis, C., Sottoriva, A., Zhao, terns of Immune Infiltration in Breast Cancer and Their Clinical Implications: Y., Hirst, M., Armisen, J., Miska, E.A., et al. (2013). The shaping and functional A Gene-Expression-Based Retrospective Study. PLoS Med. 13, e1002194. consequences of the microRNA landscape in breast cancer. Nature 497, 378–382. Ali, H.R., Dariush, A., Provenzano, E., Bardwell, H., Abraham, J.E., Iddawela, M., Vallier, A.L., Hiller, L., Dunn, J.A., Bowden, S.J., et al. (2016b). Computa- Eirew, P., Steif, A., Khattra, J., Ha, G., Yap, D., Farahani, H., Gelmon, K., Chia, tional pathology of pre-treatment biopsies identifies lymphocyte density as a S., Mar, C., Wan, A., et al. (2015). Dynamics of genomic clones in breast cancer predictor of response to neoadjuvant chemotherapy in breast cancer. Breast patient xenografts at single-cell resolution. Nature 518, 422–426. Cancer Res. 18, 21. Flensburg, C., Sargeant, T., Oshlack, A., and Majewski, I. (2018). SuperFreq: Altschul, S.F., Gish, W., Miller, W., Myers, E.W., and Lipman, D.J. (1990). Basic Integrated mutation detection and clonal tracking in cancer. bioRxiv. https:// local alignment search tool. J. Mol. Biol. 215, 403–410. doi.org/10.1101/380097. Bashford-Rogers, R.J., Palser, A.L., Huntly, B.J., Rance, R., Vassiliou, G.S., Galili, T. (2015). dendextend: an R package for visualizing, adjusting and Follows, G.A., and Kellam, P. (2013). Network properties derived from deep comparing trees of hierarchical clustering. Bioinformatics 31, 3718–3720. sequencing of human B-cell receptor repertoires delineate B-cell populations. Gaujoux, R., and Seoighe, C. (2010). A flexible R package for nonnegative Genome Res. 23, 1874–1884. matrix factorization. BMC Bioinformatics 11, 367. Bronner, I.F., Quail, M.A., Turner, D.J., and Swerdlow, H. (2014). Improved Hackl, H., Charoentong, P., Finotello, F., and Trajanoski, Z. (2016). Computa- Protocols for Illumina Sequencing. Curr. Protoc. Hum. Genet. 80, 18.2.1-42. tional genomics tools for dissecting tumour-immune cell interactions. Nat. Brown, D., Smeets, D., Szekely, B., Larsimont, D., Szasz, A.M., Adnet, P.Y., Rev. Genet. 17, 441–458. Rothe, F., Rouas, G., Nagy, Z.I., Farago, Z., et al. (2017). Phylogenetic analysis Harper, K.L., Sosa, M.S., Entenberg, D., Hosseini, H., Cheung, J.F., Nobre, R., of metastatic progression in breast cancer using somatic mutations and copy Avivar-Valderas, A., Nagi, C., Girnius, N., Davis, R.J., et al. (2016). Mechanism number aberrations. Nat. Commun. 8, 14944. of early dissemination and metastasis in Her2+ mammary cancer. Nature 540, Brown, S.D., Warren, R.L., Gibb, E.A., Martin, S.D., Spinelli, J.J., Nelson, B.H., 588–592. and Holt, R.A. (2014). Neo-antigens predicted by tumor genome meta-analysis Hoadley, K.A., Siegel, M.B., Kanchi, K.L., Miller, C.A., Ding, L., Zhao, W., He, correlate with increased patient survival. Genome Res. 24, 743–750. X., Parker, J.S., Wendl, M.C., Fulton, R.S., et al. (2016). Tumor Evolution in Two Bruna, A., Rueda, O.M., Greenwood, W., Batra, A.S., Callari, M., Batra, R.N., Patients with Basal-like Breast Cancer: A Retrospective Genomics Study of Pogrebniak, K., Sandoval, J., Cassidy, J.W., Tufegdzic-Vidakovic, A., et al. Multiple Metastases. PLoS Med. 13, e1002174. (2016). A Biobank of Breast Cancer Explants with Preserved Intra-tumor Het- , M.M.S., Hoffmann, M., Harper, K.L., Sosa, M.S., Hosseini, H., Obradovic erogeneity to Screen Anticancer Compounds. Cell 167, 260–274.e22. Werner-Klein, M., Nanduri, L.K., Werno, C., Ehrl, C., Maneck, M., et al. 2706 Cell Reports 27, 2690–2708, May 28, 2019 (2016). Early dissemination seeds metastasis in breast cancer. Nature 540, Mermel, C.H., Schumacher, S.E., Hill, B., Meyerson, M.L., Beroukhim, R., and 552–558. Getz, G. (2011). GISTIC2.0 facilitates sensitive and confident localization of the Huber, W., Carey, V.J., Gentleman, R., Anders, S., Carlson, M., Carvalho, B.S., targets of focal somatic copy-number alteration in human cancers. Genome Bravo, H.C., Davis, S., Gatto, L., Girke, T., et al. (2015). Orchestrating high- Biol. 12, R41. throughput genomic analysis with Bioconductor. Nat. Methods 12, 115–121. Munson, D.J., Egelston, C.A., Chiotti, K.E., Parra, Z.E., Bruno, T.C., Moore, Hundal, J., Carreno, B.M., Petti, A.A., Linette, G.P., Griffith, O.L., Mardis, E.R., B.L., Nakano, T.A., Simons, D.L., Jimenez, G., Yim, J.H., et al. (2016). Identifi- and Griffith, M. (2016). pVAC-Seq: A genome-guided in silico approach to cation of shared TCR sequences from T cells in human breast cancer using identifying tumor neoantigens. Genome Med. 8, 11. emulsion RT-PCR. Proc. Natl. Acad. Sci. USA 113, 8272–8277. Juric, D., Castel, P., Griffith, M., Griffith, O.L., Won, H.H., Ellis, H., Ebbesen, Murtaza, M., Dawson, S.J., Pogrebniak, K., Rueda, O.M., Provenzano, E., S.H., Ainscough, B.J., Ramu, A., Iyer, G., et al. (2015). Convergent loss of Grant, J., Chin, S.F., Tsui, D.W., Marass, F., Gale, D., et al. (2015). Multifocal PTEN leads to clinical resistance to a PI(3)Ka inhibitor. Nature 518, 240–244. clonal evolution characterized using circulating tumour DNA in a case of met- astatic breast cancer. Nat. Commun. 6, 8760. Koboldt, D.C., Zhang, Q., Larson, D.E., Shen, D., McLellan, M.D., Lin, L., Miller, C.A., Mardis, E.R., Ding, L., and Wilson, R.K. (2012). VarScan 2: somatic mu- Nik-Zainal, S., Alexandrov, L.B., Wedge, D.C., Van Loo, P., Greenman, C.D., tation and copy number alteration discovery in cancer by exome sequencing. Raine, K., Jones, D., Hinton, J., Marshall, J., Stebbings, L.A., et al.; Breast Can- Genome Res. 22, 568–576. cer Working Group of the International Cancer Genome Consortium (2012a). Koressaar, T., and Remm, M. (2007). Enhancements and modifications of Mutational processes molding the genomes of 21 breast cancers. Cell 149, primer design program Primer3. Bioinformatics 23, 1289–1291. 979–993. Lambert, A.W., Pattabiraman, D.R., and Weinberg, R.A. (2017). Emerging Nik-Zainal, S., Van Loo, P., Wedge, D.C., Alexandrov, L.B., Greenman, C.D., Biological Principles of Metastasis. Cell 168, 670–691. Lau, K.W., Raine, K., Jones, D., Marshall, J., Ramakrishna, M., et al.; Breast Cancer Working Group of the International Cancer Genome Consortium Lawrence, M.S., Stojanov, P., Polak, P., Kryukov, G.V., Cibulskis, K., Siva- (2012b). The life history of 21 breast cancers. Cell 149, 994–1007. chenko, A., Carter, S.L., Stewart, C., Mermel, C.H., Roberts, S.A., et al. (2013). Mutational heterogeneity in cancer and the search for new cancer- Nik-Zainal, S., Davies, H., Staaf, J., Ramakrishna, M., Glodzik, D., Zou, X., associated genes. Nature 499, 214–218. Martincorena, I., Alexandrov, L.B., Martin, S., Wedge, D.C., et al. (2016). Land- scape of somatic mutations in 560 breast cancer whole-genome sequences. Lawrence, M.S., Stojanov, P., Mermel, C.H., Robinson, J.T., Garraway, L.A., Nature 534, 47–54. Golub, T.R., Meyerson, M., Gabriel, S.B., Lander, E.S., and Getz, G. (2014). Discovery and saturation analysis of cancer genes across 21 tumour types. Olshen, A.B., Venkatraman, E.S., Lucito, R., and Wigler, M. (2004). Circular Nature 505, 495–501. binary segmentation for the analysis of array-based DNA copy number data. Biostatistics 5, 557–572. Lefebvre, C., Bachelot, T., Filleron, T., Pedrero, M., Campone, M., Soria, J.C., Massard, C., Lévy, C., Arnedos, M., Lacroix-Triki, M., et al. (2016). Mutational Paradis, E., Claude, J., and Strimmer, K. (2004). APE: Analyses of Phyloge- Profile of Metastatic Breast Cancers: A Retrospective Analysis. PLoS Med. 13, netics and Evolution in R language. Bioinformatics 20, 289–290. e1002201. Park, J.H., Jang, M., Tarhan, Y.E., Katagiri, T., Sasa, M., Miyoshi, Y., Kalari, Lefranc, M.P. (2011). IMGT unique numbering for the variable (V), constant (C), K.R., Suman, V.J., Weinshilboum, R., Wang, L., et al. (2016). Clonal expansion and groove (G) domains of IG, TR, MH, IgSF, and MhSF. Cold Spring Harb. of antitumor T cells in breast cancer correlates with response to neoadjuvant Protoc. 2011, 633–642. chemotherapy. Int. J. Oncol. 49, 471–478. Li, H., and Durbin, R. (2009). Fast and accurate short read alignment with Bur- Parker, J.S., Mullins, M., Cheang, M.C., Leung, S., Voduc, D., Vickery, T., Da- rows-Wheeler transform. Bioinformatics 25, 1754–1760. vies, S., Fauron, C., He, X., Hu, Z., et al. (2009). Supervised risk predictor of breast cancer based on intrinsic subtypes. J. Clin. Oncol. 27, 1160–1167. Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., Abecasis, G., and Durbin, R.; 1000 Genome Project Data Processing Sub- Pereira, B., Chin, S.F., Rueda, O.M., Vollan, H.K., Provenzano, E., Bardwell, group (2009). The Sequence Alignment/Map format and SAMtools. Bioinfor- H.A., Pugh, M., Jones, L., Russell, R., Sammut, S.J., et al. (2016). The somatic matics 25, 2078–2079. mutation profiles of 2,433 breast cancers refines their genomic and transcrip- tomic landscapes. Nat. Commun. 7, 11479. Li, H., Ye, C., Ji, G., and Han, J. (2012). Determinants of public T cell re- sponses. Cell Res. 22, 33–42. Popic, V., Salari, R., Hajirasouliha, I., Kashef-Haghighi, D., West, R.B., and Batzoglou, S. (2015). Fast and scalable inference of multi-sample cancer line- Madi, A., Poran, A., Shifrut, E., Reich-Zeliger, S., Greenstein, E., Zaretsky, I., ages. Genome Biol. 16, 91. Arnon, T., Laethem, F.V., Singer, A., Lu, J., et al. (2017). T cell receptor reper- toires of mice and humans are clustered in similarity networks around Rainer, J. (2016). EnsDb.Hsapiens.v75: Ensembl based annotation package, R conserved public CDR3 sequences. eLife 6, e22057. package version 2.1.0. Martincorena, I., Raine, K.M., Gerstung, M., Dawson, K.J., Haase, K., Van Loo, Reiter, J.G., Makohon-Moore, A.P., Gerold, J.M., Bozic, I., Chatterjee, K., Ia- P., Davies, H., Stratton, M.R., and Campbell, P.J. (2017). Universal Patterns of cobuzio-Donahue, C.A., Vogelstein, B., and Nowak, M.A. (2017). Reconstruct- Selection in Cancer and Somatic Tissues. Cell 171, 1029–1041.e21. ing metastatic seeding patterns of human cancers. Nat. Commun. 8, 14114. McGranahan, N., Favero, F., de Bruin, E.C., Birkbak, N.J., Szallasi, Z., and Robinson, J.T., Thorvaldsdóttir, H., Winckler, W., Guttman, M., Lander, E.S., Swanton, C. (2015). Clonal status of actionable driver events and the timing Getz, G., and Mesirov, J.P. (2011). Integrative genomics viewer. Nat. Bio- of mutational processes in cancer evolution. Sci. Transl. Med. 7, 283ra54. technol. 29, 24–26. McGranahan, N., Furness, A.J., Rosenthal, R., Ramskov, S., Lyngaa, R., Saini, Robinson, D.R., Wu, Y.M., Vats, P., Su, F., Lonigro, R.J., Cao, X., Kalyana- S.K., Jamal-Hanjani, M., Wilson, G.A., Birkbak, N.J., Hiley, C.T., et al. (2016). Sundaram, S., Wang, R., Ning, Y., Hodges, L., et al. (2013). Activating Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune ESR1 mutations in hormone-resistant metastatic breast cancer. Nat. Genet. checkpoint blockade. Science 351, 1463–1469. 45, 1446–1451. McGranahan, N., Rosenthal, R., Hiley, C.T., Rowan, A.J., Watkins, T.B.K., Wil- Rooney, M.S., Shukla, S.A., Wu, C.J., Getz, G., and Hacohen, N. (2015). son, G.A., Birkbak, N.J., Veeriah, S., Van Loo, P., Herrero, J., et al. (2017). Molecular and genetic properties of tumors associated with local immune Allele-Specific HLA Loss and Immune Escape in Lung Cancer Evolution. Cell cytolytic activity. Cell 160, 48–61. 171, 1259–1271.e11. Rosenthal, R., McGranahan, N., Herrero, J., Taylor, B.S., and Swanton, C. McLaren, W., Gil, L., Hunt, S.E., Riat, H.S., Ritchie, G.R., Thormann, A., Flicek, (2016). DeconstructSigs: delineating mutational processes in single tumors P., and Cunningham, F. (2016). The Ensembl Variant Effect Predictor. Genome distinguishes DNA repair deficiencies and patterns of carcinoma evolution. Biol. 17, 122. Genome Biol. 17, 31. Cell Reports 27, 2690–2708, May 28, 2019 2707 Ross, E.M., and Markowetz, F. (2016). OncoNEM: inferring tumor evolution Tang, J., Shalabi, A., and Hubbard-Lucey, V.M. (2018). Comprehensive anal- from single-cell sequencing data. Genome Biol. 17, 69. ysis of the clinical immuno-oncology landscape. Ann. Oncol. 29, 84–91. Ross, E.M., Haase, K., Van Loo, P., and Markowetz, F. (2017). Allele-specific Thorsson, V., Gibbs, D.L., Brown, S.D., Wolf, D., Bortone, D.S., Ou Yang, T.H., multi-sample copy number segmentation. bioRxiv. https://doi.org/10.1101/ Porta-Pardo, E., Gao, G.F., Plaisier, C.L., Eddy, J.A., et al. (2018). The Immune 166017. Landscape of Cancer. Immunity 48, 812–830.e14. Roth, A., Khattra, J., Yap, D., Wan, A., Laks, E., Biele, J., Ha, G., Aparicio, S., Toy, W., Shen, Y., Won, H., Green, B., Sakr, R.A., Will, M., Li, Z., Gala, K., Fan- Bouchard-Co ^ té, A., and Shah, S.P. (2014). PyClone: statistical inference of ning, S., King, T.A., et al. (2013). ESR1 ligand-binding domain mutations in clonal population structure in cancer. Nat. Methods 11, 396–398. hormone-resistant breast cancer. Nat. Genet. 45, 1439–1445. Saunders, C.T., Wong, W.S., Swamy, S., Becq, J., Murray, L.J., and Chee- Untergasser, A., Cutcutache, I., Koressaar, T., Ye, J., Faircloth, B.C., Remm, tham, R.K. (2012). Strelka: accurate somatic small-variant calling from M., and Rozen, S.G. (2012). Primer3–new capabilities and interfaces. Nucleic sequenced tumor-normal sample pairs. Bioinformatics 28, 1811–1817. Acids Res. 40, e115. Savas, P., Teo, Z.L., Lefevre, C., Flensburg, C., Caramia, F., Alsop, K., Man- Van Loo, P., Nordgard, S.H., Lingjærde, O.C., Russnes, H.G., Rye, I.H., Sun, sour, M., Francis, P.A., Thorne, H.A., Silva, M.J., et al. (2016). The Subclonal W., Weigman, V.J., Marynen, P., Zetterberg, A., Naume, B., et al. (2010). Architecture of Metastatic Breast Cancer: Results from a Prospective Commu- Allele-specific copy number analysis of tumors. Proc. Natl. Acad. Sci. USA nity-Based Rapid Autopsy Program ‘‘CASCADE’’. PLoS Med. 13, e1002204. 107, 16910–16915. Scheinin, I., Sie, D., Bengtsson, H., van de Wiel, M.A., Olshen, A.B., van Thuijl, H.F., van Essen, H.F., Eijk, P.P., Rustenburg, F., Meijer, G.A., et al. (2014). DNA Vogelstein, B., Papadopoulos, N., Velculescu, V.E., Zhou, S., Diaz, L.A., Jr., copy number analysis of fresh and formalin-fixed specimens by shallow and Kinzler, K.W. (2013). Cancer genome landscapes. Science 339, 1546– whole-genome sequencing with identification and exclusion of problematic re- 1558. gions in the genome assembly. Genome Res. 24, 2022–2032. Wang, K., Li, M., and Hakonarson, H. (2010). ANNOVAR: functional annotation Schreiber, R.D., Old, L.J., and Smyth, M.J. (2011). Cancer immunoediting: of genetic variants from high-throughput sequencing data. Nucleic Acids Res. integrating immunity’s roles in cancer suppression and promotion. Science 38, e164. 331, 1565–1570. Wang, T., Wang, C., Wu, J., He, C., Zhang, W., Liu, J., Zhang, R., Lv, Y., Li, Y., Schwarz, R.F., Trinh, A., Sipos, B., Brenton, J.D., Goldman, N., and Marko- Zeng, X., et al. (2017). The Different T-cell Receptor Repertoires in Breast wetz, F. (2014). Phylogenetic quantification of intra-tumour heterogeneity. Cancer Tumors, Draining Lymph Nodes, and Adjacent Tissues. Cancer Immu- PLoS Comput. Biol. 10, e1003535. nol. Res. 5, 148–156. Shah, S.P., Morin, R.D., Khattra, J., Prentice, L., Pugh, T., Burleigh, A., Dela- Watson, S.J., Welkers, M.R., Depledge, D.P., Coulter, E., Breuer, J.M., de ney, A., Gelmon, K., Guliany, R., Senz, J., et al. (2009). Mutational evolution Jong, M.D., and Kellam, P. (2013). Viral population analysis and minority- in a lobular breast tumour profiled at single nucleotide resolution. Nature variant detection using short read next-generation sequencing. Philos. Trans. 461, 809–813. R. Soc. Lond. B Biol. Sci. 368, 20120205. Shah, S.P., Roth, A., Goya, R., Oloumi, A., Ha, G., Zhao, Y., Turashvili, G., Ding, Wolff, A.C., Hammond, M.E., Hicks, D.G., Dowsett, M., McShane, L.M., Alli- J., Tse, K., Haffari, G., et al. (2012). The clonal and mutational evolution spec- son, K.H., Allred, D.C., Bartlett, J.M., Bilous, M., Fitzgibbons, P., et al.; Amer- trum of primary triple-negative breast cancers. Nature 486, 395–399. ican Society of Clinical Oncology; College of American Pathologists (2013). Shen, Z., Qu, W., Wang, W., Lu, Y., Wu, Y., Li, Z., Hang, X., Wang, X., Zhao, D., Recommendations for human epidermal growth factor receptor 2 testing in and Zhang, C. (2010). MPprimer: a program for reliable multiplex PCR primer breast cancer: American Society of Clinical Oncology/College of American design. BMC Bioinformatics 11, 143. Pathologists clinical practice guideline update. J. Clin. Oncol. 31, 3997–4013. Shukla, S.A., Rooney, M.S., Rajasagi, M., Tiao, G., Dixon, P.M., Lawrence, Yates, L.R., Gerstung, M., Knappskog, S., Desmedt, C., Gundem, G., Van Loo, M.S., Stevens, J., Lane, W.J., Dellagatta, J.L., Steelman, S., et al. (2015). P., Aas, T., Alexandrov, L.B., Larsimont, D., Davies, H., et al. (2015). Subclonal Comprehensive analysis of cancer-associated somatic mutations in class I diversification of primary breast cancer revealed by multiregion sequencing. HLA genes. Nat. Biotechnol. 33, 1152–1158. Nat. Med. 21, 751–759. Soetaert, K., Van den Meersche, K., and van Oevelen, D. (2009). limSolve: Yates, L.R., Knappskog, S., Wedge, D., Farmery, J.H.R., Gonzalez, S., Martin- Solving Linear Inverse Models, R-package version 1.5.1. corena, I., Alexandrov, L.B., Van Loo, P., Haugland, H.K., Lilleng, P.K., et al. Stamatakis, A. (2014). RAxML version 8: a tool for phylogenetic analysis and (2017). Genomic Evolution of Breast Cancer Metastasis and Relapse. Cancer post-analysis of large phylogenies. Bioinformatics 30, 1312–1313. Cell 32, 169–184.e167. 2708 Cell Reports 27, 2690–2708, May 28, 2019 STAR+METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Antibodies CD68 Antibody Novocastra Cat# NCL-CD68; RRID: AB_563623 CD3 (Clone SP7) Thermo Scientific Cat# RM-9107-S; RRID: AB_149922 CD19 Antibody Abcam Cat# ab134114 Anti FOXP3 Antibody Abcam Cat# ab20034; RRID: AB_445284 CD8 Monoclonal Antibody Thermo Scientific Cat# RM-9116-S; RRID: AB_149960 Anti-IL3RA Atlas Cat# HPA003539; RRID: AB_1078438 Anti-IDO1 Antibody Atlas Cat# HPA027772; RRID: AB_1846222 CD4 Antibody Novocastra Cat# CD4-368-L-CE CD56 Antibody Novocastra Cat# CD56-504-L-CE CD1A Antibody Novocastra Cat# CD1A-235-L-CE Mast Cell Tryptase Antibody DAKO Cat# M7052; RRID: AB_2206478 CD45RO Antibody DAKO Cat# M0742; RRID: AB_2237910 CD38 Antibody Novocastra Cat# NCL-L-CD38-290; RRID: AB_563555 PDL1 Antibody Cell Signaling Cat# 13684; RRID: AB_2687655 Technologies ER Novocastra Cat# NCL-ER-6F11/2; RRID: AB_876939 PR Dako Cat# M3569; RRID: AB_2532076 HER2 Abbott Diagnostics Cat# 06N46-035 Raindance Source Chips Raindance Techologies Cat# 20-04410 (BioRad) TaqMan Genotyping Master Mix Thermo Fisher Cat# 4371353 SPRIselect Reagent Beckman Coulter Cat# B23318 Critical Commercial Assays DNeasy Blood and Tissue Kit QIAGEN Cat# 69506 QIAamp DNA Mini Kit QIAGEN Cat# 51306 MiRneasy mini kit QIAGEN Cat# 217004 GoTaq DNA polymerase Promega Cat# M7808 GoTaq Flexi DNA polymerase Promega Cat# M7808 SuperScript III Reverse Transcriptase ThermoFisher Scientific Cat#18080093 Illumina Nextera Rapid Capture Exome kit Illumina Cat# FC-140-1003 Quant-IT dsDNA BR Thermo Fisher Scientific Cat# Q33130 KAPA Library Quantification Kit Illumina KAPA Biosystems Cat# KK4873 DNA 1000 Kit Agilen Cat# 5067-1504 TruSeq Stranded Total RNA HT kit with Illumina Cat# RS-122-2303 Ribo-Zero Gold RNA 6000 Nano Kit Agilen Cat# 5067-1511 PhiX control Illumina Cat# FC-110-3001 SuperScript IV First-Strand Synthesis System ThermoFisher Scientific Cat#18091050 Deposited Data Aligned DNA and RNA sequencing data Deposited at European Genome Archive (EGA) with ID: EGAS00001002703 Additional supplemental figures Deposited in Mendeley Data repository. https://doi.org/10.17632/6cv77bry6m.1 (Continued on next page) Cell Reports 27, 2690–2708.e1–e10, May 28, 2019 e1 Continued REAGENT or RESOURCE SOURCE IDENTIFIER Software and Algorithms bcl2fastq2 2.17 Illumina https://support.illumina.com/sequencing/sequencing_software/ bcl2fastq-conversion-software.html R 3.2.2 R Core Team., 2017 http://www.r-project.org MATLAB version 9.2 1994-2017 The https://www.mathworks.com/products/matlab.html MathWorks, Inc. BWA Mem v0.7.15 Li and Durbin, 2009 http://bio-bwa.sourceforge.net/ GATK 3.4.46 DePristo et al., 2011 https://software.broadinstitute.org/gatk/ HaplotypeCaller HaplotypeCaller https://software.broadinstitute.org/gatk/documentation/ tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_ haplotypecaller_HaplotypeCaller.php Novosort 3.02 Novocraft http://www.novocraft.com/products/novoalign/ Novoalign 3.02 Novocraft http://www.novocraft.com/products/novoalign/ MuTect Cibulskis et al., 2013 https://software.broadinstitute.org/gatk/download Strelka 1.0.14 Saunders et al., 2012 strelka@ftp.illumina.com VEP (The Ensembl Variant Effect Predictor) McLaren et al., 2016 http://www.ensembl.org//useast.ensembl.org/info/docs/tools/ vep/index.html?redirectsrc=//www.ensembl.org%2Finfo%2 Fdocs%2Ftools%2Fvep%2Findex.html Integrative Genomics Viewer (IGV) Robinson et al., 2011 http://software.broadinstitute.org/software/igv/ Picard v2.2.1 Picard https://broadinstitute.github.io/picard/ samtools v1.3.1 Li et al., 2009 http://www.htslib.org/ ea-utils v1.1.2 Ea-utils https://github.com/ExpressionAnalysis/ea-utils Bioconductor 3.2 Huber et al., 2015 http://www.bioconductor.org Bioconductor package QDNaseq 1.2.4 Scheinin et al., 2014 http://www.bioconductor.org GISTIC2.0 Mermel et al., 2011 https://cloud.genepattern.org/gp/landingpage/index.html iC10: A Copy Number and Expression-Based Ali et al., 2014 https://rdrr.io/cran/iC10/ Classifier for Breast Tumors pam50: PAM50 classifier for identification of Parker et al., 2009 https://rdrr.io/bioc/genefu/man/pam50.html breast cancer R package deconstructSigs 1.8.0 Rosenthal et al., 2016 http://www.r-project.org/ ASCAT 2.5 Van Loo et al., 2010 https://www.crick.ac.uk/peter-van-loo/software/ASCAT PyClone 0.12.7 Roth et al., 2014 http://www.shahlab.ca EnsDb.Hsapiens.v75 Rainer, 2016 http://bioconductor.org/packages/release/data/annotation/ html/EnsDb.Hsapiens.v75.html POLYSOLVER Shukla et al., 2015 https://software.broadinstitute.org/cancer/cga/polysolver pVAC-Seq pipeline Hundal et al., 2016 https://github.com/griffithlab/pVAC-Seq Immunophenogram Charoentong et al., 2017 https://github.com/mui-icbi/Immunophenogram MEDICC (devel branch, commit da7ed4a) Schwarz et al., 2014 https://bitbucket.org/rfs/medicc superFreq 0.9.17 Flensburg et al., 2018 https://github.com/ChristofferFlensburg/cnv-caller treeomics 1.7.3 Reiter et al., 2017 https://github.com/johannesreiter/treeomics OncoNEM 1.0 Ross and Markowetz, http://bitbucket.org/edith_ross/onconem/src 2016 Tree: Raxml v8.2.1 Stamatakis, 2014 VarScan 2.4.3 Koboldt et al., 2012 http://dkoboldt.github.io/varscan/ alleleCount 3.1.1 alleleCount https://github.com/cancerit/alleleCount QUASR Watson et al., 2013 https://sourceforge.net/projects/quasr/ BLAST Altschul et al., 1990 ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ IMGT Lefranc, 2011 http://www.imgt.org/ Primer Design: mprimer (v1.9), primer3 Koressaar and Remm, N/A (v2.3.7), in silico PCR 2007; Shen et al., 2010; Untergasser et al., 2012 (Continued on next page) e2 Cell Reports 27, 2690–2708.e1–e10, May 28, 2019 Continued REAGENT or RESOURCE SOURCE IDENTIFIER ggplot2 2.2.1 ggplot2 https://ggplot2.tidyverse.org Igraph 1.0.1 Igraph https://cran.r-project.org/web/packages/igraph/index.html Ape 4.1 Paradis et al., 2004 https://cran.r-project.org/web/packages/ape/index.html Dendextend 1.5.2 Galili, 2015 https://cran.r-project.org/web/packages/dendextend/index.html nonnegative matrix factorization Gaujoux and Seoighe, https://cran.r-project.org/web/packages/NMF/index.html 2010 limSolve Soetaert et al., 2009 https://cran.r-project.org/web/packages/limSolve/index.html Other Adaptive 587 cohort data https://genomemedicine.biomedcentral.com/articles/10.1186/ s13073-015-0238-z Silhouettes http://silhouettesfree.com/download-silhouette/liver-silhouette/ CONTACT FOR REAGENT AND RESOURCE SHARING Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Carlos Caldas (Carlos.caldas@cruk.cam.ac.uk). EXPERIMENTAL MODEL AND SUBJECT DETAILS Patients Ten metastatic breast cancer patients who underwent post-mortem warm autopsies were included in this study. Nine patients were enrolled as part of the Vall d’Hebron Institute of Oncology (VHIO) Warm Autopsy Program, and one patient was enrolled at the Cambridge University Hospitals NHS Foundation Trust, Cambridge, UK as previously described (Murtaza et al., 2015). Informed consent was obtained for all patients. Research autopsies were performed under VHIO Warm Autopsy Program protocols approved by the institutional review board (IRB) of Vall d’Hebron University Hospital (Barcelona, Spain) and under a study protocol approved by the Cambridgeshire Research Ethics Committee (Cambridgeshire 3 REC 07/Q0106/63MN.A). METHOD DETAILS Sample Nomenclature Each sample identifier (ID) follows the nomenclature NNN-sample, where NNN denotes the patient ID. For whole exome sequencing (WES) and T cell receptor (TCR) sequencing, NNN-0XX indicates sample type (solid tumor or body-fluid derived DNA) (Table S1). For targeted amplicon sequencing and shallow whole genome sequencing (sWGS), the sample nomenclature is specified in Table S1. Briefly, this was derived from the primary tumor (PR), metastasis (year sample was taken (i.e., 13), followed by ID)) or body-fluid derived cell-free DNA (plasma, CSF, ascitic fluid, pleural fluid). Nucleic acid extraction The histologic evaluation of each diagnostic primary tumor or metastatic lesion from autopsy of the nine VHIO patients was confirmed on review of routine hematoxylin and eosin-stained slides. Samples were processed as previously described (De Mattos-Arruda et al., 2015). DNA and RNA were isolated from tumor tissue using commercially available kits according to manufacturer’s specifi- cations (Key Resources Table). DNA was extracted from peripheral blood mononuclear cells and body fluids using commercially available kits as per manufacturer’s specifications (Key Resources Table). RNA and DNA were quantified using the Qubit Fluorometer (Invitrogen). The collection, processing, DNA extraction and preparation of exome-sequencing libraries of DET52 patient for tumor tissues and plasma samples have been described previously (Murtaza et al., 2015). WES AND sWGS ANALYSES DNA sequencing Libraries for Illumina sequencing were prepared using Illumina Nextera Rapid Capture Exome kit (cat. FC-140-1003, Illumina). Prior to library preparation DNA concentrations for each sample were quantified using a fluorescence-based method (Quant-IT dsDNA BR, cat. Q33130, Thermo Fisher Scientific) and 50ng of genomic DNA was used for library preparation. Samples were processed following manufacturer’s instructions (part# 15037436 Rev. J, Illumina) for WES and sWGS. Prior to first hybridization all libraries were quantified using quantitative polymerase chain reaction (qPCR). KAPA Library Quantification Kit Cell Reports 27, 2690–2708.e1–e10, May 28, 2019 e3 (cat. KK4873, KAPA Biosystems) as used as per manufacturer’s recommendations. A subset of libraries was analyzed using DNA 1000 Kit (cat. 5067-1504, Agilent). Whole genome libraries and exome libraries were normalized and pooled in equal volumes to create balanced pools. Each pool was normalized to a molarity of 4nM and used for sequencing with clustering concentration 20pM with 1% spike-in of PhiX control. Sequencing was performed on an Illumina HiSeq2500 using v4 chemistry and 50 cycles single-end for s-WGS and 75 cycles paired- end for WES. Demultiplexing was performed using Illumina’s bcl2fastq2 v.2.17 software using default options. FASTQ files were used for sub- sequent data analysis. WES analyses Adaptor and low-quality base (Phred score below 20) trimming, alignment to the GRCh37 build of the human genome, and base qual- ity recalibration were performed using Novoalign v 3.02 (Novocraft). Coordinate sorting of reads and PCR-duplicate marking was performed using Novosort (v 3.02). The resulting bam files for all samples for the same patient were locally realigned using the Genome Analysis Toolkit (GATK, v 3.4.46) (DePristo et al., 2011). MuTect (version 1.1.7) was run using default parameters (Cibulskis et al., 2013). In order to decrease the false positive rate secondary to germline variants of noisy regions within the genome, a panel of normals derived from 300 normal tissue exomes was generated using MuTect’s artifact detection mode and supplied to MuTect during variant calling. Strelka (version 1.0.14) (Saunders et al., 2012) was run with recommended starting parameters for BWA and default parameters. The isSkipDepthFilters parameter was set to 1 to skip depth filtration. Only tier 1 mutations, as well as SNVs with a QSS_NT score higher than 15 and Indels with a QSI_NT score higher than 30 were retained. Mutations present in all samples for each patient were then concatenated into one VCF, and Haplotypecaller was used in GENOTYPE_GIVEN_ALLELES mode to detect these mutations across all samples. SNVs and indels that fell into any of these categories were removed: d Rejected by MuTect for a reason different than ‘‘DBSNP Site,’’ ‘‘DBSNP Site,alt_allele_in_normal’’ or ‘‘alt_allele_in_normal.’’ Mutations rejected by MuTect by any of these reasons that were present in the 1000 Genomes Project were also rejected d Read depth less than 10 d Variant allelic frequency less than 0.05 in all samples for one patient d Minimum allelic frequency of 0.02 per sample d SNVs falling in segmental duplications regions, as annotated by annovar (Wang et al., 2010) genomic superDups, were clas- sified as potential artifacts. Those that occurred in more than one patient, or in only one patient but in less than 25% of the samples from that patient were removed d SNVs falling in simple repeat regions, microsatellites or homoplymers were removed. Somatic mutations were annotated using Variant Effect Predictor (VEP, http://grch37.ensembl.org/index.html) (McLaren et al., 2016) and visualized using IGV (Robinson et al., 2011). sWGS analyses 50 bp single-read sWGS was performed to provide copy number profiles. FASTQ files were aligned to the GRCh37 assembly of the human genome using BWA (Li and Durbin, 2009) and the bam files were merged, sorted and indexed using samtools (Li et al., 2009). Duplicates were marked using Picard. Copy number profiles were obtained using the R package QDNaseq (Olshen et al., 2004), using non-overlapping 100 kb pairs windows, and correcting for mappability and GC content. Targeted amplicon sequencing Targeted amplicon sequencing (TS) was performed in 464 unique amplicons derived from WES across 10 patients (average 46.4/case, range 16-127). Targeted sequencing libraries were prepared using droplet-based PCR amplification following the man- ufacturer’s protocols for ThunderBolts Cancer Panel with specific modifications (RainDance Technologies) as previously reported (Murtaza et al., 2015). Multiplex primers sequences are shown in Table S6. Analysis of targeted sequencing data was performed as described previously (Murtaza et al., 2015). For each locus and non-reference allele of interest, we assessed the allele fraction in eight control genomic DNA samples. A mutation was considered significantly detectable if minimum coverage at that locus was 500x, and the AF in a sample was greater than or equal to 3 standard deviations higher than the mean AF in control samples, and if present in greater than or equal to 1% allelic fraction (AF). Additionally, any samples in which over 90% of all mutations had AFs below 5% were excluded. The mutation calls generated from targeted sequencing were then used to assess the quality of the exome mutation calling pipeline. Any mutations detected on both WES and TS were defined as true positives (TP), mutations detected on WES, but not TS defined as false positives (FP), mutations detected on TS but not on WES defined as false negatives (FN), and finally mutations that were not detected on TS and WES defined as true negatives (TN). Sensitivity was subsequently calculated as TP/(TP+FN), specificity calculated as TN/(TN+FP), precision calculated as TP/(TP+FP) and accuracy calculcated as (TP+TN)/(TP+FP+FN+TN). Genotyping Germline SNPs and indels were identified in all samples using GATK HaplotypeCaller, and all tumor samples were then genotyped to the matching normal tissues by computing the percentage of shared SNPs and indels between tumors and normal. Concordance of more than 90% was taken to indicate related samples. e4 Cell Reports 27, 2690–2708.e1–e10, May 28, 2019 RNA-sequencing analysis RNA sequencing RNA sequencing libraries were prepared using the TruSeq Stranded Total RNA HT kit with Ribo-Zero Gold (cat. RS-122-2303, Illumina). Prior to library preparation samples were quantified using a fluorescence based method and RNA quality was assessed using RNA 6000 Nano Kit (cat. 5067-1511, Agilent) on Bioanalyzer2100 (Agilent). Depending on availability, 400-900ng of total RNA was used for library preparation. The RNA Integrity Number (RIN) for these samples varied from 2.3 to 7.2. Samples were pro- cessed following manufacturer’s HS (High-Sample) instruction (part# 15031048 Rev. E, Illumina). Subset of 12 libraries was analyzed using DNA 1000 Kit (cat. 5067-1504, Agilent) and the average library length was determined as 280bp. All libraries were quantified using qPCR. Serial dilutions were made and 100,000x dilution was used for quantification using KAPA Library Quantification Kit Illu- mina (cat. KK4873, KAPA Biosystems). Libraries were normalized to 40nM and pooled in equal volumes to create a balanced pool. The library pool was quantified after doing serial dilutions in triplicate and 10,000x and 100,000x dilutions were used for quantifica- tion. The final library was normalized to 4nM and sequenced at a clustering concentration of 20pM and 22pM with 1% spike-in of PhiX control (cat. FC-110-3001, Illumina). Sequencing was performed on HiSeq2500 v4 chemistry single-end flow cell (Illumina) following manufacturer’s instructions. Demultiplexing was performed using bcl2fastq2 v.2.17 software (Illumina) using default options. RNA sequencing analyses FASTQ files were aligned to the GRCh37 assembly of the human genome using STAR v 2.5.2b in two-pass mode for splice-aware read alignment. Counting of reads aligned over exonic features for the purpose of differential expression was performed using the htseq-count script in the HTSeq package (v 0.6.1) in ‘Union’ overlap resolution mode using a Gene Transfer Format (GTF) file from Ensembl (http://www.ensembl.org//useast.ensembl.org/?redirectsrc=//www.ensembl.org%2F). The gene counts for all samples were then collated and FPKM calculations per gene per sample performed using the rpkm() function in the edgeR R package. De novo transcript assembly and counting of transcripts was performed using Cufflinks v2.2.1. For variant calling, STAR-aligned BAM files were processed as per the RNA-seq GATK Best Practices. Briefly, sequencing read duplicates were marked using Picard MarkDuplicates, followed by Split’N’Trim and mapping quality reassignment using GATK SplitNCigarReads (v3.6). This was then followed by local realignment across indels and base quality recalibration using GATK. Mutations detected in the corresponding DNA sequencing data were genotyped in RNA using the GENOTYPE_GIVEN_ALLELES mode in Haplotype caller. Selection of driver mutations Breast cancer driver mutations were defined as those genes identified in previous publications (Lefebvre et al., 2016; Nik-Zainal et al., 2016; Pereira et al., 2016) and non-breast cancer drivers were defined as those present at the Cancer Gene Census and non- overlapping with breast cancer driver mutations (https://cancer.sanger.ac.uk/census/) (Table S4). Mutational signatures Somatic substitutions of each metastatic sample were organized into a 96-channel vector (where the six mutation classes and their immediate flanking sequence context are taken into account), referred to hereafter as a mutational profile. Mutational signature analysis of these mutational profiles was performed in two steps: extraction and assignment. The first step in our analysis aimed to identify any signatures previously found in associated primary tumors that are present within our cohort. It consisted of applying the widely adopted (Alexandrov et al., 2013) Non-negative Matrix Factorization (NNMF) algorithm (R-CRAN package NMF - (Gaujoux and Seoighe, 2010)) to an extended dataset, where 240 additional WES metastatic breast cancer samples (Lefebvre et al., 2016) were added to our original 86 sample cohort. NNMF extraction was performed on these mutational profiles, bootstrapped 100 times, and a KL-divergence error was used to assess the accuracy of each result. The rank of the NNMF solution (i.e., the number of extracted signatures), was allowed to vary between 2-20. Across the different extractions, cosine similarity comparison with known canonical primary tumor signatures (COSMIC) revealed the presence of Signatures 1,2,13, and 17. The use of the additional samples increased the power of the NNMF, enabling a more precise mutational profile extraction for these four well-known breast-cancer-related signatures. The second step consisted of assigning the contribution of the four COSMIC signatures identified (Signature 1, 2, 13, 17) to each sample of our original 86 metastatic cohort. This was computed using a quadratic programming algorithm (R-CRAN package limSolve - (Soetaert et al., 2009)). A minimum of either 3% of the total number of mutations of the sample or at least 10 mutations was required for a COSMIC signature to be attributed to a sample. The assignment step was performed on four versions of the muta- tional profiles for each patient, including: (i) all mutations, (ii) mutations shared across all samples (stem), (iii) mutations shared across some samples (clade), and (iv) mutations uniquely present in one sample (private). IntCluster, PAM50 and stratification into breast cancers subtypes Matched samples (n = 60) with copy number and expression data were classified into one of the 10 Integrative Clusters using the ‘iC10’ R package (Ali et al., 2014; Curtis et al., 2012). The ‘iC10’ package uses copy number and expression from breast cancer data, trains a pamr classifier with the features available and predicts the iC10 group. Each sample was classified into the 10 Integra- tive Clusters and the assignment to each model was done by consensus after manual curation in specific cases. The 50-gene sub- type predictor PAM50 was also applyied to 64 metastases with expression data using the R package genefu (Parker et al., 2009). ASCAT ASCAT v2.5 (Van Loo et al., 2010) was run by integrating copy number log ratios generated from QDNaseq and SNP allelic frequencies from WES data. The gamma technology parameter was set to 1 as recommended for exome sequencing. Cell Reports 27, 2690–2708.e1–e10, May 28, 2019 e5 Phylogenetic analyses Four multi-sample methods were used to infer the metastatic breast cancer phylogenies. MEDICC (Minimum Event Distance for Intra-tumor Copy number Comparisons) is a method for phylogenetic reconstruction and heterogeneity quantification that uses allele-specific copy number profiles (Schwarz et al., 2014). OncoNEM (Ross and Markowetz, 2016) and Treeomics (Reiter et al., 2017) are tools that utilize Bayesian inference to infer phylogenetic relationships from mutation patterns of SNVs. SuperFreq is a clon- ality tracker that uses single nucleotide variants (SNV) and copy number alterations (CNA) (Flensburg et al., 2018; Savas et al., 2016). MEDICC QDNaseq (Scheinin et al., 2014) was applied to the sWGS data to obtain sequence mappability and GC content adjusted log ratios of read depth. B-allele frequencies were calculated from WES data at previously inferred germline variant sites using alleleCounte 3.1.1. Log ratios and B-allele frequencies were segmented on a per case basis using allele specific multi-sample segmentation (Ross et al., 2017) and allele-specific copy number profiles were inferred using ASCAT 2.5 (Van Loo et al., 2010). The raw ASCAT copy number profiles were compared across samples for each case and ASCAT was rerun with adjusted ploidy and/or purity estimates where necessary to obtain the final discrete copy number profiles. Samples with copy number fits of low quality were excluded. A maximum copy number cut-off of nine was applied to both major and minor copy number profiles, replacing any values exceeding this threshold to comply with MEDICC requirements. Finally, MEDICC was used to infer the phylogenies (Schwarz et al., 2014). Treeomics Treeomics 1.7.3 (Reiter et al., 2017) was used to infer phylogenies from both WES and targeted sequencing data. For each case, Treeomics was used to calculate the posterior probabilities of a variant being present based on total read depth and number of reads covering the alternative allele. To make the Treeomics analysis computationally feasible, the number of samples had to be reduced for some of the cases. In order to keep the mutation profiles as diverse as possible and to maintain a good representation of the different tumor populations, samples that had a mutation profile similar to one of the remaining samples were excluded preferentially. Additionally, all sites whose posterior probabilities were lower than 0.5 in all samples were removed, as these were likely to be false positives. Finally, Treeomics was applied to each case with subclone detection switched on and all other parameters set to default. OncoNEM Like Treeomics, OncoNEM 1.0 (Ross and Markowetz, 2016) was used to infer phylogenies from both WES and targeted sequencing data. Binary mutation profiles were obtained from the Treeomics posterior probability matrices by setting all entries with a mutation probability smaller than 0.5 to 0 and to 1 otherwise. The OncoNEM analysis was performed using error parameter optimization over a parameter range from 0.0001 to 0.1. The Bayes factor threshold epsilon was set to 2. SuperFreq SuperFreq was used to infer the clonal composition of the different samples from WES data. Pileup files of the WES data were gener- ated using samtools 1.3.1 mpileup using a maximum depth threshold of 10000, a minimum mapping quality of 1 and a minimum base quality of 15. Liberal variant calling was performed using VarScan 2.4.3 mpileu2cns with a p value filter of 0.01, no strand-bias filter and the variant flag set to only obtain variant sites. Then SuperFreq 0.9.17 was run with default parameters using the normal of all cases apart from DET52 as reference normal samples. PyClone PyClone is a Bayesian clustering method that infers the clonal population structures for each sample (Murtaza et al., 2015; Roth et al., 2014). It integrates mutation alleles, copy number calls for each sample as input to obtain cellular frequencies for each cluster in each sample. PyClone was run using a beta- binomial density, using 40000 iterations and a burn-in sample of 20000. A minimum cluster size of 3 was selected for WES data. LICHeE The previously inferred PyClone clusters were used as input for LICHeE. To remove spurious clones from the PyClone output two filtering steps were performed. First, low prevalence clones that did not exceed a cellular prevalence of 0.1 in any of the samples were removed. Second, if multiple clusters were present in all samples, all but the cluster with the highest cellular prevalence were removed. To generate binary presence-absence patterns of mutation clusters across samples, a mutation set was classified as present in a given sample if at least 40% of its mutations had a VAF larger than 0.01. Finally, LICHeE (commit 238770c) was used to infer clone trees and sample compositions, assuming a cellular prevalence estimate error of 0.3. Cytolytic activity score Cytolitic activity score was calculated as the geometric mean of the GZMA and PRF1 expression levels from RNA gene expression data (Rooney et al., 2015). Immunophenogram and immunophenoscore calculations The immunophenogram (Charoentong et al., 2017) was applied to determine the immunophenotypes of each tumor sample and to enable the calculation of an aggregated score (immunophenoscore (IPS)), based on the expression of major determinants. These factors were classified into four categories: MHC molecules (MHC), Immunomodulators (CP), Effector cells (EC), Suppressor cells (SC), and into 20 single factors (MHC molecules, immunoinhibitors, and immunostimulators) and six cell types (effector cells: acti- vated CD4+ T cells, activated CD8+ T cells, effector memory CD4+ T cells and effector memory CD8+ T cells; suppressive cells: Tregs, and MDSCs). Briefly, the algorithm generates normalized Z scores from gene expression data for a list of cancer immunity parameters (using an input list of 162 genes (Table S5)). The outer part of the wheel illustrates sample-wise (averaged) Z scores, which is calculated for each e6 Cell Reports 27, 2690–2708.e1–e10, May 28, 2019 of the individual 26 parameters. These Z scores are positively weighted according to stimulatory factors (cell types) and negatively weighted according to inhibitory factors (cell types) and averaged. The inner wheel illustrates the weighted Z scores of the factors included into the four categories. Z scores R 3 were designated as IPS10 and z-scores % 0 are designated as IPS0. Unsupervised hierarchical clustering based on the Euclidean distance matrix of the Z scores across metastases (R heatmap.3 function) were used to produce the heatmaps. HLA typing and in silico neoantigen prediction pipeline For the metastases of all patients, the 4-digit HLA type was determined using POLYSOLVER (POLYmorphic loci reSOLVER) as previously described (Shukla et al., 2015). The pVAC-Seq pipeline (Hundal et al., 2016) was used with minor modifications. All nonsynonymous point mutations identified were translated into strings of 17-21 amino acids with the mutant amino acid situated centrally. A sliding window method was used to identify amino acid substrings within the mutant 17-mer that had a predicted MHC Class I binding affinity of % 500nM to one (or more) of the patient-specific HLA alleles. Binding affinity for the mutant and corresponding wild-type nonamer were analyzed using the NetMHCPan v3 prediction tool bundled within the IEDB MHC-I binding prediction resources. Following this, for all cases excluding DET52 (which did not have any RNA expression data), candidate neoantigens were further filtered by retaining mutations that were also present in the RNA sequencing data, as well as had a gene and transcript RPKM expression of more than 1. Neoanti- gens were subsequently classified as stem, private and clade by applying the classification derived from the WES mutational dataset. When generating Figure 5A, highly similar neoantigens generated from one mutation were counted as one. Neo-antigen simulations Simulations procedure We assumed that each gene has its own mutation rate. For each gene, we used the mutations list of 4,742 WES tumor normal pairs from (Lawrence et al., 2014) (http://www.tumorportal.org) to estimate the gene relative background mutation rate by counting the amount of mutations the gene had divided by the total number of mutations. Every mutation was then randomly assigned to a new gene based on the gene’s relative rate. The position within the gene was chosen to maintain the trinucleotide context of each mutation (the 50 and 30 nearest neighbors and the mutated position) and the variant was based on the original mutation. In addi- tion, for every base we counted how many times it was sufficiently covered for mutation detection (i.e., > = 14 reads), across 7,732 TCGA tumor WES samples. The fraction of covered patients at a given base was used as a probabilistic weight when selecting the new position for a mutation. P value calculation For a single sample analysis we defined the P-value as the fraction of replications (out 100) that had a neo-epitopes ratio larger than 1. For combining cases based on their target tissue we tried to see if their average was less expected at random. For that we generated 20 replicates and and 100 additional simulations based on each of the 20 replicates. We then calculated the average neo-epitopes ratio for each of these sets and compared it to the original ratio. We then calculated a P value for the fraction of the real average being different from the simulations average. Loss of heterozygosity in human leukocyte antigen (LOHHLA) LOH over the HLA Class I locus was identified using LOHHLA on the whole exome sequencing data, and LOH called if the copy number at HLA-A/HLA-B or HLA-C locus was less than 0.5, with a p value of less than 0.05 (McGranahan et al., 2017). T cell receptor (TCR) analysis TCR-sequencing library preparation Reverse transcription primers were designed using Primer3 (Koressaar and Remm, 2007; Untergasser et al., 2012) and multiplex PCR primers using MPrimer (Shen et al., 2010). cDNA synthesis was performed using TCR-a or TCR-b constant region-specific primers carrying a molecular barcode of 12 random/degenerate nucleotides (N12, TNNNNTNNNNTNNNNT) to enable molecule-level identification (unique molecular identifier). The molecular barcode was inserted upstream (50 ) of the sequence that recognizes the constant region(s) and downstream (30 ) of an adaptor sequence complementary to first round PCR reverse primers. Reverse transcription was performed with SuperScript IV First-Strand Synthesis System (Invitrogen, ThermoFisher Scientific) using 500ng to 2 mg of total RNA, 1 mL of barcoded TCR-a or TCR-b specific reverse primer (0.1 mM final), 1 mL of dNTPs (0.5mM final) and added with nuclease-free water to a total volume of 13 ml. This was incubated at 65 C for 5min, then on ice for 2 minutes and followed by addition of 4 mL of 5X First strand buffer, 1 mL of DTT (5mM final), 1 mL of RNaseOUT Recombinant Ribonuclease Inhibitor (40 units), and 1 mL of SuperScript IV reverse transcriptase (200 units). The reaction was incubated at 56 C for 60min followed by inactivation at 80 C for 10min. 1 mL of E. coli RNase H (2 units) was added and incubated at 37 C for 20min to remove RNA from the cDNA:RNA hybrids. The first strand cDNA was cleaned up using 1.8x Agencourt Ampure XP beads (Beckman Coulter) and eluted with 21 mL of nuclease-free water (Ambion). First-round multiplex PCR amplifications were set up in a total volume of 50 ml, with 20 mL of cDNA as template, 25 mL of Q5 Hotstart High-Fidelity DNA Polymerase Master Mix 2x (New England Biolabs) and tailed TCR-a or TCR-b forward primer set pools and sam- ple-indexed reverse primers (0.2uM final concentration each). Multiplex forward primers target different TCR-a or TCR-b V-regions and sequences are shown in Table S7. The sample-indexed reverse primers used were published previously (Bronner et al., 2014). The following PCR program was used: 30 s at 98 C, 25 cycles of 20 s at 98 C, 1 min at 55 C, and 1 min at 72 C, with a final extension Cell Reports 27, 2690–2708.e1–e10, May 28, 2019 e7 step of 2 min at 72 C. The PCR product was cleaned up using 0.75x Agencourt Ampure XP beads (Beckman Coulter) and eluted with 20 mL of nuclease-free water (Ambion). First-round forward (TCR-a or TCR-b V-region-specific) PCR primers each contained a shared sequence to allow Illumina sequencing adapters to be introduced with a second round PCR. The second-round PCR amplification step was performed on the first round PCR amplicons to generate Illumina-ready sequencing libraries. 12 cycles of PCRs were performed and the product was analyzed and quantitated using Agilent Bioanalyzer DNA 1000 chips. For each batch, equal nanomoles of each sample were pooled, double SPRI size selected (0.5x and 0.7x) and stored at 20 C until sequencing. Libraries were batched and sequenced on MiSeq sequencers (300bp paired-end reads). TCR-sequencing and barcode filtering MiSeq libraries were prepared using Illumina protocols and sequenced using 300bp paired-ended MiSeq (Illumina). Raw MiSeq reads were filtered for base quality (median Phred score > 32) using the QUASR program (https://sourceforge.net/projects/quasr/) (Watson et al., 2013). MiSeq forward and reverse reads were merged together if they contained identical overlapping region of > 50bp, or other- wise discarded. Universal barcoded regions were identified in reads and orientated to read from V-primer to constant region primer. The barcoded region within each primer was identified and checked for conserved bases (i.e., the T’s in NNNNTNNNNTNNNNT). Primers and constant regions were trimmed from each sequence, and sequences were retained only if there was > 80% sequence certainty between all sequences obtained with the same barcode, otherwise discarded. The constant region allele with highest sequence similarity was identified by k-mer matching (where k = 10bp) to the reference constant region gene IMGT database (Lefranc, 2011), and sequence trimmed to give only the region of the sequence corresponding to the variable (V-D-J) regions, where constant region usage information for each TCR was retained throughout the analysis. Sequences without complete reading frames and non-TCR sequences were removed and only reads with significant similarity to reference TRBV or TRAV and J genes from the IMGT database were retained using BLAST (Altschul et al., 1990). Sequences were annotated using IMGT. Sample clustering was performed as previously described (Bashford-Rogers et al., 2013). TCR repertoire generation and network analysis The network generation algorithm and network properties were calculated as in Bashford-Rogers et al.(Bashford-Rogers et al., 2013): each vertex represents a unique sequence, where relative vertex size is proportional to the number of identical reads. A clone (cluster) refers to a group of identical related T cells, each containing TCRs with identical CDR3 regions and TCRV gene usage. Repertoire parameters that were dependent on sequencing depth were generated by subsampling each sequencing sample to a specified depth and the mean of 20 repeats of resulting parameters were calculated using the clonality measures. These measures include 1) total repertoire clonality (vertex & cluster Gini Indices), and 2) mean cluster sizes, (3) largest cluster sizes calculated as follows: 1) Total repertoire clonality, measured by vertex & cluster Gini Indices are defined in Bashford-Rogers et al. (Bashford-Rogers et al., 2013), calculated from the distribution of the number of unique RNA molecules per vertex and the distribution of the number of unique TCRs per cluster, respectively. These were all calculated per 1000 read subsample of the each total TCR repertoire. 2) Mean cluster sizes (MCS) were within each subsample for the total TCR repertoire were calculated as follows for any given sample Y: P N TCRs per cluster in subsample MCSðSample YÞ = P x100 Total number of clusters in subsample 3) Largest cluster sizes (LCS) were within each subsample for the total TCR repertoire were calculated as follows for any given sample Y: MaxðN TCRs per cluster in subsampleÞ MCSðSample YÞ = P x100 Total number of clusters in subsample Clonal overlap analysis Clonal groups defined as TCRs sharing same V and J gene usages and identical CDR3 region sequence (nucleotide). Public TCRs were defined as TCR clusters that were shared between 2 or more T cell samples from unrelated individuals within this dataset and compared to the Adaptive Biotechnology database (Dean et al., 2015). The clonal overlap coefficient, describing the overlap between samples, was calculated either considering the whole TCR datasets, or with the public TCR clones removed. Clonal overlap, Oði; jÞ; between any two samples, i and j was calculated by Ci;j Oði; jÞ = 1 2 ðCi+ Cj Þ = Where Cx is the number of clusters in sample x and Ci;j is the number of clusters shared between samples i and j. To account for differences in sequencing depth between samples, each sample was subsampled each sequencing sample to a specified depth (1000 TCRs), and the mean of 20 repeats per sample of resulting clonal overlap was calculated. e8 Cell Reports 27, 2690–2708.e1–e10, May 28, 2019 Similarity heatmaps Similarity heatmaps were produced using Jaccard index calculated between each pair of metastases, using unique amino acids TCRs (CDR3s), for the alpha and beta chains. The Jaccard index was calculated using the ecological toolkit of the vegan R package, and the heatmaps were produced using R pheatmap package. Clustering of the heatmaps was done by the standard R hclust (hierarchical clustering) method, using the ‘‘complete’’ option. The comparison of the hclusts objects was done by the cophenetic correlation, using the dendextend package (Galili, 2015). The data were reshuffled to assign the TCR sequences to ‘‘randomized’’ metastases, and then on the reshuffled repertoire clus- tering by the Jaccard index was performed. This randomization was done 100 times. In this setting where there are 100 clusterings performed on the different randomizations and these were compared then between themselves, and to the original clustering of the biological data. To compare the randomized and ‘real’ trees we used the cophenetic correlation, and the Robinson-Foulds metric. The comparison of the Jaccard clustering trees with the genetic trees was done by using the cophenetic definition for edge- weighted trees. In this version of the cophenetic, the distance between each pair of nodes is the sum of the edges weights along the path connecting these pair of nodes. Correlation of TPM of CDR3 genes and sum of TCR reads From the RNASeq we extracted the TPM values of four genes that encode for the four different parts of the CD3 complex. These genes are CD3D, CD3G, CD3E and the zeta chain CD247 (Ensembl codes: ENSG00000167286, ENSG00000160654, ENSG00000198851, ENSG00000198821, respectively). For each sample, we calculated two measurements: 1) the sum of the RNASeq TPM values for these four genes; 2) from the TCR repertoires - the sum of the number of alpha chains and beta chains for each sample. We computed the Pearson correlation between the log10 values of these two measurements across all samples. The correlation and p value were computed using R’s cor.test method. Histopathological analyses Tissue microarrays (TMA) were prepared using duplicate 1 mm cores extracted from formalin-fixed paraffin-embedded blocks con- taining material from the individual tumors and metastases. Immunohistochemistry (IHC) was conducted for CD68, CD3, CD19, FOXP3, CD8, IL3RA, IDO1, CD4, CD56, CD1A, Mast Cell Tryptase, CD45RO, CD38, PDL1, ER, PR, and HER2 proteins. Details of reagents and protocols for IHC are provided in Key Resources Table. Stains were manually quantified by counting the absolute number or the percentage of positive stained cells. ER, PR statuses were assessed based on IHC applied to TMAs and scored using the percentage of positive tumor cells and intensity of staining (AllRead score). Herceptest was performed for all samples according to ASCO/CAP Guidelines (Wolff et al., 2013). Fluorescence in situ hybridization (FISH) (HER2-to-CEP17 ratio and gene copy number) for HER2 status was performed as previously described for all samples (Wolff et al., 2013). Positive HER2 amplification was considered when FISH ratio was higher than 2.2 or HER2 gene copy greater than 6.0. Digital pathology Whole slide images (either FFPE sections or frozen sections from tissue samples used for RNA extraction) were analyzed using CellExtractor v1.0, an open-source platform developed for high throughput analysis of histopathological images. The code was writ- ten in Python and used the OpenCV, i.e., an open source computer vision and machine learning software library written in C++, and the OpenSlide library. Full-face H&E scanned images were analyzed and divided into several sub-regions. Each sub-region is pro- cessed and segmented to compute cellular features such as centroids, shape descriptors, and pixel intensities. These features were used to train a support-vector machine (SVM) based classifier to identify cancer cells, stromal cells, and lymphocytes based on a training set of objects selected by a pathologist (W.C.) of approximately 1,000 objects for each category. Finally, based on these clas- ses descriptive statistical parameters such as cellular fractions and densities were estimated. For each detected cell density was obtained based on counting the number of nearest neighbor approach, i.e., the density within the distance to the Nth nearest neighbor calculated as follow: Sigma_N (pixel^[-2]) = N/(pi* d_N^2), Where d_N is the distance to the Nth nearest neighbor within a density-defining population. A value of N = 50 was used in order to estimate the density parameter (see (Ali et al., 2016b) for a detailed description). QUANTIFICATION AND STATISTICAL ANALYSIS All statistical analyses were performed using R version 3.2.2 and associated packages (Key Resources Table). The statistical details of experiments including the exact value of n in terms of number of samples for a given patient, the experimental method and specific statistical tests employed are reported in the relevant section, Results, Figures and Figure Legends, and Supplementary tables. For a given test (i.e., Wilcoxon rank sum, test chi-square test) significance was defined if a p value was less than 0.05. DATA AND SOFTWARE AVAILABILITY Software Custom scripts to run the analyses described in the manuscript are available at https://github.com/cclab-brca/MET-breast-landscapes/ Cell Reports 27, 2690–2708.e1–e10, May 28, 2019 e9 Data Resources Sequence data has been deposited at the European Genome-phenome Archive (EGA), which is hosted by the EBI and the CRG, under accession number EGAS00001002703. Further information about EGA can be found on https://ega-archive.org. The European Genome-phenome Archive of human data consented for biomedical research (https://idp.nature.com/authorize? response_type=cookie&client_id=grover&redirect_uri=http%3A%2F%2Fwww.nature.com%2Fng%2Fjournal%2Fv47%2Fn7%2F full%2Fng.3312.html). Supplemental Information was deposited on Mendeley at https://doi.org/10.17632/6cv77bry6m.1 e10 Cell Reports 27, 2690–2708.e1–e10, May 28, 2019