An omic and multidimensional spatial atlas from serial biopsies of an evolving metastatic breast cancer.

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.

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., Unfortunately, treatments deployed according to precision medicine principles do not always elicit positive responses, and durable control is achieved for only a subset of individuals with metastatic cancer. 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),, through which all data are available in standardized formats.

Results

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, 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).

Figure 1

Workflows and analytical platforms used to generate the OMS atlas

Figure 2

Timeline of clinical treatment and response metrics

(A) Treatment schedule and biopsy timing (red stars) over four phases of treatment (green, orange, blue, and pink areas). The timeline is sectioned into 28-day months. The duration and relative dose for each drug is indicated by the extent and width of a horizontal bar. Drug continuation after the end of phase 4 is indicated by a right arrow.

(B) Clinically reported serum levels of tumor protein biomarkers. CEA values were multiplied by 5 to ease visualization.

(C) RECIST 1.1 assessment of tumor response (orange stars) indicating partial response (PR), progressive disease (PD), or stable disease (SD). Shown are longitudinal tracking and variation in the longest-axis size of 16 representative metastatic lesions measured from serial CT images. Targets of metastatic biopsies are bolded and marked with stars. Circles represent FDG-PET imaging results, colored and centered on the lines of their corresponding lesion at interpolated lesion sizes. The diameter of each circle is proportional to the background-normalized maximum standardized uptake value (SUVmax).

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. 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 (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), 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: B×3 = 6.8×, B×4 = 7.2×; YES1: B×3 = 2.0×, B×4 = 4.0×). Importantly, both biopsies were acquired after treatment with the TYMS inhibitor capecitabine (Figures 2A, 2C, and S1G).

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.

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. 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). Human leukocyte 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).,

ctDNA increased during progression and radiotherapy

Dual index degenerate adaptor sequencing (DIDA-seq) 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; Figure 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 whole-transcriptome sequencing (RNA sequencing [RNA-seq]). Classification using the PAM50 subtype gene signature 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 (TCGA-BRCA) and gene set variation analysis (GSVA) of enriched Molecular Signatures Database (MSigDB) cancer hallmarks relative to TCGA-BRCA luminal samples., 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.

Protein and phosphoprotein abundances were measured in Bx1, Bx2, and Bx4 using reverse-phase protein arrays (RPPAs), and proteomic pathway signatures were compared with TCGA-BRCA (Figures 3D and 3E).19, 20, 21, 22 Aspects of hormone signaling varied across biopsies. Estrogen receptor (ER) protein levels from RPPA and clinical immunohistochemistry (IHC) were high in all three biopsies (Table S1). The protein pathways “hormone signaling” and “hormone receptor” were higher in Bx2 (Figure 3D), whereas the GSVA RNA hallmarks “estrogen early” and “estrogen late” (Figure 3C) showed little change, an intriguing finding because protein levels of the hormone-regulated transcription factors ER, GATA3, and adrenergic receptor (AR) were increased in Bx2 relative to Bx1 after phase 1 treatments (Figure 3E). Bx4, taken after phase 4 treatment without hormone suppression, showed continued elevation of the “hormone receptor” pathway and ER and AR protein levels relative to Bx1. However, GATA3 protein levels, the “hormone signaling” protein pathway, and the “estrogen early” and “estrogen late” GSVA hallmarks were downregulated.

PI3K/AKT/mTOR pathway signaling from RPPA was generally similar across all biopsies, even though Bx2 was collected after treatment with the mTORC1 inhibitor everolimus and contained the hotspot mutation PIK3CA p.E542K (Figure 3D). Individual protein levels within these pathways varied but did not result in changes in overall signaling. For example, Bx2 showed decreased mTORC1 complex activity based on decreased S6 phosphorylation at S235/236 and S240/244 (0.7× and 0.1× versus Bx1) but increased activity downstream of mTORC2, including increased phosphorylation of AKT (S473: 2.7× versus Bx1) and its substrates GSK3A/B (S21/S9: 1.7× versus Bx1), TSC2 (T1462: 1.4× versus Bx1), and MDM2 (S166: 1.8× versus Bx1; Figure 3E). Likewise, Bx4 showed increased phosphorylation of AKT at S473 (2.7× versus Bx1, 1.0× versus Bx2) and NDRG1 (T346: 1.8× versus Bx1, 1.6× versus Bx2) but without an accompanying increase in AKT or mTORC1 substrate phosphorylation. We also used the clinical Intracellular Signaling Protein Panel (ISPP) to quantitate phosphoproteins and total proteins in Bx2–Bx4 (Figure 3F; Table S2). ISPP showed that Bx2 had the highest AKT phosphorylation and lowest S6 phosphorylation relative to the other biopsies, consistent with RPPA results, whereas p-ERK, p-cRAF, and p-MEK were elevated in Bx3 and increased relative to other biopsies.

A transcriptional regulator analysis using a molecular interactions network derived from Pathway Commons was used to infer regulator protein activity from the gene expression data. Integrative analysis of the longitudinal changes in proteomics, gene expression, and transcriptional regulator scores between Bx1 and Bx2 was also performed using CausalPath (Figure S3D). These analyses showed strong inhibition of mTOR regulator activity (Bx2 5.1× > Bx1). Activities of multiple JAK-STAT family proteins were increased, including JAK2 (Bx2 1.8× > Bx1), phospho-STAT3 (Y705: Bx2 1.6× > Bx1), and STAT5 (Bx2 3.2× > Bx1), which, together with the oncoprotein mucin 1 (MUC1; protein Bx2 27.1× > Bx1; regulator Bx2 +3.15 versus Bx1) constitute a known feedforward loop whereby MUC1 binds STAT3 to facilitate its phosphorylation by JAK1. These observations are reinforced by elevation in “IL6/JAK/STAT3” and “IL2/STAT5” signaling from GSVA (Figure 3C). This analysis also highlighted decreases in MYC and E2F regulator activity and E2F1 total protein, consistent with decreased enrichment of “MYC targets” and “E2F1 targets” in GSVA. These analyses provide a view of the interaction dynamics of cell cycle control networks, with decreases in the expression of cell cycle progression genes (CCNB1, CDK4, CDK1, CCNE2, CCND3, and PLK1) balanced by a sharp decrease in cell cycle inhibitor genes (CDKN1A, CDKN1B, and CDKN2A), leaving RB1 phosphorylation unchanged in Bx2.

Tumor immune microenvironment evolution and barriers to T cell activation

Changes in composition and functionality of lymphoid and myeloid lineage immune cells were assessed in the primary tumor and Bx1–Bx4 using multiplex IHC (mIHC),27, 28, 29 noting possible discordance in Bx3 because of its bone origin (Figures 4 and S4B–S4E; Table S2). Total immune cell infiltration, as indicated by the percentage of CD45+ cells, was lowest in Bx2 (2.1%) but comparable between the PT (7.2%), Bx1 (9.3%), and Bx4 (10.0%; Figures 4A and 4B). Myelomonocytic cells (macrophages and monocytes) comprised the dominant CD45+ leukocyte lineage subgroup in the PT (65.8%; Figure 4C, green and brown), Bx1 (48.3%), and Bx2 (82.4%) and were reduced in Bx4 (6.1%). Analysis of the myeloid lineage revealed that the fraction of immature dendritic cells was higher in Bx1 (0.2%; Figure 4D) than in Bx2 (0.05%), whereas the proportions of CD163+ and CD163− macrophages and monocytes were higher in Bx2 (51.2%) than in Bx1 (18.9%), with the largest increase in CD163+ macrophages (Figure 4D). CD163 positivity is associated with differentiation of myelomonocytic cells toward an alternatively activated or “M2”-type state, which is considered to be pro-tumorigenic within solid tumors., CD163 expression on monocytes and macrophages is induced by interleukin-10 (IL-10) and glucocorticoids and repressed by lipopolysaccharides, tumor necrosis factor alpha (TNF-α), and interferon γ (IFNγ) and is concordant with upregulation of IL-containing GSVA gene sets in Bx2 (Figure 3C; Table S2). The dominance of macrophages and monocytes and relative lack of T cells in the PT, Bx1, and Bx2 was in stark contrast to Bx4, which had many more T cells than macrophages and monocytes (PT: 65.8% macrophages/monocytes, 10.9% T cells; Bx1: 48.3%, 20.6%; Bx2: 82.4%, 5.3%; Bx4: 6.1%, 33.1%; Figures 4C and 4E, orange).

