The breast cancer data set of 862 samples and corresponding Infinium 450K DNAm data were obtained from TCGA Genomic Data Commons (GDC), which houses molecular profiles of different cancer types from the TCGA project. The data set included four different survival outcomes, including overall survival, progression-free interval, disease-specific survival, and disease-free interval. These outcomes were defined and curated using available clinical information by Liu et el. and were made available in GDC (Liu et al., 2018). BC subtypes of samples were classified based on Spearman distance between samples’ expression profiles and 50-gene centroids (representing different subtypes) constructed using the Prediction Analysis of Microarray (PAM) algorithm (Tibshirani et al., 2002). Inclusion criterion was availability of OS or PFI endpoints. Exclusion criteria were based on tumor types, stages and missingness of covariates which were included in the model. Specifically, we excluded subjects whose menopausal status (n = 80), age (n = 1), race (n = 16) and tumor stage (n = 7) information were missing. We also excluded subjects with stage four cancer (n = 16), and adjacent normal samples (n = 110) and HER2 subtype (n = 46). We randomly chose one sample from a subject if the subject had multiple samples. Subjects whose race was not white were categorized as other, and whose stage was not three was categorized as stage one or two. The menopausal status was categorized into two categories, post-menopause, or pre- or peri-menopause. For this study, we analyzed overall survival and progression-free interval outcomes.

In the end, our study included 692 breast tumor samples from female TCGA BC participants. A majority of these participants were post-menopausal (Table 1). Tumor subtypes of TCGA participants were previously classified using the 50-gene PAM50 model (Parker et al., 2009; Koboldt et al., 2012), with 415 samples classified as luminal A, 141 as luminal-B, and 136 as basal-like (Table 1). During follow-up, 87 participants died, and 94 progressed to a new tumor event. The median follow-up time was 42 months (IQR = 55) for OS endpoint and 27 months (IQR = 36) for PFI endpoint (Table 1).

We used minfi (Version 1.34.0) in R (Aryee et al., 2014) to normalize, correct for background noise, and convert the raw signal intensities into beta-values. The raw DNAm data was preprocessed using Noob (normal-exponential out-of-band) background correction with dye-bias normalization in minfi. We included 485,512 CpG loci in our analysis of the TCGA and GEO datasets.

The purity of each tumor sample was estimated using the InfiniumPurify (Version 1.3.1) R package (Qin et al., 2018). First, the top 1000 informative, differentially methylated CpG sites (iDMCs) were selected using the rank-sum test to detect significant differences in methylation levels (beta-values) between tumor and normal samples. Additionally, the variances of methylation levels of tumor samples were required to be greater than 0.005 for iDMCs. Next, iDMCs were categorized as hypermethylated if the mean beta-values of the tumor samples were greater than the normal samples, and into the hypomethylated group, otherwise. Beta-values of hypomethylated iDMCs were transformed into 1 − beta-values, and those of hypermethylated iDMCs were identity transformed. The distribution of the transformed beta-values was estimated using the Gaussian kernel, of which the mode was taken to be estimated tumor purity. The tumor type was set as breast cancer because no normal samples were used.

The cell-type proportions of each sample was estimated using the RefFreeEWAS (Version 2.2) R package (Houseman et al., 2014). The cell proportions for each sample were estimated by non-negative matrix factorization (NMF), which decomposes the beta matrix into a cell-type-specific methylation matrix and a cell proportion matrix. We used an existing cell-type-specific methylation matrix to initialize the NMF rather than estimating it de novo. This cell-type-specific methylation matrix was previously estimated using the informative probes selected based on reference DNA methylation profiles of known cell lines, and the TCGA breast cancer data set (Onuchic et al., 2016). Eight cell types were inherited from the initial cell-type specific methylation matrix, and thus no estimation of the number of cell types was conducted. The number of iterations was set to 100. The final cell-type-specific methylation matrix did not differ much from the initial cell-type-specific methylation matrix.

All statistical analyses were performed using R. The site-specific analysis was conducted using TCGA data. PH assumptions were assessed by examining log-log survival curves of OS and PFI for covariates including age, ethnicity, stage, menopausal status, and subtype. The log-log survival curves of PFI and OS for subtype overlapped noticeably across different subtypes. For this reason, Cox PH models was stratified by tumor subtype, allowing the baseline hazards functions to vary across the subtypes. The association between DNAm and our outcomes were assessed for each probe separately, with beta-values standardized across samples, such that the effect sizes corresponded to a one standard deviation increase in DNAm. All models were adjusted for age, race, stage, menopause status, tumor purity, and cell type proportion. For each model, the genome-wide significance level of 1 × 10 − 7 was established using Bonferroni correction to account for multiple testing of 485,512 probes. Model A tested the association of beta-values and the survival outcome stratified by the tumor subtypes. Model B assumed the association between beta-values and the survival outcome could depend on tumor subtype by integrating a product term between subtype and site-specific DNAm. Thus, under this model, we were able to compute the estimates of log hazards and their P -value for each of the subtype. Model C assumed the association between beta-values and the survival outcome could depend on menopausal status (post-menopausal vs. pre- and peri-menopause), by including a product term between menopausal status and DNAm. As in Model B, we were able to compute the estimates of log hazards and their P -values for each menopausal status. During the model fitting, a small number of probes produced warnings due to a lack of convergence. These probes were excluded from further analysis.

