Bone metastases (BoMs) occur in approximately 65%–75% of breast cancer patients with relapsed disease, resulting in significant comorbidities, such as fractures and chronic pain (1). Following colonization to the bone, breast cancer cells exploit the local microenvironment by activating osteoclasts, which in turn provides proliferative fuel for tumor cells (2). This process is targeted clinically using antiosteoclast agents, such as bisphosphonates and RANKL inhibitors, yet these therapies do not confer significant survival benefits (3).

Importantly, the majority of breast cancers that metastasize to bone are estrogen receptor (ER) positive and present clinically in the context of long-term endocrine therapies, such as selective ER modulators and aromatase inhibitors (4). In vivo models of BoM have unfortunately been somewhat restricted to ER-negative disease due to the more indolent characteristics of ER-positive cell lines (5). Molecular characteristics of ER-positive specimens that have recurred in an estrogen-blunted system, which represents the major burden of breast cancer BoM, are thus essential to reinforce the significant scientific contributions made using in vivo bone metastasis models (6–9). Nonetheless, data sets are currently limited, in part due to the practical difficulties of obtaining and processing human BoM specimens (10).

Because of the untapped potential of archived, decalcified BoM specimens; the burden of BoMs in breast cancer patients; and the lack of long-term endocrine-treated tumor data sets, the performance of ecRNA-seq from decade-old, degraded, and decalcified tumor samples was assessed. Following this evaluation, ecRNA-seq was then applied to a collection of 11 ER-positive patient-matched primary breast cancers and BoMs to define transcriptional evolution in breast cancer cells following metastatic colonization in the bone and years of endocrine therapy.

ecRNA-seq of aged and decalcified breast cancers. To determine the feasibility of sequencing an aged, FFPE, and decalcified tumor cohort, ecRNA-seq on two separate sample sets was performed. The first sample set included four cases of primary breast tumors that, at the time of resection, were split in two. One section was flash frozen and stored at –80°C, and the other tumor section was FFPE and stored at room temperature. Storage times ranged from 8.2 to 12.3 years. RNA-sequencing quality control (QC) analyses after alignment showed differences in GC content and insert size, yet gene body coverage and transcript diversity assignments were largely similar (Figure 1A). After quantifying and normalizing gene abundances, expression correlations between frozen and FFPE-matched samples were assessed using log2-transformed trimmed mean of M values–normalized (TMM-normalized) CPM (log2normCPM) values. Pearson r correlations ranged from 0.929 to 0.963, with an average correlation of 0.953 (Figure 1B). The same analysis was performed using a second sample set of matched FFPE-decalcified and FFPE-nondecalcified samples. Again, no concerning deviations in ecRNA-seq quality metrics were observed between the two differently processed sample groups (Figure 1C), and Pearson r expression correlations ranged from 0.936 to 0.969 (Figure 1D). Furthermore, correlation matrices of the two sample sets showed that matched tumor sample expression values were more similar to each other than expression values from tumors with equivalent processing and storage (Supplemental Figure 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.95703DS1). Full ecRNA-seq metrics from the QC analysis did reveal differences in some metrics between FFPE and flash-frozen tissue (i.e., splice junction loci number) that may be informative for other applications (Supplemental Tables 1 and 2). In summary, ecRNA-seq shows outstanding quality metrics for analysis of aged FFPE and decalcified BoMs samples.

ecRNA-seq of breast cancer BoMs. Following the validation of ecRNA-seq, a cohort of 11 ER-positive patient-matched primary tumors and BoMs was acquired through the University of Pittsburgh Health Science Tissue Bank (Table 1 and Supplemental Table 3). Abstracted clinical records showed that nearly all patients (10 of 11) were documented as having received adjuvant endocrine therapy, and bone metastasis–free survival ranged from 0 (synchronous) to greater than 5 years, with the most common site of bone metastasis being the vertebral column. ecRNA-seq was performed on the 22 samples, yielding an average read count of 59,570,288 and an average Salmon transcript-mapping rate of 92.9% (Supplemental Table 4). Consistent with the initial QC studies above, quality metrics on these samples showed consistent gene body coverage, GC content, insert sizes, and transcript diversity, regardless of decalcification status (Supplemental Figure 2 and Supplemental Table 5). Furthermore, since samples within the cohort had been surgically excised and banked many years apart, all paired specimens underwent an analysis of shared variants, which confirmed tumor pairs were patient matched (Supplemental Figure 3).