Figure 4

Monitoring response to therapy with deep in situ immune phenotyping by mIHC

(A) Primary tumor (PT) and Bx1–Bx4 were subjected to multiplex immunohistochemistry (mIHC) analyses measuring immune (CD45+) and epithelial (PanCK+) cells in tumor compartments as a percentage of total nucleated cells.

(B) Representation of tissue composition, showing density (number of cells per square millimeter of tissue analyzed) of PanCK+ (cytokeratin), CD45+, and PanCK− CD45− (other) nucleated cells.

(C) Immune composition of seven major leukocyte lineages, as a percentage of total CD45+ cells.

(D) Deeper auditing of leukocyte lineages in Bx1 and Bx2, measuring 12 immune cell populations and functional states.

(E) CD3+ T cell proportions of total CD45+ cell populations (orange, left), and CD4+ (blue) and CD8+ T cells (periwinkle) proportions within CD45+CD3+ T cells (right).

(F) PD-1+ cells as a percentage of total CD3+T cells in the CD3+CD4+ (top) and CD3+CD8+ (bottom) T cell populations.

(G) Differentiation state of CD3+CD4+ T cells, reflected by regulatory T (Treg), Th1, and Th2, Th17, and Th0/γδ subsets (left) and CD3+CD8+ T cells, as reflected by expression of PD-1 and EOMES.

(H) Differentiation state of CD3+CD4+ T cells reflected by Treg, Th17, Th1, Th2, and Th0/γδ subsets in Bx1 and Bx2.

See also Figure S4 and Table S2.

Analyses of T cell subsets and functionality showed that only a small fraction of CD3+CD4+ and CD3+CD8+ T cells in the PT, Bx1, or Bx4 expressed the programmed cell death-1 (PD-1) protein, which is typically expressed on activated T cells following T cell priming or persistent antigen exposure (PT: 1.5% CD4+ T cells, 4.1% CD8+ T cells; Bx1: 1.3%, 1.6%; Bx4: 0%, 0.7%; compared with Bx2: 10.1%, 30.5%; Bx3: 11.3%, 6.3%; Figure 4F). However, T cell status was markedly altered in Bx2 (Figure 4C; Table S2). Notably, although T cells were least abundant in Bx2 compared with Bx1 and Bx4, the largest fraction of CD3+CD4+ and CD3+CD8+ T cells expressing PD-1 was observed in Bx2 (Bx2: 10.1% CD4+ T cells, 30.5% CD8+ T cells; compared with Bx1: 1.3%, 1.6%; Bx3: 11.3%, 6.3%; Bx4: 0%, 0.7%; Figure 4F), coincident with relatively reduced FoxP3+CD4+ regulatory T (Treg) cells (Bx1: 4.0%, Bx2: 1.1%, Bx3: 9.3%, Bx4: 0.5%) and expanded Th17 CD4+ T cells (Bx1: 2.5%, Bx2: 26.1%, Bx3: 0%, Bx4: 0%; Figures 4G and 4H).

Analyses of PD-1 and eomesodermin (EOMES) expression showed that the PT contained predominately PD-1−EOMES− (71.5%) and PD-1+EOMES− CD8+ T cells (27.2%; Figure 4G), likely reflecting naive and early effector subsets, respectively. Evolution of CD8+ T cells in Bx1, Bx2, and Bx4 indicated progressive loss of late effector PD-1−EOMES+ (61.5%, 19.0%, and 11.0%) and exhausted PD-1+EOMES+ subsets (20.3%, 23.8%, and 1.7%), with replacement by likely naive PD-1−EOMES− CD8+ T cells in Bx1 and Bx2 (16.8%, 47.6%; Figure 4G).

The Bx3 bone metastasis differed from the PT and liver metastases and had the highest percentage of CD45+ leukocytes (18.0%; Figures 4A and 4B) with comparatively high percentages of granulocytes (16.7%; Figure 4C), dendritic cells (22.6%), and CD20+ B cells (6.0%). However, like Bx4, Bx3 contained a prominent granulocyte infiltrate that most likely was predominantly neutrophils (Bx4: 2.0%, Bx3: 16.7%, compare with Bx1: 0.2%, Bx2: 0.1%). Neutrophils can exert significant pro-metastatic activities, including suppressive effects on T cells, and are associated with poor prognosis in many solid tumors, including breast cancer.34, 35, 36, 37, 38

Tumor and stromal interactions defined using cyclic immunofluorescence (CycIF) and focused ion-beam-scanning electron microscopy (FIB-SEM)

Tumor and stromal compositions and organizations of Bx1–Bx4 plus control biospecimens were assessed using CycIF (Figure 5A)., This joint analysis revealed 17 tumor and stromal clusters (Figures 5B and S5A–S5E). Three of the stromal clusters (11, 13, and 14) and five of the tumor clusters (0, 1, 5, 7, and 9) comprised major subpopulations in Bx2–Bx4. The three stromal clusters were identified as fibroblast-like cells that differed in levels of vimentin (VIM; cluster 11: 3.2× other clustered means, 13: 2.6×, 14: 4.9×). Endothelial cells (CD31) and macrophages (CD68) were excluded from cluster analysis because of loss of the normal breast and tonsil tissues needed for normalization during staining; their presence was confirmed using manual gating (Figure S5F). All tumor clusters expressed CK7/CK19 but different levels of ER, EGFR, and CK8 (Figure 5B). An additional proliferative cluster, cluster 16, appeared in Bx3 and Bx4, comprised of tumor and stromal cells expressing high levels of Ki67 (31.3×) and/or PCNA (2.4×).

Figure 5

Monitoring tumor and stromal responses to therapy using CycIF and FIB-SEM

(A) Example images of antibody staining overlaid with segmentation borders, colored by cell type. Scale bar, 50 μm.

(B) Heatmap of mean Z-scored intensity of unsupervised Leiden clustering (resolution, 0.45) on single-cell mean intensity of biopsies and control tissues and cell lines, with annotations on the right. Lum, luminal; Mes, mesenchymal; Fibro, fibroblast. Colored row labels indicate which biopsy was most dominant for each cluster: Bx2 (orange), Bx3 (blue), or Bx4 (pink). Cluster 16 is evenly split between Bx3 and Bx4.

(C) Single-cell mean intensity distributions of ER and PCNA staining of cells 0–25, 25–50, and 50–75 μm from positive collagen staining. Asterisks indicate significant (p < 0.001) differences in mean intensity between distances (ANOVA).

(D) Two views of reconstructed 3D FIB-SEM data from Bx1 showing the relationship between cancer cells (red and pink), stromal cells (blue and turquoise), and collagen (green). A full-volume view (left) shows nanoscale cell-cell interactions of stromal cells surrounding a tumor nest (collagen is not rendered in this image), whereas the close-up view (right) shows a fibroblast-like cell interposed between the tumor and collagen. Scale bars, 5 μm.

See also Figure S5 and Table S5.

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 intensity at 0–25 μm from collagen I = 780, 50–75 μm = 463, p < 0.001; Bx3: 0–25 μm = 1,058, 50–75 μm = 600, p < 0.001; Bx4: 0–25 μm = 1687, 50–75 μm = 1,105, p < 0.001. Bx2: mean PCNA intensity at 0–25 μm = 745, 50–75 μm = 218, p = 0.17; Bx3: 0–25 μm = 948, 50–75 μm = 567, p < 0.001; Bx4: 0–25 μm = 713, 50–75 μm = 406, p < 0.001; Figures 5C and S5G); p values describe differences in mean intensities between distances (ANOVA).

Tumor-tumor and tumor-stromal interactions in Bx1, Bx2, and Bx4 were explored at ∼4-nm resolution using FIB-SEM. 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).

Figure 6

Inter- and intracellular compositions and interactions revealed using FIB-SEM

