{'fraction_non_alphanumeric': 0.029580137706049343, 'fraction_numerical': 0.019369814457566775, 'mean_word_length': 5.247917687408133, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An omic and multidimensional spatial atlas from serial biopsies of an evolving metastatic breast cancer

Summary Mechanisms of therapeutic resistance and vulnerability evolve in metastatic cancers as tumor cells and extrinsic microenvironmental influences change during treatment. To support the development of methods for identifying these mechanisms in individual people, here we present an omic and multidimensional spatial (OMS) atlas generated from four serial biopsies of an individual with metastatic breast cancer during 3.5 years of therapy. This resource links detailed, longitudinal clinical metadata that includes treatment times and doses, anatomic imaging, and blood-based response measurements to clinical and exploratory analyses, which includes comprehensive DNA, RNA, and protein profiles; images of multiplexed immunostaining; and 2- and 3-dimensional scanning electron micrographs. These data report aspects of heterogeneity and evolution of the cancer genome, signaling pathways, immune microenvironment, cellular composition and organization, and ultrastructure. We present illustrative examples of how integrative analyses of these data reveal potential mechanisms of response and resistance and suggest novel therapeutic vulnerabilities.


In brief
Identifying mechanisms of response and resistance to treatment in individual cancer patients is challenging but critical for improvement of precision medicine outcomes. Johnson et al. report a comprehensive atlas from a single individual with breast cancer and demonstrate how longitudinal, integrative analyses can provide new insights.

INTRODUCTION
Precision medicine has led to substantial improvements in clinical outcomes for some individuals with cancer, increasingly through use of analytical procedures that identify people with molecular characteristics associated with an increased likelihood of response. 1,2 Unfortunately, treatments deployed according to precision medicine principles do not always elicit pos-itive responses, and durable control is achieved for only a subset of individuals with metastatic cancer. 3 We posit that the failure to control individual cancers using biomarker-guided treatments stems in large part from our imperfect understanding of the multitude of resistance mechanisms that drive an individual tumor's adaptive ability to survive as they evolve under therapy. These mechanisms may involve regulatory networks intrinsic to tumor cells, chemical and mechanical influences from proximal or distal microenvironments, and/or aspects of the immune system. They may vary between individuals with similar biomarkers, across metastases within a person, or among cell subpopulations within a single lesion and may change during treatment.
To stimulate and support community-wide investigations of cancer resistance, response, and evolutionary mechanisms in individuals, we present a comprehensive omic and multidimensional spatial (OMS) atlas composed of clinical and research data, response correlates from an affected individuals, and illustrative analytical workflows from a single person with metastatic breast cancer during 3.5 years of treatment. We also illustrate how the atlas can be used to explore spatial and temporal heterogeneity, study tumor evolution, and identify candidate resistance mechanisms and therapeutic vulnerabilities. These studies were carried out in the Serial Measurements of Molecular and Architectural Responses to Therapy (SMMART) program with support from the Human Tumor Atlas Network (HTAN), 4,5 through which all data are available in standardized formats.

Longitudinal data generation from a single individual
The focus of this OMS atlas is a female individual diagnosed with hormone receptor-positive, HER2-normal, high-OncotypeDx recurrence score, 6 0.6-cm right breast ductal carcinoma at the age of 64. She underwent a lumpectomy with intra-operative radiation therapy, followed by treatment with four cycles of adjuvant docetaxel and cyclophosphamide, 2 years of anastrozole, and 5 years of exemestane. Subsequent computed tomography (CT) and fluorodeoxyglucose-positron emission tomography (FDG-PET) scans revealed widespread metastatic disease.
The individual was then enrolled in the SMMART program ( Figure 1). Management decisions were made by the treating physician based on all clinical information plus input from a multidisciplinary tumor board (STAR Methods). This led to four treatment phases over a 3.5-year period ( Figure 2A). Temporary tumor control was achieved in the first three phases, with a new phase of therapy beginning at signs of progression. Toxicity of the combination therapies was effectively managed through supportive medication and dose reduction ( Figure S2A). Standard toxicity-related blood chemistries were monitored, including absolute neutrophil and platelet counts and liver function tests (Figures S2B-S2D; Table S1).
The clinical metadata in Table S1 link detailed treatment doses and timelines ( Figure 2A) to tumor response metrics. Serum levels of the tumor protein biomarkers carcinoembryonic antigen (CEA), cancer antigen 15-3 (CA 15-3), and cancer antigen 27-29 (CA 27-29) were routinely measured to monitor treatment response ( Figure 2B; Table S1). Increasing biomarker levels are concerning for underlying progression, but National Comprehensive Cancer Network (NCCN) guidelines do not recommend changing therapies solely based on blood biomarkers. 7 Biomarker measurements were thus complemented by periodic CT and FDG-PET imaging, with response evaluated using response evaluation criteria in solid tumors (RECIST) 1.1 criteria 8 ( Figure 2C; Table S1). Representative computed tomography (CT), fluorodeoxyglucose-positron emission tomography (FDG-PET), and ultrasound images highlight disease burden at key time points ( Figure S1).
Biospecimens collected for analysis include serial blood samples, a primary breast tumor (PT), a liver biopsy taken immediately prior to phase 1 (Bx1), a biopsy of a different liver lesion taken at the end of phase 1 (Bx2), a bone lesion biopsy taken at the end of phase 3 (Bx3), and a biopsy of a third liver lesion taken at the end of phase 4 (Bx4; Figure 2A). Importantly, Bx2- Bx4 were acquired from metastatic lesions explicitly identified on serial CT and/or FDG-PET imaging as progressing near the end of each respective treatment phase (Figures 2A and S1A-S1I). These biospecimens were analyzed using 11 distinct omic and multiscale spatial imaging workflows to generate this OMS atlas ( Figure 1).
Genomic differences between metastases were substantial Targeted DNA sequencing (GeneTrails solid tumor panel, all biopsies), whole-exome sequencing (WES; primary tumor and Bx1-Bx4), and low-pass whole genome sequencing (LP-WGS; Bx3 and Bx4) were used to identify somatic genomic alterations, including single-nucleotide variants (SNVs), insertions or deletions (indels), and copy number changes ( Figures 3A and  S3A). Ubiquitous alterations included amplification of the CDK4/6 regulatory partner cyclin D1 (CCND1; Figure S3B). Other biologically and clinically relevant alterations were private to a single biopsy. For example, Bx2 contained a hotspot PIK3CA mutation (p.E542K; GenBank: NM_006218:c.1624G>A) 9,10 that was absent from other samples and was taken from a liver lesion that increased in size during treatment with drugs that target aspects of phosphatidylinositol 3-kinase (PI3K) signaling (Figures 2A and 2C). Bx3 and Bx4 both harbored similar amplified regions on chromosome 18 that were not detected in prior biopsies or the primary tumor; Bx3 had 8 copies of a 2.3-Mb region, and Bx4 had 14 copies of a 0.7-Mb region ( Figures  3B and S3B). These amplicons contained the genes for thymidylate synthetase (TYMS) and the SRC family tyrosine kinase YES1 and were accompanied by increased TYMS and YES1 RNA relative to Bx2 (TYMS: B33 = 6.83, B34 = 7.23; YES1: B33 = 2.03, B34 = 4.03). Importantly, both biopsies were acquired after treatment with the TYMS inhibitor capecitabine (Figures 2A, 2C, and S1G).
A phylogenetic analysis revealed that Bx3 diverged from the primary tumor at an earlier evolutionary stage than Bx1, Bx2, or Bx4 ( Figure 3B) but was only detected on FDG-PET imaging 1 month before the biopsy occurred (Figures S1A, S1F, and S1G).
WES of circulating tumor DNA (ctDNA) from blood collected immediately prior to Bx1 (ctDNA1) and 23 days after Bx2 (ctDNA2) showed that ctDNA1 carried mutations identified previously as private to Bx2, Bx3, or Bx4, whereas ctDNA2 had mutations that were private to each of the four biopsies ( Figure S3A). Thus, at least some of the genomic features detected in later biopsies were present before initiation of treatment.
Tumor mutational burden (TMB) was assessed for the primary tumor and Bx1-Bx4 because a TMB of 10 or more mutations per megabase (mut/Mb) has been associated with a positive response to immune checkpoint blockade. 11 The TMB was low overall (1.2-5.2 mut/Mb), but we identified 1,271 unique neoepitopes (158-687 neoepitopes per biopsy) predicted to bind to at least one major histocompatibility complex (MHC) allele with an affinity of less than 500 nM (Table S2) Article ll OPEN ACCESS antigen (HLA) subtypes were stable across all biopsies, and no loss of heterozygosity was observed. Notably, 68 neoepitopes that might serve as targets for a personalized cancer vaccine were present in the primary tumor and all four biopsies ( Figure S4A). 12,13 ctDNA increased during progression and radiotherapy Dual index degenerate adaptor sequencing (DIDA-seq) 14 was performed, using a panel of 53 SNVs present in the individual's primary tumor, Bx1, and/or Bx2 to assess ctDNA levels from serial plasma samples collected over the first 32 months of treatment (Table S2). The average variant allele frequency (VAF) of the SNV panel remained below 0.3% of total cell free DNA during this period, with the exception of two transient increases ( Figure S2E). The first occurred immediately prior to Bx2 ( Figure S2E), coincident with rising CA 15-3 and CA 27-29 levels ( Figure 2B), followed by progressive disease (PD; Fig-ure 2C). The increase in ctDNA VAF was greatest for mutations shared by the primary tumor, Bx1, and Bx2 (30% VAF) compared with those private to the metastases (Bx1 and Bx2, 3.1%; Bx1, 0.05%; Bx2, 1.3%). A second ctDNA increase occurred after palliative radiation therapy to spinal lesions at C2-C5. Interestingly, the VAFs of all SNV groups in the panel increased at this time, including those private to Bx1 and Bx2 liver lesions.
Signaling and pathway activities evolved during therapy Signaling and pathway activities were calculated from wholetranscriptome sequencing (RNA sequencing [RNA-seq]). Classification using the PAM50 subtype gene signature 15 showed liver biopsies Bx1, Bx2, and Bx4 to be luminal A, whereas the bone biopsy Bx3 was luminal B ( Figure S3C). Table S2 summarizes RNA transcript levels and pathway activity estimates for Bx1-Bx4 relative to breast cancers in The Cancer Genome Atlas Figure 3. Genomic, transcriptomic, and proteomic profiles reveal spatiotemporal heterogeneity and evolution (A) Comparison of somatic mutations. Columns represent individual, non-silent SNVs or indels identified from WES in at least one tissue sample and classified as ubiquitous (present in all samples, blue), shared (present in at least two samples, green), or private (present in only a single sample, red). Mutational status in each sample is indicated as independently called (colored), detected in at least 2 sequencing reads but not independently called (reduced opacity), or absent (white). (B) Phylogenetic tree showing the evolutionary relationship between the PT and four metastases. (C) Transcriptomic gene set variation analysis (GSVA) of cancer hallmark pathways. The boxplot represents the distribution (upper and lower quartiles and median) of GSVA scores for the TCGA luminal breast cancer cohort. Enrichment scores are shown for each of the biopsy samples: Bx1 (green), Bx2 (orange), Bx3 (blue), and Bx4 (pink). (D) RPPA protein pathway activity assessment using pathway scores. The boxplots represent the distribution of the pathway activity of the TCGA breast cancer cohort. The pathway activities of three biopsy samples are marked as in (D). (E) Total and phosphoprotein levels from RPPA normalized within the TCGA breast cancer cohort. The heatmap shows relative protein levels for three biopsies and the fold change between sample pairs. Proteins are ordered based on the fold change difference between Bx2 relative to Bx1. Selected proteins are highlighted. (F) ISPP measurements of total and phosphoprotein levels. The boxplots represent the distribution of protein levels of 57 metastatic breast cancers. The protein levels of three biopsy samples are marked as in (D). See also Figure S3 and Table S2. (TCGA-BRCA) 16 and gene set variation analysis (GSVA) of enriched Molecular Signatures Database (MSigDB) cancer hallmarks relative to TCGA-BRCA luminal samples. 17,18 Proliferation, Immune, and Signaling were the most variable (MSigDB) Hallmark Process categories across the biopsies ( Figure 3C). Notably, Bx2 harboring the PIK3CA p.E542K mutation had reduced ''PI3K/AKT/mTOR'' compared with Bx1, although that gene set was still increased relative to TCGA-BRCA samples.
Spatial analyses indicated that tumor cells formed nests surrounded by immune, fibroblast, and endothelial cells as well as collagen I and collagen IV deposits. This was observed in all biopsies but was pronounced in Bx3. Quantitative analyses of nuclear ER and PCNA expression in Bx2-Bx4 as a function of distance to collagen I-rich tumor nest boundaries showed that cells expressing higher levels of ER and PCNA were closest to these boundaries and other stromal compositions (Bx2: mean ER in- Tumor-tumor and tumor-stromal interactions in Bx1, Bx2, and Bx4 were explored at $4-nm resolution using FIB-SEM. 41 Computational renderings of 3D images of Bx1 (Videos S1 and S2) and Bx2 (Video S3) revealed a previously unappreciated lattice-like structure for fibroblast-like cells surrounding tumor cell clusters and an intricate interaction pattern between these cells, collagen bundles, and tumor cells on the nest boundaries (Figure 5D). The production of collagen by tumor-associated fibroblast-like cells is particularly apparent in the 2D SEM image of Bx4 ( Figure 6A).

Intracellular nanobiology defined by FIB-SEM
3D FIB-SEM images of cancer cells in Videos S1, S2, and S3 provide important details about intracellular structures and interactions that may influence cell function and therapeutic response. These include the following. (1) Numerous $100-nm-diameter, micrometers-long, filopodium-like protrusions (FLPs) and lamellipodia that project from tumor cells into the stromal environment ( Figures 6B-6D; Video S3). Published work and our studies in model systems show that these protrusions have actin-rich cores and are decorated with receptor tyrosine kinases that are transported along FLPs by the actin-motor protein Myosin-X. 42 Cultured SKBR3 breast cancer cells exhibit similar FLPs, and dynamic in vitro images acquired using stochastic optical reconstruction microscopy (STORM) 43 reveal that the FLPs respond to epidermal growth factor by rapidly decreasing in length, causing cell movement toward the anchored ends of the FLPs ( Figure S6B; Video S4). (2) Alignment of mitochondria along the length of an elongated cell and insinuation into nuclear folds (Video S3; Figures 6B and 6E). (3) A high abundance of lamellipodia and macropinosomes, implicating nutrient scavenging via macropinocytosis as a possible tumor survival mechanism (Videos S1 and S3; Figures 6C, 6F, and 6G). 44,45 (4) A high prevalence of densely stained vesicles that appear to be lysosomes (Videos S1 and S3; Figures 6A and 6G). Figure 6H presents a qualitative summary of the nanoscale features described in Bx1 and Bx2, made by visual analysis of large-format 2D SEM images ( Figure S6A) and informed by 3D FIB-SEM images of selected features.

DISCUSSION
This OMS atlas is a compendium of cellular, molecular, and organizational features of four biopsies along with detailed clinical response data collected over a 3.5-year period from a single individual with metastatic breast cancer. It is intended to illustrate the feasibility of generating longitudinal multiplatform analyses in the clinical setting to support investigations of mechanisms of response and resistance beyond those that are apparent from routine omics analyses. Its features include DNA, RNA, and protein 19-21,23 profiles and spatially defined analyses, including mIHC, 27-29 CycIF, 39,40 and 2D and 3D electron microscopy. 41 Other key components are (1) preservation of samples starting within 2 min of biopsy to conserve labile molecular and architectural features, (2) precise temporal linking of clinical and molecular responses with drug treatments and doses, (3) quantitative assessment of individual lesion changes by CT imaging to measure response heterogeneity, (4) workflows that enable multiplatform measurements using material from a single biopsy, and (5) well-curated data and data standards developed by the HTAN program to facilitate community analyses and integration with other datasets.
Selected analyses encompassing multiple OMS atlas datasets illustrate approaches to uncovering mechanisms of drug resistance and response that arise over the course of treatment and that could be missed by limiting analyses to one or a few analytical platforms (Figure 7; Table S3; discussed below). It is important to note that the post hoc analyses described here were not used to guide treatment. They are not meant to definitively explain why the individual progressed on any given therapy but are intended to stimulate validation in follow-up studies.
Phase 1 treatment consisted of a combination of fulvestrant, palbociclib, and everolimus, supported by findings in Bx1 of high ER protein expression, wild-type ESR1, and two intact copies of wild-type RB1. Bx2 was taken when the tumor began to progress on this treatment ( Figure 2). Interestingly, none of the four biopsies analyzed from this individual had mutations in ESR1 or loss of expression of ER protein ( Figure 3E; Table S1) even though ESR1 mutations are frequently observed in individuals progressing on endocrine therapies. 46 Mutations in PIK3CA, ERBB2, and NF1 also are observed in individuals progressing on endocrine therapies, and RB1 is frequently lost after treatment with CDK4/6 inhibitors, 47,48 but only a PIK3CA mutation and immune-related pathway activation was seen in Bx2 ( Figure 3; Table S3). Thus, we interrogated Bx2 data to identify additional bypass mechanisms. One known mechanism by which cells become resistant to everolimus and other mTORC1 inhibitors is through activation of mTORC2. 49,50 Consistent with this, phosphoprotein analyses of Bx2 revealed decreased S6 phosphorylation, which supports continued inhibition of mTORC1 by everolimus ( Figure 3E). Concurrently, Bx2 had increased phosphorylation at an mTORC2 site on AKT (S473) and of multiple AKT substrates that together are predicted to maintain oncogenic PI3K/mTOR signaling in the presence of everolimus (Figures 3E and 7A). Everolimus efficacy might also have been reduced by the PIK3CA p.E542K activating mutation unique to Bx2, which is known to activate the PI3K/AKT/mTOR pathway. 51 Indeed, this variant was among the SNVs monitored in serial blood samples by DIDA-seq ( Figure S2E) and was only detected in ctDNA significantly above background after 7 months on phase 1 therapy (0.06% VAF, p = 0.0071, Weitzman overlapping coefficient), indicating that this mutation may have emerged because of selective pressure from one or more phase 1 drugs.
Several analyses inform on mechanisms of resistance to the CDK4/6 inhibitor palbociclib in Bx2. Loss of RB1 has been shown to drive resistance in multiple clinical trials, 48 but this gene was not mutated or deleted in Bx2. It is noteworthy that RB1 phosphorylation was at pre-treatment levels in this biopsy (P-S807/S811: 1.03 versus Bx1) because this modification pro-  Figure S6 and Videos S1, S2, S3, and S4. motes cell cycle progression and should have been decreased by palbociclib. 52 Evidence from protein profiling of key cell cycle regulators revealed that RB1 might have been phosphorylated by CDK2, which also inhibits RB1 but is not a target of palbociclib ( Figure 7B). 53-55 First, protein levels of the CDK2 inhibitor p21 were decreased 23 from Bx1 to Bx2 ( Figure 3E), possibly because of activated PI3K/AKT signaling, which maintains low p21 levels in CDK4/6 inhibitor-resistant cells. 55 Second, tumor cells with high CCND1 and activated PI3K can adapt to palbociclib via non-canonical binding of CDK2 to CCND1, 53 and Bx2 had higher CCND1 protein levels (1.73 versus Bx1) and PI3K/AKT signaling than Bx1. The CDK4/6 inhibitor abemaciclib has a broader spectrum of activity that includes CDK2 56 and might be expected to be effective in cases where palbociclib escape occurs via CDK2 activation. Indeed, abemaciclib administered subsequent to the period covered by this study showed efficacy (data not shown).
Capecitabine administered in phase 3 along with pembrolizumab, enzalutamide, and fulvestrant initially resulted in a partial response (PR), followed by PD ( Figure 2B Article ll OPEN ACCESS confers resistance to capecitabine. 57 Consequently, TYMS amplification may have provided a relative fitness advantage during capecitabine treatment and might explain the temporally late emergence of a clone that branched off early in the evolutionary process ( Figure 3B). The TYMS/YES1 amplicon arose independently in Bx4, presumably because of the earlier capecitabine exposure. But although only TYMS was overexpressed in Bx3, both genes were increased in Bx4, indicating that YES1 may have provided a growth advantage to later lesions after cessation of capecitabine ( Figure 2A). YES1 is an SRC family tyrosine kinase and a target of the broad-spectrum kinase inhibitor dasatinib, so inhibition of YES1 might be considered as a possible orthogonal therapeutic strategy for individuals who become resistant to capecitabine via amplification of TYMS/ YES1. However, dasatinib was administered subsequent to the period covered by this study and did not show efficacy (data not shown), arguing against this strategy.
Comparative analyses of the PT and serial biopsies suggested several mechanisms shaping immune contexture. The most significant was associated with palbociclib treatment at the time of Bx2. mIHC analyses showed increased macrophages/monocytes and T cells and decreased Tregs in Bx2 compared with Bx1 and Bx4 (Figure 4). Th17 cells and Treg cells arise from a common precursor but have opposing functionality upon terminal differentiation, with anti-tumor immunity promoted by Th17 cells and dampened by Treg cells. 58 This suggests that the relatively high frequency of Treg cells in the PT and Bx1 may have contributed to reduced T cell activation, as detected by a lack of PD-1 expression ( Figure 4G). Conversely, the Th17 dominance over Treg cells in Bx2 may have supported T cell activation, as evidenced by increased PD-1 expression. These changes were coincident with increases in signaling by IFNg, ILs, and STATs, as revealed by gene and protein expression profiles ( Figure 3C; Table S2), and are consistent with studies in mammary tumor models showing that CDK4/6 inhibitors promote T cell-mediated tumor cell clearance by stimulating type III interferons and suppressing Treg cell proliferation. 59 The utility of an immune checkpoint inhibitor was supported by our observations relating to Bx2 and increased PD-1 expression in T cells. Phase 2 and 3 pembrolizumab treatments were associated with a decrease in the Bx1 and Bx2 lesions ( Figure 2C), but the role of pembrolizumab in the decrease in lesion size is unknown as it was given with other drugs. Indeed, a challenge in this and other studies is in deconvoluting effects of individual agents when given in combination. Pembrolizumab and other phase 3 drugs were discontinued upon PD, whereupon the immune contexture changed again, with Bx4 showing more T cells but fewer macrophages/monocytes and Th17 T cells. Although Bx4 also contained fewer Treg cells (similar to Bx2) and the highest proportions of Th1 differentiation ( Figure 4G), there was low PD-1 expression on T cells ( Figure 4F). So, although the PT, Bx1, and Bx4 showed minimal T cell responses, this may have been due to low neoantigens and antigen presentation as likely barriers to functional anti-tumor immunity in Bx4, whereas T cell-mediated suppression was predominant in the PT and Bx1. (B) Activation status for cell cycle regulatory pathways affected by palbociclib in Bx2, as inferred from total and phosphoprotein levels. Palbociclib blocks cell division in responsive cells by inhibiting CDK4/6 phosphorylation of RB1, but Bx2 had continued high levels of phospho-RB1 and cell proliferation under treatment with this drug (RB1 P-S807/811 = 0.983 versus pre-treatment Bx1). This is possibly due to degradation of the CDK2 inhibitor p21 (0.503 versus Bx1) by activated PI3K/AKT signaling (see A), which would activate canonical cyclin E/CDK2 complexes to drive cells through G1-S. Alternatively, cell division might be proceeding through the formation of non-canonical cyclin D1/CDK2 complexes because of amplified CCND1 ( Figure S3B), high levels of cyclin D1 protein (1.673 versus Bx1), and low p21. CDK2 activation can be countered with the broad-spectrum CDK inhibitor abemaciclib. Inferred activation status is based on total protein levels or phosphorylation and is designated as relative increases (red), decreases (green), or unchanged/unknown (yellow). See also Figure 3E and Table S2. CycIF and FIB-SEM analyses showed tumor cells organized into nests surrounded by stromal cells and substantial collagen I deposits (Figures 5 and 6), suggesting that the lack of neoantigens and/or antigen presentation inferred from immune profiling may be caused, at least in part, by biophysical barriers that diminish tumor-immune cell interactions. Interestingly, the 3D images suggest that fibroblast-like cells are interposed between tumor cells and collagen bundles in most cases, raising the issue of how collagen stiffening leads to more aggressive tumor behavior 60 and how stromal barriers stimulate increased expression of ER and PCNA in closely proximal tumor cells (Figure 5D). From a technical perspective, the complex cellular interactions revealed by FIB-SEM illustrate the difficulties of properly segmenting individual cells during multiplex imaging of 2D sections using mIHC or CycIF ( Figure S6C) and of dissociating tightly interacting and potentially fragile cells for single-cell analyses.
The FIB-SEM analyses reveal several ultrastructural features that may influence tumor behavior and/or therapeutic vulnerability. These include the following. (1) FLPs that project from tumor cells into the stromal microenvironment. The receptordense, dynamic nature of FLPs may mediate proximal and distal interactions with elements of the microenvironment and enable directed movement therein. This might provide the force needed to produce the elongated tumor cell shown in Figure 6B, with mitochondria aligned along its long axis and inserted into nuclear folds ( Figure 6E). FLPs have also been implicated in protein transport between cells. 61 These functions suggest the possible utility of FLP inhibitors. (2) Insinuation of mitochondria into nuclear folds (Video S3; Figures 6B and 6E). These apparently forced interactions may increase the potential for nucleus-mitochondrion signaling that would alter DNA damage repair and/or reactive oxygen species (ROS) signaling. 62,63 This might be countered therapeutically by attacking ROS or by inhibiting FLP function. (3) High abundance of lamellipodia and macropinosomes (Videos S1 and S3; Figures 6C, 6F, and 6G). Nutrient scavenging from the intercellular space and nearby dying cells is a known tumor survival mechanism. 44,45,64 Protein-conjugated drugs might convert this survival mechanism into a therapeutic vulnerability. We speculate that this mechanism may have been partly responsible for the control achieved by treatment with liposomal doxorubicin during phase 3. Macropinocytosis may also diminish communication of neoantigens to immune cells by competing with dendritic cells for exogenous antigens released from dying tumor cells. 65 (4) High-density, densely stained vesicles that appear to be lysosomes (Videos S1 and S3; Figures 6A and 6G). Lysosomes can sequester cancer drugs via a process called lysosomotropism, in which weakly basic drugs become protonated and trapped within the acidic interior of lysosomes. 66 Lysosomotropic sequestration has been implicated as a mechanism of resistance to CDK4/6 inhibitors 67,68 and is suggested in this individual by the increase in lysosomes from Bx1 to Bx2, as seen by FIB-SEM ( Figures 6G, 6H, and S6A; Videos S1 and S3). Interestingly, hydroxychloroquine, sometimes used to counter treatment-induced rashes, has been reported to be lysosomotropic and thus might reduce treatment efficacy when co-administered with any basic drug, including CDK4/6 inhibitors. 69 Recent studies suggest that lysosomotropism-mediated doxorubicin resistance can be countered by the b-AR antagonist propranolol, which acts through a b-AR-independent mechanism to increase cytoplasmic doxorubicin concentrations and decrease lysosomal accumulation. 70 Overall, this OMS atlas shows the promises and challenges of elucidating evolving resistance mechanisms and new therapeutic vulnerabilities in individuals. The study shows that multianalyte workflows can be executed routinely and safely. Analyses of the data provide insights into mechanisms of tumor response and resistance that can be explored in subsequent studies. The ready availability of the data and protocols in standardized form will encourage further analyses and method development.

Limitations of the study
The overall goal of this study is to elucidate the mechanisms of resistance and therapeutic vulnerability experienced by an individual during extended treatment of metastatic disease, using information from multiple analytical workflows. We acknowledge the difficulty of assigning specific response and resistance mechanisms to individual drugs within a multi-drug treatment regimen, especially for drug combinations targeting overlapping biological pathways, and we note the danger of ''cherry picking'' mechanisms from the vast published literature. Moreover, working with a single human subject precludes implementation of hypothesis testing, which is de rigueur in experimental cell lines, animal models, and clinical studies that would more definitively support our proposed mechanisms. These factors demonstrate the challenges of implementing this type of program and analyzing N-of-1 data in a real-world clinical setting. However, our studies do suggest mechanisms and interpretational processes that can be tested in larger studies.
Several of our methods are too complex to be widely applied as currently implemented. However, when the utility of each assay platform is established, workflows can be simplified and streamlined. Our work shows that further development of analytical methods to integrate and interpret multi-platform omics and imaging datasets is clearly needed for the clinical and research communities. The OMS atlas can serve as a resource in that effort.
Finally, we are aware that we are inferring mechanisms from single biopsies of a metastatic disease that displays remarkable intra-and intermetastatic lesion heterogeneity. 71 A biopsy of a single metastatic lesion at any single time point is unlikely to provide a comprehensive picture of the entire heterogeneous disease within an individuals or even within the biopsied lesion. This is a fundamental limitation of any biopsy-based analytical strategy. Continued advancement of assays that report on overall tumor composition across multiple lesions, such as peripheral blood assays, is one potential avenue toward understanding heterogeneous disease burden. Indeed, our observation that radiation induced a transient increase in ctDNA in peripheral blood suggests that individuals undergoing radiotherapy might have circulating tumor nucleic acids and proteins in sufficient quantities for practical diagnostic measurement and for revealing latent, low-level molecular changes in unbiopsied lesions in almost real time.

RESOURCE AVAILABILITY
Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Joe Gray (grayjo@ohsu.edu).

Materials availability
This study did not generate new unique reagents.
Data and code availability dRaw data generated by next-generation sequencing platforms have been deposited at dbGaP and are publicly available as of the date of publication. The project accession number is listed in the key resources table. Processed next-generation sequencing data, protein expression data, as well as raw and processed image data have been deposited with the HTAN Data Coordinating Center. They are publicly available through the HTAN Data Portal as part of the HTAN OHSU Atlas as of the date of publication. The case number is listed in the key resources table. Processed images are also available for web-based viewing through the HTAN Imaging Data Portal as of the date of publication. The case number is listed in the key resources table. This paper also analyzes existing, publicly available data. The sources for these datasets are listed in the key resources table. d All original code has been deposited at Zenodo and is publicly available as of the date of publication. DOIs are listed in the key resources table. d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS
Human subjects This study was approved by the Oregon Health & Science University (OHSU) Institutional Review Board (IRB). All biospecimens and data were acquired and analyzed under the OHSU IRB-approved protocols Molecular Mechanisms of Tumor Evolution and Resistance to Therapy (IRB#16113) and Reconstructing the Tumor Genome in Peripheral Blood (IRB#10163). Participant eligibility was determined by the enrolling physician and informed written consent was obtained from all subjects. This study includes a single 64 year-old female woman (HTAN participant ID HTA9_1).  HCC1143, T47D, AU565, and HCC3153 were cultured in RPMI1640 with 10% FBS. All cell lines were derived from human female breast cancers, incubated in 5% CO2 at 37 C, and grown to near confluency. Genomic DNA was extracted from each cell line using DNeasy Blood and Tissue Kits (QIAGEN) and submitted to either Labcorp Cell Line Testing (Genetica) or the MD Anderson Cytogenetics and Cell Authentication Core (CCAC) for short tandem repeat (STR) analysis. Cell line identities were confirmed by comparison with reference STR profiles from the ATCC, DSMZ STR, or CCAC database.

METHOD DETAILS
Clinical decision making All clinical decisions were the treating physicians' discretion, with reference to clinical analytics results, established treatment guidelines, appropriate data from clinical trials, and input from an IRB-approved multidisciplinary tumor board. Research Use Only (RUO) data from exploratory analytics were not used for clinical decision making. In Phase 1, palbociclib and fulvestrant were used as the standard frontline for progression on adjuvant aromatase inhibitors. Everolimus was added due to rising tumor markers and clinical concern for symptomatic progression as well as emerging data from the TRINITI-1 trial that suggested PI3K/AKT/mTOR pathway targeting along with CDK4/6 to prevent resistance. 72 In Phase 2, doxorubicin was introduced to counter progressive disease. Pembrolizumab was added to doxorubicin based on emerging data from the TONIC trial suggesting that anthracyclines (and platinum) were the best immunotherapy combination. 73 In Phase 3, doxorubicin was replaced with capecitabine to counter progressive disease. Enzalutamide was added because of increased high expression of AR. Fulvestrant was added to counter persistent ER signaling. Pembrolizumab was continued since it was well-tolerated in Phase 2. Carboplatin was introduced in Phase 4 to counter disease progression. Drugs given to moderate aspects of therapy-induced toxicities included (a) denosumab to reduce risk of skeletal related events due to bone metastases; (b) pegfilgrastim to stimulate production of neutrophils; and (c) hydroxychloroquine for suspected druginduced amyopathic dermatomyositis as recommended by dermatology.

Radiology
FDG-PET/CT imaging was performed according to the standard institutional protocol, with patients fasting for 6 hours following 24 hours of rest. Prior to the examination and FDG injection, blood glucose levels were confirmed to be less than 200 mg/dL. The patient received a dose of 18 F-FDG of 370 to 555 MBq (10-15 mCi) on the basis of body weight. After an uptake period of 90 minutes, a vertex-to-mid-thigh FDG-PET/CT scan was performed using 3 min/bed position on a CTI Biograph duo PET/CT scanner (Siemens Medical Systems, Hoffman Estates, Illinois, USA) or a CTI Biograph TruePoint 40 PET/CT scanner (Siemens Medical Systems, Knoxville, Tennessee, USA). CT imaging was performed according to the standard institutional protocol from clavicles to mid-thigh on a Phillips Brilliance CT 128slice helical scanner (Philips Medical Systems, Amsterdam, NE).

Clinical and exploratory workflows
All blood and metastatic biopsy biospecimens used in this study were prospectively collected by trained study coordinators. Preservation procedures for biopsy tissue were started within two to five minutes of removal of tissue from the patient in order to preserve the molecular and architectural features of the tumor that may quickly degrade. Clinical analytics were performed in CLIA-certified, CAP-accredited laboratories. RUO exploratory analyses were performed in academic research laboratories or core facilities. Both manual and automated abstraction from the patient's medical record were used to generate the clinical metadata, including detailed information about anticancer treatments and supportive care for integration with the analytic results.
Biospecimens were tracked and managed using a custom implementation of the LabVantage laboratory information management system. The LabKey system was used to store and visualize both clinical data and analytic results. 74 The Galaxy computational workbench was used to create and run multi-step analysis workflows that process raw omics and imaging datasets. 75 The OMERO system was used to visualize multiplex imaging and electron microscopy datasets. 76 GeneTrailsâ solid tumor panel Formalin-fixed biopsy tissue was submitted to the CLIA-certified/CAP-accredited OHSU Knight Diagnostic Laboratories for targeted DNA sequencing with the clinical GeneTrails Solid Tumor Panel assay. There, genomic DNA was extracted from macro-dissected, tumor-rich regions of FFPE. Next-generation sequencing libraries were prepared using custom QIASeq chemistry (QIAGEN) with multiplexed PCR and sequenced on an Illumina NextSeq500/550. The DNA library was generated by 9,229 custom-designed primer extension assays covering 613,343 base pairs across 225 cancer-related genes (including whole exons of 199 genes and hotspot regions of 26 genes). This panel is routinely sequenced to an average read depth of >2,000, providing high sensitivity for SNVs, short insertions/deletions, and copy number alterations. All variants identified were reported out clinically. Blood collection and DNA isolation Up to 40 mL (range 6-40 mL) of blood were collected in 5 3 6-mL or 4 3 10-mL, purple-capped EDTA tubes. Consistent with published recommendations, blood plasma was isolated within 6 hours of collection by first spinning whole blood at 1000g for 10 min, separating the top plasma layer into 1 mL aliquots, then spinning those aliquots at 15,000g for 10 min, transferring the supernatant to cryovials, and storing at À80 C. 77 DNA extraction of tumor tissue from FFPE was carried out using QIAamp DNA FFPE Tissue Kit (QIAGEN). DNA was extracted from plasma using the NucleoSnap cfDNA kit (Macherey-Nagel) and from buffy coat using the DNeasy Blood & Tissue Kit (QIAGEN). DNA isolated from both FFPE samples and buffy coat were fragmented by sonication to 150 bp using a Covaris E220 prior to library preparation.
Whole exome sequencing Sequencing libraries were prepared from 100-500 ng of cell free DNA (cfDNA) or sonicated genomic DNA using the KAPA Hyper-Prep Kit (KAPA Biosystems), enriched using the SureSelectXT Target Enrichment System (Agilent Technologies) and the SureSelectXT Human All Exon V5 capture baits (Agilent Technologies). Next generation sequencing was carried out using the Illumina NextSeq500 or HiSeq 2500 platform with 2x79-144 cycles by the OHSU Massively Parallel Sequencing Shared Resource to an average depth of 100x per library replicate. For Bx3 and Bx4 only, DNA isolated from both FFPE samples and buffy coat were submitted to Tempus Labs, Inc. for whole exome sequencing (WES) with the Tempus xE assay (Tempus Labs, Inc., Chicago, IL, USA).

Dual index degenerate adaptor sequencing
Bait Design: Single nucleotide variants (SNVs) were filtered by frequency (>5% in the tumor and <2% in the matched normal) and depth (>30x in the tumor and >15x in a patient-matched matched buffy coat normal). A set of 20-40 SNVs per tumor tissue sample were then hand-selected for bait design based on high variant allele frequency (VAF) and potential clinical relevance. In total, 55 mutations were selected with varying presence across the primary tumor and all four biopsies, to monitor longitudinal blood draws for the presence of tumor-derived circulating tumor DNA (ctDNA). Biotinylated oligonucleotides (IDT, Coralville, IA, USA) matching the 120bp UCSC GRCh37/hg19 human genome reference sequence spanning each mutation were synthesized for use as baits in hybridization capture library preparation (Table S2). Oligonucleotides for mutations in INVS and LILRA3 were eventually discarded for inconsistent coverage and high error rates.
Library Preparation: Dual Index Degenerate Adaptor Sequencing (DIDA-Seq) error-correction libraries were created using the KAPA Biosystems HyperPrep kit (KAPA Biosystems) using at least 30 ng of cell-free DNA (cfDNA) as input as previously described, 14 using a single over-night capture incubation instead of two incubations. Briefly, DIDA-Seq adaptors (Table S2) were ligated to extracted cfDNA using KAPA Hyper Prep Reagents and protocol (with 16 C overnight incubation) and then PCR amplified for 8 cycles using library amplification primers (Table S2). Approximately 250 ng of amplified libraries were pooled (4-5 libraries per pool) and overnight hybridization capture and purification was carried out using the xGen Hybridization and Wash kit (IDT, Coralville, IA, USA). Libraries were then PCR amplified for 10 cycles as before. Amplified libraries were then purified 1:1 with Agencourt AMPure XP beads (Beckman-Coulter) and assessed using the 2100 Bioanalyzer system (Agilent Technologies, Santa Clara, CA, USA) High Sensitivity dsDNA kit. Samples were then sequenced on either the Illumina HiSeq 2500, paired-end 100 bp with dual 14-bp indexing cycles (high-capacity, rapid run mode) or the Illumina NextSeq 500, paired-end 75 bp with dual 14-bp indexing cycles (high-capacity, 150-cycle kit).
Low-pass whole genome sequencing Low-Pass Whole Genome Sequencing (LP-WGS) libraries were prepared with 50 ng of sonicated tumor DNA (extracted from FFPE as described above) and patient-matched buffy coat DNA using the KAPA Hyper-Prep kit (KAPA Biosystems) with Illumina single index TruSeq adaptors (IDT, Coralville, IA, USA) and sequenced on the Novaseq S4 platform (Illumina) to 0.9X mean coverage.
Whole transcriptome sequencing Library construction and sequencing: RNA was extracted from macro-dissected, tumor-rich regions of FFPE at the OHSU Knight Diagnostic Laboratories. Sequencing libraries are constructed with the TruSeq RNA Library Prep Kit, followed by sequencing on Illumina NextSeq500. A Universal Human Reference RNA (UHR; Agilent Technologies) was sequenced with each batch of samples to allow for assessment and removal of technical artifacts (due to, e.g., library preparation).

Reverse phase protein arrays
Tumor tissue collected from core needle biopsies was flash frozen in liquid nitrogen within two to five minutes of removal from the patient and stored at À80 C. Frozen tissue was then submitted to the MD Anderson Cancer Center Functional Proteomics RPPA Core for proteomic profiling of the abundances of 450 proteins and phosphoproteins with the Reverse Phase Protein Array (RPPA) assay. 21 Article ll region of interest). 83 Typical power densities for the 647 nm and the 405 nm lasers were 1-2 kW/cm 2 and 1-20 W/cm 2 , respectively. Raw STORM images were processed and reconstructed using custom Matlab scripts. 84

QUANTIFICATION AND STATISTICAL ANALYSIS Radiology
Pre-and on-treatment FDG PET/CT studies were reviewed by an expert nuclear medicine physician with analysis performed by a body imager with 15 years of experience in oncologic imaging. Target lesions were selected by having maximum standard uptake values (SUVmax) greater than normal mediastinum average (lymph nodes) and greater than background liver SUV (liver lesions) and were recorded on the pre-and on-treatment scans at the same tumor region. Image analysis was performed using syngo.via advanced visualization software (Siemens Healthcare GmbH, Erlangen, Germany) and Horos visualization software. All lesions meeting these criteria were recorded both on FDG-PET/CT and combined with long axis measurements (e.g., liver, splenic, lung lesions) and long and short-axis measures (lymph nodes) at all time points during the care of the patient. Variability in the measurement of the long axis of each lesion was estimated to be about 20%. 85 Change in tumor burden was assessed for each phase of treatment using RECIST 1.1 criteria. 8 All SUVmax measures were normalized by subtracting the mean background SUVmax from the organ of origin (e.g., mediastinum or liver). The uncertainty in SUVmax measurements was estimated to be up to 18% based on historical testretest reproducibility. 86 Whole exome sequencing Somatic mutation calling: sequence read FASTQ data files were aligned to the UCSC GRCh37/hg19 human genome build using BWA-MEM (0.7.12, GATK, Broad Institute), followed by marking duplicate reads (Mark Duplicates, GATK) and base recalibration (BQSR, GATK). 87,88 Bam files for replicate libraries were merged and somatic variants were called using MuTect2 (4.0.4.0, GATK, Broad Institute) between tumor or cfDNA and the patient's matched normal from buffy coat. 88 A panel of normal (PON) and the gnomAD (2.0.1) germline reference resource were used to filter out technical sequencing artifacts and common polymorphisms, respectively. 89 All analysis tools were run using an OHSU Galaxy instance (v17.09). 75 Phylogenetic and Clonal Analysis: Mutect2 (GATK) and mpileup (Samtools) were used to call or detect presence of variants across all samples. 88,90 Only sequence reads with base quality greater than 20 and mapping quality greater than 30 were used for mpileup. Variants with VAF lower than 5% or depth lower than 30 reads were filtered. The R package ape was used for phylogenetic analysis. 91 A binary table of variants present across all tumor samples was generated as input. Genetic distance was estimated using the dist.gene function with the pairwise method. Minimum Evolution (ME) fit with ordinary least-squares (OLS) using the FastME function was used to reconstruct the phylogeny.
Copy Number Analysis: Copy number analysis of WES data was performed with CNVkit (v0.9.4a0) using the tumor/ctDNA aligned reads (BAM) and a pooled normal reference. 92 On-and off-target read depths from each sample were median-centered log2 normalized, followed by GC bias correction and repeat masking. Tumor copy ratios were estimated by subtracting the log2-normallized depths for each bin. Corrected copy ratio profiles were segmented using circular binary segmentation (CBS). Tumor purity estimates were then used to call each segment's absolute integer copy number.
Tumor Mutational Burden (TMB) and Neoepitope Prediction: TMB was calculated based on the number of somatic nonsynonymous mutations per megabase of the targeted regions. Neoepitopes were identified across all biopsy samples using neoepiscope v.0.5.0, 93 with HLA types predicted for all samples using OptiType v1.3.3 94 and MHC binding affinities predicted using MHCflurry v2.0. 95 All potential neoepitopes that could arise from each identified mutation were subsequently filtered based on MHC binding affinity, retaining only neoepitopes with binding affinity <500 nM for at least one patient allele.

Dual index degenerate adaptor sequencing
Error-Correction, Bait Evaluation, and Variant Analysis: The DIDA-Seq computational pipeline was implemented as previously described, 14 based on substantial modification of previous work. 96 Indexing reads containing multiplexing barcodes and degenerate unique molecular identifiers (UMI) were appended to the read header of each set of paired-end sequencing reads. Next, paired-end reads were aligned using BWA-MEM 87 and UMI families were collapsed to generate consensus sequences (requiring at least three reads and 90% agreement between reads, otherwise resulting in read omission or an ''N'' at a given consensus site), which were output as a FASTQ file. These FASTQ files were then realigned using BWA-MEM, 87 3 bases from either end were replaced with ''N''s, and overlapping reads were collapsed to avoid double-counting. Final BAM files were used for downstream VAF calculation (mutant allele read count divided by total read count at each site of interest) for each mutation of interest at each timepoint. All hybrid capture baits were also evaluated using unrelated patient cfDNA samples as negative controls. We sequenced each library to an average post-error correction depth of 4,000-15,000x for each site of interest and determined the VAF in these negative controls.
Quantification and Analysis: To determine the significance of a given VAF measurement, we compared the mutation-specific VAF in the patient's plasma sample data to the VAF of the same site in a set of pooled negative controls (sequenced to an average post error-correction depth of 100,000X per site, giving an average error rate of 1 in 30,000 reads) as previously described. 97 We used a Bayesian approach to test the null hypothesis that the sample VAF and negative control VAF were generated from the same distribution using R. To do this, A p value was generated for each site and mutational group of sites ( Figure S2E)  Article ll OPEN ACCESS coefficient of the beta distributions between the VAF in the sample and VAF in the negative controls. 98 Any site with greater than 1% VAF in the negative controls was omitted from further evaluation. Data points having a p value greater than 0.05 were considered not statistically different from the negative controls, effectively determining our lower limit of detection given the individual or aggregated sites' sequencing depth at each time point.
Low-pass whole genome sequencing FASTQ files were aligned to the UCSC GRCh37/hg19 human genome build using BWA-MEM, 87 and copy number alterations were called using ichorCNA with window set to ''50000' '. 99 Whole transcriptome sequencing Gene Quantification: Transcript quantification followed the methods described by Tatlow and Piccolo. 100 The raw sequence reads were quality trimmed using Trim Galore, followed by pseudo-alignment and transcription quantification using Kallisto with GENCODE reference transcriptome (version 24). 101,102 Transcript level expression was aggregated to gene level abundance using the R package tximport yielding expression values for 60554 Ensembl genes. 103 Batch Correction: Genes were filtered based on a minimum of 3 transcripts per million (TPM) in at least 3 of 48 samples, which included 29 ER+ metastatic breast cancer samples and 19 UHR samples. The filtered gene expression matrix (16,364 genes) was batch corrected by removing unwanted variation (RUV; RUVSeq). 104,105 RUV correction uses factor analysis to identify the factors of unwanted variation observed in the UHR batch control and corrects for them across all samples. RUVSeq was applied by removing 1 factor (k) using the 5% of genes with the lowest standard deviation. In addition to intra-cohort batch correct, the patient samples were batch adjusted for analyses comparing to TCGA BRCA. 100 TCGA BRCA gene expression matrix was filtered to samples with a Luminal (A or B) molecular subtype and joined with the RUV corrected patient sample gene expression. The combined matrix was log transformed, filtered to genes with a minimum of 3 log2 TPMs in at least 3 samples, and batched corrected using SVA/ComBat with TCGA samples set as the reference. 106 Molecular Subtype Signature: The PAM50 subtype gene signature was used to classify samples into intrinsic molecular subtypes. 15 A cohort of 20 ER+ and 20 ER-samples was used as the background for classifying the patient samples' subtypes. The gene expression matrix using these 40 samples and the patient samples was mean centered and correlated (Spearman) to the pre-defined centroids from Parker et al. 15 The samples were assigned to the molecular subtype with the highest Spearman correlation to the subtype centroid.
Pathway enrichment analysis: Gene set variation analysis (GSVA) was used to estimate pathway enrichment of the MSigDB Cancer Hallmark Pathways (50 gene sets), All MSigDB Pathways ($20K gene sets), and Reactome Pathway Database ($2K gene sets). 17,18,107 GSVA used a Gaussian kernel for estimating the cumulative density function and the enrichment statistic was calculated as the difference between the maximum positive/negative random walk deviations. This analysis was applied to the RUV/ComBat adjusted log2 gene expression matrix that included both TCGA BRCA Luminal Samples and the patient samples.

Transcriptional regulator networks
Regulatory pathway and molecular interactions network: The regulatory network used to generate enrichment signatures is derived from the aggregation of publicly available molecular interactions and biological pathway databases provided by the Pathway Commons (PC) resource. 24 The aggregated data is represented in the standard Biological Pathway Exchange (BioPAX) language and provides the most complete and rich representation of the biological network models stored in PC. These complex biochemical reactions were reduced to pairwise relationships using rules to generate a Simple Interaction Format (SIF) representation of BioPAX interactions. The reduction of BioPAX interactions to the SIF allows for the representation of pairwise molecular interactions in the context of specific binary relationships. The feature space of the SIF regulatory network was restricted to primary and secondary downstream interactions for genes within Pathway Commons. The regulatory network was then reduced to edges that are associated with the binary relationship ''controls-expression-of'', defined as any reaction where the first protein controls a conversion or a template reaction that changes the expression of the second protein.
Network weight assignment: Weights are assigned to the protein-protein edges within the graph for each regulator-target pair within the regulatory network that is represented in the expression data set. These weights are derived from the integration of an F-test statistic to capture linear dependency and the Spearman rank-order correlation coefficient for a given regulator-target pair.
Regulon enrichment: This method leverages pathway information and gene expression data to produce regulon-based protein activity scores. Our method tests for positional shifts in experimental-evidence supported networks consisting of transcription factors and their downstream signaling pathways when projected onto a rank-sorted gene-expression signature. The gene-expression signature is derived by comparing all features to the median expression level of all samples considered within the data-set. After weights have been assigned to the regulatory network, the positive and negative edges of each regulator are rank ordered. The first component of the enrichment signature, the local delta concordance signature, is derived by capturing the concordance between the magnitude of the weight assigned to a particular edge and the ranked position of that edge. The features associated with activation, positive edges within the regulatory network, are monotonically ranked from most lowly to highly expressed in the restricted feature space, where the features that are repressed are ranked by a monotonically decreasing function. This component of the signature Cell Reports Medicine 3, 100525, February 15, 2022 e13 Article ll considers positive and negative edges independently, which captures support for an enrichment signature even if one of the edge groups is underrepresented in the network graph. The second component of the enrichment signature, the local enrichment signature, captures positional shifts in the local gene ranked signature and integrates those shifts with the weights assigned to overlapping features for a given regulon and the expression data set. The last component of the enrichment signature considers the entire feature space and projects the rank-sorted local signature onto this global ranked feature space. We derive our global enrichment signature from this projection for each regulator we consider. We use the median of robust quantile-transformed ranked positions as the enrichment scores for both the local enrichment and global enrichment signatures. We then integrate these three individual signatures together, which allows us to capture differences between individual regulator signatures within the context of an individual patient as well as at a cohort level.

Reverse phase protein arrays
In order to scale the reported protein expression values, the RPPA data from the patient samples was merged within the TCGA breast cancer RPPA dataset, using the replicate-based normalization (RBN) method. 108 The protein expression values were then z-scored by using the median and standard deviation, and a heat-map was generated using Rank-Sum ordering of the proteins fold change. The heat map was produced using publicly available Cluster 3.0 and TreeView software. 109,110 Pathway Analysis: All pathway predictors have been previously described. 22 Proteins used as predictors of the different pathways are listed in Table S2. To determine a pathway score, for each sample all positively associated predictors were summed minus the predictors that are negatively associated with the pathway. The total was then divided by the numbers of predictors in the pathway. To generate the pathway scores histograms, the distribution of each TCGA samples subtype was plotted and the value of the patient pre-and post-treatment sample was added to the histograms.

Intracellular signaling protein panel
Batch correction was performed using Removal of Unwanted Variation (RUV-III) using the replicate positive controls to estimate the factors associated with batch effect. 111 RUV parameters were optimized by measuring the consistency of replicate controls and careful evaluation of outliers to ensure validity of results. The GeneTrails Intracellular Signaling Protein Panel assay reports out protein expression levels from RUV-normalized data relative to tumor-type-matched reference cohorts, in this case a panel of 57 metastatic breast cancer patient samples.

Integrative analyses
Multi-omic integrated pathway analysis: CausalPath was used for integrated pathway analysis of protein, phosphoprotein, gene abundance, and transcriptional regulator activity. 25 CausalPath is a hypothesis generating tool that uses literature-grounded interactions from Pathway Commons to produce a graphical representation of causal relationships that are consistent with patterns in a multi-omic datasets. 25 This integrative approach allows for holistic evaluation of signaling networks and pathway activity across longitudinal biopsies. The CausalPath analysis used the log fold change of total and phosphoprotein (RPPA) and gene expression from Bx1 and Bx2 with the following parameters: threshold-for-data-significance = 0.3 for RNA, protein, and phosphoprotein, value-transformation = max, calculate-network-significance = true, permutations-for-significance = 10,000, color-saturation-value = 2.5, datatype-for-expressional-targets = rna and protein, show-all-genes-with-proteomic-data = true. The resulting network was pruned to include the neighborhoods encompassing MTOR, AKT, MUC1, STAT3, MYC, and E2F1 to highlight biologically interesting patterns discussed in the text. For additional depth, the difference in transcriptional regulon enrichment activity between Bx2 and Bx1 was mapped to and overlaid on the pruned CausalPath network.
Integrated Heatmap: The gene, protein, phosphoprotein abundances, and transcriptional regulon enrichment activities were integrated into a single heatmap. Each data type was independently scaled to À1 to 1 with the exception of protein/phosphoprotein, which were scaled together. Fold change of Bx2 to Bx1 was calculated for each scaled feature and represented as a heatmap grouped by pathway categories of interest.

Multiplex immunohistochemistry
Image analysis pipeline: 112 Regions of interest (ROIs) were selected from hematoxylin-stained images in ImageScope (Leica) based on histopathological assessment. ROIs were selected to capture all possible tumor area (composed of neoplastic cells and surrounding stroma), while excluding regions of adjacent normal appearing tissue, heavy RBC infiltrate, necrosis, acellular material, or areas of tissue deformation/folding which are all known to create artifactual results. The number of ROIs selected per sample was variable and sample dependent (range 2-4, average 3). ROI size was also variable (range 0.6-7.9 mm 2 , average 3.7 mm 2 ). The tissue area and cell counts from each ROI on a given tissue were summed in order to generate global immune composition. Across the cohort, the average amount of tissue assessed per sample was 11.1 mm 2 , containing an average of 68,223 nucleated cells. Digitally acquired images were registered in MatLab (MathWorks) utilizing the SURF algorithm 113 in the Computer Vision Toolbox. Nuclear segmentation and color deconvolution were performed using FIJI 81 (ImageJ, NIH). Watershed based segmentation on hematoxylin only stained tissue was used to identify single cell objects. In short, preprocessing to isolate signal and remove background was performed, then nuclear objects were identified by watershed and standard image processing (erosion, dilation, and noise removal). AEC chromogenic signal was extracted by converting images from RGB to CMYK in ImageJ 81 using the NIH plugin RGB_to_CMYK. The contrast