Clustering and temporal expression shifts. Unsupervised hierarchical clustering of patient-matched pairs revealed that decalcification of BoMs did not produce independent clades, with 5 of 11 BoMs clustering in the same doublet clade as their matched primary tumor (denoted with an asterisk in Figure 2A). Notably, 3 of 5 doublet-clustering cases were synchronous metastases. Discrete PAM50 intrinsic subtype assignments were identical in 7 of 11 pairs. Three pairs switched from LumA to LumB in the metastasis and another was classified as normal subtype in the primary tumor and LumB in the BoM (Figure 2B). To obtain more granularity than discrete PAM50 calls, probability scores for each PAM50 subtype were assigned (Figure 2B and Supplemental Table 6). Her2 and LumB profile gains (defined as a probability gain of >10% in a matched BoM) were the most common — both being observed in 5 of 11 cases (Figure 2B). Given shifts in expression profiles of BoMs and doublet clustering of synchronous BoMs, temporal influence on transcriptional evolution was analyzed. Pearson r correlations between each patient-matched pair using log2normCPM expression values were utilized as a metric for transcriptional similarity. Expression pair similarity was significantly correlated (Pearson r = –0.864, P < 0.001) with time from primary tumor diagnosis to bone metastasis (Figure 2C).

Bone is the most common site of distant recurrence for patients with ER-positive breast cancer, yet comprehensive sequencing data sets of endocrine therapy–treated, metastatic samples are currently limited. This is in part due to the challenge of obtaining tissue and degradation of nucleic acids caused by decalcification. In this study, we found that aged FFPE and FFPE-decalcified tumors showed highly similar transcript quantification values as matched flash-frozen and FFPE-nondecalcified tumors. As a proof of concept, we then applied ecRNA-seq to a cohort of patient-matched primary and BoMs collected over a period of 5 years. We identified subtle shifts in intrinsic subtypes and found a strong temporal influence on transcriptional evolution in breast cancer recurrences. Furthermore, we created several DEG sets/signatures that are prognostic and point toward acquired RBBP8 loss and CDK/Rb/E2F and FGFR pathway gains as mediators of ER-positive breast cancer progression. Finally, we found BoMs commonly gain or lose expression in clinically actionable genes, which may be distinct from primary tumors.

ecRNA-seq is an effective method for quantifying expression on aged, FFPE, and decalcified tumor specimens. Previous work has assessed nucleic acid amplification success, DNA sequencing, and RNA integrity metrics using decalcified samples (17, 28, 29); however, a comprehensive analysis of RNA sequencing, to our knowledge, has not yet been performed. Consistent with only very minor differences between GC content, insert sizes, and other QC metrics, gene expression values between aged-matched FFPE/flash-frozen and FFPE-decalcified/FFPE-nondecalcified tumors are highly correlated (Pearson r range 0.929–0.969). This study reinforces and should encourage the use of capture-hybridization approaches to sequence RNA from retrospectively collected, low-yield, highly degraded, and decalcified archival specimens (Supplemental Table 14 and refs. 22–24). Expanding sample sets and modalities for genome-wide characterization, especially for rare specimen cohorts that may be impractical to obtain prospectively in large numbers, will accelerate translational discoveries.

Given promising results from our evaluation, we applied ecRNA-seq in a proof-of-concept effort to characterize the transcriptome of 11 archival patient-matched ER-positive primary and recurrent metastases — 3 cases having treatment-naive, synchronous BoMs and 8 recurrent cases harboring long-term endocrine-therapy treated metastases. In the recurrent cases, bone metastasis–free survival ranged from 18 to 65 months. Despite a large portion of the BoMs being decalcified, global transcriptome QC metrics showed similar features (i.e., GC content, insert sizes, gene body coverage, and transcript assignment diversity) and no outliers. Consistent with this, unsupervised hierarchical clustering showed no distinct clusters of decalcified samples, with 5 BoMs clustering in the same doublet clade as their patient-matched primary breast cancer. Interestingly, 3 of these doublet-clustering pairs were clinically synchronous, treatment-naive BoMs, implying limited transcriptional evolution from the primary tumor in synchronous metastases. This was further corroborated with a striking negative correlation between patient-matched expression similarity and time to bone metastasis, suggesting metachronous metastases that present later clinically in their treatment course are more dissimilar from their derived primary lesions. Intrinsic subtyping revealed 4 of the 11 cases changed PAM50 subtypes, with all 4 cases switching to LumB in the metastasis. Subtle Her2 and LumB profile shifts were also the most common when observing continuous PAM50 probability scores, even in samples that remained concordant in their discrete PAM50 assignments. A recent targeted expression study analyzed PAM50 assignments in 123 matched breast cancer metastases, and the authors found similar frequencies of LumB and Her2 acquisitions in ER-positive metastatic tumors (30). Given this transcriptional evolution to more LumB and Her2 profiles, a thoughtful reevaluation of therapy selection in the advanced and perhaps the adjuvant setting may be necessary — especially considering HER2-targeted therapies are generally reserved for patients with HER2-positive primary disease.