(A) 2D SEM image from Bx4 showing the relationship between tumor cell nests and stromal collagen, along with a high density of extracted lysosomes. Scale bar, 10 μm. The selected insets show these features at high magnification. Scale bars, 3 μm.

(B) A side view of an elongated tumor cell from 3D FIB-SEM of Bx2 showing FLPs (red) and alignment of the internal mitochondria (fuchsia). Scale bar, 1 μm.

(C) Additional cells from Bx2 (the same red cell as in B) showing paddle-shaped lamellipodia (green cell) and long FLPs (red and blue cells) extending into the stroma and interacting with neighboring cells. Scale bar, 500 nm.

(D) Reconstructed 3D FIB-SEM data from Bx1 showing FLPs selectively extending toward neighboring cells and extracellular debris. Scale bar, 1 μm.

(E and F) Additional detail from Bx2 (E) and Bx1 (F) of the nuclear invaginations (blue), showing the organization of mitochondria (fuchsia) and macropinosomes (yellow) with respect to nuclear folds. Scale bars, 1 μm.

(G) 3D FIB-SEM volume of Bx2 showing large electron-dense lysosomal granules (green) dispersed between macropinosomes (red). Scale bar, 900 nm.

(H) Qualitative summary of ultrastructural feature prevalence within each biopsy. Bx4 scoring of lamellipodia is not available.

See also Figure S6 and Videos S1, S2, S3, and S4.

Video S1. Utrastructure of Bx1 from 3D FIB-SEM, related to Figure 6

The 3D FIB-SEM volume collected at 4 nm/voxel resolution from Bx1 showing ultrastructural features at the nanoscale. Individual cell contours are rendered and illustrate the cell-cell and cell-stromal interactions. Aberrant nuclear morphology, clustered macropinosomes, organized mitochondria, and the presence of lysosomes are all observed in the 25 × 20 × 6 μm3 volume.

Video S2. Tumor stromal interaction from 3D FIB-SEM, related to Figure 6

A larger 3D FIB-SEM volume (60 x 40 x 18 μm3) collected at 10 nm/voxel resolution of Bx1 shows tumor nest interaction with the fibroblasts in the stroma. The fibroblasts and stromal cells in blue are interacting with the red cancer cells and wrapping themselves around the nest to form a barrier. In this case, the fibroblasts closest to the nest are observed to be blebbing. In addition, the green collagen bundles are entwined with the stromal and tumor cells.

Video S3. Utrastructure of Bx2 from 3D FIB-SEM, related to Figure 6

The 3D FIB-SEM volume collected at 4 nm/voxel resolution from Bx2 shows ultrastructural features at the nanoscale in a 25 x 20 x 10 μm3 volume. Similar to the Bx1 volume, aberrant nuclear morphology, clustered macropinosomes, organized mitochondria, and the presence of lysosomes are all observed. However, the cell-cell interactions are remarkable. The center tumor cell squeezing between the surrounding cells is observed to have micrometer-long protrusions, while its neighbors have clear lamellipodia.

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. Cultured SKBR3 breast cancer cells exhibit similar FLPs, and dynamic in vitro images acquired using stochastic optical reconstruction microscopy (STORM) 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)., (4) A high prevalence of densely stained vesicles that appear to be lysosomes (Videos S1 and S3; Figures 6A and 6G).

Video S4. SKBR3 cell response to epidermal growth factor (EGF) from STORM, related to Figure 6

Live cell stochastic optical reconstruction microscopy (STORM) imaging of HER2 in a SKBR3 (HER2+ breast cancer) cell immediately after EGF treatment (at 10 ng/mL), highlighting the dynamics of HER2-enriched FLPs during a ∼4-minute period. Note that the tips of most FLPs remained at their original locations and the cell body extended significantly in the form of lamellipodia toward the tips of the FLPs.

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 protein19, 20, 21, profiles and spatially defined analyses, including mIHC,27, 28, 29 CycIF,, and 2D and 3D electron microscopy. 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.

Figure 7

Mechanisms of therapeutic resistance and response suggested by RPPA

(A) Phosphorylation and inferred activation of the PI3K/AKT/mTOR pathway affected by everolimus in Bx2. Decreased phosphorylation of S6 downstream of mTORC1 likely resulted from everolimus inhibition, but increased phosphorylation of proteins downstream of PI3K and AKT, possibly through mutant PI3K E542K activity and/or feedback signaling to mTORC2, may have provided continued oncogenic signaling in the presence of this drug. Proteins are noted as increased activating phosphorylation (>1.4×, red), increased inhibitory phosphorylation (>1.4×, pink), decreased activating phosphorylation (<0.7×, green), or unchanged/unknown phosphorylation (yellow). Changes in phosphorylation in Bx2 versus Bx1: PDK1 = 1.45×, AKT T308 = 1.20×, AKT S473 = 2.69×, TSC2 = 1.43×, GSK3A/B = 1.71×, MDM2 = 1.75×, p70S6K1 = 0.92×, 4EBP1 = 1.37×, S6 S235/236 = 0.69×, S240/244 = 0.14×.

(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.98× versus pre-treatment Bx1). This is possibly due to degradation of the CDK2 inhibitor p21 (0.50× 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.67× 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.

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. 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,, 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., 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. 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, 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.0× versus Bx1) because this modification promotes cell cycle progression and should have been decreased by palbociclib. 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, 54, 55 First, protein levels of the CDK2 inhibitor p21 were decreased 2× from Bx1 to Bx2 (Figure 3E), possibly because of activated PI3K/AKT signaling, which maintains low p21 levels in CDK4/6 inhibitor-resistant cells. Second, tumor cells with high CCND1 and activated PI3K can adapt to palbociclib via non-canonical binding of CDK2 to CCND1, and Bx2 had higher CCND1 protein levels (1.7× versus Bx1) and PI3K/AKT signaling than Bx1. The CDK4/6 inhibitor abemaciclib has a broader spectrum of activity that includes CDK2 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) at the time of Bx3. Analysis of Bx3 revealed a focal amplification of TYMS and YES1. TYMS is inhibited by capecitabine, and its overexpression confers resistance to capecitabine. 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. 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 IFNγ, 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. 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.

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 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 receptor-dense, 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. 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., 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.,, 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. (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. Lysosomotropic sequestration has been implicated as a mechanism of resistance to CDK4/6 inhibitors, 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. Recent studies suggest that lysosomotropism-mediated doxorubicin resistance can be countered by the β-AR antagonist propranolol, which acts through a β-AR-independent mechanism to increase cytoplasmic doxorubicin concentrations and decrease lysosomal accumulation.

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. 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.

STAR★Methods

Key resources table

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.

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).

Cell lines

Cell line SKBR3 was cultured in McCoy's 5A (modified) Medium supplemented with 10% Fetal Bovine Serum (FBS). For cell lines used as controls for the Intracellular Signaling Protein Panel, MCF7, and MDAMB468 were acquired from the Characterized Cell Line Core (CCLC) while MDAMB231, BT474, and HCC1954 were acquired from the American Type Culture Collection (ATCC). MCF7, MDAMB468, and MDAMB231 were cultured in Dulbecco's Modified Eagle's Medium (DMEM) with 10% FBS. BT474 was cultured in ATCC Hybri-Care Medium with 10% FBS. HCC1954 was cultured in Roswell Park Memorial Institute 1640 (RPMI1640) with 10% FBS. For cell lines used as controls during CycIF, HCC1143, T47D, BT474, AU565, and MDAMB436 were acquired from the ATCC while HCC3153 was acquired from UT-Southwestern. MDAMB436 was cultured in DMEM with 10% FBS. BT474, 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. 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. 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 drug-induced 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 18F-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. The Galaxy computational workbench was used to create and run multi-step analysis workflows that process raw omics and imaging datasets. The OMERO system was used to visualize multiplex imaging and electron microscopy datasets.

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 × 6-mL or 4 × 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. 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, 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.,

Intracellular signaling protein panel

