Later-generation epigenetic aging clocks outperform first-generation models in predicting survival in TCGA breast cancer.

Introduction
Epigenetic aging captures the divergence between an individual’s biological state and their chronological age. DNA methylation (DNAm) based “epigenetic clocks” quantify this process by modeling age from methylation patterns across the genome [1]. First-generation clocks, exemplified by the Horvath and Hannum models, were optimized to predict chronological age from DNAm [2, 3]. They are highly accurate for dating tissues but were originally not designed to capture morbidity or mortality risk [4]. As a result, their prognostic utility in disease settings can be limited because much of their signal aligns with the known effects of chronological age. To address these limitations, later clocks incorporated biology linked to healthspan and mortality. The Levine “PhenoAge” clock integrates DNAm surrogates of clinical chemistry and hematologic markers to approximate phenotypic aging, improving associations with morbidity and survival [5]. The GrimAge family goes further by combining DNAm surrogates for plasma proteins and smoking exposure [6]; the updated GrimAge2 improves the mortality signal and calibration relative to the original GrimAge1 [7]. Related methylation metrics such as DNAm-estimated telomere length (DNAmTL) and indices tied to stem-cell division (e.g., pcgtAge) provide complementary views on cellular aging, replicative history, and systemic stress [8, 9]. Taken together, these metrics bridge age estimation and mortality risk, positioning them for rigorous assessment in oncology cohorts.
Breast cancer is heterogeneous, with outcomes shaped by hormone receptor status, human epidermal growth factor receptor-2 (HER2) amplification, tumor stage, and patient age at diagnosis [10]. The Cancer Genome Atlas (TCGA) provides one of the largest resources linking tumor tissue-based DNAm profiles to curated clinical data in breast cancer, generating a TCGA-BRCA dataset [11]. Although individual clocks have been examined in various cancers, there remains a need for a systematic, side-by-side evaluation of first-, second-, and third-generation DNAm clocks for prognostication in TCGA-BRCA, while accounting for established clinicopathologic determinants. We previously performed such a systematic analysis in colorectal cancer using TCGA-COAD and TCGA-READ, highlighting both the promise and the caveats of DNAm clocks in cancer prognosis [12, 13]. Whether similar patterns hold in breast cancer has not been clarified.
Here, we undertake a comprehensive assessment of multiple DNAm aging metrics in TCGA-BRCA. We compared first-generation clocks (Horvath, Hannum), second-generation PhenoAge (Levine), and the GrimAge family including GrimAge2, alongside DNAmTL and a stem-cell division proxy for their effects on breast cancer survival. Our goal is to determine which DNAm clocks, if any, provide prognostic information beyond standard clinical factors in breast cancer, and to delineate which aging measurement matrices of a clock (e.g., raw value, age-acceleration residual, or age difference) are most informative.

Materials and methods

Study population
We retrieved the DNAm dataset of 1076 breast cancer patients from TCGA using the R package TCGAbiolinks. Subsequently, the clinical records were retrieved from the TCGA-BRCA “Clinical Supplement (BCR Biotab)” files and merged with curated survival data and immunohistochemistry (IHC) biomarkers (i.e., estrogen receptor/progesterone receptor [ER/PR], and HER2). Patient records were then harmonized by submitter ID and de-duplicated to retain one record per individual. Using an inner join across the DNAm–matched phenotype table, the IHC dataset, and the survival file, we eventually refined the analytic cohort of 781 unique patients for downstream analyses.
We next programmatically derived clinical annotations (Table 1). Menopausal status was defined as post-menopausal (n = 496) or pre/peri-menopausal (n = 195); unknown (n = 90) was retained for descriptive summaries. Receptor subtype was assigned from IHC as Luminal A (n = 309), Luminal B (n = 76), HER2-enriched (n = 17), or Triple-negative breast cancer (TNBC) (n = 83); cases without sufficient IHC were labeled Unknown (n = 296) and excluded from subtype-specific analyses. Tumor stages were harmonized by collapsing AJCC sub-stages to Stage I (n = 126), Stage II (n = 439), Stage III (n = 194), Stage IV (n = 11), with Unknown (n = 11). Demographic variables followed TCGA fields for ethnicity—Hispanic or Latino (n = 38), Not Hispanic or Latino (n = 671), Unknown (n = 72)—and race—White (n = 570), Black or African American (n = 158), Asian (n = 37), American Indian or Alaska Native (n = 1), Unknown (n = 15). “Unknown/not reported” categories were kept in cohort counts and excluded from subgroup comparisons requiring that attribute (Table 1).