We found 207 genes to be differentially expressed between primary tumors and patient-matched BoMs. The top upregulated genes — BGLAP, RANKL, and PTH1R — belonged to osteogenic gene sets and all showed significant expression gains. This supports in vivo modeling observations of breast cancer osteomimicry and hijacking of the bone microenvironment (31). Downregulated gene sets included genes involved in broad categories, such as cellular adhesion, hemidesmosome assembly, and epithelium development, pointing toward specific biological programs lost following metastatic colonization. Moreover, when either the upregulated or downregulated genes are expressed coordinately in primary tumors, we found that they confer worse and better outcomes, respectively, in ER-positive tumors, suggesting some tumors may develop these transcriptional programs early in their evolution. Finally, a differential expression analysis between endocrine-naive primary tumors and long-term endocrine-treated BoMs identified a larger list of DEGs. Importantly, known mediators of endocrine resistance are represented in the list, including dysregulated expression of Wnt family members (32), expression gains in FGFR1 (33) and FOXC1 (34), and loss of ESR1 expression (35). Notably, many of these genes do not overlap with the differential expression analysis that included the synchronous metastases, suggesting expression alterations specific to late recurrent therapy-treated tumors. This nonoverlapping gene set included a greater than 2-fold average expression gain of ABCG2 — a multidrug resistance protein shown to be active in breast cancer (36, 37) — in therapy-exposed metastases and loss of CDKN2A. CDKN2A encodes p16, a negative regulator of CDK4/CDK6 and is located on a common somatically deleted region (9p21) in cancer (38). Given the recent success of CDK4/CDK6-inhibiting compounds (palbociclib and ribociclib) in treating ER-positive breast cancers, this recurrent, acquired, metastatic-specific loss of CDK2NA is a clinically important observation (39–41).

Following significant gene-level changes, a GSEA defined enriched and diminished pathways in breast cancer BoMs. Enriched genes included those involved in G2M checkpoint and E2F targets. Consistent with the observed LumB enrichments, our data suggest breast cancer cells may develop a more proliferative phenotype following bone colonization, and the strong enrichment of E2F signature in metastatic disease again highlights the CDK/Rb/E2F pathway as a prime target. Interestingly, another study that utilized a targeted gene expression platform found proliferative gene signatures in ER-positive metastases may be more accurate at predicting overall survival than signatures in the primary tumor (30). A survival analysis for this work was impractical given the small set of patient-matched pairs, but future meta-analyses are warranted to determine if gene expression signatures in metastases are better predictors of overall survival in the advanced setting, especially given the significant transcriptomic shifts observed in this study.

The most significant gene set enriched in bone metastasis was an experimental perturbation gene set involving the knockdown of the tumor suppressor RBBP8 (42). RBBP8 (also known as CtIP) binds directly to Rb, mediates cell cycle regulation, and helps maintain genomic stability, and loss of RBBP8 incurs tamoxifen resistance and sensitizes breast cancer cells to PARP inhibition in vitro (43–46). Concordant with the GSEA analysis, BoMs have significant expression loss of RBBP8, with 45% of cases showing a greater than 2-fold decrease in expression. We found that low RBBP8 expression in ER-positive tumors confers poorer DSS and bone metastasis–free survival outcomes. These observations point to RBBP8 loss in metastatic breast cancers as being a compelling, perhaps therapeutically relevant candidate for further preclinical investigations.