The GeneTrails® Intracellular Signaling Protein Panel is an assay validated for clinical use within the Knight Diagnostic Laboratories at OHSU. It utilizes the NanoString nCounter Vantage 3D Solid Tumor Panel for FFPE, including 23 antibodies, with 12 targeting phosphorylated proteins, specifically designed to interrogate the MAPK and PI3K/mTOR signaling pathways (Table S2). This multiplexed panel allows for the simultaneous quantification of multiple proteins from a single section of FFPE tissue. Four micrometer FFPE sections from six control cancer cell lines [MCF7, MDA-MB-468, MDA-MB-231, BT474, HCC1954, and human Epidermal Growth Factor (Sigma-Aldrich) treated MDA-MB-468] and tumor biopsy tissue were deparaffinized in Xylene, 100% ethanol, and 95% ethanol sequentially and subjected to citrate-based antigen retrieval (pH6.0) in a high-pressure cooker. Samples were blocked with blocking buffer (Buffer W) for 1 hour and incubated overnight at 4°C with the cocktail of target-specific oligonucleotide-tagged antibodies. After washing with TBS-T buffer, the oligo-tags were released by UV light (3 minutes on a UV lightbox) and hybridized overnight with NanoString’s four color-barcoded probe tagsets. Hybridized codesets (oligo target + probe tagset) were bound and immobilized to the nCounter cartridge during purification process. Target protein expression levels were semi-quantitatively measured with codeset molecule counts using the NanoString nCounter FLEX analysis system. Six cancer cell lines were selected as positive controls and were included on every run to assess antibody and control performance and to correct for batch effects.

Multiplex immunohistochemistry

Immunohistochemical staining: Glass mounted FFPE tissue sections (5 μm) were baked at 60°C for 60 minutes, deparaffinized with xylene, and rehydrated in serially graded alcohols, then place in distilled water. Slides were stained with hematoxylin (Dako) for 1 minute, mounted with 1x Tris Buffered Saline-Tween (TBST) buffer (Boston Bioproducts), coverslipped with Signature Series Cover Glass (Thermo Scientific), and subjected to whole slide digital scanning at 20x magnification using an Aperio ImageScope (Leica Biosystems). Slides were de-coverslipped with 1 min of agitation in TBST and subjected to heat-mediated antigen retrieval in 1x pH 6.0 citrate buffer (Biogenex Laboratories) for 20 min at 95°C, followed by blocking of endogenous peroxidase activity with Peroxidase and Alkaline Phosphatase Blocking Reagent (Dako, per manufacturer’s instructions). Slides were then subjected to 12 cycles of multiplex immunohistochemistry (mIHC), each cycle consisted of either 1 or 2 rounds of IHC. Each round consisted of application of primary antibody, HRP-linked secondary antibody (Histofine Simple Stain Max PO, Nacalai USA, 414311F, 414144F, 414134F), and HRP-mediated development of AEC chromogen (Vector Laboratories, SK4200), followed by whole slide scanning. Citrate antigen retrieval was used between cycles to remove primary antibodies, and HRP inactivation was used between rounds (Dako, S2003, per manufacturer’s instructions) to eliminate HRP carry-over as described previously., Several antibody panels (and variations thereof) were utilized for the current study (Table S2). Each mIHC antibody panel required one FFPE tissue specimen, thus, in instances where more than one antibody panel was used on a single timepoint (e.g., PT, Bx1, and Bx2), serial sections were used. In all other instances (e.g., Bx3 and Bx4), only one FFPE tissue section was used for one antibody panel (Discovery). Where IHC and chromogenic staining did not pass QC, they were not included in analysis: e.g., PD-L1 and CSF1R on the myeloid panel, and CD68 and ICOS on the functional panel(s). Several antibodies were not common across all or some panels, thus not included in results: IDO on functional panel (Bx1), Tryptase on myeloid panel (Primary, Bx1 and Bx2), RORyT and GATA3 on the lymphoid panel (primary, Bx1, and Bx2), and CCR2, HLA class-I, CD169, CD11b, and CD11c on the discovery panel (23 antibodies) (Bx3 and Bx4).

Cyclic immunofluorescence

Immunofluorescence analysis of tumor tissue: FFPE biopsy tissues and control biospecimens prepared from normal breast tissue, tonsil, and six cell lines representing basal-like (HCC1143, HCC3153), claudin-low (MDAMB436), luminal (T47D), and HER2 positive (BT474, AU565) breast cancers were sectioned at 4 μm and mounted on adhesive microscope slides (Mercedes Scientific). The slides were baked overnight in an oven at 55°C (Robbin Scientific, Model 1000) and an additional 30 minutes at 65°C (Clinical Scientific Equipment, NO. 100). Tissues were deparaffinized with xylene and rehydrated with graded ethanol baths. Two step antigen retrieval was performed in the Decloaking Chamber (Biocare Medical) using the default settings. After completion of the first step in 10 mM citrate buffer pH 6 (Sigma-Aldrich) in the chamber, slides were further incubated in pre-boiled Target Retrieval Solution, pH 9 (Agilent) for 15 minutes. Slides were then washed in two brief changes of deionized water (diH2O) for ∼2 seconds and once for 5 minutes in 1x phosphate buffered saline (PBS), pH 7.4 (Fisher Scientific). Sections were blocked in 10% normal goat serum (NGS, Vector Laboratories), 1% bovine serum albumin (BSA, Sigma-Aldrich) in PBS for 30 minutes at 20°C in a humid chamber, followed by PBS washes. Direct labeled primary antibodies (Table S2, key resources table) were diluted in 5% NGS, 1% BSA in 1x PBS and applied overnight at 4°C in a humid chamber, covered with plastic coverslips (IHC World). Following overnight incubation, tissues were washed 3 × 10 min in 1x PBS. Rectangular Cover Glass (Corning) were mounted in SlowFade Gold Antifade Mountant plus DAPI mounting media (Thermo Fisher Scientific).

Fluorescence Microscopy: Fluorescently stained slides were scanned on the Zeiss AxioScan.Z1 (Zeiss, Germany) with a Colibri 7 light source (Zeiss). The filter cubes used for image collection were DAPI (Zeiss 96 HE), Alexa Fluor 488 (AF488, Zeiss 38 HE), AF555 (Zeiss 43 HE), AF647 (Zeiss 50), and AF750 (Chroma 49007 ET Cy7). The exposure time was determined individually for each slide and stain to ensure good dynamic range but not saturation. Full tissue scans were taken with the 20x objective (Plan-Apochromat 0.8NA WD = 0.55, Zeiss), and stitching was performed in Zen Blue image acquisition software (Zeiss).

Quenching Fluorescence Signal: After successful scanning, slides were soaked in 1x PBS for 10–30 minutes in a glass Coplin jar, waiting until the glass coverslip slid off without agitation. Quenching solution containing 20 mM sodium hydroxide (NaOH) and 3% hydrogen peroxide (H2O2) in 1x PBS was freshly prepared from stock solutions of 5 M NaOH and 30% H2O2, and each slide placed in 10 ml quenching solution. Slides were quenched under incandescent light, for 30 minutes for FFPE tissue slides. Slides were then removed from the chamber with forceps and washed times for 2 min in 1x PBS. The next round of primary antibodies was applied, diluted in blocking buffer as previously described, and imaging and quenching were repeated over ten rounds for FFPE tissue slides.

Scanning electron microscopy

Sample Fixation:, Tumor tissue for scanning electron microscopy (SEM) was collected at the time of biopsy and placed into SEM-specific fixative (2.5% paraformaldehyde, 2.5% glutaraldehyde in 0.1M sodium cacodylate buffer) as rapidly as possible to preserve tissue ultrastructure. Tissues were then stored in fixative at 4°C indefinitely until processing could take place. No tissue was available for SEM from Bx3.

Sample Preparation:, Tissue samples were prepared for SEM by post-fixation heavy metal infiltration followed by epoxy-resin embedding with the EMbed 812 Embedding Kit. Heavy metal staining using osmium tetroxide, uranyl acetate, and lead aspartate provided contrast for imaging by dissociating the metals and allowing them to bind to lipids and proteins within cellular membranes and organelles. After staining and resin embedding, polymerized blocks were mounted directly to SEM pin-style stubs (Ted Pella, Inc.) and trimmed to create a flat surface using a Ultramicrotome Leica EM UC7 (Leica Microsystems) equipped with trim 90 diamond knives (DiATOME). Mounted blocks were conductively coated with 8-nm carbon using an EM ACE600 High Vacuum Sputter Coater (Leica Microsystems).