We used the functional class scoring (FCS) robust rank aggregation approach implemented in the methylGSA (Version 1.6.1) R package (Ren and Kuan, 2019) for the pathway analysis. This approach uses the site-specific EWAS inference to identify gene sets enriched among probes that are associated with survival outcomes. Two types of biases should be corrected for functional pathway analysis with DNAm data. One is the length bias arising from differing lengths of pre-defined gene sets based on molecular and cellular functions, which is addressed by FCS via Gene Set Enrichment Analysis (Subramanian et al., 2005). The other is the probe bias arising from different numbers of probes mapped to gene sets of equal size, which is addressed by the robust rank aggregation (RRA). RRA computes combined p-value based on the order statistic of p-values of probes (from site-specific analyses) mapped to a gene, assuming independent identical uniform distribution of p-values. We restricted our evaluation of enrichment to gene sets containing between 100 and 500 genes. The adjusted p-values of gene sets (enriched pathways) account for multiple hypothesis testing and the size of gene sets. We defer readers to (Subramanian et al., 2005) for details on calculation of adjusted p-values.

The ipDMR approach implemented in the Enmix (Version 1.25.1) R package (Xu et al., 2020) was used to identify differentially methylated regions (DMRs). The initial DMRs are identified as a region bound by a pair of adjacent probes across the methylome. These regions are filtered based on their sizes and combined P -values, each of which is computed as some function of p-values of probes (from site-specific analyses) within a region. We defer readers to (Xu et al., 2020) for details on the form of this function that computes the combined p-values. The initial DMRs go through another round of filtering and merging regions based on their sizes, distances, and adjusted combined P -values. A DMR was deemed significant if its adjusted combined P -value was less than 0.05.

We used the Illumina 450K array data of 295 BC tumor tissue samples available in GEO database (GSE72308) to validate site-specific analysis results. The samples were obtained from retrospectively selected frozen samples and prospective cohorts who were treated with adjuvant and neoadjuvant therapies at Jules Bordet Institute from 1995 to 2009 (Jeschke et al., 2017). IHC-based tumor subtypes available in the dataset were taken as the subtypes of the samples. To be consistent with site-specific analysis of TCGA data, we removed samples of HER2 subtype. OS endpoint was the only available outcome that was in common with survival outcomes available in TCGA data. Age and tumor subtype were the only available covariates in the data set. Tumor purity and cell type proportion were estimated using identical methods in the site-specific analysis. The beta-values of each probe were normalized across samples before conducting EWAS. Since menopausal status information was not available, Model C was not fitted. Model A and Model B were identical to the ones fitted using TCGA data except for the covariates such as race, stage, and menopausal status, which were not available in the GEO data set. EWAS results from each model were compared to the results from TCGA site-specific analyses. In particular, log hazards of the probes were compared to examine whether the directions of the association between beta-values to survival outcomes were similar between TCGA and GEO EWAS results.

The association between site-specific DNAm and each clinical outcome was assessed across methylome. Time to event was modeled using Cox proportional hazards (PH) models, stratified by tumor subtype, and adjusted for age, race, stage, menopausal status, DNAm-estimated tumor purity and cell proportions (8 cell types). Based on this model, one probe was associated with OS at genome-wide significance (p < 1 × 10 − 7 ; Figure 1; Table 2, Supplementary Figure S1). Effect measure modification by tumor subtype or menopausal status was evaluated by further incorporating product terms between these characteristics and DNAm. Three additional luminal A subtype-specific probes met genome-wide significance. These probes were significantly associated with OS only after accounting for effect measure modification by tumor subtype. For all four probes, an increase in DNAm was associated with a lower hazard of all-cause mortality. All four probes were located in open sea regions and none of the probes were located within promoter regions (Table 2). Kaplan-Meier (KM) plots for the probes presented in Table 2 are presented in Figure 2, stratified by both tumor subtype and menopause status, and with site-specific DNAm dichotomized by the median. These results were consistent with the direction of associations estimated from the site-specific Cox PH models.

Similar models were used to identify site-specific associations with PFI. A majority of genome-wide significant sites were found to be luminal A, luminal B and/or post-menopausal specific (Figure 1; Figure 2; Table 2, Supplementary Figure S1). There were 10 probes that were significantly associated with PFI only after accounting for effect measure modification by tumor subtype or menopausal status. Higher DNAm was associated with a reduced hazard of disease progression for all but two probes, cg17735983 and cg10678486, both of which were specific to the luminal A subtype (Table 2). While most of the probes were in open sea regions and gene bodies, cg17735983 was located in an CpG island, and cg10678486 was located within 200 bp of the transcription start site (TSS200) of ELAC1. Probe cg09926728 was associated with PFI overall, but this association was strongest among participants with the luminal A subtype and postmenopausal women (Supplementary Figure S5I). The estimated hazard ratios (HRs) for this probe were similar across different models ranging between 0.57–0.59. Similarly, higher methylation at cg16976520 was associated with a reduced hazard of progression overall, which was strongest among the post-menopausal women (Figure 2F; Supplementary Figure S5M). Among post-menopausal women, the rates of progression were higher among women with the luminal A and luminal B subtypes with lower methylation at cg16976520, relative to those with the basal-like subtype (Figure 2F).