Finally, considering that we have previously shown that brain metastases acquire highly recurrent gains in clinically actionable genes following colonization (13), particularly in HER2, we analyzed an expanded set of these genes in BoMs. All tumors harbored significant gains and losses, some of which were highly recurrent. PIK3C2G, a relatively uncharacterized gene in the PI3K pathway, was the most recurrent gene expression loss. Other notable losses included ESR1, CDKN2A, and GATA3. Intriguingly, GATA3 is one of the most recurrently mutated genes in breast cancer, being particularly enriched in ER-positive disease (47). Moreover, GATA3 inhibits breast cancer metastasis in various model systems, and given that losses of GATA3 in ER-positive BoMs are common, further evaluation of GATA3 as a potentially targetable breast cancer metastasis suppressor gene should be encouraged (34, 48, 49). Metastatic expression gains were found in FGFR family members (FGFR3, FGFR4, FGFR1), ALK, and KDR — all protein products having small molecules currently in clinical trials. Some highly recurrent expression gains (i.e., EPHA3, PTPRD, PDGFRA, PTCH1) were exclusive to long-term endocrine-treated BoMs, suggesting them as clinically actionable candidate mediators of therapy resistance. Collectively, these observations provide further evidence of acquired transcriptional remodeling in metastatic lesions and suggest that precision care in advanced cancers may benefit from defining longitudinal changes in tumor transcriptomes. Further research into these longitudinal transcriptional remodeling events is demanded — especially given their high recurrence rates — including identifying events that may be specific or more likely to occur in metastases from certain anatomic sites, such as HER2 gains in brain metastases (13).

Although this study points toward ecRNA-seq as being a viable option to characterize the transcriptome of archived, decalcified specimens, there are limitations. First, multiple methods are used for decalcification with varying effects on nucleic acids, and we were unaware of this information for the profiled specimens, as it is rarely recorded in clinical notes (17). Second, in metastasis-versus–primary tumor expression studies, it is difficult to deconvolute expression contributions from tumor cells and cells within the altered microenvironment of the distant organ site. To limit these artifacts in this study, regions of high tumor cellularity in the bone metastasis were cored by a trained molecular pathologist for RNA extraction, which was corroborated by ecRNA-seq–derived tumor purity estimates — as no significant tumor purity differences between primary and metastatic tumors (Supplemental Table 15) were observed (50). Nonetheless, single-cell sequencing approaches will be crucial to bring cell-level resolution to identifying transcriptional differences between primary and metastatic cells. Novel computational methods that deconvolute heterogeneous sample sets, until single-cell sequencing becomes more widely adopted, will also be essential (51–53). All of this withstanding, features of this data set are encouraging, such as patient-matched tumors clustering together, intuitive PAM50 assignments, corroboration of other groups’ findings, and treatment-specific gains and losses. A third limitation is performing an analysis on already colonized metastatic lesions, as this likely masks some of the intermediate steps involved in metastasis, such as epithelial mesenchymal transition and tumor-initiating programs. Finally, another limitation of this study is the small sample size. Hopefully, these results will encourage the use of ecRNA-seq to transcriptionally profile other highly degraded samples and begin a collection of genomic data from metastatic or rare tissues for integration. Importantly, deidentified clinical data should be provided alongside the sequencing, as in this study, to allow more fluid merging of data sets and inspire clinical phenotype-driven analyses.

In summary, this study both validates the use of ecRNA-seq to transcriptionally profile highly degraded RNA from decade-old and decalcified tumor specimens and defines multiple acquired and lost transcriptional programs in ER-positive BoMs. We highlight acquired changes in the CDK/Rb/E2F and FGFR pathways, particularly relevant given the recent clinical use of CDK4/6 inhibitors, and point toward RBBP8 as a particularly compelling candidate in breast cancer progression. We also found significant gains in clinically actionable genes that may have not been appreciated in primary tumors, reinforcing the need for longitudinal characterizations of tumor transcriptomes to guide clinical care.

Sample acquisition. Eleven sets of FFPE primary breast tumors and patient-matched BoMs (total of 22 samples) were obtained from the Health Sciences Tissue Bank, a certified honest broker facility at the University of Pittsburgh that maintains an IRB-approved protocol for collecting excess tissue and biological materials. A molecular pathologist reviewed hematoxylin and eosin slides from each sample and then subsequently cut one to two 0.6- to 1-mm cores from the paraffin block exclusively from regions of high tumor cell purity for RNA extraction.

Tissue processing and RNA extraction. Tissues were digested overnight with shaking at 300 rpm at 56°C in PKD buffer with the addition of proteinase K (Qiagen). RNA extraction was then performed with the FFPE RNeasy kit (Qiagen, catalog 73504) according to the manufacturer’s instructions under sterile RNase/DNase-free conditions. RNA concentration was determined with the Qubit 3.0 Fluorometer (ThermoFisher Scientific). Quality RNA integrity number scores, and fragment sizes (DV200 metrics) were obtained utilizing either the Agilent 2100 Bioanalyzer or the Agilent 4200 TapeStation.