Preparation of DNAm data: global level of DNAm array
Global level DNAm data was obtained from TCGA. DNA samples derived from breast tissues were analyzed via the Illumina Infinium450K array and normalized via normal-exponential out-of-band (Noob) background correction by minfi [14]. By using RefFreeEWAS V.2.2 [15], cell-type proportions (cancer and normal epithelial, stromal, and immune cells) of each tumor sample were estimated and adjusted for the heterogeneities in estimating DNAm age. All DNAm clock scores were then computed using the original published coefficients; no model parameters were re-trained, recalibrated or optimized using TCGA-BRCA survival outcomes.

Survival analysis with demographic factors and cancer characteristics
Overall survival (OS) was evaluated with the survival and survminer R packages (functions survfit and ggsurvplot). The event indicator was derived from TCGA vital status (Alive = 0; Deceased = 1). To limit undue influence of extreme follow-up, OS time was truncated at 4000 days (OS.time_4000 = pmin(OS.time, 4000)); individuals whose follow-up exceeded 4000 days were censored at 4000 days (event = 0 at that time point). Kaplan–Meier curves were plotted with two-sided log-rank p-values, 95% confidence (95% CI) bands, and risk tables [16].
Groupings for Kaplan–Meier comparisons were done according to the programmatic annotations described in Study Population in Table 1. Specifically, we compared: (i) menopause (pre/peri vs post), (ii) receptor subtypes in prespecified pairwise contrasts (TNBC vs Luminal A; TNBC vs Luminal B; Luminal A vs Luminal B; and HER2-enriched vs each luminal subtype) and collapsed contrasts (Luminal A&B vs TNBC; and Luminal A&B vs HER2-enriched), (iii) tumor stage contrasts (Stage 4 vs others; Stage 1 vs others; Stage 2 vs Stage 3&4; and Stage 3 vs Stage 4), and (iv) demographic groups (by race: Asian vs White; Black or African American vs White; by ethnicity: Hispanic or Latino vs Not Hispanic or Latino). For race-by-ethnicity KM plots, analyses were repeated within each race stratum (White only; Black or African American only; and Asian only). Entries labeled “Unknown” or “not reported” for a given grouping were treated as missing and excluded from that specific comparison.

Cox regression analysis with demographic factors and cancer characteristics
We modeled OS with Cox proportional hazards using the survival package (coxph). The event‐time object was Surv(OS.time_4000, deceased_4000), where OS was truncated at 4000 days and participants followed longer were censored at day 4000.
For univariate analyses, we fit separate models for each predictor. Age at diagnosis was modeled as a continuous covariate, and hazard ratios (HRs) were interpreted per 1-year increase. Menopausal status was coded as post-menopausal compared with pre/peri-menopausal. For race, White served as the reference and was compared with Black or African American and Asian. For ethnicity, Not Hispanic or Latino served as the reference and was compared with Hispanic or Latino. For receptor subtype, Luminal A served as the reference and was compared with Luminal B, HER2-enriched, and TNBC. For tumor stage, Stage I served as the reference and was compared with Stages II–IV. Categorical variables were explicitly re-leveled to enforce these reference categories. For each model, we reported HRs, 95% CIs, and Wald p values, and displayed the results as log-scale forest plots with a reference line at HR 1.
For the multivariable analysis, we fit a combined model including age at diagnosis, menopausal status, race, and ethnicity, receptor subtype, and tumor stage. Prior to modeling, we excluded records with unknown receptor status, unreported tumor stage, race recorded as not reported or American Indian or Alaska Native, and ethnicity recorded as not reported; unused factor levels were dropped. HRs with 95% CIs and two-sided p values (α = 0.05) were summarized and presented in a forest plot.