Imaging:, Two-dimensional large-format SEM maps were collected on trimmed block faces using a Helios NanoLab G3 DualBeam™ (FEI) focused ion-beam-scanning electron microscope (FIB-SEM) equipped with the Thermo Scientific Maps Software package. Using this software for automation, hundreds of tiled images were collected over the entire block surface and stitched together, creating a pyramidal viewing architecture that provides observations starting at the millimeter-scale and zooms all the way down to 4-nm/pixel spatial resolution. Imaging conditions were 3 keV, 200–400 pA, 4-mm working distance, and 3 μs dwell time using the concentric backscatter detector (CBS). A custom script converts these large maps from TIFF into OME-TIFF format (maps2ometiff, part of the ometiff_converters library) for web-based viewing and sharing via OMERO.

Regions of interest for three-dimensional electron microscopy (3DEM) were selected from the high-resolution maps. Three separate 3DEM datasets collected using FIB-SEM technology were generated using vendor-specific automated serial-sectioning software: two high-resolution, small volumes (4-nm/voxel, 25 × 20 × 6–10 μm3) on each respective biopsy and one lower resolution, larger volume (10-nm/voxel, 48 × 48 × 17 μm3). The high-resolution image stacks were collected using the aforementioned Helios FIB-SEM with the same electron beam conditions and the In-Column Detector (ICD). The large volume was collected from the pre-treatment biopsy using a Gemini 550 Crossbeam FIB-SEM (ZEISS International), using 1.5 keV, 1.0 nA, 5-mm working distance, 1.6 μs dwell time, and the Energy-Selective Backscatter (EsB) detector.

Segmentation: 2D SEM maps were manually reviewed using ImageJ and IrfanView. Segmentation of image stacks was performed manually with the assistance of a CloudFactory managed workforce using Microscopy Image Browser. Deep learning models were utilized for nucleus and nucleoli segmentation on the high-resolution image stacks (em_segmentation; 10.1101/2021.05.27.446019). 3D reconstruction and movies were created using Amira Software and Dragonfly.

Stochastic optical reconstruction microscopy

Alexa Fluor 647 conjugated trastuzumab was prepared using Alexa Fluor 647 NHS Ester (Thermo Fisher Scientific) and purified according to manufacturer recommended procedures; the final dye to antibody conjugation ratio was measured to be around 2:1 using a UV-Vis spectrometer.

SKBR3 cells were cultured in McCoy’s 5A (Modified) Medium (Thermo Fisher Scientific) supplemented with 10% Fetal Bovine Serum (Thermo Fisher Scientific). For Stochastic Optical Reconstruction Microscopy (STORM) experiments, the cells were plated in Nunc Lab-Tek II Chambered Coverglass (Thermo Fisher Scientific) for 36 to 48 hours before labeling and imaging. To prepare for imaging, the cells were first serum starved overnight (∼16 h); on the day of imaging, the cells were treated with 100 nM Alexa Fluor 647 conjugated trastuzumab for ∼15 min, washed with pre-warmed blank medium, and placed on the microscope stage for imaging. Next, fresh STORM imaging buffer was added at 1:1,000 v/v dilution to the medium; the buffer is Phosphate Buffered Saline (PBS) supplemented with 0.5 mg/mL glucose oxidase (Sigma-Aldrich), 40 μg/mL catalase (Sigma-Aldrich), and 10% D-Glucose (Fisher Scientific); this was followed by addition of 10 mM (final concentration) cysteamine (Sigma-Aldrich). The sample was then explored at low 647 nm laser power (∼100 W/cm2; this avoids unnecessary loss of AF647 due to photobleaching) to identify regions of interest. Human Epidermal Growth Factor (EGF; Cell Signaling Technology) was diluted from a 1 mg/mL stock in PBS to 10 mg/mL and then added to the cell culture at 1:100 v/v dilution to yield a final concentration of 10 ng/mL. Image acquisition was initiated right after adding EGF, as described below. Throughout the imaging process, the cells were kept in an on-stage incubator (TokaiHit) at 37°C with 5% CO2.

The STORM microscope setup was the same as described previously. A custom single-molecule fluorescence imaging setup was built on a Nikon Ti-U microscope frame, with other essential components including an objective lens with high numerical aperture (Nikon 60x oil, TIRF, NA = 1.49), a 647 nm laser (Coherent OBIS, max output = 140 mW; for exciting and converting AF647 into a dark state), a 405 nm laser (Coherent CUBE; for converting AF647 to fluorescent on-state), and an EM-CCD (Evolve 512 Delta, Photometrics), as well as other components including dichroic mirrors and emission filters. Image acquisition was performed using microManager with an EM-CCD gain setting typically set at 300 and the frame acquisition time 8 ms (possible by selecting a small region of interest). Typical power densities for the 647 nm and the 405 nm lasers were 1–2 kW/cm2 and 1–20 W/cm2, respectively. Raw STORM images were processed and reconstructed using custom Matlab scripts.

Quantification and statistical analysis

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%. Change in tumor burden was assessed for each phase of treatment using RECIST 1.1 criteria. 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 test-retest reproducibility.

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)., 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. 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. All analysis tools were run using an OHSU Galaxy instance (v17.09).

Phylogenetic and Clonal Analysis: Mutect2 (GATK) and mpileup (Samtools) were used to call or detect presence of variants across all samples., 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. 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. 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, with HLA types predicted for all samples using OptiType v1.3.3 and MHC binding affinities predicted using MHCflurry v2.0. 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.

Error-Correction, Bait Evaluation, and Variant Analysis: The DIDA-Seq computational pipeline was implemented as previously described, based on substantial modification of previous work. 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 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, 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. 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) using the overlap coefficient of the beta distributions between the VAF in the sample and VAF in the negative controls. 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.

FASTQ files were aligned to the UCSC GRCh37/hg19 human genome build using BWA-MEM, and copy number alterations were called using ichorCNA with window set to “50000”.

Gene Quantification: Transcript quantification followed the methods described by Tatlow and Piccolo. 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)., Transcript level expression was aggregated to gene level abundance using the R package tximport yielding expression values for 60554 Ensembl genes.

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)., 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. 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.

Molecular Subtype Signature: The PAM50 subtype gene signature was used to classify samples into intrinsic molecular subtypes. 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. 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).,, 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. 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 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.

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. 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.,

Pathway Analysis: All pathway predictors have been previously described. 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.

Batch correction was performed using Removal of Unwanted Variation (RUV-III) using the replicate positive controls to estimate the factors associated with batch effect. 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. 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. 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, data-type-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.

