Therapeutics for breast cancer based on molecular classification

Abstract

Systems and methods for assessment of breast cancer are provided. The assessment can stratify breast cancer individuals based on genomic archetypes of the cancer. The assessment can be utilized as a diagnostic to determine a treatment regimen that would benefit a patient based on the genomic archetypes. The treatment regimen can include a targeted treatment that targets the molecular characteristics that are associated with the genomic archetype. A targeted treatment can be optionally combined with chemotherapy, immunotherapy, hormonal therapy, or other therapies for treating breast cancer.

Description

THERAPEUTICS FOR BREAST CANCER BASED ON MOLECULAR CLASSIFICATIONCROSS REFERENCE TO RELATED APPLICATIONS

[0001] This application claims the benefit of U.S. Provisional Patent Appl. No. 63 / 631 ,397, entitled “Therapeutics for Breast Cancer Based on Molecular Classification,” filed April 8, 2024, the disclosure of which is hereby incorporated by reference in its entirety.STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT

[0002] This invention was made with Government support under contracts NIH Pioneer Award and U54 MetNet awarded by the National Institutes of Health. The Government has certain rights in the invention.TECHNICAL FIELD

[0003] The disclosure is generally directed to systems and methods for assessing breast cancer.BACKGROUND

[0004] Breast cancer is the most common malignancy in women, accounting for more than 15% of new cancer cases in the USA annually. Clinically, breast tumours are stratified into three immunohistochemistry subtypes— ER+HER2", HER2 and triplenegative breast cancer (TNBC)— on the basis of the expression of ER, progesterone receptor and HER2. Although heterogeneity in gene expression, especially measures of proliferation, within these subtypes correlates with prognosis and patterns of relapse, and is used to guide therapy, ultimately the paradigm of three major subtypes dictates our understanding of and approach to the disease.SUMMARY

[0005] Various embodiments are directed to assessment of breast cancer. In many embodiments, a breast cancer of a patient is assessed to determine a genomic archetype.In several embodiments, a genomic archetype is map is utilized to determine where the genomic archetype of the breast cancer. In many embodiments, to project the genomic archetype on the map, a sequencing result of the breast cancer is assessed for copy number aberrations (and optionally structural variants) via a computational assessment using matrix factorization or a predictive computational model. In several embodiments, the map projection is a genomic archetype enrichment score (or vector) that is utilized to determine a treatment regimen and / or further assessments to be performed. In various embodiments, one or more of the following additional assessments is performed: assessment of integrative cluster, assessment of replication stress, assessment of homologous recombination deficiency, and assessment of genetic mechanisms of immune escape. In many embodiments, the treatment regimen is also based on the additional assessment.BRIEF DESCRIPTION OF THE DRAWINGS

[0006] The description and claims will be more fully understood with reference to the following figures and data graphs, which are presented as exemplary embodiments of the invention and should not be construed as a complete recitation of the scope of the invention.

[0007] Fig. 1 provides a flow chart of a clinical method for assessing a breast cancer.

[0008] Fig. 2A provides a schematic that provides an example of determining a genomic archetype of a breast cancer.

[0009] Fig. 2B provides a visualization of a genomic archetype enrichment score that can be utilized to determine a treatment regimen.

[0010] Fig. 2C provides treatment regimens based on archetype enrichment score.

[0011] Fig. 2D provides a flow chart of possible assessments that can be performed in conjunction with an archetype enrichment score.

[0012] Fig. 2E provides a schematic of a methodology to classify a breast cancer into an integrative cluster subtype.

[0013] Figs. 3A-3F provide data and schematics of computational method ENICIust and identification of the IC subtypes. Fig. 3A, Schematic of the study design. BC, breast cancer; HTAN, Human Tumor Atlas Network; sWGS, shallow WGS; TCGA, The CancerGenome Atlas; PCAWG, Pan-Cancer Analysis of Whole Genomes. Fig. 3B, Schematic of the ENiClust IC classifier. WES, whole-exome sequencing. Fig. 3C, Kaplan-Meier curves of distant relapse-free (DRF) survival of the ER+typical-risk and ER+high-risk classes detected by the four IC subtype classifiers. Shaded area represents 95% confidence interval. HR, hazard ratio. Fig. 3D, Difference in distant relapse-free survival probability (top) or delta in Cox proportional hazard ratio (bottom) between ER+typicalrisk and ER+high-risk classes detected by the four different IC classifiers. Error bars represent the difference in 95% confidence intervals between ER+typical-risk and ER+high-risk in each model. Fig. 3E, Differential pattern of relapse across the ICs, illustrated by the cumulative (black) and annual (red) risk of relapse over time. Fig. 3F, IC subgroup (left) and subtype (right) distributions across disease stages. P, primary; M, metastatic.

[0014] Figs. 4A-4I provide data and schematics of quality assessment and Ensemble Integrative Clustering (ENiClust) predictive model performance. Scatterplots showing the consistency between (Fig. 4A) single-nucleotide variant (SNV), (Fig. 4B) structural variant (SV), and (Fig. 4C) fraction of genomic alteration (FGA) burden reported in the present study and in the PCAWG study 11 . Fig. 4D, Heatmaps indicate to the performance metrics (P: Precision, R: Recall, F1 : F1 -score, and PR AUC: PrecisionRecall Area Under Curve) of the DNA-based predictive models iC10 DNA-only and ENiClust compared to iC10 DNA+RNA or ENiClust compared to iC10 DNA-only on the Testing dataset (n=187 TCGA-WES, left), the with-held Validation datasets: ICGC (n=189 WGS, middle) and METABRIC (n=986 arrays, right). Crossed tiles represented cases when less than 4 samples were counted for the given class in either the prediction or the reference model. Fig. 4E, Kaplan-Meier curves of distant relapse free (DRF) survival of the ER+ Typicalrisk reclassified as ER+ High-risk (top) and the ER+ High-risk reclassified as ER+ Typicalrisk (bottom) by ENiClust (purple) vs. iC10 DNA-only (blue) (left) or (right) by ENiClust (purple) compared with the reference (iC10 DNA+RNA, red). Fig. 4F, Forest plot of the hazard ratios of the ER+ Typical vs. ER+ High-risk classification from the four IC subgroups classifiers (ENiClust, iC10 DNA+RNA, iC10 DNA-only, and iC10 RNA-only) for DRF survival. Squares correspond to estimated hazard ratios and the segments correspond to their 95% confidence intervals. Fig. 4G, Barplots illustrate the proportionand count of IC subtype (left) and subgroup (right) tumors in sensitive, intermediate and resistant cases from the clinical trial (NCT00651976). In this trial, tumors were classified based on their post-treatment Ki67 score as sensitive (Ki67 < 2.7%), intermediate (Ki67 2.8-7.3%), or resistant (Ki67 > 7.3%) to endocrine therapy. Fig. 4H, Boxplots illustrate the PAM50 proliferation score (Martin et al.48, encompassing 10 genes) from mRNA in ER+ High-risk (left) and ER+ Typical-risk cases (right) at pre-treatment (Pre-tx) and early on- treatment (Early-tx) time points. Each dot represents one sample and each line connects samples from the same patient from Pretx to Early-tx. The displayed values on top of the boxplot are the P-values from Wilcoxon tests comparing the two IC subgroups’ proliferation mRNA levels at the two time points. Fig. 4I, Boxplot shows purity (y-axis) for each subgroup (x-axis). HR: hazard ratio. Survival analysis was performed using Cox’s proportional hazard models corrected for clinical covariates.

[0015] Figs. 5A-5J provide data and schematics of complex rearrangements fuel ER+and HER2+breast tumors. Fig. 5A, IC subgroup (left) and subtype (right) across stages of progression in ER+ samples. Fig. 5B, Inferred ancestry (primary samples, left or metastatic samples, right) across IC subgroups. Fig. 5C, IC subgroup (left) and subtype (right) across inferred ancestry in primary (top) and metastatic (bottom) stages. Fig. 5D, IC subgroup (left) and subtype (right) across inferred ancestry in primary (top) and metastatic (bottom) stages in ER+ samples. Fig. 5E, IC subtypes (left) and subgroups (right) across paired primary and metastatic samples with WES data. Fig. 5F, Graphical network representing primary / primary or primary / metastatic pairs; dots corresponding to a tumor biopsy colored by IC subgroup, the edge between two dots indicates whether the classification is stable (black) or changes (red) through metastasis. The surrounding color represents the PAM50 subtype (gray indicates missing data). Fig. 5G IC subtype or subgroup consistency with PAM50 in pre-invasive DCIS (left), primary (middle), and metastatic (right) samples. Fig. 5H, ER early transcriptional signature according to subgroup. Fig. 5I, ER early signaling transcriptional signature in primary and metastatic tumors. Fig. 5J, Schematic overview of IC-specific amplification peaks and associated genes. ES, effect size. DCIS, ductal carcinoma in situ; LumA, luminal A; LumB, luminal B.

[0016] Figs. 6A to 6C provide data and schematics for SVs define three distinct genomic archetypes. Fig. 6A, IC group-level CNA profile (shaded area; dark denotes amplification, light denotes deletion) with SV burden (line) as overlay and total alteration burden in primary and metastatic samples. Fig. 6B, Pareto front projection on ternary plot of CNA and SV signature profiles from primary (left) and metastatic (right) tumors independently, resulting in three genomic archetypes. Each plotted circle represents a tumor. Fig. 6C, Lollipop plots illustrating the correlation between mutational features and the distance to each archetype, amp., amplification; BFB, breakage-fusion-bridge; TIC, templated insertion chain; LOH, loss of heterozygosity; WGD, whole-genome doubling; FGA, fraction of genome altered.

[0017] Figs. 7A to 7J provide data and schematics for genomic features of primary IC subgroups. Fig. 7A, IC group-level copy number profile with SV burden overlay in DCIS. Fig. 7B, Fraction of genome altered by subgroup across pre-invasive, primary invasive and metastatic tumors. Boxplot represents median, 0.25 and 0.75 quantiles with whiskers at 1.5x interquartile range. Fig. 7C, Alteration burden in metastatic tumors split based on treatment prior to biopsy. The sample size is indicated at the top of each bar. Fig. 7D, Fraction genome altered, fraction LOH, and number of damaging SVs in IC10 and IC4ER- subtypes. Fig. 7E, Proportion of IC10 and IC4ER- tumors with alterations in genes involved in three key pathways: cell cycle, DNA damage response (DDR), and ubiquitination. Fig. 7F, Alteration burden distribution in metastatic samples across metastatic sites. The sample size for each group is at the top of each bar. Figs. 7G-7H, Activity of each of the six rearrangement signatures across the IC subgroups (Fig. 7G) or the ER+High-risk subtypes (Fig. 7H) in primary tumors. Fig. 7I, Copy number and SV profiles of primary (left) and metastatic (right) samples, each representative of either a TNBC -enriched, ER+Typical -enriched, ER+High / HER2+-enriched, or mixed profile in the center of the Pareto front. Fig. 7 J , Proportion that each complex SV event contributes to the total complex SV burden stratified by subgroup. DEL, deletion; LOH, loss-of- heterozygosity; FDR, false discovery rate; BFB, bridge-fusion breakage; CPXDM, complex double minute; DM, double minute; INVDUP, inverted-duplication; TIC, templated insertion chain; TRA, translocation.

[0018] Figs. 8A-8G provide data and schematics for mutational processes underpinning primary breast cancer. Fig. 8A, Representation of the six SV signatures discovered. SVs were categorized based on size, type of rearrangement and proximity to neighboring SVs. Barplots show the proportion of each category attributed to the six SV signatures. Fig. 8B, Scatterplot shows stability (silhouette index) and similarity (cosine distance) scores, y-axis, used to select the optimal number of SV signatures, x-axis, in the primary tumors. Gray shading indicates the number selected. Fig. 8C, Heatmap shows the Spearman’s correlation between de novo SV signature composition and Nik- Zainal et al. signature composition 8. Fig. 8D, Heatmap shows the Spearman’s correlation between the SV signatures and 24 COSMIC CNA signatures. Covariate on the right indicates the category of the CNA signature as previously defined by Steele et al.21 . Fig. 8E, Violinplots show activity of the nine most recurrent CN signatures across the subgroups in primary. Fig. 8F Scatterplot showing absolute correlation of SV types burden with RS4 and RS6 activity across all subgroups (left) and per subgroup (right). Fig. 8G Boxplot showing RS4 (top) and RS6 (bottom) signature activity in tumors with or without alterations in specific SV types, with the number of tumors indicated, p-values from linear regression model for SV type burden correcting for RS4 and RS6 activity are shown. Del, deletion; Tds, tandem duplications; Inv, inversions; T, translocations; HRD, homologous recombination deficiency; LOH, loss-of-heterozygosity; cLOH, chromosomal-scale LOH; fLOH, focal LOH; SV, structural variant; CNA, copy number aberrations; WGD, whole genome doubling.

[0019] Figs. 9A-9J provide data and schematics of validation of Pareto Front archetypes. Figs. 9A to 9B, The proportion (Fig. 9A) or raw activity (Fig. 9B) of the six SV signatures and 24 CNA signatures active in each primary tumor were projected onto a two-dimensional space using PCA. Scatterplot represents each tumor as a dot along the first and second principal components and colored by its IC subgroup. Fig. 9C, Arrow color, direction and length indicate the contribution of each RS and CN signature to principal components 1 and 2. Fig. 9D, Barplot indicates the magnitude of the contribution of the top 10 RS and CN signatures to the principal components 1 and 2. Fig. 9E, Pareto front model fits robustness statistics (explained variance and t-ratios) for each number of PCs / number of Archetypes combination for each architectural PCAs. Fig. 9F, Pareto frontmodel fits robustness statistics (explained variance and t-ratios) across the architectural PCAs with k=3 archetypes within PC1 -PC2 spaces. Fig. 9G, Lollipop plots illustrating the correlation between various mutational features and the distance to each archetype. Fig. 9H, Proportion of TNBC-enriched Archetype in BRCA1-like, BRCA2-like, and non-HRD- like tumors within ER+and TNBC tumors. Fig. 9I, Signature contribution of SBS signatures (left) and activity of SV signatures (right) of BRCA1 -like, BRCA2-like, and non- HRD-like ER+ and TNBC tumors. P-values are from Mann-Whitney Rank Sum test. Fig. 9J, Barplot shows the proportion (top) or number (bottom) of samples predicted to be BRCA1-like (dark) or BRCA2-like (light) in ER+High risk IC subtypes. In h and i, p-values are from Mann-Whitney Rank Sum testecDNA, extrachromosomal DNA (called by JaBba); inv, inversion; invdup, inverted-duplication.