DNAm-based age analysis in association with OS
We evaluated OS in TCGA-BRCA using DNAm–based aging metrics. Categorical covariates were encoded identically to those used in the survival and Cox analyses. OS time and event were defined by OS.time_4000 and deceased_4000.
For each clock metric, in addition to aging scale, we estimated deviation from chronological age, such as DNAmDiff by subtracting chronological age from DNAmAge and age acceleration (AgeAccel) by taking the residual from regressing DNAm age on chronologic age. For Horvath and Hannum clocks, we further derived intrinsic epigenetic age acceleration (IEAA), which represents estimates from DNAm age adjusted for variation in blood cell-type composition and reflects cell-intrinsic aging effects.
Specifically, we analyzed Horvath DNAmAge together with IEAA and DNAmDiff (Horvath DNAmAge minus chronological age); Hannum DNAmAge together with IEAA.Hannum (intrinsic epigenetic age acceleration for the Hannum clock) and DNAmDiff.Hannum (DNAmAgeHannum minus chronological age); and Levine DNAmPhenoAge together with AgeAccelPheno (the residual from regressing DNAmPhenoAge on chronological age) and DNAmDiff.Pheno (DNAmPhenoAge minus chronological age). For GrimAge (both predicted-age GrimAge1 and real-age GrimAge2), we likewise derived an age-acceleration residual and an age-difference measure. DNAm teleomere length (DNAmTL) and the stem-cell division proxy pcgtAge were used on their native scale.
We generated Kaplan–Meier curves under two grouping schemes: a median split (into two groups) and tertiles (into three groups). Curves were compared with log-rank tests; 95% CI bands and risk tables. Univariate Cox proportional hazards models were fit separately for each DNAm aging metric to estimate HRs per unit increase of aging scale, with 95% CIs and Wald p-values. Multivariable Cox models were constructed separately for each DNAm clock family. For the Horvath clock, we fit models with DNAmAge, the age-acceleration residual IEAA, and the age-difference metric DNAmDiff. For the Hannum clock, we used DNAmAgeHannum, the age-acceleration residual IEAA.Hannum, and the age-difference metric DNAmDiff.Hannum. For the Levine clock, we used DNAmPhenoAge, the age-acceleration residual AgeAccelPheno, and the age-difference metric DNAmDiff.Pheno. For GrimAge2, we evaluated the predicted-age and real-age versions, together with their age-adjusted residuals and age-difference metrics. All models were adjusted for menopausal status; age at diagnosis; receptor subtype with Luminal A as the reference; tumor stage comparing Stages II–IV with Stage I while excluding cases coded as Others;race comparing Black or African American and Asian with White while excluding those who did not reported and American Indian or Alaska Native; and ethnicity comparing Hispanic or Latino with Not Hispanic or Latino. Each of the three representations for a given clock was fit in a separate multivariable model.

Results

Kaplan–Meier analysis: demographic features associated with patients’ survival
To establish baseline survival differences in the TCGA-BRCA cohort, we evaluated demographic and clinical covariates using Kaplan–Meier analyses. Post-menopausal patients exhibited significantly shorter overall survival than pre/peri-menopausal patients (log-rank p = 0.0037; Fig. 1a). Ethnicity showed a modest association, with non-Hispanic/Latino patients experiencing worse survival than Hispanic/Latino patients (p = 0.046; Fig. 1b). By contrast, race alone (Asian and Black/African American in comparison with White) was not associated with survival (p ≥ 0.46; Fig. 1c and 1d). Stratification by combined race and ethnicity categories did not reveal additional differences: White Hispanic/Latino versus White non-Hispanic/Latino showed only a trend (p = 0.071), and comparisons within Black/African American and Asian groups were non-significant (p ≥ 0.60; Fig. S1). Together, these results indicate that menopausal status is a primary demographic determinant of survival in this cohort, whereas race and ethnicity show limited or no independent associations.