We evaluated gene ontology (GO) enrichment among the site-specific associations with BC outcomes, taking into account the length bias of a gene and the number of probes per gene (Supplementary Figures S2 and Supplementary Figure S3). For each endpoint, GO enrichment was assessed based on the overall model, as well as models incorporating effect modification by either subtype or menopausal status. The number of significantly (adjusted p-value ≤ 0.05 ; Materials and methods section) enriched pathways ranged from 7 to 155 for OS, and from 0 to 217 for PFI. While the ten most significantly enriched GOs differed between each model, there were some overlap between significantly enriched GOs by different prognosis outcomes and different strata. The GOs related to axon development were associated with OS among both post-menopausal, and pre- and peri-menopausal women. These GOs were also associated with PFI for both luminal A subtypes and post-menopausal women. The GOs related to insulin regulation (GO:0050796 and GO:0030073) were significantly associated (adjusted p-value ≤ 0.05 ) with OS and suggestively associated (adjusted p-value ≈ 0.1 ) with PFI among pre- and peri-menopausal women. The GOs related to regulation of neuron death were associated with PFI for the luminal B subtype.

To identify DMRs associated with our outcomes, we applied an approach that integrates site-specific inference between proximal loci. Significant DMRs differed by tumor subtype and menopausal status (Table3; Table 4; Figure 3, Supplementary Table S1, Supplementary Figure S4). These included a DMR within the promoter of ELAC1 that was associated with PFI among the luminal A tumor samples, but not among the luminal B and basal-like samples (Figure 3A). In particular, this DMR was the most significant among all DMRs containing at least four probes, and included a genome-wide significant CpG probe, cg10678486, which was identified in the site-specific analysis to be associated with PFI for women with luminal A tumor subtype. We identified 271 luminal A specific DMRs that were associated with OS (Table 3). Only one out of three genome-wide significant luminal A subtype-specific probes associated with OS appeared within these significant DMRs. Similarly, only one of the three probes significantly associated with PFI was within one of the 42 DMRs associated with PFI (Table 3). These results indicated a majority of DMRs were not driven by one strong site-specific association, but several weaker associations in the same direction among proximal loci (Table 4). Several DMRs were not close to known genes. For example, an identified DMR specific to luminal A subtype did not reside near any known gene, but in a region indicated by chromatin state of promoter upstream TSS (Supplementary Table S1, Supplementary Figure S4F).

Validation of results of TCGA site-specific analyses was carried out on Illumina 450K array data profiled from BC tissue of 180 patients available in GEO database (GSE72308) (Jeschke et al., 2017). We were only able to analyze associations with OS since GEO dataset did not include any other common type of survival outcome. These validation models did not adjust for race and stage, or evaluate effect modification by menopausal status as these characteristics were missing from the GEO data. Tumor purity and cell proportion of samples were estimated and included in the models. The age of a majority of subjects was less than 60 years old in both TCGA and GEO. At alpha level of 0.05, there were no significant age differences between the two datasets (p-value = 0.148) based on binary categorization of the age variable. More than 50 percent of tumor subtypes of TCGA samples was luminal A, whereas the subtypes of GEO samples were relatively balanced across luminal A, luminal B and basal-like subtypes (Table 1). The median follow-up times for all-cause mortality were not significantly different (p-value = 0.789) between the two datasets at alpha level of 0.05. We compared the top 20 most significant (arranged by increasing order of p-values) associations with OS identified from TCGA site-specific analyses to those estimated among the GEO dataset. While the directions of associations were similar, most did not reach statistical significance, potentially due to the smaller sample size of the GEO dataset. Associations with similar magnitudes between the two datasets include the overall association with DNAm at cg05497253 (Figure 4A), as well as the luminal A subtype-specific associations at cg05860556 and cg13328771 (Figure 4B), and the basal-like-subtype specific associations at cg04706201 and cg14253517 (Figure 4D).

Further, we compared p-values for the 20 most significantly enriched GOs identified from the TCGA site-specific analysis results to the p-values of matching GOs identified from the GEO site-specific analysis results (Figure 5). Several GOs seemed to be consistent across datasets. Specifically, for the luminal A subtype-specific associations, GO:0031012 (extracellular matrix) and GO:0042383 (sarcolemma) showed consistent results between TCGA and GEO (Figure 5B). For the luminal B subtype-specific associations, GO:0007389 (pattern specification process), GO:0001216 (DNA-binding transcription activator activity), and GO:0001228 (DNA-binding transcription activator activity) showed consistent results between TCGA and GEO (Figure 5C). For the basal-like subtype-specific associations, GO:0007269 (neurotransmitter secretion) and GO:0099643 (signal release from synapse) showed consistent results across datasets (Figure 5D).