[0020] Figs. 10A-10O provide data and schematics for genomic features are conserved though elevated through metastasis. Fig. 10A, Pareto front projection with tumors colored by presence of co-amplification of two or more amplifications in the following cytobands: 17q23 (IC1 ), 11q13 (IC2), 17q12 (IC5 / HER2+), 8p12 (IC6) or 8q24 (IC9). Fig. 10B, Pareto front projection with tumors colored by HRD and ER status. Fig. 10C, Barplot shows the proportion and number of samples predicted to be BRCA1 -like or BRCA2-like across the subgroups. Fig. 10D, Proportion of various SV events in BRCA1 - like, BRCA2-like or non-HRD tumors across the subgroups. Fig. 10E, Replication of the Pareto front projection using the GEL (primary) cohort. Each dot represents the architecture profile of each tumor colored by IC subgroup. Fig. 10F, Activity of six SV signatures across the IC subgroups in metastatic tumors. Fig. 10G, Distribution of primary and metastatic tumors on Pareto front. Fig. 10H, Comparison of SV signatures in primary and metastatic tumors across the IC subgroups. Barplot shows the log fold change of each rearrangement signature between primary and metastatic tumors across the IC subgroups. Fig. 101, Transition vector corresponding to the difference in position on the Pareto fronts from (Fig. 10G) between the centroid of primary samples and the centroid of metastatic samples in each IC group. Fig. 10J, Replication of the Pareto front projection using the METABRIC (primary) cohort. Each dot represents the architecture profile of each tumour colored by IC subgroup. Fig. 10K, Forest plot shows the association between the proportion of archetypes and distant relapse free (DRF) survival, correcting for ERand HER2. Dots correspond to estimated hazard ratios and segments to 95% confidence intervals. Fig. 10L, Association between recurrence in the ER+Typical-risk samples and the distance to each archetype (spanning from 0 to 1 , linear regression, top), transcriptom ic proliferative and HRD LOH scores (linear regression, bottom) and histological type (IDC or ILC, fisher’s exact test, bottom) in the METABRIC dataset. Significance: P < 0.05 (*), P < 0.01 (**), and P < 0.001 (***). Fig. 10M, Differential pattern of relapse across ER+IC subgroups and by histology (IDC: invasive ductal carcinoma and ILC: invasive lobular carcinoma), illustrated by the cumulative (black) and annual (red) risk of relapse. Fig. 10N, ER+Typical IDC and ILC distribution on the Primary-Discovery Pareto front. Fig. 100, ER+Typical IDC and ILC distribution on the METABRIC Pareto front. SV, structural variant; WGD, whole-genome doubling; HRD, homologous repair deficiency; LOH, loss-of-heterozygosity; IDC, invasive ductal carcinoma; ILC, invasive lobular carcinoma.

[0021] Figs. 11A-11 M provide data and schematics for mutational processes underpinning metastatic breast cancer mirror primary breast cancer. Fig. 11 A, Scatterplot shows stability (silhouette index) and similarity (cosine distance) scores, y-axis, used to select the optimal number of SV signatures, x-axis, in the metastatic tumors. Gray shading indicates the number selected. Fig. 11 B, Heatmap shows the Spearman’s correlation between primary and metastatic SV signatures. Covariates on the top and right indicate whether the signature was identified in primary or metastatic tumors. Fig. 11 C, The proportion of the six SV signatures and 24 CNA signatures active in each metastatic tumor were projected onto a two-dimensional plane using PCA. Scatterplot represents each tumor as a dot along the first and second principal components, colored by IC subtype. The ternary plot closely resembles that for primary tumors. Fig. 11 D, SV and CNA signatures of primary and metastatic tumors projected onto two-space dimensional space using PCA. Fig. 11 E, Lollipop plots illustrating the correlation between various mutational features and the distance to each archetype with archetypes defined by Pareto (top) or NMF (bottom). Fig. 11 F, Heatmap shows concordance of archetypes defined by Pareto and NMF in primary, metastatic and primary+metastatic samples combined. Fig. 11 G, Forest plot shows the association (x-axis) between the three genomic archetypes and the breast cancer subgroups (y-axis) with archetypes definedby Pareto (orange) and NMF (purple). Fig. 11 H, Barplot shows proportion of samples with whole genome doubling (WGD) per subgroup in primary and metastatic tumors. Fig. 111, Violinplots show activity of the nine most recurrent CN signatures across the subgroups in metastatic tumors. Fig. 11 J, Schematic overview of how archetypes were identified in METABRIC in the absence of SV detection. Briefly, a random forest model was trained to map the CNA signatures in the primary discovery cohort to the first and second principal components generated considering both SV and CNA signatures together. The random forest model was then used to predict the first and second principal components in METABRIC using the same CNA signatures. Fig. 11 K, Scatterplot illustrates the principal components predicted in the discovery cohort using the CNA-only random forest model. Fig. 11 L, Scatterplot illustrates the principal components predicted in METABRIC using the CNA-only random forest model. Fig. 11 M, IC subtype proportions in invasive lobular carcinoma (ILC) and invasive ductal carcinoma (IDC) cases in METABRIC. Del, deletion; Tds, tandem duplications; Inv, inversions; T, translocations; HRD, homologous recombination deficiency; LOH, loss-of-heterozygosity; cLOH, chromosomal-scale LOH; fLOH, focal LOH; SV, structural variant; CNA, copy number aberrations; WGD, whole genome doubling.

[0022] Figs. 12A-12I provide data and schematics for cyclic amplifications are early mutational processes in ER+high-risk and HER2+breast tumors. Fig. 12A, Proportion (top) and number (bottom) of samples with at least one cyclic or complex non-cyclic amplification in primary or metastatic tumors. Fig. 12B, The density of SNVs occurring before amplification in primary (top) and metastatic (bottom) tumors. Boxplot represents median, 0.25 and 0.75 quantiles with whiskers at 1.5* the interquartile range. Fig.12C, Illustration showing copy number (CN) and SVs linking together disjoint segments in ecDNA (top), ratio of read depth in the tumor versus normal sample (middle) and location of oncogenes in ecDNA (bottom) in a representative primary IC2 tumor. Fig. 12D, Ratio of sequencing coverage in digested versus parental UCD65 (IC2) cell line in the predicted ecDNA region (dashed red line) compared to 1 ,000 null regions. Fig. 12E, Proportion of tumors within each IC subtype that harbor cyclic, complex non-cyclic or linear amplification in IC-specific oncogenes. Fig. 12F, Schematic for the genesis of cyclic amplifications. TC-NER, transcription-coupled nucleotide-excision repair. Fig. 12G, Thedensity of ER-induced R-loops in cyclic versus complex non-cyclic amplifications, h, The percentage of breakpoints that overlap ER-induced R-loops with (+) or without (-) E2 treatment. Error bars represent the standard deviation across three replicates. Fig. 121, The distance of each oncogene to the nearest ER-induced R-loop.

[0023] Figs. 13A-13L provide data and schematics for cyclic amplifications preferentially amplify IC-specific oncogenes. Fig. 13A, Proportion or number of samples with at least one cyclic or complex non-cyclic amplification in primary GEL tumors. Fig. 13B, Proportion or number of samples with at least one cyclic or complex non-cyclic amplification in HER2+ primary tumors stratified by ER positivity. Fig. 13C, Proportion or number of primary samples with at least one cyclic amplification according to JaBbA. Fig. 13D, Proportion of samples where AmpliconArchitect called a cyclic amplification but JaBbA called an alternative type of alteration. Colors indicate which alteration JaBbA called. Fig. 13E, Proportion of HER2+ primary tumours that harbor cyclic or linear amplification in ER+ High-risk-specific oncogenes (left), and the SV types (right). Fig. 13F, Proportion or number of samples with at least one cyclic or complex non-cyclic amplification in DCIS lesions. Left panel: DCIS cohort stratified by subgroup, right panel: DCIS cohort and additional samples from GEL stratified by sequencing method. Fig. 13G, Proportion of cyclic amplifications, stratified by subgroup, that amplify IC-specific or alternative oncogenes in primary tumors, both the discovery and replication (GEL) cohorts. The number of amplifications in each category are included on each bar. Fig. 13H, Number ecDNA involving more than one IC-specific oncogene. Fig. 131, Number of oncogenes per megabase involved in ecDNA in each subgroup. Boxplot represents median, 0.25 and 0.75 quantiles with whiskers at 1.5x interquartile range. Fig. 13J, Ratio of oncogenes amplified on ecDNA compared vs. oncogenes in the IC-specific cytoband per megabase. Fig. 13K, Proportion of ER+ Typical-risk ecDNA that incorporate each oncogene. Fig. 13L, Proportion of each archetype in ER+ High-risk, ER+ Typical and ER+ Typical containing ecDNA tumors.

[0024] Figs. 14A-14K provide data and schematics complex amplifications are enriched in ER+ High-risk and HER2+ breast cancer. Fig. 14A, Barplot shows the proportion of samples with cyclic, complex non-cyclic, linear or no amplification as identified in the current study vs. as identified previously. Fig. 14B, Barplot shows theproportion of cyclic and complex non-cyclic amplifications that are coincident with SV clusters, considering various minimum thresholds (x-axis), identified using orthogonal methods. Fig. 14C, Barplots illustrate the proportion of HER2+ metastatic tumors that harbor cyclic or linear amplification in ER+ High-risk-specific oncogenes (left), and the SV types (right). One HER2+ tumor possessed 11 q13 linear amplification. Fig. 14D, Barplots illustrate the proportion of ecDNA in BRCA1 -like, BRCA2-like or non-HRD tumors stratified by subgroup as called using Am pl icon Architect (top) and JaBba (bottom). Fig. 14E, Barplots illustrate the proportion of tumors in which ecDNA was called from both downsampled and full coverage bams (green); full coverage WGS only (blue); downsampled bams only (brown); and neither. Fig. 14F, Scatterplot shows correlation between number of clock-like mutations (y-axis) with age at diagnosis (x-axis). Tumors are colored by subgroup. Fig. 14G, Tracks depict ecDNA: copy number and structural variants linking together disjoint segments (top), ratio of read depth in the tumor vs. normal sample (middle) and location of oncogenes in ecDNA (bottom) in a representative primary IC6 tumor. Fig. 14H, Ratio of sequencing coverage in digested vs. parental LICD12 (IC6) cell line in the predicted ecDNA region (dotted red line) compared to 1 ,000 null regions. Fig. 141, Barplot showing proportion of samples with MYC and PVT1 amplifications in primary and metastatic breast tumors. Figs. 14J-14K, Scatterplot shows the correlation between Iog10 MYC mRNA abundance and MYC copy number in primary (Fig. 14J) and metastatic (Fig. 14K) tumors. Pink dots represent IC9 tumors and correlation is Spearman’s p. DEL, deletion; DUP, duplication; INV, inversion; TRA, traversion; FDR, false discovery rate.

[0025] Figs. 15A-15F provide cyclic amplifications are maintained in metastatic tumors. Figs. 15A-15B, Number of recurrent oncogenes per Mbp for each subgroup in primary (Fig. 15A) and metastatic (Fig. 15B) tumors. Figs. 15C-15D Proportion of cyclic amplifications, stratified by subgroup, that amplify IC-specific or alternative oncogenes in metastatic (Fig. 15C) and DCIS (Fig. 15D) lesions. The number of amplifications in each category is included on each bar. Fig. 15E Proportion of metastatic tumours within each IC subtype that harbor cyclic, complex non-cyclic or linear amplification in the IC-specific oncogenes. The number of tumors within each subtype are indicated at the top of eachsubpanel. Fig. 15F Two representative examples of ER+ High-risk DCIS lesions harboring ecDNA containing at 8p11 (IC6). AMP, amplification.

[0026] Figs. 16A-16H provide data and schematics for elevated replication stress in TNBC, ER+ High-risk and HER2+ tumors. Figs. 16A-16B, Replication stress signature stratified by IC subgroup (left) or by histology within ER+ Typical-risk subgroup IDC and ILC (right) in the TCGA (Fig. 16A) and METABRIC (Fig. 16B) datasets. FDR adjusted p- values are reported. Figs. 16C-16D, Replication stress signature stratified by IC subtypes in TCGA (Fig. 16C) and METABRIC (Fig. 16D) datasets. Fig. 16E, Pareto projection of METABRIC tumors, colored by replication stress. Fig. 16F, Replication stress signature in HER2+ tumors; IC1 , IC6, and IC9 subtypes of ER+ High-risk; and TNBC IC10 stratified by presence of ecDNA. ER+ High IC2 was excluded due to lack of sample size (n = 2 for ecDNA+ IC2). Fig. 16G-16H, cGAS / STING signature stratified by IC subgroup (left) or by histology within ER+ Typical-risk subgroup IDC and ILC (right) in the TCGA (Fig. 16G) and METABRIC (Fig. 16H) datasets. FDR adjusted p-values are reported. In a-d, effect sizes (ES) and FDR-adjusted p-values from Mann-Whitney Rank Sum test are shown. In f, ES and p-values from linear regression correcting for cohort are shown. Additionally, the amplicon copy number was corrected for amplicon-driven HER2+ and ER+ Typical tumors.

[0027] Figs. 17A-17F provide data and schematics for mechanism of ecDNA genesis in ER+ breast tumors. Fig. 17A, Type-I IFN signature in HER2+ tumors ; IC1 , IC6, and IC9 subtypes of ER+ High-risk; and TNBC IC10 stratified by presence of ecDNA. ER+ High IC2 was excluded due to lack of sample size (n=2 for ecDNA+ IC2). Effect Size (ES) and p-values from linear regression correcting for cohort are shown. Additionally, the amplicon copy number was corrected for amplicon-driven HER2+ and ER+ Typical tumors. Fig. 17B, Proportion of SVs in cyclic and non-cyclic amplifications that are translocations in ER+ High-risk and IC10 primary and metastatic tumors. Boxplot represents median, 0.25 and 0.75 quantiles with whiskers at 1 ,5x interquartile range. Fig. 17C, Density of APOBEC3B and ESR1 ChlP-Seq peaks within cyclic and complex non- cyclic amplifications in metastatic tumors. Fig. 17D, Difference in number of R-loops between A3B knockout (KO) wildtype (WT) MCF10A cell lines overlapping cyclic or non- cyclic amplifications at baseline or after A3B activation (PMA treatment) considering onlyER+ High-risk or HER2+ tumors. Fig. 17E, Barplot highlights the proportion of ecDNA+ or ecDNA- tumors, as called using AmpliconArchitect (top) and JaBba (bottom), that harbor at least one large duplication (>100kbp) stratified by subgroup. Fig. 17F, Boxplot shows replication timing weighted average in cyclic vs non-cyclic amplifications stratified by subgroup. ES and p-values from Mann-Whitney Rank Sum test.

[0028] Figs. 18A-18L provide data and schematics of model for ER-Induced R-Looks in ecDNA genesis. Fig. 18A, Simplified schematic illustrating model. Blue letters correspond to figures in Figs. 18B-18L. Fig. 18B, Number of translocations in cyclic vs. non-cyclic amplifications across subgroups. Boxplot represents median, 0.25 and 0.75 quantiles with whiskers at 1.5x interquartile range. Fig. 18C, ESR1 mRNA abundance in cyclic amplification-positive vs. -negative (top) and non-cyclic amplification-positive vs. - negative (bottom), TOGA or metastatic tumors stratified by the IC subgroups, considering ER+ High-risk and HER2+ subgroups. Odds ratio from logistic regression correcting for tumor purity and error bars represent 95% confidence intervals. Fig. 18D, Density of APOBEC3B and ER ChlP-Seq peaks within cyclic and complex non-cyclic amplifications in primary tumors. Fig. 18E, ER early signaling transcriptional signature in DCIS and primary ER+ Typical vs. ER+ High-risk tumors. Figs. 18F-18G, Density of ER-induced R- loops in cyclic and complex non-cyclic amplifications stratified by IC subgroups in primary (Fig. 18F) and metastatic (Fig. 18G) tumors. Fig. 18H, Density of all R-loops in cyclic vs. non-cyclic amplifications in primary and metastatic tumors. Fig. 181, Difference in number of R-loops between A3B knockout (KO) wildtype (WT) MCF10A cell lines overlapping cyclic or non-cyclic amplifications at baseline or after A3B activation (PMA treatment). Figs. 18J-18K, Median distance between a translocation and its closest ER-induced R- loop considering translocations within or outside cyclic amplifications in primary (Fig. 18J) and metastatic (Fig. 18K) tumors. Fig. 18L Percent of breakpoints that overlap any R-loop with (+) or without (-) E2 treatment. Error bars represent the standard deviation across three replicates, m) Proportion of samples with or without ecDNA stratified by inferred APOBEC3B germline copy number. The total number of samples is included at the top of each bar. In b and j-l fold change (FC) and p-values or false discovery rates (FDR) are from Mann-Whitney Rank Sum test. In Figs. 18D-18I, effect sizes (ES) are the differencein medians and p-values are from Mann-Whitney Rank Sum test. BER, base-excision repair; TC-NER, transcription-coupled nucleotide excision repair; E2, estrogen.