Kaplan–Meier analysis: receptor status associated with patients’ survival
We next evaluated baseline receptor-defined subtypes in relation to OS using Kaplan–Meier curves and log-rank tests (Fig. 2, Fig. S2). Luminal A showed markedly better survival than TNBC (p = 0.0019; Fig. 2a), and the combined luminal group (A + B) also outperformed TNBC (p = 0.006; Fig. 2c), whereas Luminal B alone did not differ from TNBC (p = 0.42; Fig. 2b). Luminal A had superior survival compared with all other subtypes combined (p = 0.0027; Fig. 2d), while Luminal B versus others was not significant (p = 0.47; Fig. 2e). Consistent with these contrasts, TNBC had worse outcomes than the remaining subtypes (p = 0.0091; Fig. 2f). Pairwise luminal comparisons were nonsignificant (Luminal A vs B, p = 0.12; Fig. S2a). HER2-enriched cases showed no detectable difference versus any comparator (all p > 0.05; Fig. S2b–f), noting the small HER2-enriched sample. Collectively, these results position Luminal A as the most favorable and TNBC as the least favorable survival phenotype, with Luminal B intermediate.

Kaplan–Meier analysis: cancer Clinical stage associated with patients’ survival
Using Kaplan–Meier curves with log-rank tests, we observed a clear stage–survival gradient in TCGA-BRCA. Stage 1 disease showed superior overall survival compared with Stages 2–4 combined (p = 0.014; Fig. 3a). Stage 2 also outperformed aggregated advanced disease (Stages 3 and4 combined; p = 2.9 × 10⁻4; Fig. 3b). Within the advanced disease state, Stage 3 had better survival than Stage 4 (p < 0.0001; Fig. 3c). Consistently, Stage 4 had the poorest prognosis relative to Stages 1–3 combined (p < 0.0001; Fig. 3d). Together, the curves indicate progressively worse outcomes with increasing stage, with the steepest decline in Stage 4. Of note, times were truncated at 4000 days; cases with missing/indeterminate stages were excluded from stage-specific contrasts.

Results from inferential analyses for demographic and clinical factors in association with survival
To complement the Kaplan–Meier analyses, we fit Cox proportional-hazards models (univariate and multivariable) using pre/peri-menopause, Luminal A, Stage I, White race, and Not Hispanic or Latino as reference categories. Age at diagnosis was modeled as a continuous variable. Forest plots with 95% CIs are shown in Fig. 4a, b. In the univariate analysis, post-menopausal status associated with worse OS (HR 2.26, 95% CI 1.3–4.0, p = 0.0048), but the effect disappeared after adjustment (HR 1.02, 95% CI 0.29–3.6, p = 0.98), consistent with confounding by tumor features. Relative to Luminal A, TNBC showed a robust increase in hazard in both models (univariate HR 2.55, 95% CI 1.3–4.9, p = 0.0043; multivariable HR 3.06, 95% CI 1.2–7.5, p = 0.015), whereas Luminal B and HER2-enriched were not significant. Increases in stage strongly tracked with mortality: Stage III and Stage IV remained significant after adjustment (multivariable HR 4.20, 95% CI 1.5–12, p = 0.007; HR 120.25, 95% CI 16–910, p = 3.6 × 10−6, relatively), while Stage II was not. Race and ethnicity showed no association with OS after adjustment (p > 0.1). Older age at diagnosis was consistently detrimental, with an ~ 8% higher hazard per year in the adjusted model (HR 1.08, 95% CI 1.0–1.1, p = 4.7 × 10⁻5). Together, these results reinforce the KM findings and indicate that TNBC status, advancing stage, and older diagnostic age are the strongest correlates of poorer overall survival in this cohort although the very wide CI for Stage IV reflects a small sample size and warrants cautious interpretation.

Kaplan–Meier analysis: survival analysis of DNAm clocks
We next used Kaplan–Meier analyses to assess the prognostic value of different DNAm clocks in TCGA-BRCA. Patients were stratified into a binary group by the median value of each clock (Fig. 5; Fig. S3a–c) or into three groups by tertiles (Fig. S3d–j). In addition, models based on DNAmAge residuals and age acceleration measures derived from Horvath, Hannum, Levine, GrimAge1 and GrimAge2 clocks were further examined using median-based dichotomization in Kaplan–Meier analyses (Fig. S3k-t). First-generation clocks (Horvath, Hannum) showed no difference in overall survival between high and low DNAmAge groups (log-rank p = 0.98 and p = 0.55; Fig. 5a,b). By contrast, higher values from the second generation clocks were associated with worse outcomes: DNAmPhenoAge (Levine) and GrimAge2 each stratified patients into distinct risk groups (Fig. 5c, Levine: p = 0.017 and Fig. 4d, GimAge2: p = 3.9 × 10⁻4). GrimAge1, DNAmTL, and the stem-cell division proxy (pcgtAge) did not yield significant separations (Fig. S3a–c).

