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)

Quick start

For the impatient: install (Section @ref(getting-started)), 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.

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.

clock_info("Horvath")

#> <clock> Horvath  [chronological / gen 1 / linear]
#>   resource: coefHorvath
#>   platform: 27K/450K
#>   status  : ok

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.

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

methylclock(mvals)          # M-values: converted to betas automatically
methylclock(grset)          # a minfi GenomicRatioSet
methylclock(se)             # a SummarizedExperiment

data(methylclock_betas)
demo_clocks <- c("Horvath", "Levine", "Wu")
truth <- as.data.frame(methylclock(methylclock_betas, clocks = demo_clocks))

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

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 450K

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

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 @ref(quality-control-does-a-clock-fit-your-data)).

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.

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)

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 @ref(quality-control-does-a-clock-fit-your-data) makes this quantitative).

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

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.

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)

Distribution of each clock’s estimates.

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)

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)

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.

cc <- cellCounts(betas)          # samples x blood cell types

cc <- cellCounts(betas, reference = "FlowSorted.BloodExtended.EPIC")

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

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

res <- DNAmAge(betas, clocks = "Horvath")
ieaa <- IEAA(res, age = pheno$age, betas = betas)

res <- DNAmAge(betas, clocks = "Hannum")
eeaa <- EEAA(res, age = pheno$age, betas = betas)

chen <- ageAccelerationChen(res$Hannum, age = pheno$age,
                            cells = immune_proportions)   # the Chen subsets

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 318

plotClockDensities(demo, group = demo$sex,
                   clocks = c("Horvath", "Hannum", "BLUP", "EN"))

plotDNAmAge(demo, age = demo$age, group = demo$sex,
            clocks = c("Horvath", "Hannum"), ncol = 2)

plotGroupDifference(demo, group = demo$sex, age = demo$age)

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 @ref(seeing-and-interpreting-the-results)) 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 CpGs

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

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.

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 @ref(which-clocks-are-available-and-what-they-are-for) 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.15

The 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 @ref(which-clocks-are-available-and-what-they-are-for)) are there so you can choose deliberately.

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

Beyond age: a trait predictor

Not every clock estimates an age. The trait and exposure predictors (Section @ref(trait-and-exposure-predictors)) 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.532425

library(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")

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.

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.

Session info

sessionInfo()

#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_4.0.3      methylclock_1.99.1 BiocStyle_2.41.0  
#> 
#> loaded via a namespace (and not attached):
#>  [1] KEGGREST_1.53.6       gtable_0.3.6          xfun_0.60            
#>  [4] bslib_0.12.0          httr2_1.3.0           lattice_0.23-1       
#>  [7] Biobase_2.73.2        bitops_1.1-0          vctrs_0.7.3          
#> [10] tools_4.6.1           generics_0.1.4        stats4_4.6.1         
#> [13] curl_8.0.0            tibble_3.3.1          AnnotationDbi_1.75.2 
#> [16] RSQLite_3.53.3        blob_1.3.0            BiocBaseUtils_1.15.1 
#> [19] pkgconfig_2.0.3       Matrix_1.7-6          data.table_1.18.6.1  
#> [22] dbplyr_2.6.0          RColorBrewer_1.1-3    S7_0.2.2             
#> [25] S4Vectors_0.51.9      lifecycle_1.0.5       compiler_4.6.1       
#> [28] farver_2.1.2          stringr_1.6.0         textshaping_1.0.5    
#> [31] Biostrings_2.81.9     Seqinfo_1.3.2         htmltools_0.5.9      
#> [34] sys_3.4.3             buildtools_1.0.0      sass_0.4.10          
#> [37] RCurl_1.98-1.20       yaml_2.3.12           crayon_1.5.3         
#> [40] pillar_1.11.1         jquerylib_0.1.4       cachem_1.1.0         
#> [43] nlme_3.1-171          ExperimentHub_3.3.2   AnnotationHub_4.3.2  
#> [46] tidyselect_1.2.1      digest_0.6.39         stringi_1.8.9        
#> [49] purrr_1.2.2           dplyr_1.2.1           BiocVersion_3.24.0   
#> [52] labeling_0.4.3        splines_4.6.1         maketools_1.3.2      
#> [55] grid_4.6.1            fastmap_1.2.0         cli_3.6.6            
#> [58] magrittr_2.0.5        withr_3.0.3           filelock_1.0.3       
#> [61] scales_1.4.0          rappdirs_0.3.4        bit64_4.8.6          
#> [64] XVector_0.53.0        rmarkdown_2.32        httr_1.4.9           
#> [67] bit_4.6.0             otel_0.2.0            png_0.1-9            
#> [70] kableExtra_1.4.1      memoise_2.0.1         evaluate_1.0.5       
#> [73] knitr_1.52            IRanges_2.47.5        BiocFileCache_3.3.0  
#> [76] viridisLite_0.4.3     mgcv_1.9-4            rlang_1.3.0          
#> [79] BigDataStatMeth_2.0.4 Rcpp_1.1.2            glue_1.8.1           
#> [82] DBI_1.3.0             BiocManager_1.30.27   xml2_1.6.0           
#> [85] BiocGenerics_0.59.12  svglite_2.2.2         rstudioapi_0.19.0    
#> [88] jsonlite_2.0.0        R6_2.6.1              systemfonts_1.3.2