Estimating DNA methylation age with methylclock
A user’s guide: from a methylation matrix to interpretable clock estimates
2026-09-10
Package
methylclock 1.99.1
Contents
- 1 Introduction
- 2 Quick start
- 3 Which clocks are available, and what they are for
- 4 Getting started
- 5 Running the clocks and reading the output
- 6 Seeing and interpreting the results
- 7 Age acceleration
- 8 Comparing groups
- 9 Quality control: does a clock fit your data?
- 10 Reference ranges: your samples against a named cohort
- 11 Choosing the right clock for your data
- 12 Longitudinal data: trajectories and slopes
- 13 Beyond age: a trait predictor
- 14 EPICv2 arrays
- 15 Session info
- References
1 Introduction
DNA methylation — a chemical mark on the genome — changes with age in a remarkably orderly way. Some CpG sites gain methylation as we get older, others lose it, and the pattern is regular enough that a model reading anywhere from a handful to hundreds of thousands of these sites can estimate a person’s age to within a few years. Such a model is an epigenetic clock.
Since the first clocks in 2013 [1,2], many have been published. They do not all measure the same thing: some estimate calendar (chronological) age, others estimate a biological age tied to health and mortality, others the age of a fetus, the number of divisions a cell lineage has gone through, or the rate at which someone is ageing. Because they answer different questions, it is common to compute several and compare them.
methylclock does exactly that. You give it a methylation matrix; it returns each clock’s estimate and a small toolbox to plot, compare and interpret the results. This guide walks through the whole path: which clocks are available and what each is for, what data to provide, how to run them, and — the part that matters most — how to read what comes back.
library(methylclock)2 Quick start
For the impatient: install (Section 4), hand
methylclock() a beta matrix, and plot.
library(methylclock)
res <- methylclock(betas) # betas: CpGs x samples, values in [0, 1]
plotDNAmAge(res, age = pheno$age)Everything below unpacks that: which clocks exist and how to choose one, what the input must look like, and how to read what comes back. The figures in this guide run on bundled example estimates, so you can reproduce them before touching your own data.
3 Which clocks are available, and what they are for
Before running anything, it helps to know what you can ask for. The package groups
clocks by what they measure. The tables below are generated from the package
itself, so they always match the version you have installed. You can reproduce
any of them, or filter differently, with clock_list().
Each table lists, for every clock: its name (what you pass to
methylclock(..., clocks =)), what it estimates and in which units, the
tissue it was trained on (a clock reads best the tissue it came from:
a buccal-trained clock applied to blood answers a question nobody asked), the
arrays it supports (a clock listed as EPIC only will not apply to 450K
data), the number of CpGs it uses, and a numbered reference to the paper
the clock comes from — follow the number to the reference list at the end of
the vignette. All of it is also available programmatically from clock_list(),
clock_info() and clock_catalog.
3.1 Chronological-age clocks
These estimate calendar age in years, and are the workhorses: use them to
predict age, or to derive age acceleration (Section 7).
Horvath’s clock works across almost any tissue; Hannum’s, BLUP and elastic net
(EN) are for blood; the skin & blood clock suits skin and blood samples;
AltumAge and BNN are neural-network clocks; Wu and PedBE are tuned for
children, where the ageing signal is steepest; and Garma was trained
only on probes shared by the 450K, EPICv1 and EPICv2 arrays, so it is the one
clock here that runs natively on EPICv2 data (Section
14).
| Clock | Estimates | Tissue | Arrays | CpGs | Ref |
|---|---|---|---|---|---|
| AltumAge | Age (years) | multi-tissue | 27K/450K/EPIC | 20318 | [3] |
| BLUP | Age (years) | whole blood | 450K | 318395 | [4] |
| BNN | Age (years) | whole blood | 450K | 353 | [5] |
| EN | Age (years) | whole blood | 450K | 514 | [4] |
| Garma | Age (years) | whole blood | 450K/EPIC/EPICv2 | 4962 | [6] |
| Hannum | Age (years) | whole blood | 450K | 69 | [2] |
| Horvath | Age (years) | multi-tissue | 27K/450K | 353 | [1] |
| Lin | Age (years) | whole blood | 450K | 99 | [7] |
| NEOaPMA450K | Age (years) | buccal (preterm infants) | 450K | 408 | [8] |
| NEOaPMAEPIC | Age (years) | buccal (preterm infants) | EPIC | 520 | [8] |
| NEOaPNA450K | Age (years) | buccal (preterm infants) | 450K | 302 | [8] |
| NEOaPNAEPIC | Age (years) | buccal (preterm infants) | EPIC | 508 | [8] |
| PedBE | Age (years) | buccal epithelium (children) | 450K | 94 | [9] |
| skinHorvath | Age (years) | skin and blood | 450K/EPIC | 391 | [10] |
| VidalBralo | Age (years) | whole blood | 27K/450K | 8 | [11] |
| Weidner | Age (years) | whole blood | 27K/450K | 3 | [12] |
| Wu | Age (years) | whole blood (children) | 27K/450K | 111 | [13] |
3.2 Biological-age clocks
Where a chronological clock asks how old are you, a biological-age clock asks
how old does your body look. Levine (PhenoAge) is trained on clinical
biomarkers of ageing, so a sample reading older than its calendar age is flagged
as biologically aged — the basis of age acceleration. TL estimates telomere
length, which shortens with age.
3.3 Gestational-age clocks
For samples taken before or at birth (placenta, cord blood), these estimate gestational age in weeks rather than years.
| Clock | Estimates | Tissue | Arrays | CpGs | Ref |
|---|---|---|---|---|---|
| Bohlin | Gestational age (weeks) | cord blood | 450K | 96 | [16] |
| EPIC | Gestational age (weeks) | cord blood | EPIC | 176 | [17] |
| Knight | Gestational age (weeks) | cord blood | 450K | 148 | [18] |
| Lee.CPC | Gestational age (weeks) | placenta | 450K/EPIC | 1125 | [19] |
| Lee.refRPC | Gestational age (weeks) | placenta | 450K/EPIC | 1125 | [19] |
| Lee.RPC | Gestational age (weeks) | placenta | 450K/EPIC | 1125 | [19] |
| Mayne | Gestational age (weeks) | placenta | 450K | 62 | [20] |
3.4 Mitotic (cell-division) clocks
These do not estimate age but the cumulative number of stem-cell divisions a tissue has undergone — a proxy for proliferative history used mainly in cancer research.
| Clock | Estimates | Tissue | Arrays | CpGs | Ref |
|---|---|---|---|---|---|
| epiTOC1 | Cell divisions (count) | multi-tissue (proliferative) | 450K/EPIC | 385 | [21] |
| epiTOC2 | Cell divisions (count) | multi-tissue (proliferative) | 450K/EPIC | 163 | [22] |
| HypoClock | Cell divisions (count) | multi-tissue (proliferative) | 450K/EPIC | 678 | [22] |
| RepliTali | Cell divisions (count) | cultured cells | 450K/EPIC | 87 | [23] |
| stemTOC | Cell divisions (count) | multi-tissue (proliferative) | 450K/EPIC | 371 | [24] |
3.5 Pace of ageing, and other measures
DunedinPACE estimates the rate of biological ageing (biological years per
calendar year) rather than an age. A set of causal clocks separate damaging from
adaptive ageing signals.
3.6 Trait and exposure predictors
Beyond age, the package also carries methylation predictors of lifestyle and
clinical traits — body mass index, smoking, alcohol, cholesterol and others.
These are not ages, and are opt-in: they are not computed by default, so
request them by name (e.g. methylclock(betas, clocks = "McCartney.BMI")).
| Clock | Estimates | Tissue | Arrays | CpGs | Ref |
|---|---|---|---|---|---|
| McCartney.Alcohol | Trait score | whole blood | 450K/EPIC | 450 | [27] |
| McCartney.BMI | Trait score | whole blood | 450K/EPIC | 1109 | [27] |
| McCartney.BodyFat | Trait score | whole blood | 450K/EPIC | 968 | [27] |
| McCartney.Education | Trait score | whole blood | 450K/EPIC | 373 | [27] |
| McCartney.HDL | Trait score | whole blood | 450K/EPIC | 737 | [27] |
| McCartney.LDL | Trait score | whole blood | 450K/EPIC | 233 | [27] |
| McCartney.Smoking | Trait score | whole blood | 450K/EPIC | 233 | [27] |
| McCartney.TotalChol | Trait score | whole blood | 450K/EPIC | 204 | [27] |
| mCigarette | Trait score | whole blood | EPIC | 1255 | [28] |
The Ref column of every table above links to the original paper for that clock
in the reference list at the end. The same information is available
programmatically (clock_catalog$reference, or clock_info(name)).
To inspect any single clock in full — its engine, platform and generation —
use clock_info():
clock_info("Horvath")#> <clock> Horvath [chronological / gen 1 / linear]
#> resource: coefHorvath
#> platform: 27K/450K
#> status : ok4 Getting started
Install from Bioconductor:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("methylclock")Between Bioconductor releases, development accumulates on GitHub; to try the latest development version instead:
if (!requireNamespace("remotes", quietly = TRUE))
install.packages("remotes")
remotes::install_github("isglobal-brge/methylclock")Installation pulls in the packages methylclock builds on — you do not install
them by hand. Most are lightweight (preprocessCore, ggplot2, Rcpp and a few
others); the one worth knowing about is BigDataStatMeth (version 2.0.4 or
later), which does the heavy numerical work and lets the clocks run on
large, disk-backed (HDF5) data without loading the whole matrix into memory
— if your library carries an older copy, update it first.
4.1 What data to provide
The estimators take a methylation beta matrix: one row per CpG (row names are
the CpG identifiers, like cg00000029), one column per sample, and values between
0 and 1. A data frame whose first column holds the CpG identifiers is also
accepted.
betas[1:3, 1:3]
#> Sample_1 Sample_2 Sample_3
#> cg00000029 0.52 0.48 0.55
#> cg00000108 0.91 0.89 0.92
#> cg00000165 0.12 0.10 0.14You do not have to convert your data first. If you pass M-values (the
log-ratio scale, outside the 0–1 range), they are detected and converted to
betas automatically. And besides a matrix or data frame, methylclock() accepts
the objects you typically hold after processing IDATs — a minfi object
(GenomicRatioSet, MethylSet), a SummarizedExperiment or an ExpressionSet
— pulling the methylation matrix out of them for you:
methylclock(mvals) # M-values: converted to betas automatically
methylclock(grset) # a minfi GenomicRatioSet
methylclock(se) # a SummarizedExperimentYou do not need to worry about a clock’s CpGs being missing from your array
either: a clock is computed only when enough of its CpGs are present, and
otherwise returns NA for every sample rather than a misleading number. Missing
values among the CpGs that are present are filled in automatically.
For very large cohorts, the same call also reads betas straight from an on-disk (HDF5) store and keeps the results there, so memory use stays bounded instead of growing with the array.
4.2 Missing values and imputation
Real arrays have gaps — probes that failed QC on some samples. A clock cannot
skip a CpG for one sample and not another, so the missing values among the CpGs a
clock uses are filled in before it runs. How they are filled is up to you, through
impute:
"mean"(the default) fills each gap with that CpG’s average across your samples — simple and predictable."knn"fills it from the CpGs that behave most like it (its nearest neighbours across the samples), so co-methylation is used rather than a flat average. It runs entirely in compiled code and, on a disk-backed array, block by block, so it stays memory-safe even for whole-array clocks."reference"uses a stored reference level for each CpG, and"none"does not impute at all — any sample missing a used CpG is returned asNA.
Does the choice matter? methylclock_betas is a small, complete real beta matrix,
so we can blank out some values on purpose and see how close each method gets back
to the truth. First, the estimates on the complete data:
data(methylclock_betas)
demo_clocks <- c("Horvath", "Levine", "Wu")
truth <- as.data.frame(methylclock(methylclock_betas, clocks = demo_clocks))Now blank out 8% of the values at random and recompute, once with each method:
set.seed(42)
blanked <- methylclock_betas
blanked[sample(length(blanked), round(0.08 * length(blanked)))] <- NA
est <- function(method) {
as.data.frame(methylclock(blanked, clocks = demo_clocks, impute = method))
}
by_mean <- est("mean")
by_knn <- est("knn")
mae <- function(e) {
vapply(demo_clocks, function(cl) mean(abs(e[[cl]] - truth[[cl]])), numeric(1))
}
data.frame(clock = demo_clocks,
`MAE, mean impute` = round(mae(by_mean), 2),
`MAE, KNN impute` = round(mae(by_knn), 2),
check.names = FALSE, row.names = NULL)#> clock MAE, mean impute MAE, KNN impute
#> 1 Horvath 1.35 0.82
#> 2 Levine 1.44 1.18
#> 3 Wu 0.14 0.13The error is in years of clock age. For Horvath and Levine, filling gaps from
co-methylating CpGs ("knn") recovers the true estimate noticeably better than a
flat CpG average; for Wu the two are almost identical. The lesson is practical:
with little missing data the default is fine, but if a meaningful fraction of a
clock’s CpGs are missing, "knn" recovered them better in the example
below — at some extra compute.
5 Running the clocks and reading the output
methylclock() is the single entry point. Called on its own it computes every
clock your data support:
res <- methylclock(betas)Restrict the set by name, or let the package filter by platform or generation — useful when you know your array:
methylclock(betas, clocks = c("Horvath", "Levine"))
methylclock(betas, platform = "450K") # every clock native to 450KDNAmAge() and DNAmGA() are shortcuts that preselect the non-gestational and
gestational clocks, respectively.
The result is a methylclock object. Turn it into a plain data frame — one row
per sample, one column per clock — whenever you want to work with the numbers:
res <- methylclock(betas)
head(as.data.frame(res))For cohorts too large to hold in memory, the same call can read the betas from an on-disk (HDF5) store and keep the results there too, so the workflow does not change as your data grow.
6 Seeing and interpreting the results
The rest of this guide uses a bundled example so every plot runs as you read.
methylclock_demo holds the clock estimates for a public whole-blood dataset with
known ages (656 adults, 19–101 years); only the estimates, age and sex are
stored, not the methylation matrix. One thing to keep in mind while reading:
this cohort was part of the training of several of the clocks shown, so it
is used here to demonstrate the visualizations — the out-of-sample
evaluation comes later, on an independent cohort (Section
9).
data(methylclock_demo)
age_clocks <- c("Horvath", "Hannum", "Levine", "skinHorvath",
"BLUP", "EN", "AltumAge", "Lin")Every plotting function returns a ggplot object with a clean default look, so
you can drop it straight into a report or keep tailoring it — add a title, change
the theme, tweak the scales. The default theme is exported as theme_methylclock()
if you want to match it on your own figures.
6.1 Predicted versus real age
The first question is whether a clock recovers age. plotDNAmAge() plots each
clock’s estimate against the real age, one panel per clock. The dashed line is
predicted = real: the closer the points hug it, the more accurate the clock, and
each panel is labelled with its R^2 (how much of the age variation the clock
explains).
plotDNAmAge(methylclock_demo, age = methylclock_demo$age, clocks = age_clocks)
Figure 1: Estimated vs. real age
Points on the dashed line are perfect predictions.
How to read it. A tight, on-the-line cloud means the clock both follows age and reports it on the right scale. A cloud that is tilted or shifted off the line still tracks age but is miscalibrated for this data — worth knowing before you trust the absolute numbers (Section 9 makes this quantitative).
6.2 Do the clocks agree with each other?
plotClockCorrelation() shows the correlation between every pair of clocks.
Age clocks form a block of high positive correlations; a marker like telomere
length (TL) sits at the opposite sign, because it shortens as the others rise.
plotClockCorrelation(methylclock_demo, clocks = c(age_clocks, "TL"))
Figure 2: Correlation between clocks
Warm = agree; TL anchors the opposite sign.
How to read it. Strong agreement across age clocks is expected and reassuring. A clock that correlates weakly with the rest is a flag — either it measures something different (by design) or it is struggling on your data.
6.3 How the estimates are distributed
plotClockDistributions() compares the spread of each clock, which is handy for
spotting a clock on a different scale or with outliers.
plotClockDistributions(methylclock_demo, clocks = age_clocks)
Figure 3: Distribution of each clock’s estimates
7 Age acceleration
Two people the same age can have blood that looks epigenetically older or younger. Age acceleration captures that difference: it is the part of a clock’s estimate not explained by chronological age (the residual after regressing one on the other). A positive value means a sample looks older than its years. This — not the raw estimate — is what most studies correlate with exposures, disease or lifestyle.
plotAgeAcceleration(methylclock_demo, age = methylclock_demo$age,
clocks = age_clocks)
Figure 4: Age acceleration per clock (residual in years, centred at zero)
Split it by a sample variable to look for associations. Here, by sex, men sit slightly above women across most clocks — the commonly reported pattern of faster epigenetic ageing in men.
plotAccelerationByGroup(methylclock_demo, age = methylclock_demo$age,
group = methylclock_demo$sex, clocks = age_clocks)
Figure 5: Acceleration by sex
A group consistently above zero looks epigenetically older.
How to read it. Look for a group whose boxes sit consistently above (or below) zero across clocks: that is a candidate association worth testing formally. A shift seen across several clocks is a useful robustness pattern — though not independent evidence, since the clocks are highly correlated with each other.
7.1 Cell composition and immune-aware acceleration
Blood is a mix of cell types whose proportions themselves change with age, so some
of what a clock reads is really a shift in cell composition rather than ageing
within cells. To separate the two, first estimate the cell proportions from the
methylation matrix with cellCounts() (reference-based deconvolution):
cc <- cellCounts(betas) # samples x blood cell typesSeveral reference panels are available (listCellReferences()). For
higher-resolution immune profiling you can also use the FlowSorted panels —
"FlowSorted.Blood.EPIC" (6 cells) and "FlowSorted.BloodExtended.EPIC" (12
immune subtypes) — which are not shipped with the package: install the
corresponding package yourself, under its own licence, and cellCounts() uses it
if present.
cc <- cellCounts(betas, reference = "FlowSorted.BloodExtended.EPIC")ageAcceleration() then returns, per clock, the raw difference, the classic
age-acceleration residual (residual, predicted regressed on
chronological age) and — given the cell proportions — the cell-adjusted
residual (residualCells, further adjusted for cell composition):
methylclock_cells bundles clock estimates together with the estimated cell
proportions for 795 adult blood samples, so this runs live:
data(methylclock_cells)
cells <- c("Bcell", "CD4T", "CD8T", "Eos", "Mono", "Neu", "NK")
aa <- ageAcceleration(methylclock_cells, age = methylclock_cells$age,
cell_counts = as.matrix(methylclock_cells[, cells]))
head(aa)#> id clock ageAcc residual residualCells
#> 1 200928190033_R01C01 Horvath 14.813071 7.974972 6.2796193
#> 2 200928190033_R02C01 Horvath 17.840316 7.816071 7.2924363
#> 3 200928190033_R03C01 Horvath 11.086294 3.679240 0.8910541
#> 4 200928190033_R04C01 Horvath 9.624874 3.014356 2.7308029
#> 5 200928190033_R05C01 Horvath 2.296949 -4.199778 -3.8488830
#> 6 200928190033_R06C01 Horvath 6.200526 -3.482347 -2.4712663The two residuals are close but not identical — adjusting for cell composition shifts each sample’s value by a bit:
h <- aa[aa$clock == "Horvath", ]
c(correlation = cor(h$residual, h$residualCells),
mean_abs_diff = mean(abs(h$residual - h$residualCells)))#> correlation mean_abs_diff
#> 0.9767419 0.8132283Use the plain residual when the immune ageing captured by cell-composition
shifts is part of what you want to measure, and residualCells when you
want ageing independent of it — for example to compare tissues or to
control for cell-type confounding. Without cell_counts, only ageAcc and
residual are returned.
Two terminology notes for readers of the literature. The canonical IEAA
is this same adjustment idea on the Horvath clock, but defined with a
specific set of seven immune covariates [29]: naive and exhausted
CD8+ T cells, plasmablasts, CD4+ T cells, natural killer cells, monocytes
and granulocytes. IEAA() computes it as defined — the three rare
subsets from the methylation data itself, using the published CpG
estimators (see EEAA() below), and the four common ones from
cellCounts() (or a matrix you pass). Because it needs the beta matrix,
the example is shown but not run; for any other clock, or any other cell
panel, residualCells is the honest label:
res <- DNAmAge(betas, clocks = "Horvath")
ieaa <- IEAA(res, age = pheno$age, betas = betas)And the published EEAA is not the plain residual but a specific
cell-weighted blend of the Hannum clock with three immune cell estimators
(plasmablasts, exhausted and naive CD8+ T cells), its weights fixed once in
the Women’s Health Initiative [29]. EEAA() computes exactly that
measure: the three cell estimators are derived from the methylation data
itself as published linear combinations of CpGs — the same route the
original software takes — and the fixed numeric parameters are the ones
published in European patent EP 3 494 210 B1 [30]. It therefore needs the beta
matrix (blood, 450K/EPIC), not just the clock estimates, so the example is
shown but not run — the bundled demo data do not carry the estimator CpGs:
res <- DNAmAge(betas, clocks = "Hannum")
eeaa <- EEAA(res, age = pheno$age, betas = betas)ageAccelerationChen() offers a related but deliberately different measure,
in the spirit of Chen et al. [29]: a Klemera–Doubal weighted blend
of the Hannum clock with immune cell measures you supply, with the blend
weights re-estimated from your data rather than fixed. (The patent
documents this re-estimated flavour too, as its “dynamic” variant, and
reports it highly correlated with the fixed-weight measure.) Because both
the weights and the cell inputs differ from the canonical recipe, its
output column is named chenAcc, not EEAA. It expects the specific immune
subsets the measure was defined on (naive/exhausted CD8, plasmablasts),
which the reference panels bundled here do not resolve — they report
total CD8+ T cells. For that reason this example too is shown but not
run: substituting the closest available columns would compute a quantity
that looks like Chen’s measure but is biologically a different one, and
this guide prefers no number to a misleading one.
chen <- ageAccelerationChen(res$Hannum, age = pheno$age,
cells = immune_proportions) # the Chen subsetsIf you need those subsets, the
proportions can be obtained elsewhere and passed straight to cells: from the
advanced analysis of the online epigenetic-age calculator [1],
which reports these subsets directly. An extended immune deconvolution such
as the EPIC IDOL-Ext library [31] is sometimes used instead, but
note what it does and does not resolve: naive and memory CD8+ T cells and B
cells, yes; exhausted CD8+ T cells and plasmablasts, no — so proportions
from it are stand-ins, and the result is a measure in the spirit of Chen’s
rather than one computed on the original subsets. That library comes under
its own licence from the authors’ institution; install it yourself and cite
the source you used. Bear in mind that the rarer subsets are reported by
some algorithms as scores linearly related to abundance rather than as
cell-type proportions.
8 Comparing groups
Most studies do not ask “how old is this sample” but “do these two groups differ” — cases and controls, exposed and unexposed, one population and another. The package gives that question three complementary views: the shape of each group’s values, the trend of each group along age, and a one-look summary of the differences with their uncertainty.
The bundled methylclock_demo cohort records each sample’s reported sex,
which makes a convenient two-group example:
data(methylclock_demo)
demo <- methylclock_demo
table(demo$sex)#>
#> F M
#> 338 3188.1 Shapes: densities by group
plotClockDensities() draws each clock’s density curve, and with group =
overlays one curve per group, so a shift, a widening or an extra mode is
visible directly:
plotClockDensities(demo, group = demo$sex,
clocks = c("Horvath", "Hannum", "BLUP", "EN"))
There is a reading to be careful with here. If the groups differ in
chronological age, their raw estimates differ with it, and the curves
separate for that reason alone. The more comparable scale for a group
contrast is the age acceleration, which removes that trend — and only
that trend:
what = "acceleration" (with age =) draws the same figure on it. Two more
layouts are available when they say it better: style = "mirror" (two
groups drawn as facing halves) and style = "ridge" (every clock stacked in
one column).
8.2 Trends: fitted lines by group
plotDNAmAge() accepts the same group =: points and fitted lines are
drawn per group over the identity line. A group whose line runs above the
other’s is predicted older at the same chronological age; lines with
different slopes mean that gap changes along the age range:
plotDNAmAge(demo, age = demo$age, group = demo$sex,
clocks = c("Horvath", "Hannum"), ncol = 2)
8.3 Summary: the difference, clock by clock
plotGroupDifference() condenses the comparison to one row per clock: the
difference between the two groups’ mean acceleration, with its confidence
interval. The dashed line at zero means “no difference”: an interval that
reaches it is compatible with none, and points whose interval leaves zero
out are drawn filled.
plotGroupDifference(demo, group = demo$sex, age = demo$age)
When both levels matter — one group being high tells a different story
than the other being low, even with the same gap — style = "dumbbell"
draws the two group means joined by a segment instead. And when no ages are
known, what = "estimate" compares the raw estimates; any difference in
chronological age between the groups is then part of the number, so read it
as a description of the groups as they are, not as an age-adjusted contrast.
A boundary worth stating plainly: all of these are descriptive, exploratory summaries. Removing the age trend adjusts for nothing else — not cell composition, batch, cohort or any other covariate — and the intervals describe each contrast in isolation. When a group difference matters to your conclusions, fit a model that carries the covariates your design needs; these figures are for seeing, not for deciding.
9 Quality control: does a clock fit your data?
Two checks belong here regardless of whether you know the samples’ real ages.
How many of each clock’s CpGs your array carries. A clock computed from only a
fraction of its CpGs is on thin ice (and below a threshold it is not computed at
all). clockCoverage() reports this from the methylation matrix:
clockCoverage(betas)
#> clock n_cpgs present pct_present
#> Horvath 353 349 98.9
#> Hannum 71 65 91.5
#> ...A note on ages. The clocks estimate an age (or score) without needing the
true age — so you can always compute them, compare them (Section
6) and check their coverage. You only
need the real chronological age to measure accuracy (R2, MAE) or age
acceleration; without it, those are simply not reported.
When you do have ages, qcReport() gathers everything — descriptives, coverage,
and the accuracy of each age clock, optionally broken down by a variable such as
sex, case/control status or ethnicity — in one call:
data(methylclock_validation)
grady <- subset(methylclock_validation, cohort == "adult (GSE132203)")
qcReport(grady, age = grady$age, by = grady$sex,
clocks = c("Horvath", "Hannum", "skinHorvath", "BLUP", "EN"))#> methylclock QC report
#> - samples: 795
#> - age (years): Min. 18 1st Qu. 32 Median 45 Mean 42.3 3rd Qu. 52 Max. 76
#> - groups: F=571 M=224
#> - clocks computed on all samples: 5 of 5
#> - accuracy vs chronological age:
#> group clock n r R2 MAE
#> M Horvath 224 0.880 0.774 7.68
#> M Hannum 224 0.895 0.802 5.61
#> M skinHorvath 224 0.932 0.869 2.54
#> M BLUP 224 0.955 0.911 2.33
#> M EN 224 0.952 0.906 2.96
#> F Horvath 571 0.898 0.807 7.40
#> F Hannum 571 0.915 0.837 5.75
#> F skinHorvath 571 0.950 0.903 2.43
#> F BLUP 571 0.963 0.927 2.06
#> F EN 571 0.961 0.923 2.85
#>
#> Reading of each section
#>
#> Overview
#> This report covers 795 samples with chronological ages from 18 to 76
#> years (median 45). All 5 clocks returned a value for every sample.
#>
#> Sample variables
#> Sample variables supplied: age and a grouping variable (F, M), so the
#> report breaks the accuracy down by group.
#>
#> Agreement
#> Clocks correlate with one another because they read the same ageing
#> signal from different sets of CpGs. Here the median correlation
#> between pairs of clocks is 0.94. The pair that tracks each other most
#> closely is BLUP and EN (0.98); the pair furthest apart is Hannum and
#> Horvath (0.91). Clocks built for different targets --- calendar age,
#> cell divisions, a trait --- measure different quantities, so a low
#> correlation between two of them follows from their design and not
#> from the data.
#>
#> Distributions
#> This shows where each clock places your samples and how widely it
#> spreads them. Median estimates run from 38.2 (Hannum) to 51.7
#> (Horvath). Clocks reported in different units do not share a scale
#> --- years of age, weeks of gestation and cumulative cell divisions
#> are not comparable numbers --- so a clock sitting apart from the rest
#> is often on its own scale rather than disagreeing about the samples.
#>
#> Accuracy against age
#> r is how tightly a clock's estimate follows chronological age, R2 the
#> share of the age variation it accounts for, and MAE the typical gap
#> between estimate and age, in years. Across these clocks r runs from
#> 0.88 (Horvath [M]) to 0.96 (BLUP [F]), and the median absolute error
#> from 2.1 years (BLUP [F]) to 7.7 (Horvath [M]). The two columns say
#> different things, and a clock can score high on one and low on the
#> other: a high r with a large MAE is the signature of a clock that
#> orders the samples by age correctly while sitting on a different
#> scale for this dataset, which is what happens when a predictor is
#> applied to a population or tissue other than the one it was fitted
#> on.
#>
#> Bias
#> The bias is the average signed gap between estimate and chronological
#> age: above zero the clock reads older than the calendar, below zero
#> younger. The dashed lines hold 95% of the samples. It runs from -5.7
#> years (Hannum) to +7.6 years (Horvath). No clock's gap drifts
#> markedly across the age range: the offsets are roughly constant from
#> the youngest samples to the oldest.
#>
#> Per-sample quality
#> A sample is flagged when its estimates sit far from the rest of the
#> dataset --- a mean absolute z-score above 3 across the clocks. 0 of
#> 795 samples are flagged. The flag is relative to this dataset: it
#> marks a sample as unusual among these samples, which can come from
#> the assay or from the sample being genuinely atypical.This section deliberately switches cohorts. methylclock_demo (GSE40279) is
the training cohort of Hannum, BLUP and EN, so accuracy computed on it
describes training fit, not performance on new data — it remains a fine
cohort for illustrating the functions, which is how the earlier sections
use it. The cohort here (GSE132203) is independent of every clock shown, so
its numbers read as out-of-sample accuracy.
The accuracy has two columns that say different things. A high r/R2 means the
clock follows age; a low MAE means it reports it on the right scale. A clock
can have high correlation but a large error — it tracks age with a constant
offset. Reading them per group is how you catch a clock that works on average but
drifts in one subgroup. For just the table, clockAccuracy() returns that piece
alone.
One more angle worth a look is the samples themselves: sampleQC() flags any that
stand out in clock space, and — given the methylation matrix — how much of each
sample’s data had to be imputed, so you can spot the samples whose estimates are
closer to an approximation than a measurement.
sampleQC(res, betas = betas) # per-sample outlier score and % missing CpGsIts discordance column measures something the outlier score does not: how
much the clocks contradict each other on a sample. The two are different
questions — a sample can be extreme with every clock in agreement, or
unremarkable on average while the clocks disagree. plotSampleDiscordance()
draws both at once: each sample placed by the clocks’ consensus
(horizontally) and their disagreement (vertically), with the most discordant
samples labelled by id. Most samples form a low band; one far above it is a
sample the clocks genuinely disagree about — a technical problem, or a
biology worth a second look:
plotSampleDiscordance(methylclock_demo)
To share the check, output = "html" writes the whole thing — the summary, the
coverage, accuracy and per-sample tables, and the key plots (agreement,
distributions, predicted-vs-real and Bland-Altman) — to a self-contained HTML
file. Pass betas = to include the coverage and missingness, and impute = to
record which imputation method produced the estimates; the report then also lists,
per clock, how many CpG values had to be filled in:
qcReport(betas = betas, x = res, age = pheno$age, by = pheno$sex,
impute = "knn", output = "html", file = "qc.html")Each section of the report is written out as well as plotted: what the section
shows, what your own numbers read as, and what produces such a pattern, with a
glossary of the terms at the end. The same text is on the object, under
notes, so it is also there when you work at the console rather than in HTML.
10 Reference ranges: your samples against a named cohort
The growth-chart idea, applied to epigenetic age: plotReferenceRange()
draws, for each clock, the percentile bands of a reference cohort along
chronological age (5th–95th light, 25th–75th darker, median line), with
your samples on top. A sample inside the central band sits where most of
that reference does at its age; one outside the light band sits where fewer
than one in ten of them do. Where the reference has too few samples near an
age, no band is drawn: there is no rule to measure against there.
The bands describe the named cohort, nothing more — a reference measures
well what resembles it, in population, tissue, array and processing. The
most meaningful reference is therefore your own control samples, passed
as a data frame with an age column. When none is given, the function falls
back to a cohort shipped with the package, says so in a message, and the
subtitle always states which reference was used and its size; the shipped
options and what each one is (a population-based Swedish cohort, two African
American cohorts with their recruitment caveats, the two ethnicity strata of
GSE40279, a paediatric set, and more) are listed in ?plotReferenceRange.
None of them stands for people in general.
pred <- methylclock_demo[1:12, ]
plotReferenceRange(pred, age = pred$age, reference = "GSE87571",
clocks = c("Horvath", "Hannum", "skinHorvath", "BLUP"))
How much the choice of reference matters is easy to see with the two
population strata of GSE40279: the same samples drawn over
"GSE40279-Hispanic" and over "GSE40279-Caucasian" land differently,
because each rule was made from a different population. That is the reason
the reference is always declared — and why, if your study has controls,
they are the rule to prefer.
11 Choosing the right clock for your data
The single most useful habit is to pick a clock by what it was built for — its
tissue and age range — not by its reputation. The catalogue tables of Section
3 carry a Tissue
column for exactly this, and the registry exposes it as clock_info(x)$tissue:
most clocks here read whole blood, but the package also ships clocks trained
on buccal epithelium, cord blood, placenta, skin-and-blood and multi-tissue
panels, and their answers are only calibrated on their own ground. A second bundled dataset,
methylclock_validation, makes the point: it holds estimates for two public blood
cohorts, 795 adults (Grady Trauma Project, GEO accession GSE132203) and
134 children (1–17 years) [32].
data(methylclock_validation)
adult <- subset(methylclock_validation, cohort == "adult (GSE132203)")
kids <- subset(methylclock_validation, cohort == "pediatric (GSE36054)")clockAccuracy() on the children shows the lesson clearly:
clockAccuracy(kids, age = kids$age,
clocks = c("Horvath", "skinHorvath", "Wu", "BLUP", "Levine"))#> clock n r R2 MAE
#> 1 Horvath 134 0.934 0.873 0.59
#> 2 skinHorvath 134 0.967 0.936 0.42
#> 3 Wu 134 0.977 0.955 0.66
#> 4 BLUP 134 0.962 0.925 7.16
#> 5 Levine 134 0.845 0.714 24.15The contrast carries the lesson. BLUP, an adult-trained clock, still
correlates with age (high r) but is badly miscalibrated (large MAE):
it reads children’s ages with a big offset, because nothing in its training
looked like a toddler. Levine, which measures biological rather than
chronological age, should not be read as a calendar age here at all. The
clocks built to include the young — Wu, the skin & blood clock,
Horvath — land far closer. (One candour note: this very cohort was among
Wu’s training datasets, so Wu’s sub-year error here reads as training fit
rather than independent validation; the miscalibration lesson rests on the
adult-trained clocks, which never saw these children.)
The practical rule: correlation tells you a clock follows age; only the right
clock for your tissue and age range reports it correctly. The registry’s
platform and target fields (Section 3)
are there so you can choose deliberately.
12 Longitudinal data: trajectories and slopes
When the same subjects were sampled at several ages, a new question opens:
not how old each sample reads, but how each person’s estimate moves over
time. The package ships a real example of that shape —
methylclock_longitudinal, 68 children sampled repeatedly from birth to
about age three — and three small tools to work with it.
clockTrajectories() pairs the samples into per-subject trajectories (one
row per subject, visit and clock, ordered in time). Technical replicates at
the same age are averaged, and subjects seen only once are excluded, each
with a message — one point has no trajectory. The table it returns is also
exactly what a mixed-model package expects, for example
lme4::lmer(estimate ~ age + (age | subject), data = tr): this package
shapes and displays longitudinal data; the formal inference belongs to those
tools.
plotTrajectories() is the longitudinal view itself: one thin line per
subject, the group’s binned median in bold, and the identity line dashed —
a subject moving parallel to it gains one epigenetic year per calendar year:
data(methylclock_longitudinal)
d <- methylclock_longitudinal
plotTrajectories(d, subject = d$subject, age = d$age,
clocks = c("Wu", "PedBE", "Horvath"))
The figure carries its own lessons: the paediatric clocks (Wu, PedBE)
climb steeply at these ages — a mixture of genuinely rapid developmental
change in the methylome and of each clock’s own calibration — while
Horvath, calibrated for adults, barely moves here. The wobble of the thin
lines is faithful too: repeated clock estimates are noisy, and the plot does
not smooth it away.
trajectoryRates() turns each subject’s visits into a slope: how many
epigenetic years the estimate advances per calendar year, 1 meaning it
tracks the calendar. This is a property of the estimate’s trajectory, not by
itself a measure of biological ageing pace — of the clocks here, only
DunedinPACE was designed and validated to estimate such a pace. With two
visits the slope is a difference quotient and carries no per-subject
uncertainty (se is NA: read it as orientative); with three or more it is
a fitted slope with its standard error. Read individual slopes with care: the
technical replication error of many clocks spans several years
[33], and the difference of two measurements carries the noise
of both, so over a short follow-up a single slope can be mostly measurement.
The steadier signal is in the aggregate — the distribution of slopes in a
group, or their group difference, for which
plotGroupDifference() accepts these rates directly (one row per subject).
head(trajectoryRates(d, subject = d$subject, age = d$age, clocks = "Wu"))#> subject clock n_visits span_years slope se
#> 1 SV1 Wu 4 2.035616 0.06136641 2.6275210
#> 2 SV11 Wu 5 3.035616 2.87504626 0.2551833
#> 3 SV12 Wu 5 3.000000 2.72067011 0.1446377
#> 4 SV15 Wu 4 3.005479 2.41810487 0.8114296
#> 5 SV16 Wu 5 3.205479 2.12923399 0.2214229
#> 6 SV18 Wu 4 3.000000 2.65676918 0.236085913 Beyond age: a trait predictor
Not every clock estimates an age. The trait and exposure predictors (Section 3.6) are methylation scores for lifestyle and clinical traits — smoking, body mass index, alcohol and others. They are opt-in, so you request them by name. A quick example with a smoking predictor, on a public cohort with recorded smoking status:
methylclock(betas, clocks = "McCartney.Smoking")methylclock_smoking bundles that score for 464 blood samples labelled never,
former or current smokers (GEO accession GSE50660 [34]),
so we can see whether it tracks the exposure:
data(methylclock_smoking)
tapply(methylclock_smoking$McCartney.Smoking, methylclock_smoking$smoking,
median)#> never former current
#> 3.129504 4.146654 5.532425library(ggplot2)
ggplot(methylclock_smoking,
aes(smoking, McCartney.Smoking, fill = smoking)) +
geom_boxplot(alpha = 0.85, outlier.size = 0.6) +
scale_fill_manual(values = c("#CC6677", "#332288", "#44AA99")) +
labs(x = NULL, y = "McCartney smoking score") +
theme_methylclock() + theme(legend.position = "none")
Figure 6: Methylation smoking score by smoking status
The score climbs from never- through former- to current-smokers, exactly the
ordering you would hope for. Two caveats to read it correctly: a trait score is a
relative value, not a calibrated quantity like age — it is meaningful within
a dataset, compared across groups or against a reference, not as an absolute
number. And this is an illustration of how to compute and use such a clock, not
a validation of the predictor. As always, the plot is a plain ggplot, so you can
restyle it for your own report.
14 EPICv2 arrays
The MethylationEPIC v2 array (2023) renamed its probes (suffixes such as
_BC21), removed some that older arrays carried, measures some CpGs with
several replicate probes, and redesigned others. The clocks in this package
were built on 450K/EPICv1 names, so on EPICv2 data they fail visibly: the
probe names do not match, each clock’s coverage collapses, and the result is
NA with a warning rather than a number.
One clock is the exception by construction: Garma [6] was trained
only on probes present on the 450K, EPICv1 and EPICv2 arrays, so it runs
natively on all three with no conversion.
A caution for the impatient: stripping the version suffixes from EPICv2 probe names by hand makes the names match again — and quietly ignores the removed, replicated and redesigned probes, which can produce plausible but wrong estimates. Support for running the older clocks on EPICv2 through a proper probe mapping is planned as future work.
15 Session info
sessionInfo()#> R version 4.6.1 Patched (2026-06-24 r90190)
#> Platform: x86_64-apple-darwin20
#> Running under: macOS Ventura 13.7.8
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
#>
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] ggplot2_4.0.3 methylclockData_1.21.3 futile.logger_1.4.9
#> [4] methylclock_1.99.1 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] DBI_1.3.0 bitops_1.1-0 RBGL_1.89.0
#> [4] httr2_1.3.0 formatR_1.14 rlang_1.3.0
#> [7] magrittr_2.0.5 otel_0.2.0 compiler_4.6.1
#> [10] RSQLite_3.53.3 mgcv_1.9-4 GenomicFeatures_1.65.0
#> [13] png_0.1-9 systemfonts_1.3.2 vctrs_0.7.3
#> [16] stringr_1.6.0 pkgconfig_2.0.3 crayon_1.5.3
#> [19] fastmap_1.2.0 magick_2.9.1 dbplyr_2.6.0
#> [22] XVector_0.53.0 labeling_0.4.3 biocViews_1.81.1
#> [25] Rsamtools_2.29.0 rmarkdown_2.32 graph_1.91.0
#> [28] UCSC.utils_1.9.0 tinytex_0.60 purrr_1.2.2
#> [31] bit_4.6.0 xfun_0.60 AnnotationHubData_1.43.0
#> [34] cachem_1.1.0 cigarillo_1.3.1 GenomeInfoDb_1.49.1
#> [37] jsonlite_2.0.0 blob_1.3.0 BiocParallel_1.47.0
#> [40] parallel_4.6.1 R6_2.6.1 bslib_0.12.0
#> [43] stringi_1.8.9 RColorBrewer_1.1-3 BigDataStatMeth_2.0.4
#> [46] rtracklayer_1.73.0 GenomicRanges_1.65.4 jquerylib_0.1.4
#> [49] Rcpp_1.1.2 Seqinfo_1.3.2 bookdown_0.48
#> [52] knitr_1.52 IRanges_2.47.5 BiocBaseUtils_1.15.1
#> [55] Matrix_1.7-6 splines_4.6.1 tidyselect_1.2.1
#> [58] stringdist_0.9.17 rstudioapi_0.19.0 dichromat_2.0-1
#> [61] yaml_2.3.12 codetools_0.2-20 RUnit_0.4.33.1
#> [64] curl_8.0.0 lattice_0.23-1 tibble_3.3.1
#> [67] S7_0.2.2 Biobase_2.73.2 withr_3.0.3
#> [70] KEGGREST_1.53.6 evaluate_1.0.5 AnnotationForge_1.55.0
#> [73] lambda.r_1.2.4 BiocFileCache_3.3.0 xml2_1.6.0
#> [76] ExperimentHub_3.3.2 Biostrings_2.81.9 pillar_1.11.1
#> [79] BiocManager_1.30.27 filelock_1.0.3 stats4_4.6.1
#> [82] OrganismDbi_1.55.1 generics_0.1.4 RCurl_1.98-1.20
#> [85] BiocVersion_3.24.0 S4Vectors_0.51.9 BiocCheck_1.49.31
#> [88] scales_1.4.0 glue_1.8.1 tools_4.6.1
#> [91] ExperimentHubData_1.39.1 AnnotationHub_4.3.2 BiocIO_1.23.3
#> [94] data.table_1.18.6.1 GenomicAlignments_1.49.2 XML_3.99-0.24
#> [97] grid_4.6.1 AnnotationDbi_1.75.2 nlme_3.1-171
#> [100] restfulr_0.0.17 cli_3.6.6 rappdirs_0.3.4
#> [103] kableExtra_1.4.1 textshaping_1.0.5 futile.options_1.0.1
#> [106] viridisLite_0.4.3 svglite_2.2.2 dplyr_1.2.1
#> [109] gtable_0.3.6 sass_0.4.10 digest_0.6.39
#> [112] BiocGenerics_0.59.12 rjson_0.2.23 farver_2.1.2
#> [115] memoise_2.0.1 htmltools_0.5.9 lifecycle_1.0.5
#> [118] httr_1.4.9 bit64_4.8.6