Bio Applied Epigenetic Clocks

Epigenetic Clocks and Aging Analysis

When to Use

• Predicting chronological or biological age from Illumina 450K/EPIC methylation beta values
• Building or re-implementing an elastic-net DNAm clock (Horvath, Hannum, PhenoAge, GrimAge style)
• Computing epigenetic age acceleration (EAA = DNAmAge − chronological age) and adjusting for cell composition
• Testing EAA associations with exposures (smoking, BMI, disease status) or mortality/morbidity outcomes
• Comparing first- vs second-generation clocks for downstream epidemiological analysis

Version Compatibility

scikit-learn ≥1.3, pandas ≥2.0, numpy ≥1.24, scipy ≥1.11, Python ≥3.10. Concepts map directly onto the Bioconductor methylclock R package (≥1.6) if working from IDAT/GEO methylation matrices in R.

Prerequisites

pip install scikit-learn pandas numpy scipy matplotlib. Assumes beta values (0–1 methylation fraction) already QC'd, normalized (Noob/BMIQ), and matrix-aligned to a CpG reference panel — see bio-methylation-analysis-methylation-calling and bio-methylation-analysis-dmr-detection for upstream processing.

Goal: Predict age from a CpG methylation matrix using an elastic-net regression (the Horvath-clock strategy: age-transform the response, then L1/L2-penalized regression on thousands of correlated CpGs).

Approach: Apply Horvath's log/linear age transform to compress childhood dynamics, fit ElasticNetCV with cross-validated predictions to avoid overfitting-driven MAE, then invert the transform to report years.

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.linear_model import ElasticNetCV
from sklearn.model_selection import cross_val_predict

def transform_age(age, adult_age=20):
    """Horvath (2013) age transform: compresses rapid childhood methylation
    dynamics onto the same scale as slower adult change."""
    age = np.asarray(age, dtype=float)
    return np.where(age < adult_age, np.log(age + 1) - np.log(adult_age + 1),
                     (age - adult_age) / (adult_age + 1))

def inv_transform_age(t, adult_age=20):
    """Inverse of transform_age(); maps model output back to years."""
    t = np.asarray(t, dtype=float)
    return np.where(t < 0, np.exp(t + np.log(adult_age + 1)) - 1,
                     t * (adult_age + 1) + adult_age)

def fit_epigenetic_clock(beta_matrix, chronological_ages, l1_ratio=0.5, cv=5):
    """Fit an elastic-net DNAm clock and return cross-validated age predictions.

    beta_matrix: (n_samples, n_cpgs) array/DataFrame of methylation beta values.
    chronological_ages: (n_samples,) true ages in years.
    Returns (predicted_age_years, mae, pearson_r, fitted_model).
    """
    y = transform_age(chronological_ages)
    model = ElasticNetCV(l1_ratio=l1_ratio, cv=cv, n_alphas=50, max_iter=5000)
    # cross_val_predict avoids optimistic in-sample MAE from the same folds used to pick alpha
    y_pred_cv = cross_val_predict(model, beta_matrix, y, cv=cv)
    pred_age = inv_transform_age(y_pred_cv)
    mae = np.abs(pred_age - chronological_ages).mean()
    r, _ = stats.pearsonr(chronological_ages, pred_age)
    model.fit(beta_matrix, y)  # final model on all data, for scoring new samples
    n_nonzero = np.sum(model.coef_ != 0)
    print(f"DNAm clock: {n_nonzero} non-zero CpGs, MAE={mae:.2f} yr, r={r:.3f}")
    return pred_age, mae, r, model

Goal: Turn clock predictions into epigenetic age acceleration (EAA) and test whether it associates with an exposure such as smoking, correcting for the confound that raw EAA still carries a chronological-age trend.

Approach: Regress DNAmAge on chronological age (optionally plus cell-type proportions for intrinsic EAA), take the residual, then run standard linear association tests.

import numpy as np
import pandas as pd
from scipy.stats import linregress, pearsonr, ttest_ind