ecRNA-seq. Sequencing library preparation was performed using a minimum of 25 ng RNA according to Illumina’s TruSeq RNA Access Library Preparation protocol. Indexed, pooled libraries were then sequenced on the Illumina NextSeq 500 platform with high-output flow cell–producing stranded, paired-end reads (2 × 75 bp). A target count of 50 million reads per sample was used to plan indexing and sequencing runs.

RNA sequencing expression quantification and normalization. RNA transcripts from paired-end FASTQ files were mapped and quantified using k-mer–based quasi-mapping with seqBias and gcBias corrections (Salmon v0.7.2, 31-kmer index built from GRCh38 Ensembl v82 transcript annotations) (54). Transcript-level abundance estimates were collapsed to gene-level estimates using tximport (55). To filter out nonexpressed genes or genes with low expression, only genes harboring a TPM value of more than 0.5 in at least 10% of samples were considered. Gene-level counts or log2normCPM values were implemented for subsequent analyses (56, 57).

Expression correlations and ecRNA-seq quality assessment. Exome-capture ecRNA-seq was performed on two cohorts. A set of four aged (ranging from 8–12 years) primary breast cancer specimens that, at the time of surgical resection, were split in half and either immediately embedded in optimal cutting temperature compound and flash frozen for storage at –80°C or FFPE and stored at room temperature. A second cohort consisted of three breast cancer BoMs that, at the time of resection, were split in half and either decalcified or nondecalcified and processed to FFPE. These data sets were quantified and normalized as described above. Pearson r correlations between all samples were determined using log2normCPM values. Reads and mapping rates were obtained from Salmon. More detailed ecRNA-seq metrics were calculated and plotted using QoRTs (v1.1.8) following two-pass read alignment with STAR (v2.4.2a) for the 11 patient-matched cases (58, 59).

tumorMatch patient-matched sample identifier. To confirm samples were patient-matched, variants from ecRNA-seq were called using GATK’s Best Practices for variant calling on ecRNA-seq (60). Output.vcf files were then provided to tumorMatch, a custom R script that analyzes a pool of.vcf files and calculates the proportion of shared variants (POSV) between each sample. These proportion values were visualized using corrplot in R (61).