Consistent results were observed with tertile stratification: only GrimAge2 produced a graded, statistically significant decline in survival across increasing tertiles (Fig. S3h), whereas other clocks did not show significant differences (Fig. S3d–g,i,j). These findings indicate that later-generation clocks, particularly GrimAge2, provide stronger prognostic discrimination in breast cancer than earlier clocks.
Interestingly, analyses based on age‐acceleration metrics revealed the opposite trend in several models. Specifically, patients in the high acceleration groups defined by the Horvath, Hannum, and GrimAge1 clocks displayed significantly improved overall survival compared with their low-acceleration counterparts (Fig. S3l, n, r). This observation is conceptually consistent with the notion that tumors with more advanced epigenetic aging may be more susceptible to therapeutic eradication, thereby yielding more favorable outcomes.

Results from inferential analyses for DNAm clocks in association with survival
We further evaluated associations between each DNAm clock and OS using Cox proportional hazards models. For clocks except DNAmTL and stem cell clocks, we independently analyzed three measurement metrics—epigenetic clock value, age‐adjusted residual, and age difference (Table 2). Multivariable models adjusted for menopause status, age at diagnosis, receptor subtype, tumor stage, race, and ethnicity, showing results in forest plots (Fig. S4). In univariate analyses, several clocks as continuous variable showed significant associations with survival. DNAmTL was inversely associated with risk (HR 0.54, 95% CI 0.33–0.89, p = 0.0147). Second‐generation clocks were positively associated with risk: Levine (HR 1.01, 95% CI 1.00–1.01, p = 0.037), GrimAge1 (HR 1.01, 95% CI 1.00–1.03, p = 0.0365), GrimAge2 (HR 1.04, 95% CI 1.02–1.05, p = 2.62 × 10−6), and the GrimAge2 residual (HR 1.04, 95% CI 1.01–1.06, p = 0.0032). Horvath and Hannum clocks, as well as their residual and age‐difference metrics, were not significant. After multivariable adjustment, only GrimAge1 retained independent prognostic value (raw: HR 1.03, 95% CI 1.01–1.06, p = 0.0159; residual: HR 1.03, 95% CI 1.00–1.06, p = 0.0257; age difference: HR 1.03, 95% CI 1.01–1.06, p = 0.0159). Associations for GrimAge2 and Levine attenuated and were no longer significant. In addition, DNAmTL kept showing a negative association with survival, but the trend did not reach significance (HR 0.48, 95% CI 0.20–1.13, p = 0.0922). Collectively, while second‐generation clocks capture prognostic signal in unadjusted analyses, only GrimAge1 provides risk information that remains statistically significant, independent of clinicopathologic covariates in TCGA-BRCA.