[0029] Figs. 19A-19D provide data and schematics for complex alterations contribute to IC-Specific immune escape. Fig. 19A, Schematic of TME subtypes and select immune escape pathways. Fig. 19B, Comparison of TME subtypes by IC subgroup. The number of tumors in each subgroup is indicated on the top of each bar. Fig. 19C, Top: proportion of primary and metastatic samples in each IC subgroup with GIE. Bottom: proportion of samples with alterations in each pathway stratified by IC subgroup and stage of progression. Fig. 19D, Proportion of alteration types in primary and metastatic samples for each of the immune escape pathways.

[0030] Figs. 20A-20G provide data and schematics for IC subgroups harbor distinct TMEs. Fig. 20A, Schematic illustrating additional transcriptom ic profiles and overlap with genomic profiles induced in Fig. 3A. Fig. 20B, Mean proportion of different cell types from IMC data by TME subtypes. The Wilcoxon test significance was reported above each comparison as follow: ns: not significant, P < 0.05 (*), P < 0.01 (**), P < 0.001 (***), and P < 0.0001 (****). Fig. 20C, Proportion of TME subtypes in primary and metastatic samples for the ER+ High-risk ICs and IC5 (HER2+) by ER status. Fig. 20D, Proportion of TME subtypes for primary samples (METABRIC) in ER+ Typical invasive IDC and ER+ Typical ILC. Fig. 20E, Mean proportion of fibroblasts and T cells in TNBC samples with IMC proteomic data obtained from bootstrapping (n = 1000). Fig. 20F Proportion of TME subtypes for primary and metastatic samples stratified by ER status. Fig. 20G, Proportion of TME subtypes for primary samples and liver metastases by groups. IMC, imaging mass cytometry; SMA, smooth muscle actin.

[0031] Figs. 21A-21 G provide data and schematics for validation of TME subtypes and genetic mechanisms of immune escape. Fig. 21A, Y-axis represents the proportion of selected cell types obtained from in silico deconvolution using CIBERSORTx and a custom reference matrix from a single cell breast cancer transcriptom ic dataset. Comparing proportions across TME subtypes for primary samples from TCGA, primary samples from METABRIC, and metastatic samples. The Wilcoxon test significance was reported above each comparison as follow: ns: not significant, P < 0.05 (*), P < 0.01 (**), P < 0.001 (***), and P < 0.0001 (****). Fig. 21 B, Proportion of TME subtypes for primarysamples from METABRIC by groups. Fig. 21 C, Proportion of immune phenotypes for primary samples from TCGA by subgroups. Fig. 21 D, Proportion of immune phenotypes for ER+ High-risk primary samples from TCGA by IC subtypes. Fig. 21 E, CIBERSORTx estimates for cancer associated fibroblasts (CAFs) and T cells in triple negative tumors in TCGA and METABRIC. Fig. 21 F, Proportion of immune-depleted samples (D-desert or F-fibrotic), where primary samples were downsampled to the number of metastatic samples; 1000 permutations were performed. Fig. 21 G, Forest plot shows the association between neoantigens and GIE, or neoantigens and GIE + log(SNV). The points correspond to estimated coefficients for GIE and the segments correspond to their 95% confidence intervals. CAF, cancer-associated fibroblasts; SNV, single nucleotide variant.

[0032] Figs. 22A-22F provide data and schematics for genetic mechanisms of immune escape in IC subgroups. Fig. 22A, Proportion of primary and metastatic samples in each ER+ High-risk subtype with genetic immune escape (GIE) alterations, where values correspond to the number of pathways altered (top). Proportion of samples with alterations in each pathway stratified by IC subtype and disease stage (bottom). Fig. 22B, Proportion of primary ER+ High-risk samples with co-amplification of IDO1 with FGFR1 or ZNF703 by IC subgroup. Fig. 22C, Proportion of IC6 tumors with immune enriched (IE or IE / F) or immune depleted (D or F) TME subtypes stratified by the co-amplification of IDO1 with FGFR1 or ZNF703 in METABRIC and TCGA. Fig. 22d, Proportion of ER+ Typical IDC and ILC with GIE (left). Odds ratio and p-value from Fisher’s exact test. Proportion of pathways altered in IDC and ILC with GIE (right). Fig. 22E, Number of alteration in immune escape pathways for primary and metastatic samples, normalized by number of samples with alterations. Fig. 22F, Odds ratio for the frequency of GIE pathway alterations, comparing metastatic to primary samples. Background shading indicates FDR adjusted p-values (Fisher’s exact test). The color of the dot represents the direction and magnitude of the odds ratio while the dot size indicates the number of samples with a GIE in each pathway (y-axis). LOH, loss-of-heterozygosity.

[0033] Figs. 23A-23B provide data and schematics for genomic and microenvironmental evolution of breast cancer subgroups. Fig. 23A, Schematic summary of the genomic and microenvironmental characteristics of the three dominant genomic archetypes in breast cancer. Fig. 23B, Temporal changes in genomic stability, ERsignaling and immune enrichment from pre-invasive, primary invasive to metastatic disease across subgroups.DETAILED DESCRIPTION

[0034] The disclosure generally describes systems and methods for assessing breast cancer, which can be utilized as a diagnostic for determining a treatment regimen. Integrative clustering subtyping of breast cancer is a way to cluster breast cancers based on genetic aberrations (see, C. Curtis, et al., Nature. 2012 Apr 18;486(7403):346-52; and 0. M. Rueda, et al., Nature. 2019 Mar;567(7748):399-404; the disclosures of which are hereby incorporated by reference). The clustering method found that breast cancer falls within one of eleven integrative cluster (IC) subtypes. The eleven ICs can also be grouped into one of four categories: HER2+, ER+ high risk, ER+ typical risk, and triple negative breast cancer (TNBC). Experimentation and analysis was performed by assessing globally assessing genetic features and signatures of thousands of breastcancers to comprehensively understand the genomic archetypes. In sum, three extremes of genomic archetypes were found: Archetype 1 (ER+ cancer with typical risk), Archetype 2 (ER+ with high risk cancer and HER2+ cancer), and Archetype 3 (TNBC). All breast cancer fall within the triad of these archetype extremes. As such, a breast cancer can be mapped within a ternary map with each archetype at a vertex. Based on where a breast cancer is mapped, and in accordance with several embodiments of the disclosure, a breast cancer can be diagnosed. And in many embodiments, a treatment regiment is determined based on a genomic archetype enrichment score.

[0035] Several embodiments are directed to determining an archetype enrichment score. Generally, genetic data of a breast cancer is computationally assessed to yield an archetype enrichment score. In some embodiments, the score is used to determine a treatment regimen. In some embodiments, the score in used to determine whether a further follow-on assessment is to be performed, which can further dictate the treatment regimen. In some embodiments, a co-assessment is performed, which can further dictate the treatment regimen.Assessment of breast cancer and treatment regimens

[0036] Provided in Fig. 1 is a method to assess a cancer, which can be utilized as a diagnostic. The method can be utilized to determine a genomic archetype of the breast cancer. The method can also be utilized to determine a treatment regimen. And further, the method can optionally administer a treatment to a patient in accordance with the treatment regimen.

[0037] Method 100 can begin by assessing (101 ) a breast cancer to yield its genomic archetype. As mentioned, a genomic archetype is a comprehensive understanding the genetic features that give rise to the cancer. To assess the breast cancer, genetic material can be collected from a breast cancer. Sources of genetic material include (but are not limited to) DNA extracted from cancer cells, RNA extracted from cancer cells, cell-free DNA collected from a non-cellular source of a patient, or cell-free RNA extra collected from a non-cellular source of a patient. Cancer cells can be collected from a biopsy of the patient, which can be a primary tumor, a metastatic tumor, a recurrent tumor, etc. Examples of non-cellular sources for collecting cell-free nucleic acids include blood, plasma, lymph, urine, stool, saliva, mucous, etc. The genetic material can assessed via microarray / SNP array based copy number inference, microarray based gene expression, whole genome sequencing (WES / WGS), targeted (panel) sequencing, RNA-sequencing, targeted (capture) RNA-sequencing, exome sequencing, etc.

[0038] Using genetic data, an archetype can be determined. In some implementations, an archetype map is utilized, which is generated using a discovery cohort of breast cancers. As is described in greater detail in reference herein, in some embodiments, nonnegative matrix factorization (NNMF) considering several features of genetic data from each breast cancer is utilized to generate a map (an example is provided in Fig. 2B). In some embodiments, the features utilized for assessing the genetic data of each breast cancer include: single nucleotide variants, genetic amplifications, and structural variants of genomes. A list of features is shown in Fig. 6C.

[0039] In several embodiments, using a generated archetype map, the genetic data of the breast cancer is assessed via NNMF to yield a projection on the archetype map. The map projection is an archetype enrichment score that can be utilized to provide a quantitative assessment of the cancer.

[0040] Alternative to non-negative matrix factorization, a prediction computational model using machine-learning can be utilized to predict an archetype enrichment score. Machine-learning prediction modeling can facilitate the determination of an archetype enrichment score by reducing the amount of genetic data and / or features necessary to yield the comprehensive archetype enrichment score. By requiring less genetic data needed, the clinical utility of an archetype enrichment score is expanded. Simpler molecular protocols (e.g., exome sequencing, microarray, copy number aberration assessment, gene expression assessment, etc.) can be utilized as a surrogate to yield a full, comprehensive genomic assessment. The machine-learning model can be trained using breast cancer samples of a cohort or patients, in which the model is trained using reduced genetic data input, reduced features, and the comprehensive archetype enrichment score determined by NNMF.

[0041] In several embodiments, using a trained computational machine-learning model, the genetic data of the breast cancer is assessed via the model to yield an archetype enrichment score that can be utilized to provide a quantitative assessment of the cancer.

[0042] Method 100 can also determine (103) a treatment regimen based on the genomic archetype of the breast cancer. Various treatment regimens that would be beneficial are based on genetic vulnerabilities associated with genomic archetypes (Fig. 2C). Upon determining an archetype enrichment score of the breast cancer, the treatment regimen can be determined. In some implementations, the proximity of a breast cancer archetype enrichment score to an archetype extreme (e.g., a vertex of the archetype map) is utilized to determine the treatment regimen. In some implementations, when a breast cancer archetype enrichment score is between two or three extreme archetypes, one or more of treatment regimens associated with the extreme archetype is administered. In some implementations, when a breast cancer archetype enrichment score is between two or three extreme archetypes, two or more of treatment regimens associated with the extreme archetype is administered. In some implementations, when a breast cancer archetype enrichment score is between three extreme archetypes, the three treatment regimens associated with each extreme archetype is administered. Thresholds and / orranges of archetype enrichment scores can be utilized to further determine treatment regimens.

[0043] Further assessment (or co-assessment) of the breast cancer genetic data can also be performed. It has been determined that some breast cancers can be further assessed for determining an optimal treatment regimen by performing further analysis. Optional further assessments can include (but are not limited to): classification of the breast cancer into an integrative cluster, assessing the breast cancer for replication stress, assessing the cancer for homologous recombination deficiency, and assessing the breast cancer for genetic mechanisms of immune escape. The optional further assessments can further stratify a breast cancer to determine targeted treatment options.

[0044] Method 100 optionally administers (105) the treatment regimen to the patient. The treatment regimen to be administered is determined based on the archetype enrichment score and optionally another genetic assessment that is performed. Options for treatments include (but are not limited to) chemotherapy, surgery, endocrine therapy, Ataxia Telangiectasia and Rad3-related (ATR) protein kinase inhibitors, WEE1 inhibitors, checkpoint kinase (CHK) inhibitors, APOBEC inhibitors, Uracil-DNA glycosylase (UNG) inhibitors, topoisomerase II inhibitors (TOP2) inhibitors, poly [ADP-ribose] polymerase (PARP) inhibitors, immune checkpoint inhibitors, targeted antibody treatments, and immunotherapies.

[0045] Provided in Fig. 2A are examples of methodologies to yield a genomic archetype enrichment score for assessing breast cancer. The method can be utilized to assess the breast cancer of patient to determine a treatment regimen. As shown in Fig. 2A, genetic data that is utilized as input are structural variants and copy number variations to yield a genomic architecture map. To yield this map, the copy number variations and structural variants are determined from genetic data of a cohort of breast cancer patients. The genetic data can be genomic sequencing or other similar data that provide can provide comprehensive analysis of the copy number variations and structural variants (e.g., genomic sequencing). Feature data is extracted from the copy number variation data and structural variant data and utilized as input in an NNMF. The output is used to yield an archetype map, which a projection of all the breast cancers of the cohort (Fig.2B). Projections can be yielded by any appropriate method, such as principal component analysis (PCA), t-SNE, UMAP, MDS, etc.

[0046] When plotting breast cancer archetypes, the archetype map yielded a triangular shape. Further, when looking at where particular molecular subtypes of breast cancer were projected, it was noted that each vertex was heavily represented by one of the subtypes. Specifically, at a first vertex was represented by ER+ cancers with typical risk, a second vertex was represented by HER+ cancers and ER+ cancers with high risk, and a third vertex was represented by TNBC cancers. Accordingly, archetype extremes can be labeled as: archetype 1 (ER+ cancers with typical risk), archetype 2 (HER+ cancers and ER+ cancers with high risk), and archetype 3 (TNBC). As can be readily appreciated, the numbering and ordering of vertices and archetypes is merely nominal and any numbering or ordering can be utilized.

[0047] In many embodiments, a genomic archetype map is utilized to yield an archetype enrichment scoring system. Provided in Fig. 2B is an example of a ternary map with archetype extremes at each of the three vertices. Based on a projection within the map, a score that is representative of the enrichment of an archetype is computed. In the example in Fig. 2B, scores are normalized between 0 and 1 dependent on location between two vertices.

[0048] Fig. 2A also provides an alternative assessment for archetype enrichment score in which a computational machine-learning model was trained to be a surrogate methodology to predict archetype enrichment score. In this model, copy number aberrations (CNA) from exome sequencing were utilized to train a random forest model to predict archetypes. Data features from CNA signatures are extracted and utilized as input in the model to yield an output of an archetype enrichment score. The use of CNA signatures is useful because exome sequencing (which is cheaper and more often performed) instead of whole genome sequencing can be utilized to predict archetypes. CNA features can also be extracted or derived from other genetic data, such as transcriptome data, microarray data, cfDNA data, cfRNA data, etc.