Unsupervised hierarchical clustering and intrinsic subtyping. Hierarchical clustering was performed using the heatmap.3 function (https://raw.githubusercontent.com/obigriffith/biostar-tutorials/master/Heatmaps/heatmap.3.R) in R on log2normCPM values of the top 5% most variable genes (defined by IQR) with 1 minus Pearson correlations as distance measurements and the “average” agglomeration method. PAM50 calls were generated using the molecular.subtyping function in genefu (62). A separate cohort of exome-capture RNA-sequencing expression data from primary tumors (n = 12 ER negative, 9 ER positive) was merged with the bone metastasis cohort to help account for test set bias and increase the stability of the PAM50 assignments (63). To call PAM50 subtypes, for each sample in the bone metastasis cohort, a random subset of primary tumor expression data were added to enforce a balanced distribution of ER-positive and ER-negative tumors. This was repeated 20 times, and the discrete PAM50 subtype was designated as the mode of this 20-fold PAM50 assignment test, while the final probability score was an average of all 20 probability scores from genefu.

Differential gene expression.Salmon gene-level counts with effective lengths of target transcripts were used to call DEGs between primary tumors and BoMs using DESeq (64). Given that samples were patient matched, a multifactor design was implemented (~patient + tumor [i.e., metastasis vs. primary]). Genes with an FDR-adjusted P value of less than 0.10 were assigned as differentially expressed. An unclustered heatmap using log2normCPM values from the 207 DEGs, first segregated by metastatic log2 fold change gains and losses and then sorted by DESeq2-adjusted P values, was created in R using heatmap.3. DEGs within the MsigDB database that were gained or lost in BoMs were separately interrogated for gene ontology (GO: Biological Process) enrichment by computing significant (top 10 gene sets) gene overlaps using the MsigDB online tool (27).

ssGSEA signatures and METABRIC survival analyses. Microarray expression along with DSS data were obtained from the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) through Synapse (https://www.synapse.org/, Synapse ID: syn1688369), following IRB approval for data access from the University of Pittsburgh (26). Normalized expression values from IHC-confirmed ER-positive tumors were used to develop a ssGSEA for strongly DEGs (adjusted P < 0.05) between primary tumors and BoMs (25). 48 genes that carried positive log2 fold change values and had a corresponding gene expression value in METABRIC were assigned to the “boneMetSigUp” signature; 74 genes with negative log2 fold change values were assigned to the “boneMetSigDown” signature. A ssGSEA score for each sample from both gene sets was calculated using the ssGSEA method implemented in the GSVA R package (65). Binary dichotomization of samples (low vs. high) based on ssGSEA signature score strata (10th, 25th, 50th, 75th, 90th percentiles) and log-rank testing were used to assess significant differences in DSS (66). The strata with the most significant log-rank P values were plotted using survminer from CRAN (67).

RBBP8 survival analysis.RBBP8 expression was further interrogated and plotted using log2normCPM values from patient-matched tumors. RBBP8 expression influence on DSS in METABRIC ER-positive patients was interrogated as described above. RBBP8 expression influence on bone metastasis–free survival was assessed by querying a GCRMA-normalized microarray expression data set (GSE12276) from 204 primary tumors and associated survival data as described above (68).

Gains and losses in clinically actionable genes. The clinically actionable gene set was obtained using the Drug Gene Interaction Database (DGBIdB 2.0) (69). Considering that metastatic fold change distributions calculated from log2normCPM values for all genes were slightly different for each case, stringent case-specific fold change thresholds were used to transform continuous fold change values into discrete “expression alterations.” More specifically, if the fold change value for a clinically actionable Gene_X was greater than the 95th percentile of all gene fold change values, in that case, Gene_X would be designated as having a significant, case-specific expression gain. If the fold change value for Gene_Y was lower than the 5th percentile, Gene_Y was designated as having a significant, case-specific expression loss (Supplemental Figure 6 and Supplemental Table 13). After assigning discrete expression alteration calls to clinically actionable genes, data were visualized using the oncoprint function in ComplexHeatmap (70).

Data and code. Raw expression values for all samples, as well as R code related to this study, are deposited in GitHub (https://github.com/npriedig/).

Statistics. To determine DEGs between patient-matched primary tumors and BoMs, DESeq2 was used. DESeq2 is designed for ecRNA-seq gene-based count abundance estimates and assigns differential expression P values based on a negative binomial distribution. For Kaplan-Meier curves, the log-rank test was used to determine statistically significant differences in event probabilities (i.e., death or time to metastasis) based on binary expression or signature strata. For single-gene queries, paired Wilcoxon signed-rank tests on log2normCPM values were used. A P value of less than 0.05 was considered significant. If error bars are shown, they represent mean ± SD. All statistical tests are 2 tailed, unless otherwise specified.

Study approval. A protocol for this study was reviewed and approved by the University of Pittsburgh IRB Office. Tissue and associated data were obtained from the Health Sciences Tissue Bank, a certified honest broker facility at the University of Pittsburgh that maintains an IRB-approved protocol for collecting excess tissue and biological materials. Requirement for informed consent was waived, considering all samples were deidentified, there was no more than minimal risk to human subjects, and all tissue was obtained as part of routine clinical care.

This project used the University of Pittsburgh Health Sciences Core Research Facilities Genomics Research Core and Health Sciences Tissue Bank and the UPMC Hillman Cancer Center Tissue and Research Pathology Services, supported in part by award NIH grant P30CA047904. The authors would like to thank the patients who contributed samples to this study as well as Lori Miller (University of Pittsburgh), Alma E. Heyl (UPMC), and Jorge A. Rios (UPMC) for their efforts in collecting tissues. This research was also supported in part by the University of Pittsburgh Center for Research Computing and the University of Pittsburgh Center for Simulation and Modeling through the resources provided. Research funding for this project was provided in part by Susan G. Komen Scholar awards (SAC110021 to AVL and SAC160073 to SO), the Breast Cancer Research Foundation (to AVL and SO), the Fashion Footwear Association of New York, the Magee-Womens Research Institute and Foundation, the Nicole Meloche Foundation, the Penguins Foundation, the Penguins Alumni Foundation, and the Mario Lemieux Foundation as well as through a postdoctoral fellowship awarded to RJW from the Department of Defense (BC123242). NP was supported by a training grant from the NIH/National Institute of General Medical Sciences (2T32GM008424-21) and an individual fellowship from the NIH/National Cancer Institute (5F30CA203095).