def compute_age_acceleration(dnam_age, chronological_age, cell_props=None):
    """Residualize DNAmAge on chronological age (extrinsic EAA), or additionally
    on cell-type proportions (intrinsic EAA) if cell_props is given.

    cell_props: optional (n_samples, n_celltypes) array, e.g. CD4T/CD8T/NK/B/Mono/Gran.
    Returns EAA residuals (years); positive = biologically older than expected.
    """
    covariates = {"chron_age": chronological_age}
    if cell_props is not None:
        for i in range(np.asarray(cell_props).shape[1]):
            covariates[f"celltype_{i}"] = np.asarray(cell_props)[:, i]
    X = pd.DataFrame(covariates)
    X.insert(0, "intercept", 1.0)
    beta, *_ = np.linalg.lstsq(X.values, dnam_age, rcond=None)
    fitted = X.values @ beta
    return dnam_age - fitted

def test_eaa_association(eaa, exposure, group=None):
    """Test EAA vs a continuous exposure (linear regression) or a binary
    group label (Welch t-test). Returns a dict of effect size + p-value."""
    if group is not None:
        eaa = np.asarray(eaa)
        mask = np.asarray(group).astype(bool)
        t, p = ttest_ind(eaa[mask], eaa[~mask], equal_var=False)
        return {"test": "t-test", "effect": eaa[mask].mean() - eaa[~mask].mean(), "p": p}
    slope, intercept, r, p, se = linregress(exposure, eaa)
    return {"test": "linregress", "slope": slope, "r": r, "p": p, "se": se}

Goal: Compare first-generation (Horvath, Hannum) vs second-generation (PhenoAge, GrimAge) clocks by accuracy against chronological age and by predictive power for a mortality outcome.

Approach: Score each clock's Pearson r / MAE against chronological age, then bin the mortality-tuned clock's acceleration into quintiles and compare observed death rates — mirrors how GrimAge's superiority over Horvath is demonstrated in the literature.

import numpy as np
from scipy import stats

def clock_mortality_quintiles(grimage_eaa, died, n_bins=5):
    """Bin EAA into quintiles and return the observed event rate per bin.
    Use to visualize/report a mortality-tuned clock's discriminative power."""
    edges = np.percentile(grimage_eaa, np.linspace(0, 100, n_bins + 1))
    rates = []
    for lo, hi in zip(edges[:-1], edges[1:]):
        mask = (grimage_eaa >= lo) & (grimage_eaa < hi)
        rates.append(died[mask].mean() if mask.sum() > 0 else np.nan)
    return np.array(rates), edges

def compare_clocks(clock_estimates, chronological_age):
    """clock_estimates: dict of {clock_name: predicted_age_array}.
    Returns a DataFrame of MAE and Pearson r vs chronological age per clock."""
    rows = []
    for name, est in clock_estimates.items():
        r, _ = stats.pearsonr(chronological_age, est)
        mae = np.abs(chronological_age - est).mean()
        rows.append({"clock": name, "MAE_years": round(mae, 2), "pearson_r": round(r, 3)})
    import pandas as pd
    return pd.DataFrame(rows).sort_values("MAE_years")

Pitfalls

• Cell-type composition confound: bulk blood methylation mixes cell types; report intrinsic EAA (adjusted for Houseman/reference-based deconvolution proportions), not raw extrinsic EAA, when claiming "cell-intrinsic aging"
• Missing clock CpGs: 450K vs EPIC array differences or QC-failed probes mean some canonical clock CpGs are absent — impute (mean or KNN across samples) before scoring, never silently drop and renormalize weights
• Batch effects: chip, row position, and processing date all shift beta values; normalize (Noob + BMIQ) and apply ComBat before comparing DNAmAge across batches
• Overfitting MAE: report cross-validated (not in-sample) predictions when fitting your own elastic net — in-sample MAE on thousands of CpGs with hundreds of samples is optimistic
• Causality: EAA-disease associations are correlational; don't imply an intervention reverses aging without longitudinal or Mendelian-randomization evidence
• Tissue specificity: a clock trained on blood does not transfer cleanly to other tissues — check clock performance in the target tissue or use a tissue-matched clock

See Also

• bio-methylation-analysis-methylation-calling — generate beta/M-values from bisulfite data before clock scoring
• bio-methylation-analysis-dmr-detection — differential methylation regions, a complementary analysis to age prediction
• bio-causal-genomics-mendelian-randomization — testing causal direction between EAA and disease outcomes
• bio-machine-learning-biomarker-discovery — general elastic-net/regularized regression workflow for omics biomarkers