[0049] Although particular methods are depicted in Fig. 2A, it should be understood that alternative methodologies that yield a genomic archetype of breast cancer or a surrogate assessment of genomic archetype can be utilized. For example, the right sideof Fig. 2A depicts a computational machine-learning model that is trained to predict an archetype enrichment score as a surrogate assessment. In this example, the computational machine learning model using copy number aberrations (CNA) as input, but any genetic information can be utilized to train a machine-learning model to predict an archetype enrichment score in which the trained model provides adequate prediction. Trained models can be selected by various means, such as prediction sensitivity, prediction accuracy, ease of collecting / generating data as input, computation power needed, etc.

[0050] Several embodiments are directed towards the use of an archetype enrichment score to determine a treatment regimen. The archetype extremes are associated with genetic vulnerabilities, and as such treatment regimens that target these vulnerabilities can be determined based on an archetype enrichment score. Scores proximate to vertex can be provided the treatment regimen associated with that extreme vertex. Provided in Fig. 2C are treatment regimens associated with archetype 1 (ER+ cancers with typical risk), archetype 2 (HER+ cancers and ER+ cancers with high risk), and archetype 3 (TNBC).

[0051] Treatment regimens for archetype type enrichment scores proximate to archetype 1 can include endocrine therapy and chemotherapy. Treatment regimens for archetype 2 are associated with greater deficiencies with control of replication and cell division. Treatment regimens for arche type are associated with genomic instabilities.

[0052] Provided in Fig. 2D are optional assessments that can be performed. In some embodiments, an assessment is a further assessment in which the assessment is performed based on archetype enrichment score result. Alternatively, an assessment is a co-assessment, in which the assessment was not performed but based on archetype enrichment score result, but is can still be combined with an enrichment score to further stratify the treatment regimen.

[0053] Method 200 optionally classifies a breast cancer into a molecular class. In some embodiments, the molecular class is an integrative cluster. In some embodiments, the molecular class is ER+ with typical risk, ER+ with high risk, HER2+, or TNBC. One method to classify a breast cancer is the ENiClust computational model (Fig. 2E), which is further described in the Experimental Data and Results section. Molecular classes are aspectrum across the archetype enrichment map, an archetype enrichment score may not perfectly represent a molecular class. Further, archetype includes ER+ and HER2+ breast cancer, which are treated with separate therapies (e.g., endocrine therapy vs. HER2- targeted therapy). Therefore, there is benefit in classifying into a molecular class to stratify a patient into determining therapy.

[0054] Method 200 also optionally assesses (203) for a replication stress signature, which can be treated with ATR inhibitors, WEE1 inhibitors, CHK inhibitors, APOBEC inhibitors, UNG inhibitors, TOP2 inhibitors, or a combination there. A replication stress signature can assess expression of one or more of the following genes: NAT10, DDX27, ZNF48, C8ORF33, MOCS3, and MPP6.

[0055] Method 200 also optionally assesses (205) homologous recombination deficiency (HRD), which can be treated with PARP inhibitors. HRD can be predicted with a pan-cancer classifier that can discriminate between BRCA1 -like and BRCA2-like HRD.

[0056] Method 200 also optionall assess (207) for genetic mechanism of immune escape, which means a treatment regimen would not benefit from an immune checkpoint inhibitor, an antibody-based therapy that relies on an immunological response, or an immunotherapy. Disruptions of genes (e.g., as determined by presence of a mutation that yields a missense, nonsense, splicing disruption, etc) involved in immune response likely will not benefit from immune-based treatments. Fig. 19A has list of immune genes that when disrupted lead to GIE.

[0057] Dosing and therapeutic regimens can be administered appropriate to the cancer to be treated, and can be determined by preclinical and clinical studies.

[0058] In some embodiments, drugs are administered in a therapeutically effective amount as part of a course of treatment. As used in this context, to "treat" means to ameliorate at least one symptom of the disorder to be treated or to provide a beneficial physiological effect. For example, one such amelioration of a symptom could be reduction of tumor size.

[0059] A therapeutically effective amount can be an amount sufficient to prevent reduce, ameliorate or eliminate the symptoms of breast cancer. In some embodiments, a therapeutically effective amount is an amount sufficient to reduce the growth and / or metastasis of a cancer.

[0060] While specific examples of methods for assessing mutational burden and determining a treatment regimen are described above, one of ordinary skill in the art can appreciate that various steps of the method can be performed in different orders and that certain steps may be optional according to some embodiments of the disclosure. As such, it should be clear that the various steps of the process could be used as appropriate to the requirements of specific applications.EXPERIMENTAL DATA AND RESULTS

[0061] The embodiments of the disclosure will be better understood with the several examples provided. Experimental data and results support the various embodiments of treating breast cancer based on its archetype.Complex rearrangements fuel ER+ and HER2+ breast tumours

[0062] We previously defined eleven subtypes of breast cancer on the basis of integrative clustering (IC) of genomic and transcriptional profiles, and demonstrated their distinct prognosis and relapse trajectories (See, C. Curtis, et al., Nature. 2012 Apr 18;486(7403):346-52; and O. M. Rueda, et al., Nature. 2019 Mar;567(7748):399-404; the disclosures of which are hereby incorporated by reference). Among patients with ER+ cancer (80% of cases), one-quarter had a 45% chance of distant recurrence two decades post-diagnosis. This ER+ ‘high-risk’ subgroup, corresponding to IC1 , IC2, IC6 and IC9 subtypes, is enriched for luminal B tumours harbouring focal oncogene amplification and overexpression, similar to ERBB2-amplified tumours (IC5, 10-15%). Moreover, genes within these amplicons mediate resistance to hormonal therapy. TNBC comprises genome-unstable basal-like IC10 and IC4ER- tumours, the latter with relapse risk that persists beyond 5 years.

[0063] While IC subgroups improve relapse prediction and define new drivers, their origins, evolution, and tumor-immune microenvironments (TME) remain unknown. To investigate, we assessed the genomic architecture and microenvironmental composition of breast tumors from a meta-cohort of 1 ,828 tumors spanning pre-invasive DCIS, primary and metastatic lesions, profiled using whole genome sequencing (WGS) and transcriptome sequencing. We further implemented a machine learning framework todetermine IC subtypes from DNA-based profiles alone. Our analyses reveal three primary genomic archetypes of breast cancer: (i) TNBC: ICs (IC10, IC4ER-); (ii) typical-risk ER+ / HER2- (IC3, IC4ER+, IC7, IC8); and (iii) high-risk ER+ / HER2- (IC1 , IC2, IC6, IC9) and HER2+ (IC5) (referred to as ER+ High-risk / HER2+). The latter group is characterized by early, recurrent amplifications, including ecDNA owing to APOBEC3B-editing at ER- induced R-loops. These genomic patterns, accompanied by variable TMEs, implicate complex rearrangements as a major driver of immune escape and highlight new therapeutic vulnerabilities in aggressive subgroups.ResultsEvolution of the IC subgroups

[0064] The mutational processes underlying breast cancer initiation and progression are incompletely understood. Herein we uniformly processed 1 ,828 samples from DCIS (n = 406), primary (n = 702) and metastatic (n = 720) lesions using a harmonized, state- of-the-art bioinformatics pipeline to identify single nucleotide variants (SNVs), copy number aberrations (CNAs), structural variants (SVs), ecDNA and mutational signatures (Fig. 3A, Figs. 4A-4C). Owing to shallow coverage of the archival DCIS cohort, SNVs and SVs were not called. Additionally, we used the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) cohort of primary invasive tumors (n = 1 ,894) with both RNA and DNA profiles and about 20 years of clinical follow-up. To our knowledge, this represents the largest collection of uniformly processed breast tumors spanning all disease stages.

[0065] Although the ICs predict distant relapse and delineate genomic drivers current methods fail to accurately capture them using DNA profiles alone. Accordingly, we developed Ensemble Integrative Clustering (ENiClust), which reliably infers IC subtypes from whole-exome (WES) or WGS, across all stages of disease (Fig. 3B). The final ensemble model yields a nine-class prediction (Fig. 3B), which is further split into ten based on ER status of IC4 (i.e. IC4ER+ and IC4ER-). These ten classes comprise four clinically distinct IC subgroups - TNBC (IC10, IC4ER-), HER2+ (IC5), ER+ Typical-risk (IC3 / 7, IC4ER+, IC8), and ER+ High-risk (IC1 , IC2, IC6, IC9). Throughout we refer to HER2+ tumors as those classified as IC5, enriching for HER2 amplification. ENiClustoutperformed iC10 DNA-only (Fig. 4D) and improves patient stratification, with High-risk tumors exhibiting worse distant recurrence-free survival (METABRIC; Figs. 3C to 3E, Figs. 4D to 4F). Thus, ENiClust identifies clinically meaningful subgroups with distinct biology.

[0066] Using ENiClust, we interrogated the distribution of ICs across disease stages. DCIS was enriched for IC5 tumors (Fisher’s exact test P = 2.98x10-6; Fig. 3F). ER+ High- risk ICs were enriched amongst metastatic tumors, consistent with their elevated relapse risk (Fig. 3D, Fig. 5A). IC10 basal-like tumors were depleted in the metastatic cohort, potentially due to differences in ancestry (Figs. 5B to 5D). The ICs were largely stable from primary to metastasis (concordance = 71 .8%; Figs. 5E to 5F).

[0067] There was an increased proportion of Luminal B versus Luminal A from pre- invasive to primary (ALumB / LumA+B=+11 %) and primary to metastasis (ALumB / LumA+B=+29%; Fig. 5G). Amongst primary tumors, ER signaling in ER+ High- risk was more akin to HER2+ / ER+ tumors and significantly lower than ER+ Typical-risk tumors (Fig. 5H), with no difference between primary and metastatic tumors (Fig. 5I). Compared to ER+ Typical-risk, ER+ High-risk was enriched amongst endocrine therapyresistant patients (OR > 5.58, P < 0.03; Fig. 4G). In a clinical trial (NCT00651976) in early- stage ER+ breast cancer, High-risk tumors had decreased proliferation score with letrozole treatment but remained significantly higher than Typical-risk tumors (P < 0.02; Fig. 4H). Thus, ER+ High-risk tumors may experience persistent proliferation despite endocrine treatment. New therapies (SERDs, PROTACs) that more fully suppress proliferation might particularly benefit this subgroup.Early IC-Specific SVs Fuel Progression

[0068] The IC subtypes have distinct CNA landscapes (Fig. 5J), but their SV landscape and evolution have not been investigated. Leveraging ENiClust, we found that the IC subgroup-specific genomic landscape of breast cancer is consistent throughout disease progression despite an increased burden of alterations (Fig. 6A; Fig. 7A to 7B). Both HER2+ and ER+ High-risk primary and metastatic tumors display characteristic sharp increases in SV burden at their respective recurrently amplified loci (IC5: 17q12; IC6: 8p11 ; IC2: 11q13; IC1 : 17q23). The peak of SV burden at 17q12 (ERBB2 / HER2),suggests that ERBB2 amplification is fueled by complex alterations, such as ecDNA. The mutational burden in primary ER+ Typical-risk tumors was minimal (Fig. 41) but increased in metastatic disease (Fig. 6A), in part due to treatment (Fig. 7C). IC10 and IC4ER- tumors exhibit diffuse genome-wide instability with an elevated SV burden, although the latter show an attenuated pattern and harbor fewer pathogenic SVs and alterations in DNA repair pathways, confirming previous reports (Fig. 7D to 7E). Across metastatic sites, the cumulative burden of alterations was higher in lung and subcutaneous metastases and lower in soft-tissue and in-breast recurrences (Fig. 7F). These subgroup-specific alterations were seen in DCIS (Fig. 7A), emphasizing early oncogene addiction and mechanisms of malignant transformation.

[0069] Next, we characterized CNA and SV signatures in 702 primary breast tumors, replicating the 24 CN and 6 rearrangement signatures (RS) previously reported (Figs. 8A to 8C). RS3 and RS5 (associated with HRD; Fig. 8D), and CN17, were enriched in IC10 tumors, while RS4, RS6 (associated with complex amplifications) and CN7 were enriched in ER+ High-risk and HER2+ tumors (Figs. 7G to 7H; Figs. 8E to 8G). ER+ Typical-risk tumors were enriched for CN1 (associated with diploid genomes; Figs. 8D to 8E).

[0070] Projected on a two-dimensional plane (Fig. 9A to 9B), the architectural profiles follow a continuum and form a polyhedron reminiscent of Pareto optimum theory, which illustrates trade-offs between biological tasks. Primary breast cancers map onto three dominant genomic archetypes (Fig. 9C to 9F): TNBC-enriched, ER+ Typical-enriched and ER+ High / HER2+-enriched. Tumors dominated by a single mutational process are proximal to a vertex, while those characterized by multiple processes cluster at the center (Fig. 7I; Fig. 6B). The TNBC-enriched archetype was positively correlated with genomic instability, HRD and APOBEC-editing SNVs (Fig. 6C; Fig. 9G). Compared to ER+ High- risk, HER2+ tumors were enriched for tyfonas (Fig. 7J). ER+ High / HER2+-enriched archetype was positively correlated with complex amplifications, reactive oxygen species and APOBEC-associated SNVs harboring co-amplification of multiple cytobands (Fig. 10A). In contrast, the ER+ Typical-enriched archetype negatively correlated with most genomic features.

[0071] Tumors predicted to be BRCA-like based on germline or somatic genomic features map to the TNBC-enriched archetype (Fig. 10B). Indeed, both BRCA1-like andBRCA2-like ER+ and ER- tumors demonstrated significantly higher TNBC-archetype scores than non-HRD tumors, and HRD-like ER+ High-risk tumors were closer to TNBC- enriched archetype than their non-HRD-like counterparts (OR = 5.09; P = 6.5x10-4). Additionally, the mutational patterns of BRCA1 -like and BRCA2-like ER- and ER+ tumors were highly concordant (Fig. 9H to 9I). Notably, while 43.6% of TNBC tumors were HRD- like, 13.2% of ER+ High-risk tumors were also predicted to be HRD-like, with the majority being ER+ High-risk IC1 or IC9 (OR = 4.43; P = 0.03; Fig. 10C, Fig. 9J). Indeed, while foldback inversions (FBI) and pyrgos were enriched in TNBC (FBI: 17.3%, P = 2.00x10- 3; pyrgos: 18.8%, P = 9.33x10-4), these mutational events were also observed in ER+ tumors (5.1 % and 4.1 % respectively; Fig. 10D). These data reinforce multiple mechanisms of genome instability in TNBC 25 that also impact a subset of ER+ tumors.

[0072] The three genomic archetypes replicated in an independent cohort of 2,229 primary tumors from Genomics England (GEL) (Fig. 10E). Overall, the genomic landscape of primary breast tumors falls along a continuum with mutational patterns captured by three main genomic archetypes, namely, genome-stable, diploid genomes (ER+ Typical-risk-enriched), genome-wide instability (TNBC-enriched), and focal, complex amplifications (ER+ High / HER2+-enriched).

[0073] Metastatic lesions exhibit elevated SNV and SV burdens compared to unpaired primary tumors, likely due to therapy. Using the above approach, we identified six de novo SV signatures in metastases which correlated with those in primary tumors (Fig. 11A to 11 B) and showed similar subgroup-specific enrichment patterns (Fig. 10F). Two- dimensional projection again revealed three dominant archetypes (Fig. 11 C) that overlap with those in primary tumors (Figs. 6B to 6C, Fig. 10G; Fig. 11 D). Our results were robust to choice of dimensionality reduction algorithm (Figs. 11 E to 11 F). Thus, the three genomic archetypes of breast cancer are conserved in metastatic disease.