Image analysis pipeline: 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 mm2, average 3.7 mm2). 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 mm2, containing an average of 68,223 nucleated cells. Digitally acquired images were registered in MatLab (MathWorks) utilizing the SURF algorithm in the Computer Vision Toolbox. Nuclear segmentation and color deconvolution were performed using FIJI (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 using the NIH plugin RGB_to_CMYK. The contrast of AEC chromogen intensities on a 0–255 scale in the yellow channel, as compared to RGB or the built-in AEC deconvolution vector, utilizes the full range of intensity without a threshold. For single cell quantification, each channel was normalized by dividing all pixels in each image by the max intensity of that image to rescale intensity values to a range of 0–1. Next, mean intensity from each stain was quantified for every indexed nuclear object in Cell Profiler 3.1.5 (Broad Institute). Image cytometry hierarchical gating was performed in FCS Express Image Cytometry RUO 6.1.4 (DeNovo Software) to quantify distinct populations of cells.

Quantification and Analysis: Each image acquired was registered based on DAPI features acquired from each round of staining. Cellpose, a generalist algorithm for cellular segmentation, was used to generate nuclear and cell segmentation masks with a pre-trained neural network classifier. Extracted single-cell features included centroids and mean intensity of each marker from its biologically relevant segmentation mask (e.g., Ecad_Cytoplasm, Ki67_Nuclei). The last round DAPI image was used to filter out cells lost during each round of cyclic immunofluorescence staining.

For cell type determination and composition analysis, single cell mean intensities from each biopsy were batch corrected using the ComBat algorithm. ComBat was used to adjust the mean and variance of fluorescence intensity on control tissue-microarrays (TMAs) that were stained with each biopsy, and the same adjustments were applied to the corresponding biopsies. Eighteen markers were selected for clustering; some markers were excluded due to tissue loss in the TMA controls. Principal component analysis was performed with scanpy to reduce dimensionality, and Umap (https://arxiv.org/abs/1802.03426) was run on the top 17 principal components to calculate a nearest neighbor graph based on the 30 nearest neighbors. Leiden clustering was performed on the nearest neighbor graph to define clustering-based cell types. The Leiden clustering resolution of 0.5 was selected based on appropriate clustering of technical replicates in the control TMAs.

Immune, endothelial, and stromal cells were identified by manual thresholding and gating. Endothelial cells were defined as CD31+, immune cells were either CD45+ or CD68+ and CD31-, and stromal cells were cytokeratin-, E-cadherin-, CD31-, CD45-, and CD68-. Tumor was defined as cytokeratin+, and proliferating cells were Ki67+. Cell segmentation borders of manually defined cell types were visualized on the images using napari.

To calculate distance to extracellular matrix proteins, a threshold was applied to create a pixel mask of positive staining. The distance from each nuclear centroid to the nearest mask pixel was measured. Cells were grouped into bins of 0–25 microns, 25-50 microns, and 50-75 microns from the mask, and the intensity distributions of cells (n = 2 to 32324 cells) were compared using ANOVA implemented in SciPy. Significance was assigned to Bonferroni-corrected p values < 0.001.

Additional resources

Serial Measurements of Molecular and Architectural Responses to Therapy (SMMART): https://www.ohsu.edu/knight-cancer-institute/serial-measurements-molecular-and-architectural-responses-therapy

Human Tumor Atlas Network (HTAN): https://humantumoratlas.org/

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Mouse monoclonal anti-PD-1	Abcam	Cat# ab52587; RRID: AB_881954
Rabbit monoclonal anti-CD3	Thermo Fisher Scientific	Cat# MA1-90582; RRID: AB_1956722
Mouse monoclonal anti-RORgT	EMD Millipore Sigma	Cat# MABF81; RRID: AB_11205416
Mouse monoclonal anti-NKp46	R&D Systems	Cat# MAB1850; RRID: AB_2149153
Mouse monoclonal anti-CD8a	Thermo Fisher Scientific	Cat# MA5-13473; RRID: AB_11000353
Rabbit monoclonal anti-T-bet	Cell Signaling Technology	Cat# 13232; RRID: AB_2616022
Mouse monoclonal anti-GATA-3	BioCare Medical	Cat# CM 405 A; RRID: AB_10895444
Mouse monoclonal anti-FoxP3	Thermo Fisher Scientific	Cat# 14-4777-82; RRID: AB_467556
Rabbit monoclonal anti-PD-L1	Cell Signaling Technology	Cat# 13684; RRID: AB_2687655
Mouse monoclonal anti-CD20	Abcam	Cat# ab9475; RRID: AB_307267
Mouse monoclonal anti-CD20	Santa Cruz Biotechnology	Cat# sc-70582; RRID: AB_1120279
Mouse monoclonal anti-CD45	Thermo Fisher Scientific	Cat# 14-0459-82; RRID: AB_467274
Mouse monoclonal anti-Tryptase	Abcam	Cat# ab2378; RRID: AB_303023
Mouse monoclonal anti-CD68	Abcam	Cat# ab783; RRID: AB_306119
Rabbit monoclonal anti-CSF1R	Abcam	Cat# ab183316; RRID: AB_2885197
Mouse monoclonal anti-DC-SIGN (DC28)	Santa Cruz Biotechnology	Cat# sc-65740; RRID: AB_1121347
Mouse monoclonal anti-CD66b	BD Bioscience	Cat# 555723; RRID: AB_396066
Rat monoclonal anti-DC-LAMP	Novus	Cat# DDX0191P-100; RRID: AB_2827532
Mouse monoclonal anti-HLA-DPB1	Abcam	Cat# ab157210; RRID: AB_2827533
Mouse monoclonal anti-CD163	Thermo Fisher Scientific	Cat# MA5-11458; RRID: AB_10982556
Mouse monoclonal anti-CD4	Thermo Fisher Scientific	Cat# MA5-12259; RRID: AB_10989399
Mouse monoclonal anti-CD56	Thermo Fisher Scientific	Cat# 07-5603; RRID: AB_2532931
Mouse monoclonal anti-pan cytokeratin	Abcam	Cat# ab27988; RRID: AB_777047
Rabbit monoclonal anti-Ki67	Abcam	Cat# ab16667; RRID: AB_302459
Rabbit polyclonal anti-EOMES (Tbr2)	EMD Millipore Sigma	Cat# AB2283; RRID: AB_10806889
Mouse monoclonal anti-IDO	EMD Millipore Sigma	Cat# MAB10009; RRID: AB_1977068
Rabbit monoclonal anti-Granzyme-B	EMD Millipore Sigma	Cat# 262R-1; RRID: AB_2889344
Mouse monoclonal anti-IL-10	LifeSpan Bio	Cat# LS-B7411-500; RRID: AB_11233179
Rabbit monoclonal anti-ICOS/CD278	LifeSpan Bio	Cat# LS-C210350; RRID: AB_2827535
Rabbit monoclonal anti-CD4	Abcam	Cat# ab213215; RRID: AB_2861280
Rabbit monoclonal anti-CD8	Abcam	Cat# 4207-1; RRID: AB_764503
Mouse monoclonal anti-Siglec-1/CD169	Novus	Cat# NB 600-534; RRID: AB_526814
Rabbit monoclonal anti-CD11b	Abcam	Cat# ab133357; RRID: AB_2650514
Rabbit monoclonal anti-MHC class I (HLA A+ HLA B)	Abcam	Cat# 2307-1; RRID: AB_1267243
Rabbit monoclonal anti-CD11c	Abcam	Cat# ab52632; RRID: AB_2129793
Rabbit monoclonal anti-T-bet/Tbx21	Abcam	Cat# ab150440; RRID: AB_2889209
Mouse monoclonal anti-CCR2	R&D Systems	Cat# MAB150; RRID: AB_2247178
Mouse monoclonal anti-MHC class II HLA-DP/DR/DQ	LifeSpan Bio	Cat# LS-C58506-200; RRID: AB_1511620
Histofine Simple Stain Human MAX PO (Rat) (for Rat primary antibody)	Nacalai USA	Cat# 414311F
Histofine Simple Stain Human MAX PO (R) (for Rabbit primary antibody)	Nacalai USA	Cat# 414144F
Histofine Simple Stain Human MAX PO (M) (for Mouse primary antibody)	Nacalai USA	Cat# 414134F
Mouse monoclonal anti-α-SMA	Santa Cruz	Cat# sc-32251; RRID: AB_262054
Rabbit monoclonal anti-CD3	Abcam	Cat# ab213608; RRID: AB_764498
Rabbit monoclonal anti-CD31	Abcam	Cat# ab218582; RRID: AB_2857973
Rabbit monoclonal anti-CD4	Abcam	Cat# ab196147
Rabbit monoclonal anti-CD44	Abcam	Cat# ab216647; RRID: AB_764499
Rabbit monoclonal anti-CD45	Abcam	Cat# ab200317; RRID: AB_726545
Rabbit monoclonal anti-CD45	Abcam	Cat# ab214437; RRID: AB_726545
Mouse monoclonal anti-CD68	Biolegend	Cat# 916104; RRID: AB_2616797
Mouse monoclonal anti-CD8	Abcam	Cat# ab213017; RRID: N/A
Mouse monoclonal anti-CK14	Abcam	Cat# ab7800; RRID: AB_306091
Rabbit monoclonal anti-CK17	Abcam	Cat# ab185032; RRID: AB_2889195
Mouse monoclonal anti-CK19	Biolegend	Cat# 628502; RRID: AB_439773
Rabbit monoclonal anti-CK5	Abcam	Cat# ab193894
Rabbit monoclonal anti-CK7	Abcam	Cat# ab203434
Rabbit monoclonal anti-Ecad	Abcam	Cat# ab201499
Rabbit monoclonal anti-ER	Abcam	Cat# ab205851; RRID: AB_2728817
Mouse monoclonal anti-HER2	Santa Cruz	Cat# sc-33684; RRID: AB_627996
Rabbit monoclonal anti-Ki67	Cell Signaling Technology	Cat# 12075; RRID: AB_2728830
Mouse monoclonal anti-PCNA	Cell Signaling Technology	Cat# 8580; RRID: AB_11178664
Rabbit monoclonal anti-Phospho-Histone H3	Cell Signaling Technology	Cat# 3465; RRID: AB_10695860
Rabbit monoclonal anti-Vimentin	Cell Signaling Technology	Cat# 9854; RRID: AB_10829352
Rabbit monoclonal anti-Collagen I	Abcam	Cat# ab215969
Mouse monoclonal anti-Collagen IV	Thermo Fisher Scientific	Cat# 51-9871-82; RRID: AB_10853027
Rabbit monoclonal anti-EGFR	Cell Signaling Technology	Cat# 5108; RRID: AB_10694337
Rabbit monoclonal anti-CK8	Abcam	Cat# ab192467; RRID: AB_2864346
Trastuzumab	Genentech	N/A
RPPA Antibodies	MD Anderson Cancer Center Functional Proteomics RPPA Core Facility	https://www.mdanderson.org/research/research-resources/core-facilities/functional-proteomics-rppa-core/antibody-information-and-protocols.html
nCounter Vantage 3D Protein Solid Tumor Panel for FFPE (D) Antibodies	Nanostring	https://www.nanostring.com/wp-content/uploads/2021/01/LBL-10372-03_nCounter_Vantage_3D_Protein_Solid_Tumor_Panel_D_Probe_List_FFPE.xlsx
Biological samples
Universal Human Reference RNA	Agilent Technologies	Cat# 740000
Chemicals, peptides, and recombinant proteins
Human Epidermal Growth Factor	Sigma-Aldrich	Cat# SRP3027
Human Epidermal Growth Factor	Cell Signaling Technology	Cat# 8916
Tris Buffered Saline-Tween (with 0.05% Tween-20, pH 7.4)	Boston Bioproducts	Cat# IBB-181R
McCoy’s 5A (Modified) Medium	Thermo Fisher Scientific	Cat# 16600082
Fetal Bovine Serum	Thermo Fisher Scientific	Cat# 10082147
Alexa Fluor 647 NHS Ester (Succinimidy Ester)	Thermo Fisher Scientific	Cat# A37566
Glucose Oxidase	Sigma-Aldrich	Cat# G2133–50KU
Catalase	Sigma-Aldrich	Cat# C100-50MG
Dextrose (D-Glucose), Anhydrous	Fisher Scientific	Cat# D16–500
Cysteamine	Sigma-Aldrich	Cat# 30070
Hematoxylin, ready-to-use	Dako/Agilent	Cat# S330130-2
Antigen Retrieval Citra Plus Solution (10x Concentrated)	BioGenex	Cat# HK0809K
Peroxidase and Alkaline Phosphatase Blocking Reagent (Dual Endogenous Enzyme-Block)	Dako/Agilent	Cat# S2003
AEC Substrate Kit, Peroxidase (HRP), (3-amino-9-ethylcarbazole)	Vector Laboratories	Cat# SK-4200
Mol Bio Grad Ethanol (200 proof)	Sigma-Aldrich	Cat# E7023-500ML
Citrate Monohydrate	Sigma-Aldrich	Cat# C1909; CAS: 5949-29-1
Target Retrieval Solution, pH 9 (10X)	Agilent	Cat# S236784-2
Phosphate Buffered Saline (10X)	Fisher Scientific	Cat# BP39920
Normal Goat Serum	Vector Laboratories	Cat# S-1000
Bovine Serum Albumin	Sigma-Aldrich	Cat# A7906
SlowFade Gold Antifade Mountant with DAPI	Thermo Fisher Scientific	Cat# S36938
Sodium Hydroxide (Pellets)	Fisher Scientific	Cat# S318-500; CAS: 1310-73-2
Hydrogen peroxide solution, 30% (w/w)	Sigma-Aldrich	Cat# H1009-500ML; CAS: 7722-84-1
Paraformaldehyde	Electron Microscopy Sciences	Cat# 15714; CAS: 30525-89-4
Glutaraldehyde	Electron Microscopy Sciences	Cat# 16120; CAS: 11-30-8
Sodium cacodylate	Electron Microscopy Sciences	Cat# 12300; CAS: 124-65-2
Sucrose	J.T. Baker	Cat# 4072; CAS: 57-50-1
Osmium tetroxide	Ted Pella, Inc.	Cat# 18463; CAS: 20816-12-0
Potassium ferricyanide	Sigma-Aldrich	Cat# 702587; CAS: 13746-66-2
Thiocarbohydrazide	Sigma-Aldrich	Cat# 223220; CAS: 2231-57-4
Uranyl acetate	Electron Microscopy Sciences	Cat# 22400; CAS: #541-09-3
Lead nitrate	Electron Microscopy Sciences	Cat# 17900; CAS: 10099-74-8
Aspartic acid	Sigma-Aldrich	Cat# 11195; CAS: 323194-76-9
Acetone	Electron Microscopy Sciences	Cat# 10014; CAS: 67-64-1
EMbed 812 Embedding Kit with BDMA	Electron Microscopy Sciences	Cat# 14121
Critical commercial assays
nCounter FLEX analysis system	Nanostring	https://www.nanostring.com/products/ncounter-analysis-system/flex-system/
nCounter Vantage 3D Protein Solid Tumor Panel for FFPE (D)	Nanostring	Cat# VPODC-SPKP-HSTF-12
Agencourt AMPure XP PCR Purification Beads	Beckman-Coulter	Cat# A63880
2100 Bioanalyzer High Sensitivity DNA Kit	Agilent	Cat# 5067-4626
KAPA HyperPrep Kit	Roche	Cat# KR8500
KAPA HiFi Hotstart PCR master mix	Roche	Cat# KK8500
Qubit 3 HS dsDNA Assay Kit	Thermo Scientific	Cat# Q33231
xGen Hybridization and Wash Kit	IDT	Cat# 1080584
TruSeq RNA Library Prep Kit	Illumina	Cat# 20020189
Tempus xE	Tempus Labs, Inc.	https://www.tempus.com/genomic-profiling/
GeneTrails Solid Tumor Panel	OHSU Knight Diagnostic Laboratories	https://knightdxlabs.ohsu.edu/home/test-details?id=GeneTrails+Comprehensive+Solid+Tumor+Panel
Reverse Phase Protein Array	MD Anderson Cancer Center Functional Proteomics RPPA Core Facility	https://www.mdanderson.org/research/research-resources/core-facilities/functional-proteomics-rppa-core.html
QIAamp DNA FFPE Tissue Kit	QIAGEN	Cat# 56404
DNeasy Blood & Tissue Kit	QIAGEN	Cat# 69504
NucleoSnap cfDNA Kit	Macherey-Nagel	Cat# 740300.10
SureSelectXT Reagent Kit, HiSeq Platform	Agilent Technologies	Cat# G9611B
SureSelectXT Human All Exon V5, 96	Agilent Technologies	Cat# 5190-6209
Deposited data
gnomAD	89	https://gnomad.broadinstitute.org/
Human reference genome UCSC hg19 (GRCh37/hg19)	Genome Sequencing Consortium	https://genome.ucsc.edu/cgi-bin/hgGateway
GENCODE	102	https://www.gencodegenes.org/
Raw next-generation sequencing data	This work.	https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs002371.v1.p1
TCGA BRCA gene expression	100	https://osf.io/gqrz9
MSigDB	18	https://www.gsea-msigdb.org/gsea/msigdb/
Processed next-generation sequencing data	This work.	Participant ID: HTA9_1; https://htan-portal-nextjs.vercel.app/explore
Protein expression data from RPPA and the Intracellular Signaling Protein Panel	This work.	Participant ID: HTA9_1; https://htan-portal-nextjs.vercel.app/explore
Raw and processed image data from CycIF, mIHC, and EM	This work.	Participant ID: HTA9_1; https://htan-portal-nextjs.vercel.app/explore
Processed images (web viewing)	This work.	Participant ID: HTA9_1; https://idp.tissue-atlas.org/
Experimental models: Cell lines
Human MCF7	Characterized Cell Line Core (MDACC)	https://www.atcc.org/products/all/HTB-22.aspx; RRID: CVCL_0031
Human MDA-MB-468	Characterized Cell Line Core (MDACC)	https://www.atcc.org/products/all/HTB-132.aspx; RRID: CVCL_0419
Human MDA-MB-231	ATCC	https://www.atcc.org/products/all/HTB-26.aspx#generalinformation; RRID: CVCL_0062
Human BT474	ATCC	https://www.atcc.org/products/all/HTB-20.aspx; RRID: CVCL_0179
Human HCC1954	ATCC	https://www.atcc.org/products/crl-2338; RRID: CVCL_1259
Human SKBR3	ATCC	RRID: CVCL_0033; https://www.atcc.org/products/htb-30
HCC1143	ATCC	RRID: CVCL_1245; https://www.atcc.org/products/crl-2321
HCC3153	UT-Southwestern	RRID: CVCL_3377
T47D	ATCC	RRID: CVCL_0553; https://www.atcc.org/products/htb-133
AU565	ATCC	RRID: CVCL_1074; https://www.atcc.org/products/crl-2351
MDAMB436	ATCC	RRID: CVCL_0623; https://www.atcc.org/products/htb-130
Oligonucleotides
See Table S2	This paper	N/A
Software and algorithms
LabVantage	LabVantage Solutions Inc.	https://www.labvantage.com/
LabKey	74	https://www.labkey.com/
Removal of Unwanted Variation (RUV-III)	111	https://cran.r-project.org/package=ruv
Removal of Unwanted Variation (RUVSeq)	105	https://bioconductor.org/packages/release/bioc/html/RUVSeq.html
MATLAB	MathWorks	https://www.mathworks.com/products/matlab.html
syngo.via	Siemens Healthcare GmbH	https://www.siemens-healthineers.com/en-us/molecular-imaging/pet-ct/syngo-via
Horos	Nimble Co LLC	https://horosproject.org/download-horos/
Cluster 3.0	110	http://bonsai.hgc.jp/∼mdehoon/software/cluster/software.htm
TreeView	109	http://jtreeview.sourceforge.net/
MATLAB, Computer vision toolbox	MathWorks	https://www.mathworks.com/help/vision/
SURF algorithm (for MATLAB)	113	https://www.mathworks.com/help/vision/ref/detectsurffeatures.html
ImageJ	81	https://imagej.nih.gov/ij/
RGB to CMYK (FIJI Plugin)	81	https://imagej.net/tutorials/rgb-to-cmyk
CellProfiler 3.0	114	https://cellprofiler.org/
FCS Express Image Cytometry RUO	DeNovo Software	https://denovosoftware.com/
mpileup (Samtools)	90	http://samtools.sourceforge.net/
BWA-MEM	87	http://bio-bwa.sourceforge.net/
GATK	88	https://gatk.broadinstitute.org/hc/en-us
Galaxy	75	https://galaxyproject.org/
ape	91	http://ape-package.ird.fr/
CNVkit	92	https://github.com/etal/cnvkit
Trim Galore	Babraham Institute	https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/
Kallisto	101	https://github.com/pachterlab/kallisto
tximport	103	https://bioconductor.org/packages/release/bioc/html/tximport.html
SVA/Combat	106	https://bioconductor.org/packages/release/bioc/html/sva.html
GSVA	17	https://www.bioconductor.org/packages/release/bioc/html/GSVA.html
Pathway Commons	24	https://www.pathwaycommons.org/
Transcriptional Regulator Analysis	This work.	https://dx.doi.org/10.5281/zenodo.5608590
CausalPath	25	https://causalpath.org/; https://github.com/PathwayAndDataAnalysis/causalpath
Neoepiscope	93	https://github.com/pdxgx/neoepiscope
OptiType	94	https://github.com/FRED-2/OptiType
MHCflurry	95	https://github.com/openvax/mhcflurry
DIDA-Seq Error Correction Scripts	14	https://github.com/ohsu-cedar-comp-hub/DIDA-Seq
R 3.6	R Foundation for Statistical Computing	https://www.R-project.org/
DIDA-Seq Bayesian Overlap Script	97	https://github.com/mheskett/SMMART
Zen 2.3 Slidescan	Zeiss	https://www.zeiss.com/microscopy/us/products/microscope-software.html
Cellpose	117	https://github.com/MouseLand/cellpose
ComBat	118	https://doi.org/10.1093/biostatistics/kxj037
scanpy	119	https://github.com/theislab/scanpy
SciPy	122	https://github.com/scipy/scipy
Umap	https://arxiv.org/abs/1802.03426	https://github.com/lmcinnes/umap
napari	121	https://zenodo.org/record/3555620
OMERO	76	https://www.openmicroscopy.org/omero/
Cyclic Immunofluorescence Analysis	This work.	https://dx.doi.org/10.5281/zenodo.5637447
maps2ometiff	This work.	https://dx.doi.org/10.5281/zenodo.5608828
em_segmentation	https://doi.org/10.1101/2021.05.27.446019	https://github.com/archana2890/em_segmentation
Microscopy Image Browser	82	http://mib.helsinki.fi/
Amira Software	Thermo Fisher Scientific	https://www.thermofisher.com/us/en/home/electron-microscopy/products/software-em-3d-vis/amira-software.html
Dragonfly	Object Research Systems	https://www.theobjects.com/dragonfly/
IrfanView	Irfan Skiljan	https://www.irfanview.com/
Maps Software	Thermo Fisher Scientific	https://www.thermofisher.com/us/en/home/electron-microscopy/products/software-em-3d-vis/maps-software.html
ichorCNA	99	https://github.com/broadinstitute/ichorCNA
microManager	83	https://micro-manager.org/Download_Micro-Manager_Latest_Release
Matlab packages for raw PALM/STORM image data processing	84	https://www.ohsu.edu/school-of-medicine/nan-lab/resources
Other
Signature Series Cover Glass	Thermo Fisher Scientific	Cat# 12460S
Microscope Slides	Mercedes Scientific	Cat# TNR WHT45AD
Plastic Coverslips	IHC World	Cat# IW-2601
Rectangular Cover Glass	Corning	Cat# 2980-243, 2980-245
Nunc Lab-Tek II Chambered Coverglass	Thermo Fisher Scientific	Cat# 155409
SEM pin stub	Ted Pella, Inc.	Cat# 16145
Leitsilber conductive paint	Ted Pella, Inc.	Cat# 16035
H20E Epo-TEK Silver conductive epoxy	Ted Pella, Inc.	Cat# 16014
Fisherbrand Transfer pipettes	Fisher Scientific	Cat# 13-711-7M
Axygen Centrifuge microtubes	Fisher Scientific	Cat# 14-222-180
Cytivia Whatman filter paper	Fisher Scientific	Cat# 1001-090
Flat embedding mold	Ted Pella, Inc.	Cat# 10504
Helios NanoLab G3 DualBeam FIB-SEM	FEI Company (now Thermo Fisher Scientific)	https://www.thermofisher.com/us/en/home/electron-microscopy/products/dualbeam-fib-sem-microscopes.html
Gemini 550 Crossbeam FIB-SEM	ZEISS International	https://www.zeiss.com/microscopy/us/products/fib-sem-instruments/crossbeam.html
Ultramicrotome Leica EM UC7	Leica Microsystems	https://www.leica-microsystems.com/products/sample-preparation-for-electron-microscopy/p/leica-em-uc7/
Trim90 diamond knife	DiATOME	https://www.diatomeknives.com/product.aspx?pid=416
EM ACE600 High Vacuum Sputter Coater	Leica Microsystems	https://www.leica-microsystems.com/products/sample-preparation-for-electron-microscopy/p/leica-em-ace600/
Managed workforce (800 hours)	CloudFactory	https://www.cloudfactory.com/