Discussion
In this study, leveraging the TCGA-BRCA cohort, we systematically evaluated seven epigenetic aging metrics: the Horvath, Hannum, and Levine clocks; GrimAge1 and GrimAge2; a DNAm-based stem-cell division age; and DNAmTL. A key strength of TCGA is the availability of tumor tissue-level DNAm profiles linked to curated clinical annotations, providing a complementary perspective to population cohorts that primarily analyze blood.
Our first examination confirmed expected clinicopathologic survival patterns in TCGA-BRCA. Menopausal status was associated with overall survival in unadjusted KM analyses, with postmenopausal patients showing higher hazard ratios and lower survival probabilities. Race was not associated with survival, whereas ethnicity was—Hispanic/Latino patients experienced significantly better survival than non-Hispanic/Latino patients. Among receptor-defined subtypes, TNBC had the worst overall survival, whereas Luminal A had the most favorable outcomes. Stage-based analyses aligned with clinical expectations, showing progressively poorer outcomes at later stages. Age at diagnosis was inversely associated with survival in Cox models. These observations are consistent with prior reports [17–19] and indicate that our analyses reproduced the expected clinicopathologic survival patterns in TCGA-BRCA, supporting the robustness of our pipeline.
Using DNA methylation data from the TCGA-BRCA cohort, Horvath’s clock is the earliest DNAm-based epigenetic clock and was trained across multiple tissues. Hannum’s clock was developed in blood. Together they are considered first-generation models. In our survival analyses, using continuous values, median splits or tertiles, neither clock measures significantly separated overall survival. Interestingly, in univariate survival analysis, the Horvath age difference was associated with a lower hazard (HR 0.87, p = 0.0015), suggesting that in this cohort, higher Horvath’s aging acceleration phenotype corresponded to slightly better outcomes, warranting further validation studies.
Second-generation clocks anchor DNAm CpG features to clinical indicators of morbidity and mortality. The Levine PhenoAge clock integrates DNAm surrogates of composite clinical phenotypic measures and shows improved prognostic performance. GrimAge1 uses a two-stage approach by first deriving DNAm-based surrogates for plasma proteins and smoking exposure, then regressing time-to-death on these surrogates to produce a mortality predictor. GrimAge2 updates the model by adding DNAm estimators for high-sensitivity C-reactive protein and hemoglobin A1c. In TCGA-BRCA, PhenoAge significantly separated survival between high and low groups, and higher PhenoAge was associated with increased mortality in univariate Cox models. Likewise, GrimAge2 stratified survival by median split and showed a strong univariate association with mortality. Together, these findings indicate that second-generation clocks perform better prognostic discrimination than first-generation age-prediction clocks.
Beyond the traditional clocks, Yang et al. proposed a DNAm-based metric approximating a mitotic clock in normal and cancer tissues (epiTOC or pcgtAge in our study) [9]. In the TCGA-BRCA cohort, tertile stratification by epiTOC age revealed a significant survival difference, whereas a median split did not. Although Cox regression analyzing it as a continuous variable did not reach statistical significance with wide confidence intervals, the point estimate suggested a markedly elevated hazard (HR > 4), implying a substantial role in survival prediction in breast cancer patients.
Telomere length is a canonical marker of biological aging and its shortening has been linked to aging and increased cancer risk [20]. In TCGA-BRCA, tertile stratification of DNAmTL showed a significant inverse association with overall survival in which shorter DNAmTL corresponded to poorer outcomes. Our findings are consistent with prior breast cancer studies reporting that short telomeres are enriched in aggressive clinicopathologic subsets and correlate with adverse features (21). In our univariate Cox models, each unit increase in DNAmTL was associated with lower mortality risk, indicating that shorter DNAmTL predicts worse survival.
Although this analysis has the notably strength by leveraging tumor tissue-delivered DNAm profiles from the TCGA dataset, which provides a disease relevant methylome rather than relying on peripheral blood, several limitations should be acknowledged. First, minority racial/ethnic subgroups are under-represented in TCGA-BRCA, reducing statistical power for these subgroup comparisons. Second, clinical and lifestyle epidemiologic annotation is limited and partly dated (e.g., a lack of information about treatment, comorbidity, lifestyle, and socioeconomic variables), constraining multivariable adjustment and leaving room for residual confounding. Finally, this is a single, retrospective dataset with missingness and potential batch effects, warranting external validation. These limitations may attenuate or inflate effect sizes; thus, findings should be cautiously interpreted.
In conclusion, we assessed multiple DNAm-based aging metrics for prognostication in breast cancer in a comprehensive fashion within the TCGA-BRCA data. Second-generation clocks demonstrated stronger prognostic signal than first-generation models in unadjusted analyses. Among them, GrimAge1 retained independent prognostic value after multivariable adjustment. DNAm-based telomere and mitotic clocks showed potential for discriminating survival risk. Prospective validation and integration with richer clinical covariates will be important next steps.

Supplementary Information
Below is the link to the electronic supplementary material.