[0074] SV signatures were generally conserved, though elevated, in metastatic tumors except for RS4 and RS6 in ER+ High-risk and HER2+ tumors, respectively, which were stable (Fig. 10H). These data support the early occurrence of complex rearrangements and their persistence through metastasis. Although the distribution of CN signatures mirrored primary tumors, the Pareto front revealed increased alteration burden and more intermixed profiles in metastasis, consistent with increased WGD and genomic instability( Fig. 10l, Figs. 11 H to 111). Thus, metastatic tumors retain the scars of subgroup-specific mutational processes operative in early-stage disease.

[0075] While ER+ Typical-risk tumors have favorable prognosis, 29% of patients experience distant relapse. We investigated whether the genomic archetypes improve risk stratification. Mapping METABRIC onto the Pareto front (Methods; Fig. 10J; Figs. 11 J to 11 L), the position of ER+ Typical-risk tumors was predictive of relapse, with recurrent tumors mapping closer to the ER+ High / HER2+ archetype (Figs. 10K to 10L) accompanied by higher HRD LOH score, invasive lobular carcinoma (ILC) histology and increased proliferation.

[0076] In METABRIC, ILCs were enriched in ER+ Typical-risk tumors (OR = 2.20, P = 2.27x1 O’3, Fisher’s exact test, Fig. 11 M). Within ER+ High-risk tumors, ILCs exhibited a higher 5-year recurrence risk (39% vs. 30%) and cumulative recurrence risk (62% vs. 54% at 20 years, Fig. 10M). This difference was more striking amongst ER+ Typical-risk tumors (55% vs. 37% at 20 years). ILCs were closer to the ER+ Typical-risk archetype than IDC counterparts (P = 2.10x10-5, Figs. 10N to 10O) given their lower levels of WGD, ploidy, and FGA. Thus, given comparable genomic architectures, lobular histology remains a high-risk feature.ER-Induced R-Loops Fuel ecDNA Genesis

[0077] ER+ High-risk and HER2+ breast tumors were enriched for complex amplifications in two independent cohorts (OR>10.1 ; P<2.2x10-16 Fig. 12A; Fig. 13A) motivating further exploration of their origin and nature (Figs. 14A to 14B). There was no difference in cyclic amplifications in HER2+ / ER- primary tumors compared to HER2+ / ER+ primary tumors (Fig. 13B). Leveraging two independent ecDNA inference methods, 43- 67% of primary ER+ High-risk and HER2+ cases were predicted to harbor ecDNA (Figs. 13C to 13D). A proportion of HER2+ primary tumors (25.7%) harbored amplifications in loci specific to the ER+ High-risk subgroup (Fig. 13E, Fig. 14C), with 8.57% predicted to be on ecDNA. Additionally, we observed a modest enrichment of inversions at 11q13 locus in primary tumors. HRD and ecDNA were mutually exclusive in primary ER+ High- risk and IC10 tumors (OR = 0.21-0.29; FDR<0.02; Fig. 14D). We interrogated complex amplifications in 406 pre-invasive DCIS profiled via shallow WGS (5X median coverage).We predicted 35 cyclic and 205 complex non-cyclic amplifications, enriched in ER+ High- risk / HER2+ tumors (OR = 4.21 ; P = 2.48x10-4; Fig. 13F; Fig. 14E). This pattern replicated in 12 DCIS samples from GEL (92.8X). Leveraging the clock-like accumulation of mutations, SNV density informs the timing of cyclic amplifications (Methods). Compared to IC10 tumors, cyclic amplifications in ER+ High-risk and HER2+ tumors had lower SNV density prior to amplification, suggesting an earlier origin (Fig. 12B; Fig. 14F). Median time of cyclic amplification in ER+ High-risk and HER2+ tumors occurs decades earlier than in IC10 tumors, respectively, implicating cyclic amplifications as early events.

[0078] Most cyclic amplifications in ER+ High-risk (88%) and HER2+ (96%) tumors overlapped at least one COSMIC-defined oncogene (Fig. 13G). Of these, 79-92% involved oncogenes in IC-associated cytobands (Fig. 5J) and 15% involved two or more cytobands (Fig. 13H). In cell line models of IC2 (UCD65) and IC6 (UCD12) before and after linear DNA digestion, significantly higher sequencing coverage occurred at regions predicted to encode ecDNA, corroborating our computational predictions (Figs. 12C to 12D; Figs. 14G to 14H). Oncogene incorporation varied across subtypes, with HER2+ tumors harboring the largest number per megabase (Figs. 131 to 13J). 82% of IC2, 59% of IC5 (HER2+), 48% of IC6 and 32.5% of IC1 tumors had predicted cyclic amplifications at subgroup-defining cytobands, while 3% of IC1 / IC9 tumors harbored cyclic amplifications at 20q13, spanning the NCOA3 oncogene (Fig. 12E). Overall 42% of IC9 tumors harbor ecDNA, but these ecDNA are diffuse along the genome and don’t include MYC. In support, focal SV peaks were not observed at 8q24 spanning the MYC oncogene in IC9 primary or metastatic tumors. Rather, a broader region is subject to enhancer hijacking by the long-non-coding RNA, PVT1. PVT1 co-amplifies with MYC in ~90% of tumors (Fig. 141). Frequent enhancer hijacking at MYC may explain the weak correlation between MYC CN and mRNA abundance (Figs. 14J to 14K).

[0079] The subset of ER+ Typical risk tumors harboring ecDNA fell along the ER+ Typical vs. High-risk archetype continuum (Figs. 13Kto 13L). In contrast, ER- tumors with ecDNA had limited structural conservation (Fig. 15A). Across all subgroups, similar patterns were observed in metastatic and pre-invasive tumors (Figs. 15B to 15F).

[0080] Elevated replication stress has been associated with response to checkpoint and DNA-repair inhibitors and hence is a therapeutic vulnerability in TNBC. Assessingreplication stress across the IC subgroups, we found elevated levels of oncogene-induced replication stress in ER+ High-risk and HER2+ tumors compared to ER Typical-risk, IC10 and IC4ER- tumors (FDR<0.026; Figs. 16A to 16B). The replication stress signature was positively correlated with TNBC-enriched and ER+ High / HER2+-enriched genomic archetypes (ES>0.154, P<4.98x10-15; Figs. 16C to 16E). Within ER+ Typical-risk, ILC had higher replication stress than IDC (FDR=4.08x10-3). Meta-analysis suggests a positive association between ecDNA and replication stress in HER2+, IC1 and IC6 tumors (Fig. 16F) and higher levels of Type-I IFN signature in ecDNA+ tumors (Fig. 17A). Finally, ER+ High-risk and HER2+ tumors demonstrated elevated cGAS / STING activity (Figs. 16G to 16H), a possible therapeutic target linked to chromosomal instability and replication stress.

[0081] We found that cyclic amplifications were significantly enriched for translocations compared to complex non-cyclic amplifications in ER+ primary tumors (Fig. 12F; Figs. 18A to 18B; Fig. 17B). These cyclic-amplified ER+ High-risk tumors had higher ESR1 mRNA abundance ([3 = 1.27; P = 6.90x10’3; Fig. 18C) and enriched ER binding within the amplified region (Fig. 18D; Fig. 17C). Nonetheless, ER signaling was lower in ER+ High-risk compared to Typical-risk tumors (Fig. 5H). Given the evidence for ecDNA in pre-malignant lesions, we hypothesized that ER signaling is elevated in DCIS lesions that classify as ER+ High-risk and subsequently decreases in invasive disease. Leveraging 18 paired ER+ DCIS and primary tumors with transcriptome sequencing, we observed decreased ER signaling in ER+ High-risk tumors (ES = 0.33; P = 0.03; Fig. 18E). These data support the role of ER in ecDNA genesis via translocations and emphasize their early origin.

[0082] The mechanism by which ER activation induces translocations remains unknown. ER recruitment of APOBEC3B (A3B) promotes double stranded breaks (DSB) at ER binding sites (Fig. 12F; Fig. 18A). Increased ER-induced transcription leads to the formation of R-loops producing ssDNA, a substrate for A3B editing. A3B deaminates cytosine to uracil which can be repaired by base excision repair (BER). Single strand nicks induced by BER coupled with transcription-coupled nucleotide excision repair processing of the R-loop can result in DSBs. Taken together, A3B can exacerbate chromosomal instability in the pre-invasive setting. We hypothesized that ER-induced R-loops initiate translocation-bridge amplifications via A3B-editing and confirmed that A3B binding in ER+ cell lines was enriched in cyclic vs. non-cyclic amplifications (Fig. 18D; Fig. 17C). Treatment with estradiol (E2) in MCF7 cell lines induced R-loops (nR-loops = 212) in the same regions where cyclic amplifications were observed in patient tumors (Fig. 12G; Figs. 18F to 18G). This finding was specific to ER-induced R-loops (nR-loops = 13,965; Fig. 18H). Unresolved R-loops due to A3B knockout in MCF10A cells were preferentially enriched in regions of cyclic amplifications in primary breast tumors (Fig. 181; Fig. 17D). Tumors containing ecDNA were also enriched for transcription-replication collision-associated large tandem duplications (>100kbp), indicative of impaired R-loop resolution (Fig. 17E). Translocations within cyclic amplifications were significantly closer to ER-induced R-loops than those outside cyclic amplifications (Figs. 18J to 18K). These data support a role for A3B in R-loop resolution, contributing to ecDNA formation.

[0083] Accordingly, we hypothesized that estrogen-induced SV breakpoints would be enriched at ER-induced R-loops. Comparing SV patterns in E2-treated MCF7 cells via High-Throughput Genome-Wide Translocation Sequencing (HTGTS) of double stranded breaks forming translocations induced by CRISPR-Cas9, we confirmed the enrichment for E2-induced breakpoints (Fig. 12H) compared to all R-loops (Fig. 181). There was no difference in replication timing between cyclic and non-cyclic amplifications (Fig. 17F). ER-induced R-loops were enriched closer to IC-specific oncogenes PAK1 (IC2), ZNF703 (IC6), and MYC (IC9) than all other COSMIC-defined oncogenes, including ERBB2 (Fig. 121). There was no enrichment of ER-induced R-loops near IC1 oncogenes.

[0084] Germline copy number polymorphisms in A3B have been associated with APOBEC-dependent mutations and immune activation in breast cancer. Despite limited power, we found a modest but non-significant decrease in ecDNA prevalence in ER+ High and Typical-risk samples with the homozygous deletion allele (n=5; Fig. 18M). Taken together, these data indicate that ER activity promotes cyclic amplifications via R-loop formation and A3B-editing.The ICs Harbor Distinct TMEs

[0085] Tumor clonal composition and genomic features are sculpted by immune pressures, while oncogenic alterations promote both pro-tumor and anti-tumor immuneresponses. Using transcriptom ic profiles, we characterized the TME in primary (nTCGA = 1 ,015; nMETABRIC = 1 ,894) and metastatic (n = 360) tumors focusing on four subtypes defined by immune infiltration and stromal composition: immune-enriched fibrotic (IE / F); immune-enriched non-fibrotic (IE); fibrotic (F); and depleted (D) (Fig. 19A; Fig. 20A). The reproducibility of the TME subtypes is supported by single-cell spatial proteomic profiling (n = 384; Fig. 20B) and cell type proportions estimated from bulk transcriptom ics (Fig. 21A).

[0086] We then quantified microenvironmental differences across the IC subgroups. Primary IC10 and IC4ER- were enriched for immune rich (IE and I E / F) TMEs (OR = 3.004, P = 5.17x10'11, Fisher’s exact test; Fig. 19B; Fig. 21 B). ER+ High-risk and HER2+ primary tumors harbored immune-depleted TMEs (OR = 3.09, P = 1 .06x10-15, Fisher’s exact test), while genome stable ER+ Typical-risk and IC4ER- primary tumors were enriched for fibrotic signatures (subtypes F and IE / F; OR = 5.619, P < 2.2x10-16, Fisher’s exact test). These observations were replicated using a second transcriptional immune score (Figs. 21 C to 21 D). Within ER+ High-risk tumors, immune enrichment did not differ across subgroups (Fig. 20C). Amongst ER+ Typical tumors, ILCs were enriched for the IE / F subtype compared to IDCs (OR=2.18, P=1.17x10-3; Fig. 20D).

[0087] IC4ER- tumors have a more favorable prognosis but longer-term risk of recurrence than IC10 tumors despite similar genomic landscapes (Fig. 6A). To investigate differences in their TME, we leveraged single cell spatial proteomic data and discovered an increased proportion of fibroblasts and T cells in IC4ER- compared to IC10 (Fig. 20E; Fig. 21 E). In support, previous work has linked increased T-cell infiltration with improved overall survival in TNBC. Compared to primary tumors, ER- metastatic tumors were depleted of IE and IE / F features, (OR = 3.01 ; P = 2x10-4) (Fig. 20F). In contrast, HER2+ and ER+ tumors exhibited stable TMEs through metastasis (Fig. 19B; Fig. 20G; Fig. 21 F), consistent with prior reports that ER promotes immune-suppression and immunoediting in pre-invasive lesions.

[0088] We found that 43.86% of primary and 47.67% metastatic tumors exhibited genetic immune escape (GIE), most of which occurred in a single pathway with varying prevalence across IC subgroups (Fig. 19C). IC2 and IC6 tumors were more immune- depleted than IC1 and IC9 (Fig. 20C) but harbored fewer GIE (Fig. 22A). Instead, 60% ofprimary IC6 tumors amplified indoleamine 2,3-dioxygenase (IDO1 ), a heme-containing enzyme located within 8p11 .21 that metabolizes tryptophan involved in immune tolerance (Fig. 22B to 22C). ER+ Typical-risk ILCs exhibit fewer GIE alterations than ER+ Typicalrisk IDC tumors (Fig. 22D) and GIE was not associated with antigen burden (Fig. 21 G).

[0089] Complex alterations and SVs have been overlooked when evaluating GIE. We found that -20% of primary and metastatic tumors with GIE harbored SVs or complex amplifications (Fig. 19D; Fig. 22E). HER2+ tumors demonstrated the largest increase in GIE between primary and metastatic disease, potentially due to greater pressure to evade anti-HER2 therapies (OR = 2.23, FDR = 0.19, Fisher’s exact test; Fig. 22F). These data illuminate the role of complex alterations in immune escape and tumor-immune coevolution during disease progression.Discussion

[0090] Here, we identify three dominant genomic archetypes of breast cancer driven by distinct mutational processes, describing a continuum of genomic profiles and providing a mechanistic basis for these patterns (Fig. 23A). These three archetypes overlap with the major clinical breast cancer subgroups with a notable difference. For a sizable proportion of ER+ tumors (43.2%), the ER+ High-risk / HER2+ archetype dominates and the mutational processes are indistinguishable from HER2+ tumors. Rather than amplifying ERBB2 / HER2, these ER+ High-risk tumors harbor focal amplifications of other oncogenes (Fig. 5J) and have an elevated risk of recurrence akin to HER2+ tumors prior to the introduction of anti-HER2 therapies. These ER+ High-risk tumors may similarly benefit from agents directed at their amplified oncogenic drivers and / or shared vulnerabilities.

[0091] A defining feature of the ER+ High / HER2+ archetype is the generation of focally amplified ecDNA via ER-induced R-loops and A3B editing. ER-induced R-loops create ssDNA, which serves as a substrate for A3B editing. DSBs arising from BER and NER are resolved in the form of inter-chromosomal translocations. Dicentric chromosomes can form chromosome bridges during mitosis and breakage of these bridges can generate ecDNA 30. ecDNA formation preferentially occurs at loci that define the four ER+ High- risk subgroups and HER2+ disease. While ecDNA genesis depends on ER, circularamplification may reduce reliance on ER by increasing a particular oncogene’s copy number and rewiring its regulatory network. This is supported by reduced ER signaling in ER+ High-risk tumors from DCIS to invasive disease. Since ER can cause doublestranded breaks, ecDNA formation may balance increased oncogenic signaling with protection against further ER-induced genomic instability (Fig. 23B), and hence reflects an evolutionary trade off, consistent with mutual exclusivity between complex amplifications and diffuse genome instability.

[0092] Beyond tumor subtype, the mutational processes captured by our architectural map may be indicative of distinct therapeutic vulnerabilities. For example, HRD-like tumors are sensitive to PARP inhibition and this has become a mainstay of therapy for TNBC. We find that 44% of TNBC have HRD-like profiles based on WGS, while 13% of ER+ High-risk tumors exhibit BRCA2-like patterns. While HRD as measured from sequencing data is not confirmed to correlate with PARP inhibitor sensitivity, this result implies that additional patients may benefit from these agents. Further, we find that focal ly amplified ER+ High-risk tumors exhibit elevated replication stress pathway activities, suggesting potential sensitivity to novel agents targeting this pathway. Additionally, while APOBEC3 mutagenesis can occur early during tumorigenesis, given its impact on ER activity, A3B represents a potential target in ER+ High-risk subgroup for which inhibitors are in development.

[0093] The mutational processes that generate and propagate genomic instability both sculpt oncogenic signaling and mediate interactions between tumor cells and the TME. More specifically, SVs contribute to GIE in 9% of breast tumors, but have been overlooked, owing to the need for WGS. Basal-like IC10 tumors, which harbor both high genomic instability and immune infiltrates, likely adapt to this immune pressure via GIE. In contrast, ER+ tumors, both Typical and High-risk, are more immune-depleted at the onset with fewer GIE events, suggesting non-GIE mechanisms. This is noteworthy given the evolving utility of immunotherapy in breast cancer 46. Despite high immune infiltration, up to 62% of TNBC are resistant to current immunotherapies, potentially due to GIE, while 38% of ER+ tumors have immune-enriched TMEs, making them candidates for such agents. Our findings illuminate multiple potential strategies for personalizing breast cancer treatment, which will be the focus of ongoing pre-clinical and translational studies.MethodsData Collection

[0094] We uniformly processed whole-genome sequenced tumor and normal pairs for 1 ,422 breast tumors encompassing primary invasive disease (TCGA / ICGC / Hartwig; n=702) and metastatic lesions (Hartwig / TCGA; n=720). Of these, 1 ,028 TCGA, 189 ICGC and 469 Hartwig primary and metastatic lesions additionally had transcriptom ic profiling. Further, we processed 406 ductal carcinoma in situ lesions profiled with shallow whole genome sequencing (HTAN; median coverage = 5 reads); leveraged 1 ,894 primary tumors (METABRIC) with array-based genomic and transcriptomic profiles along with 20 years of clinical follow-up; and 3042 tumors from Genomics England (GEL) (v18) with whole genome sequencing data, which were processed within the Genomics England Research Environment. Finally, we considered whole exome sequencing data from 54 breast cancer patients where both primary and metastatic tumors were profiled (nMetastaec = 82).Data Analysis Using Isabl Platform and Containerized Applications

[0095] To ensure the computational reproducibility, harmoni.zing and processing the genomic data from these large cohorts at scale, we have deployed Isabl platform56 locally and developed containerized and version control applications. We integrated tools for identification of Single Nucleotide Variants (SNVs), Copy Number Alterations (CNAs), Structural Variants (SVs), Extrachromosomal DNA (ecDNA), HLA typing, Neoantigen prediction, IC-subtype prediction and RNA quantification that are described below.DNA Sequencing Data

[0096] Alignment Starting from FASTQ files, we aligned whole genome sequenced tumors (primary nTcGA / / ccc = 639, nHartwi9= 63; metastatic nrcGMCGc = 2, nHartwi9= 718; paired primary and metastatic r?=136) to hs37d5 using BWA-MEM as implemented in the TCGA-ICGC-PanCancer PCAP-core docker.Quality Controls

[0097] Whole genome sequencing. We assessed data quality and excluded samples if median tumor coverage was < 25 or median normal coverage < 10.Copy Number Analysis

[0098] Whole genome sequencing. We detected copy number aberrations in all primary and metastatic WGS samples with allele-specific CNA caller. First, we generated snp pileup files using snp-pileupwrapper.R using dbSNP (v138) and default parameters. Next, we ran run-facets-wrapper.R testing two combinations of the purity-cval and eval parameters: -purity-cval 2000 -eval 1000 or -purity-cval 1000 -eval 500 with three different -normal-depth {25,30,35} for a total of six runs. The optimal parameter selection for each sample was determined by an aggregated score based on six possible extremes: extreme diplogR defined as abs(dipLogR) > 0.8, number of whole chromosome losses, number of whole chromosome loss of heterozygosity, extreme fraction of homozygous deletions defined as > 0.8, extreme loss of heterozygosity defined as > 0.5 of the genome, extreme ploidy defined as >6 or <1.8. We selected the parameter choice that minimized these extreme cases.

[0099] Whole exome sequencing (TCGA). We detected copy number aberrations in all primary TCGA samples with allele-specific CNA caller. First, we generated snp pileup files using snp-pileup-wrapper.R using dbSNP (v146) and default parameters. Next, we ran run-facets-wrapper.R using the following parameters: -purity-cval 600 -eval 300 - normal-depth 25.

[0100] Whole exome sequencing (paired primary and metastatic samples): Copy number aberrations were detected with allele-specific CNA caller testing two combinations of the purity-cval and eval parameters: -purity-cval 1000 -eval 500 or - purity-cval 500 -eval 250 with three different normal-depth {25,30,35} for a total of six runs.

[0101] Shallow whole genome sequencing (DCIS). We detected copy number aberrations using QDNASeq with bin size set to 50 kbp. First, we calculated genomewide read counts in 50kbp bins with a minimum mapping quality of 37. Next, we filtered reads based on overlapping blacklisted regions and mappability using the applyFiltersfunction. Third, we corrected read counts for GC content using the estimatecorrection and correctBins functions. We normalized outliers and smoothed using the normalizeBins and smoothOutlierBins functions before inferring copy number aberrations using the segmentBins and normalizeSegmentedBins functions. Finally, we estimated purity and ploidy using ACE and copy number segments recalled, correcting for estimated purity using the callBins function and cellularity = estimated cellularity.Structural Variant Analysis

[0102] Structural variant detection. We detected structural variants in both primary (n=702) and metastatic (n=720) tumors using a consensus approach leveraging four independent callers: MANTA (v1.6.0), DELLY (vO.9.1 ), SvABA (v1.1.3) and GRIDSS (v2.13.2). We ran MANTA on tumor / normal pairs using default settings and DELLY on tumor / normal pairs excluding unmappable regions annotated in the DELLY blocklist. We ran SvABA on tumor / normal pairs using default parameters and excluding ENCODE blocklist regions (--blacklist hg19-blacklist-nochr.v2.bed) and annotated with dbSNP (v138) (-dbsnp-vcf dbsnp_138.b37.vcf). Fourth, we ran GRIDSS on tumor / normal pairs using the following parameters: --blacklist hg19-blacklist-nochr.v2.bed picardoptions VALIDATION_STRINGENCY=LENIENT. We filtered somatic SVs called by GRIDSS with the gridss_somatic_filter function using default parameters. Finally, we retained SV calls if they were identified by GRIDSS and at least one other SV caller, i.e. MANTA, DELLY or SvABA.

[0103] Damaging structural variants. We detected damaging structural variants by overlapping the genomic location of breakpoints to gene sequences. Briefly, we used the GENCODE Release 3967 reference to map the breakpoints using GenomicRanges R package (v1 .50.2) and assigned the consequences of each event on the long non-coding RNA (IncRNA) and protein-coding gene sequences to remove structural variants that lead to harmless changes for the sequence of the gene such as large in-frame deletions in a single intron. Single nucleotide variant (SNV) detection Similar to SV detection, we detected SNVs in both primary (n=702) and metastatic (n=720) tumors using a consensus approach leveraging two independent callers: Mutect2 (v4.1.7.0) and Strelka2 (v2.9.10). We ran Mutect2 on tumor / normal pairs as part of the nf-core / sarek Nextflow (v20.12.0)pipeline with the following parameters: -r 2.7.1 -step variantcalling -skip_qc all -tools Mutect2 -profile singularity. We ran Strelka2 on tumor / normal pairs with default parameters and leveraging indels previously called by MANTA (see Structural variant detection) using the parameters -indelCandidates. A consensus variant call set was obtained by combining Mutect2 and Strelka2 variants using ‘VariantFilter’.

[0104] Inferred genetic ancestry. In Hartwig, we genotyped 128 ancestry informative markers to produce gvcfs, followed by joint genotyping and variant recalibration. Using the 128 ancestry informative markers, we estimated the genetic ancestry of each case via principal component comparison with reference populations.

[0105] Germline variants in homologous repair genes. We determined germline variants in 14 homologous repair genes (BRCA1 , BRCA2, ATM, PALB2, CHEK2, CHEK1 , BARD1 , BRIP1 , FANCL, RAD54L, RAD51 B, RAD51 C, RAD51 D, CDK12). We merged gVCFs followed by joint genotyping and variant recalibration with -truth-sensitivity-filter- level 99. Finally, we determined pathogenic and likely pathogenic variants.

[0106] Prediction of homologous recombination deficiency (HRD). We downloaded HRD labels for the PCAWG and Nik-Zainal (BRCA-EU) cohorts, as predicted by a pan-cancer classifier that can discriminate between BRCA1-like and BRCA2-like HRD. The pan-cancer classifier outputs the probability of BRCA1-like or BRCA2-like HRD, the overall probability of HRD is given by the sum of these two probabilities and samples are classified as HRD if this sum is > 0.5. For samples called as HRD, we assigned BRCA1 -like or BRCA2-like labels based on the larger probability. Samples where HRD status could not be determined were excluded from analyses, resulting in HRD calls from 571 samples in our primary cohort and 562 samples in our metastatic cohort.

[0107] Extrachromosomal DNA (ecDNA) detection. First, we detected candidate amplifications with the following parameters: -diagram -scatter -annotate refFlat.txt - access access-hs37d5-cnvkit.bed -method wgs. We used amplifications (copy number > 4) to seed by first trimming segments using seed_trimmer_v2.py within the PrepareAA module setting -minsize 50000. Next, we ran PrepareAA.py on the trimmed seeds using the following parameters: -cngain 4.0 -cnsize_min 50000 -downsample 10. Finally, we ran AmpliconClassifier on the resulting output with parameter -min_size 5000. Wevisualized predicted ecDNA structures. As quality control, we compared the amplicon classifications for a subset of primary samples to previously reported classifications (Fig. 14A). Further, we calculated the proportion of complex amplicons that overlapped at least one predicted SV cluster based on our consensus SV calls (Fig. 14B). We defined SV clusters as regions of the genome harboring SV with breakpoints significantly closer to one another than expected by chance alone. Finally, to assess the accuracy of ecDNA calling in shallow whole genome sequencing, 57 primary breast cancer tumors from TCGA were downsampled to match the sequencing depth of the DCIS samples (~5X). Seed amplifications were called using the same protocol as the DCIS samples (see Copy Number Analysis: Shallow whole genome sequencing (DCIS)). Complex amplifications were called and the same protocol outlined above. Using down-sampled WGS, we estimate a false positive rate of 3.5% and false negative rate of 22.8% in ecDNA detection from sWGS (Fig. 14E) indicating 8.6% as a lower bound estimate of ecDNA prevalence in DCIS.Classification of SV Events using Junction Balance Analysis (JaBbA)

[0108] We first calculated the genome wide coverage for all tumors and matched normal samples at 1 kbp resolution. We used the consensus SV calls (see Structural variant detection) as the junction input along with purity and ploidy estimates. Finally, we ran JaBbA with default parameters and predicted classes and visualized rearrangements. In our cohort, JaBba predicted 13 distinct SV events: deletion (del), duplication (dup), inversion (inv), templated insertion chains (tic), inverted duplications (invdup), double minutes (dm), complex double minutes (cpxdm), breakage fusion bridge cycles (bfb), chromoplexy, chromothripsis, tyfonas, pygro and rigma. On a per sample basis, we found 77% concordance between samples that were predicted to harbor at least one ecDNA according to JaBbA.Signature Analysis

[0109] Copy number and structural variant signatures. We inferred the 24 COSMIC copy number signatures for both primary (n=702) and metastatic (n=720) tumors separately. We also inferred the 24 COSMIC copy number signatures in theMETABRIC cohort and the same parameters as the primary discovery cohort. We discovered de novo SV signatures from ensemble SV calls using DELLY, MANTA, SVABA and GRIDSS (see Structural variant detection). We generated input data from individual level SV calls in bedpe format and discovered signatures with default parameters on primary and metastatic tumors separately. We selected the optimal signature number by maximizing stability (Silhouette score) and minimizing distance (Cosine distance) between signatures inferred from 100 NMF iterations and from reproducibility in the metastatic cohort. We assessed the level of correlation of highly cooccurent SV signatures RS6 and RS4 and found that RS6 demonstrated stronger associations with complex amplifications than RS4 (Figs. 8F to 8G). Finally, we integrated the copy number and SV signatures with principal component analysis (PCA) using the prcomp implementation in R and visualized across the top two principal components (PC1 and PC2).

[0110] Single nucleotide variant signatures. We performed mutational signature analysis for SNVs. We used known mutational signatures of breast cancer as input signatures (SBS1 , 2, 5, 8, 9,10a, 10b, 13, 15, 17a, 18, 19, 21 , 26, 29, 30, 34, 37, 38, 39, 40, 41 , and 44) and generated mutational profiles of 96 contexts with the mut_matrix function. We leveraged the fit_to_signature function to fit known mutational signatures and calculated the cosine similarity by comparing the original mutation matrix with the reconstructed matrix using the cos_sim function.Pareto Front Projection

[0111] We first used an unsupervised approach without any assumption and performed a PCA to visualize the distribution of the architectural profiles. We observed that this distribution represents a continuum rather than distinct clusters over the PC1 - PC2 space (Fig. 9A). Additionally, we noticed this very particular polyhedron shape suggesting that the Pareto front theory might be applied with three archetypes at play. The Pareto optimum theory has demonstrated that organisms that perform multiple tasks evolve phenotypes along a continuum in the shape of a polyhedron where the vertices represent organisms highly specialized in one specific task. Akin to this concept, the triangular PCA projection of breast cancer tumors based on their rearrangementsignatures suggested that breast tumors may exist along a mutational continuum with vertices defined by the dominance of one distinct mutational process. To address this question, we projected the distribution of (i) primary-only; (ii) metastatic-only, and (iii) primary-metastatic samples within the above-mentioned genomic PCA individually on a Pareto front with the assumption of three different genomic archetypes.

[0112] In more details, the algorithm tries to find polyhedra by testing successively 1 to n axes, adding them one after another in decreasing order of variance explained. For each number n of axes used, ParetoTI identifies the position of the n + 1 = k vertices (archetypes) in the molecular map defined, and we used 200 bootstraps, each taking 75% of the data to measure the variability in archetype position and infer archetype positions robust to outliers (function fit_pch_bootstrap with the parameters bootstrap = T and bootstrap_N = 200). In addition, the ParetoTI framework calculates the likelihood of observing data shaped as polyhedron given no relationship between variables by disrupting the relationships between them and finding a polytope that best describes the data (function random ise_fit_pch with the parameters bootstrap = T and bootstrap_N = 200). The ParetoTI framework uses two main metrics to measure the uncertainty of the Pareto front model fit: the explained variance and the t-ratio. The explained variance represents the proportion of variation explained by the model, ranging from 0 (all points far outside the polyhedron) to 1 (all points inside the polyhedron). The t-ratio represents the volume ratio between the tested polyhedron to the convex hull of data. Of note, the explained variance increases monotonically with k, therefore the optimum value corresponds to (an ill-defined) elbow, while maximizing the t-ratio.

[0113] For replication purposes, we also projected the METABRIC primary samples onto a Pareto front by first learning the embeddings of the CNA signatures onto the archetypes in the absence of SV signatures in the discovery cohort. To this end, leveraging the primary discovery cohort (n=702), we built two Random Forest models (mtry = 4) using the randomForest function as implemented in the randomForest R package (v4.6-14; seed = 1234) to predict genomic PC1 and PC2 from the CNA signatures alone (PC ~ CN signature proportions and PC2 ~ CN signature proportions). The Random Forest predictions based on CNA alone captured the same three archetypes as the SV and CNA signatures together. The Random Forest models were then used topredict the respective genomic PC1 and PC2 in the METABRIC cohort from CNA signatures alone using the predict function. The predicted METABRIC PCs were then projected to the Pareto front using the protocol described above. The same procedure was repeated on the Genomics England cohort.

[0114] We reported the robustness statistics of all the Pareto fits performed on the four above-mentioned PCAs in Fig. 9E.Nonnegative Matrix Factorization

[0115] In order to assess the robustness of the main findings captured by the Pareto front analysis described above, we chose to use a Nonnegative Matrix Factorization (NMF) model on the three signature matrices (i) Primary-only, (ii) Metastatic-only, and (iii) Primary+Metastatic, separately. We implemented a multiple run method with nrun= 200 to select for the best NMF fit. We found three ranks that correlate strongly with the described archetypes (Figs.11 E to 11 G).Replication Stress Signatures

[0116] We scored and evaluated the signature using the Iog2 transformed TCGA expression data and Iog2 intensities from the METABRIC expression array. For each cohort, we considered the six genes that make up the signature (NAT10, DDX27, ZNF48, C8ORF33, MOCS3, and MPP6). We calculated the mean expression value for each sample, and linearly transformed these values to the range of 0 and 1 within each cohort. We compared the signature scored in this study with the original publication for the 817 TCGA samples in both sets (Spearman’s p = 0.82). To compare the replication stress signature by ecDNA status, we used ecDNA calls as detailed above (see Extrachromosomal DNA (ecDNA) detection). We considered tumors ecDNApositive when one or more ecDNA events were observed. We performed logistic regression, correcting for cohort and amplicon copy number for the amplicon-driven ER+ High and HER2+ subgroups. Specifically, the copy number of ERBB2 was considered for HER2+ / IC5; PPM1 D and MSI2 for IC1 ; FGFR1 for IC6; and MYC for IC9.cGAS-STING, ESR1, and Type-IFN Signatures

[0117] Due to the link between replication stress and cGAS-STING, we calculated the enrichment of a cGAS-STING signature for TCGA and METABRIC expression data. For TCGA, we first applied the variance stabilizing transformation (VST), while for METABRIC we used the Iog2 intensities from the expression array. We calculated a cGAS-STING signature based on seven genes (C6orf150, CCL5, CXCL10, IRF3, TBK1, TMEM173, and STAT1). TMEM173 was not found in the METABRIC expression data, so the signature was scored for the remaining 6 genes.

[0118] For both signatures, we used the gsva to calculate the enrichment for each sample.Ensemble Integrative Clustering (ENiClust) Classifier

[0119] Reference preparation. With both RNA and DNA-sequencing data (iC10 DNA+RNA) to predict the reference IC subtypes in both the TCGA-WES and the Hartwig datasets. Briefly, the iC10 method leverages copy number and / or expression profiles to predict IC subtypes. Predictions are based on a nearest shrunken centroid classification, estimating the closest class-based centroid from the training dataset to assign the testing set class. Here, we downloaded the normalized expression data for TCGA and reprocessed the Hartwig RNA-sequencing data. For DNA-sequencing data, we extracted total copy number (CN) (see Feature preparation section below) and the copy number log ratio (median_cnlr_seg column) per gene from the gene-level FACETS outputs.

[0120] Feature preparation. ENiClust classification relies on a set of 624 features, including 611 gene-level CN, four ICspecific cytoband CN levels, five genome-level CN features, expression status for ER and HER2 / ERBB2, and the sequencing technology as additional features, as described below.

[0121] Copy number. The list of 639 selected genes for gene-level CN features was generated as (i) being part of the features for n=614 genes or (ii) manual curation of genes residing in amplified cytobands in IC1 , IC2, IC6, IC9, and IC8 (n=25 genes). In total, the model has been trained on 611 out of 639 genes since 27 genes were absent from the TCGA-WES and Hartwig gene-level CN datasets and LY6L was mostly missing in the TCGA-WES dataset (missing CN value in 1001 / 1013 samples). The selected genes,annotations, and corresponding aliases present in the FACETS gene annotation files (GENCODE V29 for both hg19 and hg38). Of note, in TCGA-WES, there were 144 / 611 genes for which the CN level was missing in at least one sample, 63 / 611 in Hartwig, and 11 / 611 in METABRIC. We ensured values were available in at least 1 / 4 of samples in each dataset.

[0122] For Hartwig and TCGA-WES, we extracted the copy number aberrations. We ran the Hartwig WGS tumors with the following parameters: -puritycval 1800 -eval 1000 -normal-depth 20. For the majority of TCGA-WES tumors we used the following parameters: -purity-cval 600 -eval 300 -normal-depth 25 and we manually refitted eighty- three tumors after visual inspection. Of note, we removed two TCGA-WES samples (TCGA-AR-A0TQ and TCGA-C8-A135) estimated as outliers in their CN profile, especially for their ploidy. For METABRIC, we extracted the copy number data. To ensure consistency between datasets, we obtained the gene-level CN profile of the samples from the published ASCAT segment files with a modified version of the corresponding FACETS function (gene_level_changes).

[0123] From these FACETS (TCGA-WES and Hartwig) or ASCAT (METABRIC) CN calls, we extracted the gene-level copy number profile of the 611 selected genes. Then, we calculated the profiles of the four cytoband-level CN (17q24, 11 q13, 8p11 , and 8q24) by aggregating the transformed CN value of every gene falling in each of the four IC- specific chromosomal cytobands (see Transformation and model development - Feature transformation section). We also extracted the five genome-level features (ploidy, wholegenome doubling (WGD), average breakpoints per chromosome, fraction of CNA, HRD LOH score) from FACETS except for the “average breakpoints per chromosome”. This last feature was calculated from the average CN breakpoint burden for each mega bp using the sizes of chromosomes. To ensure consistency between datasets, we generated the genome-level features in METABRIC by calling the corresponding FACETS functions (calculate_fraction_cna and calculate_hrdloh) on the segment-level data. Importantly, a bug in the FACETS HRD LOH calculation incorrectly returns a score of 0 for a subset of samples, so all samples with a score of 0 were manually re-run with a modified version of the function (calculate_hrdloh_fixed) and we corrected the corresponding values for 9 samples in TCGA-WES, 16 samples in Hartwig, and 2 samples in METABRIC.

[0124] Receptor status. We utilize ER and ERBB2 / HER2 receptor status of METABRIC and Hartwig samples. For TCGA, we report the inferred ER and HER2 status from RNA-seq, where ER and HER2 status were classified based on their empirical expression distributions using a 2-component Gaussian mixture model. Of note, both ER and HER2 status were unknown for 21 / 297 in Hartwig.

[0125] Additional features. Sequencing technology, namely WES or WGS, was included as a feature.Transformation and Model Development

[0126] Feature transformation. Missing values were set to 0 before further transformation. To minimize the bias that a WGD event could generate at the gene- and cytoband-level CN values, we first doubled the gene-level CN profile of non-WGD samples. Then, to avoid overfitting in the use of gene-level CN profile, we transformed the total CN values into bins as the following: 0 = [0;6[; 1 = [6;8[; 2 = [8; 12[; 3 = [12; 15[; 4 = [15;20[; 5 = [20;60[ and 6 = [60;500[ except for 63 genes belonging to the IC1 cytoband (17q24 and adjacent cytobands 17q23 and 17q25) for which the bins are: 0 = [0;4[; 1 = [4;8[; 2 = [8;12[; 3 = [12; 15[; 4 = [15;20[; 5 = [20;60[ and 6 = [60;500[. Indeed, the IC1- associated amplification is globally lower than the events occurring in IC2, IC6, and IC9 cases. Then, we calculated the profiles of the four cytoband-level CN (17q24, 11 q13, 8p11 , and 8q24) with the averaged binned CN value of all genes falling in each of the four IC-specific chromosomal cytobands. Finally, we weighed the remaining features as follows: fraction of CNA has been multiplied by factor 2; WGD, HRD LOH, and receptor status features have been multiplied by 4; and average breakpoints per chromosome has been multiplied by 10. Of note, ploidy, WES, and WGS are the only features that were not transformed.

[0127] Model development strategy. We trained, tested, and validated the ENiClust classifier on a set of breast cancer specimens with WES data (TCGA-WES, n=1 ,013), mostly primary (n=1 ,009 / 4 primary / metastatic lesions) and on a set of samples with WGS data (Hartwig, n=297), mostly metastatic (n=30 / 267 primary / metastatic lesions). We inferred the reference classification by using both DNA and RNA-sequencing data. Of note, because of their close CN profiles and similar clinical outcomes, we merged IC3and IC7 as a single class ( IC3 / IC7), thus we trained ENiClust to predict the nine following classes: IC1 , IC2, IC3 / IC7, IC4, IC5, IC6, IC8, IC9, and IC10. We first split the initial dataset (with TCGA-WES and Hartwig-WGS) following an 80:20 ratio in training / testing:validation datasets, with 1 ,047 and 263 samples, respectively. We later used the validation set as “unseen” (holdout dataset) data for performance evaluation. We further split the remaining training / testing dataset again across a 5-fold crossvalidation approach resulting in training and testing sets, with ~n=838 and ~n=209 samples, respectively. To assure the balanced IC distribution at each split, we used the stratified split methods. To compensate for class imbalance in our datasets (due to the variable distribution of IC subtypes), we used the SMOTE methods (SMOTE function from imbalanced-learn python package v.0.5) on the training set at each fold.

[0128] We developed the ENiClust model based on three complementary learning algorithms: (i) random forest (RF) with 111 trees with maximum depth set at 70, minimum of 2 samples at the leaf node, and no sample bootstrapping when building trees; (ii) Elastic-Net (ENET) using the SAGA solver algorithm for a relatively large dataset and multiclass problem in a balanced mode and with the maximum number of iterations at 1000 for the algorithm, mixing parameter of 0.9, the inverse of regularization strength at 0.2; and (iii) gradient booster (GBT) with the learning rate set up at 0.2 boosting states with maximum depth set at 3, and minimum of 4 samples at the leaf node. The remaining parameters were set at their default value. We performed the tunning of these learning algorithms’ parameters (RF: n_estimators, max_depth, and min_samples_leaf; ENET: I1_ratio and C; GBT: learning_rate, max_depth, and min_samples_leaf) through the 5- fold cross-validation using Random izedSearchCV (from scikit-learn) through 100 iterations and for RF, ENET, and GBT performance optimization individually.

[0129] We selected hyperparameters based on a customized performance score evaluated on each fold: avg( 0.75 x minimum(Precision-Recall (PR) AUCim, icz, ice, ice), 0.25 x minimum(F1-scoreici, ic2, ice, ics, ER+ Typical, ics, TNBC)) where ER+ Typical = IC8, IC3 / IC7, IC4ER+ and TNBC = IC10, IC4ER-.PR AUC is indicated for evaluating the classification performance for minority classes as the ER+ High-risk ICs (< 30%). We included the F1 -score across IC subgroups to optimize ENiClust global performance across all subgroups. The ER+ Typical ICs (IC3 / IC7, IC4ER+, IC8) were grouped giventheir similar molecular profiles and clinical outcomes as were TNBC ICs (IC10, IC4ER-) that have proven challenging to distinguish. We weighted PR AUCs in the ER+ High-risk ICs 3 times more than the F1 -score across all subgroups because of the importance and challenge in distinguishing these classes of aggressive disease. The performance score was calculated for TCGA tumors samples.

[0130] Therefore, we built the ENiClust ensemble model following the voting method (voting function from scikit-leam python package v.0.21.3, soft cut), trained on the full dataset (training / testing, n=1 ,047) with the above selected hyperparameters. We performed the final performance evaluation of ENiClust on two independent validation sets: the validation set defined above (TCGA n=187) and the independent dataset with primary breast cancer RNA- and WGS-sequenced samples (n=189). Of note, we reprocessed the RNA-sequenced ICGC set as described in the RNA sequencing section.

[0131] ENiClust yielded high balanced accuracy for the four ER+ High-risk ICs (IC1 , IC2, IC6, and IC9) subtypes (TCGA-WES Validation: 0.77-0.91 ; Nik-Zainal-WGS Validation: 0.72-0.99; METABRICASCAT: 0.74-0.87), as well as IC5 (HER2+, TCGA- WES Validation: 0.95; Nik-Zainal-WGS Validation: 0.63; METABRIC-ASCAT: 0.91 ) and IC10 (TNBC, TCGA-WES Testing: 0.88; NikZainal-WGS Validation: 0.89; METABRICASCAT: 0.84) (Precision, Recall, F1 -score, and Precision / Recall AUC displayed in Fig. 4D).Pair Primary and Metastatic Tumor Subtypes

[0132] We used ENiClust to classify the 54 primary and 82 metastatic tumors from the paired dataset. We obtained receptor status from the respective original publications. For 21 patients, receptors were profiled by sample, while for the rest they were profiled only for the primary sample.Estimating ecDNA Timing

[0133] The number of pre-amplified mutations in each amplified region was used to estimate the timing of the amplification. After selecting SNVs located in each amplicon, the mutated copy number (n ) for each SNV was calculated using the previously described formula below.

[0134] where v is variant-allele frequency, p is tumor purity, and n and n are tumor and normal absolute copy number of the region. SNVs were classified as pre- or postamplified mutation using the n values as described previously102. If the n value was higher than 75% of the major allele copy number of the segment, the SNV was classified as a pre-amplified one. Since kataegis do not have the clock-like property, kataegis was excluded from the mutation count. Kataegis was defined if three or more consecutive mutations were located within 3,000bps. To exclude the effect of non-clock-like mutations, mutational signature analysis was performed for merged preamplification mutations from each subgroup of tumors (ER+ High-risk, ER+ Typical-risk, HER2+, and IC10). The final number of pre-amplified mutations was corrected by the proportion of clocklike mutational signatures (SBS1 , 5, and 40).

[0135] Assessing ecDNA presence in cell line models. We adapted the scEC&T- seq (single cell extrachromosomal circular DNA and transcriptom ic sequencing) protocol to detect ecDNA in representative IC-subtype matched cell lines, namely UCD12 (IC6) and LICD65 (IC2)104. The protocol involves several key steps: (1 ) DNA extraction using HMW DNA extraction kits; (2) Ampure Beads cleanup; (3) Exonuclease digestion of linear DNA followed by a 3-day incubation at 37°C and subsequent heat inactivation; (4) Addition of Plasmid-Safe ATP-dependent DNAse every 24 hours with ATP refreshment; (5) PEG buffer treatment and magnetic bead separation; (6) Incubation with buffer D2 at 65°C and addition of STOP solution; (7) Rolling circle amplification using REPLI-g sc DNA polymerase for 8 hours at 30°C; (8) Inactivation of REPLI-g sc DNA polymerase at 65°C for 3 minutes; (9) Ampure bead cleanup; (10) Elution in EB buffer and splitting into two tubes for Nanopore long-read and shotgun sequencing. The method was capable of capturing both large and small circular DNAs in the specimens. DNA input samples were documented for reference, and the amplified exonuclease-digested DNA quantified for subsequent sequencing analysis.

[0136] Assessing coverage in predicted ecDNA. For each sample, i.e. parental and digested, we binned raw read counts per 50 kB genomic bin across the genome. We calculated the average coverage on each of these bins by multiplying the read count by150 (the length of the sequencing reads) and dividing by the bin size (50 kB). We calculated average coverage on the predicted ecDNA segments (or linear amplicons in the case of the negative control), for both the digested and parental samples, and normalized this value by their respective median average coverage across all bins. We generated null “ecDNA regions” (i.e. simulated ecDNA composed of randomly generated genomic regions of the same size, number and gene density as the original ecDNA) only from mappable regions on the genome to remove unmappable regions. To do so, we randomly selected a chromosome followed by a genomic position within its mappable regions. If a segment of the desired length could be obtained from this starting position, we stored it as a segment of the null ecDNA. If the mappable region ended before the full length of the predicted ecDNA segment could be matched, we randomly sampled a new start position. Once all length-matched segments were generated, we returned the null “ecDNA regions”. We ran the same procedure as described previously for 1000 iterations. For each sample, we compared the average coveraged in the predicted ecDNA region to this null distribution and calculated a p-value by subtracting the proportion of null log ratios that were less than the predicted log ratio from 1 .

[0137] R-loops in MCF10A cell lines. We defined all R-loops as peaks with adjusted p-value < 0.05 and Iog2 fold change > 2 at baseline (TO) compared to input. We defined ER-induced R-loops as peaks with adjusted p-value < 0.05 and Iog2 fold change > 2 at 24 hours after E2 treatment compared to baseline. We calculated the number of each type of R-loop in cyclic and non-cyclic amplifications and calculated the R-loop density by dividing the number of R-loops by the size of the amplicon (number of R-loops per Mbp). One limitation of defining R-loops in MCF7 cells is that MCF7 cells harbor ecDNA on chromosome 17 and therefore may not reflect the regulatory elements present when ecDNA is formed in IC1 or HER2+ tumors for which chr17 amplicons have been found in the three stages of progression. Next, we downloaded DRIP-Seq data in wildtype and APOBEC3B knockdown MCF10A cells treated with DMSO or phorbol 12-myristate 13- acetate (PMA; 6 hours). We compared the change in R-loop density (number of R-loops per Mbp) in wildtype vs A3B knockout cells in cyclic vs non-cyclic amplifications at baseline and after PMA treatment. Finally, we calculated the distance of the IC-specificoncogenes recurrently incorporated in ecDNA to the nearest ER-induced R-loop, and compared them the distance of each gene as the background distribution.

[0138] Replication timing in MCF7 cell lines. Replication timing weighted average (WA) scores were calculated. Briefly, the PNDV value represents the percentage of replication occurring at a given 1 kb locus within a timing fraction. PNDVs were then converted to a single replication timing WA score per 1 kb loci. The following formula was used: WA = (0.917*G1 ) + (0.750*S1 ) + (0.583*S2) + (0.417*S3) + (0.250*S4) + (0*G2). We compared the median WA in cyclic vs complex non-cyclic regions using Mann- Whitney Rank Sum test.

[0139] ER-induced structural variants. Leveraged high-throughput genome-wide chromosomal translocation sequencing (HTGTS) to characterize translocations with and without estradiol treatment in MCF7 breast cancer cell lines initiated by double stranded DNA breaks in SHANK2 or RARA. We calculated the proportion of breakpoints in MCF7 cells within and without E2 treatment that overlapped all R-loops or only ER-induced Rloops in MCF7 cells (defined by GSE81851 ; see R-loops in MCF7 and MCF10A cell lines).

[0140] APOBEC3B deletion polymorphism. To predict the presence of APOBEC3B germline deletion polymorphism, we calculated the average coverage depth for genes APOBEC3B, APOBEC3C, APOBEC3D, APOBEC3F, APOBEC3G and APOBEC3H. We compared the average depth of APOBEC3B (d) with the average depth of the remaining genes (d2) and performed a maximum likelihood test considering three hypotheses:• Assuming no APOBEC3B deletion polymorphism, expected depth is Pois(d2)• Assuming heterozygous deletion allele, expected depth is Po / s(di), where di = d2 / 2• Assuming homozygous deletion allele, expected depth is estimated as number of misreads, set as do = d2 / 20

[0141] Average depths d, d2, di and do were rounded to the nearest integer to fit the Poisson distribution. The polymorphism was called for primary samples with subgroup and ecDNA calls, considering only one sample per donor (n = 2659).Survival Analysis

[0142] We performed survival analysis and corrected for key clinical covariates (age, grade, tumor size, lymph node, and ER status). We report ER status from immunohistochemistry (IHC) based clinical annotation when available, and RNA expression otherwise. We generated Kaplan-Meier and forest plots. The probability of relapse over time within each IC subtypes and histo-IC subgroups is reported.Treatment Response

[0143] We assessed the predictive value of each SV signature leveraging the treatment-response annotations in the Hartwig dataset for which metastatic biopsies were collected from patients with advanced cancer not curable by local treatment options. We compared three response categories: Progressive Disease, Stable Disease, and Complete / Partial Response. We used linear regression models testing the relationship between each SV signature and treatment response within each treatment group: (i) chemotherapy (n=291 ), (ii) endocrine therapy (n=256), and (iii) targeted therapy (n=224); and corrected for other additional treatment and the four IC subgroups (ER+ Typical, ER+ High, HER2+, and TNBC).

[0144] We also assessed the association of the ER+ IC subtypes (High-risk and Typical-risk) with endocrine therapy response taking advantage of two datasets named Giltnane and Selli. For the Giltnane dataset, we normalized the expected read counts from the RNA-seq of n=50 samples to library size and transformed to approximate homoscedasticity. The patients’ ER+ tumors have been classified as sensitive, intermediate and resistant cases from the clinical trial. In this trial, tumors were classified based on their post-treatment Ki67 score as sensitive (Ki67 < 2.7%), intermediate (Ki67 2.8-7.3%), or resistant (Ki67 > 7.3%) to endocrine therapy. We used this RNA-seq data and the iC10 R package classifier to call the IC subtypes as described in the Ensemble Integrative Clustering (ENiClust) classifier - Reference preparation section. For the Selli dataset, we downloaded the non-normalized RNA microarray data from NCBI GEO (n=101 samples from GSE111563, n=75 samples from GSE59515, and n=36 samples from GSE55374) and we normalized the expression data before calling the IC subtypes using the iC10 R package classifier to call the IC subtypes as described in the EnsembleIntegrative Clustering (ENiClust) classifier - Reference preparation section). The mRNA proliferation score has been calculated as the average expression of the 10 / 11 PAM50 proliferation score genes available for this dataset: BIRC5, CCNB1, CDC20, NUF2, CEP55, MKI67, PTTG1, RRM2, TYMS, and UBE2C.TME Subtypes

[0145] To classify the additional primary breast cancer samples from the METABRIC cohort and metastatic samples from the Hartwig cohort into the TME subtypes. We scored the 29 signatures defining the TME subtypes, performed median normalization of the resulting scores, and classified them considering the pan-cancer TCGA samples. We performed median normalization of primary samples across the cohort, while for metastatic samples, we normalized by metastatic site, considering only those sites with at least 10 samples.Single-Cell Data

[0146] To define the proportion of immune and stromal cells in 384 METABRIC samples and compare these with the TME subtypes called from bulk RNA (see TME subtypes). The results from single-cell proteomic data supported our TME subtyping strategy: significantly higher proportions of lymphocytes (B-cells and T-cells) were present in immune-enriched phenotypes, while significantly higher proportions of fibroblasts in fibrotic phenotypes.Cell Type Deconvolution

[0147] We applied the Cell Fraction analysis module to the TCGA RNA-seq expression data, Iog2 intensities from METABRIC expression arrays, and the expected counts from the Hartwig cohort. We built a custom reference matrix. We used S-mode batch correction, 0 as the fraction parameter for the construction of the reference matrix, 100 permutations for p-value calculation, and all other default parameters.Genetic Immune Alterations

[0148] To define the alterations associated with genetic mechanisms of immune escape (GIE), we considered 29 genes across 7 pathways associated with this phenomenon. Additionally, we considered the HLA-II pathway, made up of the following genes: HLA-DRA, HLA-DRB1, HLA-DMA, HLA-DPA1, HLA-DPB1, HLADMB, HLA-DQB1 and HLA-DQA1. We considered alterations that would result in immune escape, which varies depending on the pathway. For the Epigenetic and PD-L1 pathways, we only considered amplifications as the damaging alterations. For the rest of the pathways, we considered single nucleotide variants (SNVs), damaging structural variants (SVs; see Damaging structural variant), and homozygous deletions. Additionally, for the HLA-I and HLA-II pathways, we considered HLA loss of heterozygosity (HLA-LOH) as another damaging alteration.

[0149] Amplifications were called (see Extrachromosomal DNA (ecDNA) detection). From FACETS, we considered a sample to have an amplification in a specific gene if the CN state was called as follows: AMP, AMP (BALANCED), AMP (LOH) and AMP (many states). We additionally considered complex cyclic and non-cyclic amplifications as called by Amplicon Architect as additional alteration classes. If an amplification was called from both programs, only the complex amplification from Amplicon Architect was retained.

[0150] Homozygous deletions were called, considering the copy number state HOMDEL. HLA loss of heterozygosity was called, defined as having a minor CN of zero for genes in either the HLA-I or HLA-II pathways. For HLA-I, if a sample was originally homozygous for the specific allele, HLA LOH was not called. For the samples in our cohort, we compared this approach with the published calls from LILAC by examining the number of samples with HLA LOH in any of the HLAI alleles, and observe strong concordance in primary (97.7%, n = 217), and metastatic (98.2%, n = 679) samples.

[0151] For SNVs, we considered variants as described above (see Single nucleotide variant (SNV) detection), and restricted our analyses to protein coding variants corresponding to the following classes: Missense_Mutation, Splice_Site, Frame_Shift_Del, Splice_Region, Nonsense_Mutation, Frame_Shift_lns, ln_Frame_lns, ln_Frame_Del, Translation_Start_Site, and Nonstop_Mutation. Finally, structural variantswere called from the consensus approach, and restricted our analyses to damaging SVs as described above (see Structural variant detection).HLA Typing

[0152] We extracted reads mapping to chromosome 6 and remapped them to HLA reference fasta. We then back-extracted mapped reads fastq and used them as input to OptiTypePipeline.py using default parameters.Neoantigen Calls

[0153] Taking in our previously called SNVs and HLA (see Single nucleotide variant (SNV) detection and HLA typing), we perform antigen prediction on each sample. Briefly, we estimated three parameters for each putative epitope-HLA pair: 1 ) the binding strength between epitope and HLA; 2) foreignness score; 3) agretopicity score. We then filtered the initial list of putative neoantigens based on default thresholds of foreignness score and agretopicity score, as well as less than 10OOnM for binding strength. We called clonal neoantigens. Briefly, we inferred the cancer cell fraction (CCF) values for SNVs from the corresponding neoantigen based on their VAF, local ploidy, and purity. We categorized neoantigens as clonal if the corresponding SNV was annotated as clonal, and subclonal if otherwise. Finally, we assessed the association between neoantigens and SNV burden using a linear regression model neoantigen ~ GIE or neoantigen ~ GIE + log(SNV). For metastatic samples, we only considered tumors with purity greater than 0.6.

Claims

WHAT IS CLAIMED IS:1 . A method to treat an individual having breast cancer, comprising: classifying a breast cancer into a high-risk ER+ molecular group or into a HER2+ molecular group; and administering the individual an antagonist of cGAS-STING.

2. The method as in claim 1 , wherein the breast cancer is classified into a high- risk ER+ molecular group selected from IC1 , IC2, IC6, or IC9.

3. The method as in claim 1 , wherein the breast cancer is classified into the HER2+ molecular group IC5.

4. A method to treat an individual having breast cancer, comprising: classifying a breast cancer into a high-risk ER+ molecular group or into a HER2+ molecular group; administering the individual an antagonist of PARP.

5. The method as in claim 4, wherein the breast cancer is classified into a high- risk ER+ molecular group selected from IC1 , IC2, IC6, or IC9.

6. The method as in claim 4, wherein the breast cancer is classified into the HER2+ molecular group IC5.

7. A method to treat an individual having breast cancer, comprising: classifying a breast cancer into a high-risk ER+ molecular group or into a HER2+ molecular group; administering the individual an antagonist of BET.

8. The method as in claim 7, wherein the breast cancer is classified into a high- risk ER+ molecular group selected from IC1 , IC2, IC6, or IC9.

9. The method as in claim 7, wherein the breast cancer is classified into theHER2+ molecular group IC5.

10. A method to treat an individual having breast cancer, comprising: classifying a breast cancer into a TNBC molecular group; and administering the individual an antagonist of cGAS-STING.

11. The method as in claim 10, wherein the breast cancer is classified into a TNBC molecular group selected from IC4ER- or IC10.

12. A method to treat an individual having breast cancer, comprising: classifying a breast cancer into a high-risk ER+ molecular group or into a HER2+ molecular group; and administering the individual an antagonist of APOBEC editing.

13. The method as in claim 12, wherein the breast cancer is classified into a high-risk ER+ molecular group selected from IC1 , IC2, IC6, or IC9.

14. The method as in claim 12, wherein the breast cancer is classified into the HER2+ molecular group IC5.

15. A method to treat an individual having breast cancer, comprising: classifying a breast cancer into a high-risk ER+ molecular group or into a HER2+ molecular group; and administering the individual an antagonist of TopB1 .

16. The method as in claim 15, wherein the breast cancer is classified into a high-risk ER+ molecular group selected from IC1 , IC2, IC6, or IC9.

17. The method as in claim 15, wherein the breast cancer is classified into the HER2+ molecular group IC5.

18. A method to treat an individual having breast cancer, comprising: classifying a breast cancer into a high-risk ER+ molecular group or into a HER2+ molecular group; and administering the individual a modulator of G-quadruplexes.

19. The method as in claim 18, wherein the breast cancer is classified into a high-risk ER+ molecular group selected from IC1 , IC2, IC6, or IC9.

20. The method as in claim 18, wherein the breast cancer is classified into the HER2+ molecular group IC5.

Images