{ "cells": [ { "cell_type": "markdown", "id": "7d820e72", "metadata": {}, "source": [ "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/MeiraDavidson/biology_aging/blob/main/notebooks/build-your-own-epigenetic-clock.ipynb)\n\n", "# Build Your Own Epigenetic Clock 🧬⏱️\n", "### A hands-on companion to Chapter 2\n", "\n", "In the book you learned that an **epigenetic clock** is really a machine-learning model β€” a *penalized regression* that reads your age from a few hundred chemical marks on your DNA. Here you'll build one yourself, watch it predict age to within a few years… and then watch it get **fooled**. That last part is the most important lesson in the chapter, and it lands far harder when you make it happen with your own hands.\n", "\n", "**This notebook gives you the same clock two ways:**\n", "\n", "- **Part 1 β€” Simulated data.** We *invent* a tiny methylome. Runs in two seconds, needs no download, gives identical results for everyone. Perfect for understanding the idea.\n", "- **Part 2 β€” Real human data.** We download a slice of a *real published study* β€” the actual data behind one of the clocks the book talks about β€” and build the same clock on real people. Takes about half a minute and shows the lesson isn't a trick of our toy.\n", "\n", "Do Part 1 first. Then run Part 2 to see it on the real thing. Just press β–Ά on each cell, top to bottom.\n", "\n", "---" ] }, { "cell_type": "markdown", "id": "ea6de0ca", "metadata": {}, "source": [ "# Part 1 β€” A clock you can build in two seconds *(simulated data)*\n", "\n", "We start with invented data on purpose: it's instant, it's identical for everyone, and it lets us see the *machinery* clearly before we wrestle with a real dataset." ] }, { "cell_type": "markdown", "id": "3c27d183", "metadata": {}, "source": [ "## Step 0 β€” Grab our tools\n", "Standard scientific-Python libraries. All of these come pre-installed in Google Colab." ] }, { "cell_type": "code", "execution_count": null, "id": "9ebe78b5", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from sklearn.linear_model import ElasticNetCV, ElasticNet, LinearRegression\n", "from sklearn.model_selection import train_test_split\n", "from sklearn.metrics import mean_absolute_error\n", "\n", "rng = np.random.default_rng(42) # fixed seed -> everyone gets the same results\n", "print(\"Ready.\")" ] }, { "cell_type": "markdown", "id": "362daf2f", "metadata": {}, "source": [ "## Step 1 β€” Make a toy \"methylome\"\n", "\n", "**1,000 people**, each with **600 methylation sites**. A methylation value is a *beta-value* between 0 (never methylated) and 1 (always methylated).\n", "\n", "The crucial trick β€” and the whole point β€” is that **only 20 of the 600 sites actually drift with age.** The other 580 are pure noise. A good clock has to *find* the 20 needles in the 600-site haystack. (The real Horvath clock finds ~353 needles in a haystack of 450,000.)" ] }, { "cell_type": "code", "execution_count": null, "id": "9d2f096e", "metadata": {}, "outputs": [], "source": [ "N_PEOPLE, N_SITES, N_SIGNAL = 1000, 600, 20\n", "age = rng.uniform(20, 90, N_PEOPLE) # true chronological age\n", "\n", "signal_idx = rng.choice(N_SITES, N_SIGNAL, replace=False) # which sites track age\n", "slopes = rng.uniform(0.10, 0.28, N_SIGNAL) * rng.choice([-1, 1], N_SIGNAL)\n", "baseline = rng.uniform(0.30, 0.70, N_SITES)\n", "\n", "M = np.tile(baseline, (N_PEOPLE, 1)) + rng.normal(0, 0.11, (N_PEOPLE, N_SITES)) # noise\n", "for k, j in enumerate(signal_idx):\n", " M[:, j] += slopes[k] * (age - 55) / 35.0 # the 20 age-linked sites\n", "M = np.clip(M, 0, 1)\n", "\n", "print(f\"Methylome: {M.shape[0]} people x {M.shape[1]} sites\")\n", "print(f\"Only {N_SIGNAL} of {N_SITES} sites carry a real age signal. The rest is noise.\")" ] }, { "cell_type": "markdown", "id": "88ce499b", "metadata": {}, "source": [ "## Step 2 β€” Never grade yourself on data you trained on\n", "\n", "The book's *\"compared to what?\"* discipline at its most basic: fit the clock on a **training set**, then judge it on a **test set** it has never seen. A model scored only on data it already memorized is grading its own homework." ] }, { "cell_type": "code", "execution_count": null, "id": "3d34a057", "metadata": {}, "outputs": [], "source": [ "X_train, X_test, y_train, y_test = train_test_split(M, age, test_size=0.3, random_state=0)\n", "print(f\"Training on {len(y_train)} people, testing on {len(y_test)} the clock has never seen.\")" ] }, { "cell_type": "markdown", "id": "39823918", "metadata": {}, "source": [ "## Step 3 β€” Train the clock (an elastic-net regression)\n", "\n", "An **elastic net** is a linear regression with a penalty on large coefficients. Its L1 (\"lasso\") part does something magical: it drives most coefficients to **exactly zero**, so the model doesn't just weight the sites β€” it *selects* them. `ElasticNetCV` even tunes its own penalty by cross-validation. Watch how many of the 600 sites survive." ] }, { "cell_type": "code", "execution_count": null, "id": "10df995f", "metadata": {}, "outputs": [], "source": [ "clock = ElasticNetCV(l1_ratio=[.5, .7, .9, .95, 1.0], cv=5, max_iter=50000, random_state=0)\n", "clock.fit(X_train, y_train)\n", "\n", "pred = clock.predict(X_test)\n", "mae = mean_absolute_error(y_test, pred)\n", "print(f\"Test-set error (MAE): {mae:.1f} years\")\n", "print(f\"Sites the clock actually uses: {int(np.sum(clock.coef_!=0))} of {N_SITES} (true signal was {N_SIGNAL})\")\n", "\n", "plt.figure(figsize=(5.2,5.2))\n", "plt.scatter(y_test, pred, alpha=0.5, s=18)\n", "plt.plot([15,95],[15,95],'k--',lw=1,label=\"perfect prediction\")\n", "plt.xlabel(\"true age (years)\"); plt.ylabel(\"clock's predicted age (years)\")\n", "plt.title(f\"Your epigenetic clock - off by ~{mae:.1f} years\"); plt.legend(); plt.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "id": "81ba36a7", "metadata": {}, "source": [ "From 600 sites β€” 580 of them pure noise β€” the elastic net found the handful that matter and reads a stranger's age to within a few years. That's a real epigenetic clock in miniature.\n", "\n", "> *Carry-away #1, in code:* a clock is a strikingly **economical** model. Most of its 600 weights are zero." ] }, { "cell_type": "markdown", "id": "248d1751", "metadata": {}, "source": [ "## Step 4 β€” Would a bigger model do better? (No.)\n", "\n", "Tempting thought: *use every site and let a more powerful model sort it out.* Let's fit an ordinary regression on **all 600 sites** with no penalty, and compare its error on data it trained on versus new data." ] }, { "cell_type": "code", "execution_count": null, "id": "d2205368", "metadata": {}, "outputs": [], "source": [ "ols = LinearRegression().fit(X_train, y_train)\n", "ols_train = mean_absolute_error(y_train, ols.predict(X_train))\n", "ols_test = mean_absolute_error(y_test, ols.predict(X_test))\n", "print(f\"Big unpenalized model - train error {ols_train:.1f} yrs, test error {ols_test:.1f} yrs\")\n", "print(f\"Disciplined elastic net - test error {mae:.1f} yrs\")\n", "\n", "plt.figure(figsize=(5.5,3.6))\n", "plt.bar([\"Elastic net\\n(test)\",\"All-sites\\n(train)\",\"All-sites\\n(test)\"],\n", " [mae, ols_train, ols_test], color=[\"#2a9d8f\",\"#bbbbbb\",\"#e76f51\"])\n", "plt.ylabel(\"error (years)\"); plt.title(\"Memorizing noise is not learning age\")\n", "plt.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "id": "52e0df03", "metadata": {}, "source": [ "The big model nails the people it trained on and **falls apart on new ones** β€” it memorized noise. *More parameters do not create more reality.* (In the real field, elaborate deep-learning clocks perform only *comparably* to the humble elastic net, for exactly this reason.)" ] }, { "cell_type": "markdown", "id": "4295057d", "metadata": {}, "source": [ "## Step 5 β€” The punchline: a thermometer is not a thermostat\n", "\n", "The most important cell in the notebook. Your clock reads age off a specific set of sites. Imagine a \"treatment\" that changes **none** of a person's biology β€” it just nudges those exact sites toward youthful values. How little do we have to nudge to fake years of age?" ] }, { "cell_type": "code", "execution_count": null, "id": "9bff1abc", "metadata": {}, "outputs": [], "source": [ "def fake_younger(clock, X, years):\n", " \"\"\"Nudge only the sites the clock uses, toward 'younger', just enough to fake `years`.\"\"\"\n", " used = np.where(clock.coef_ != 0)[0]\n", " leverage = np.sum(np.abs(clock.coef_[used])) # 'years' of clock movement per unit nudge\n", " nudge = years / leverage # how far to shift each used site\n", " X_gamed = X.copy()\n", " X_gamed[:, used] = np.clip(X[:, used] - nudge*np.sign(clock.coef_[used]), 0, 1)\n", " return X_gamed, nudge, len(used)\n", "\n", "FAKE_YEARS = 12 # <-- try changing this and re-running\n", "\n", "X_gamed, nudge, n_sites = fake_younger(clock, X_test, FAKE_YEARS)\n", "pred_gamed = clock.predict(X_gamed)\n", "\n", "print(f\"We nudged each of the clock's {n_sites} sites by only {nudge:.3f} on the 0-1 scale.\")\n", "print(f\"Mean predicted age BEFORE: {pred.mean():.1f} years\")\n", "print(f\"Mean predicted age AFTER : {pred_gamed.mean():.1f} years\")\n", "print(f\"Apparent rejuvenation: {(pred-pred_gamed).mean():.1f} years - and nobody's true age moved a day.\")\n", "\n", "plt.figure(figsize=(5.4,5.4))\n", "plt.scatter(y_test, pred, alpha=0.5, s=18, color=\"#264653\", label=\"honest reading\")\n", "plt.scatter(y_test, pred_gamed, alpha=0.5, s=18, color=\"#e76f51\", label=\"after gaming\")\n", "plt.plot([15,95],[15,95],'k--',lw=1)\n", "plt.xlabel(\"true age (years) - UNCHANGED\"); plt.ylabel(\"clock's predicted age (years)\")\n", "plt.title(\"Same people. Same true age. 'Younger' clock.\"); plt.legend(); plt.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "id": "910d39bf", "metadata": {}, "source": [ "You just \"reversed aging\" by years β€” on paper β€” with a nudge so small it would be invisible to any biology. You taught to the test.\n", "\n", "This is the exact danger the chapter warns about: a drug screen rewarding *\"lowers the clock\"* will hand you the compounds best at **gaming the clock**, which is not the same as the ones that slow aging β€” *until the clock is proven a valid surrogate*. Set `FAKE_YEARS` higher and re-run: the clock will fake as many years as you ask, for almost nothing. That should worry you β€” and now you understand, in your fingers, why." ] }, { "cell_type": "markdown", "id": "30d40258", "metadata": {}, "source": [ "---\n", "# Part 2 β€” The same clock on *real human data*\n", "\n", "Part 1 ran on numbers we invented. Now let's do it for real, so you can see the lesson is not a trick of our toy.\n", "\n", "### Where this data comes from\n", "\n", "We'll use a famous public dataset called **GSE40279**, hosted on the **Gene Expression Omnibus (GEO)** β€” a giant free archive run by the US National Center for Biotechnology Information (NCBI) where scientists deposit their raw data so anyone can check and reuse it.\n", "\n", "This particular dataset is the one behind **Hannum et al. (2013)** β€” the study that built one of the first epigenetic clocks, and one the book cites in Chapter 2. It contains:\n", "\n", "- **656 real people**, aged **19 to 101**,\n", "- DNA methylation measured in their **whole blood**,\n", "- on the **Illumina Infinium 450K** chip, which reads **~450,000 CpG sites** per person.\n", "\n", "These are real, de-identified human samples, released publicly so that anyone β€” including you, right now β€” can study them. A *methylation beta-value* is exactly what we simulated: a number from 0 (this CpG is never methylated) to 1 (always methylated).\n", "\n", "### One catch: real data is *big*\n", "\n", "The full file of beta-values is about **1.2 gigabytes** (compressed!). We don't need all 450,000 sites to learn the lesson, so we'll **stream just the first 10,000** of them β€” enough to build a working clock in under a minute. (The ages live in a separate, tiny description file.)" ] }, { "cell_type": "markdown", "id": "6d29f0bf", "metadata": {}, "source": [ "**Step A β€” Get everyone's age** from the dataset's small metadata header. (This downloads only the top of the description file, so it's instant.)" ] }, { "cell_type": "code", "execution_count": null, "id": "b8f94fb4", "metadata": {}, "outputs": [], "source": [ "import urllib.request, gzip\n", "import pandas as pd\n", "\n", "BASE = \"https://ftp.ncbi.nlm.nih.gov/geo/series/GSE40nnn/GSE40279\"\n", "\n", "ids = ages = None\n", "with urllib.request.urlopen(f\"{BASE}/matrix/GSE40279_series_matrix.txt.gz\") as resp:\n", " for raw in gzip.GzipFile(fileobj=resp):\n", " line = raw.decode(\"utf-8\", \"replace\")\n", " if line.startswith(\"!Sample_source_name_ch1\"):\n", " ids = [c.strip().strip('\"') for c in line.split(\"\\t\")[1:]] # sample IDs\n", " elif line.startswith(\"!Sample_characteristics_ch1\") and \"age (y)\" in line:\n", " ages = [float(c.split(\":\")[1].strip().strip('\"')) for c in line.split(\"\\t\")[1:]]\n", " if line.startswith(\"!series_matrix_table_begin\"):\n", " break\n", "age_by_id = dict(zip(ids, ages))\n", "print(f\"Got ages for {len(age_by_id)} real people, ranging {min(ages):.0f} to {max(ages):.0f} years old.\")" ] }, { "cell_type": "markdown", "id": "afd84420", "metadata": {}, "source": [ "**Step B β€” Stream the first 10,000 methylation sites** for those people, and take a look at the real numbers." ] }, { "cell_type": "code", "execution_count": null, "id": "2fba3445", "metadata": {}, "outputs": [], "source": [ "N_SITES_REAL = 10000 # of ~450,000; we stream just a slice to keep it quick\n", "\n", "beta = pd.read_csv(f\"{BASE}/suppl/GSE40279_average_beta.txt.gz\",\n", " sep=\"\\t\", index_col=0, nrows=N_SITES_REAL, compression=\"gzip\")\n", "print(f\"Downloaded a {beta.shape[0]} sites x {beta.shape[1]} people block of REAL methylation data.\\n\")\n", "print(\"A peek (rows = CpG sites, columns = people, values = methylation 0-1):\")\n", "display(beta.iloc[:4, :5])\n", "\n", "X = beta.T # rows = people, columns = sites\n", "X = X.fillna(X.mean()) # a few missing values -> column average\n", "y = np.array([age_by_id[s] for s in X.index]) # line each person up with their real age\n", "\n", "plt.figure(figsize=(5.5,3.2)); plt.hist(y, bins=25, color=\"#457b9d\")\n", "plt.xlabel(\"real age (years)\"); plt.ylabel(\"number of people\")\n", "plt.title(f\"{len(y)} real people, ages {y.min():.0f}-{y.max():.0f}\"); plt.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "id": "f1038178", "metadata": {}, "source": [ "**Step C β€” Build the clock on real people.** Same idea as Part 1: split into train/test, fit a penalized regression, judge it on people it never saw. (Real data has age signal spread across many sites, so we use a single fast elastic-net fit rather than the slower self-tuning one.)" ] }, { "cell_type": "code", "execution_count": null, "id": "ff3ae8d0", "metadata": {}, "outputs": [], "source": [ "Xtr, Xte, ytr, yte = train_test_split(X.values, y, test_size=0.3, random_state=0)\n", "\n", "real_clock = ElasticNet(alpha=0.01, l1_ratio=0.5, max_iter=10000).fit(Xtr, ytr)\n", "pred_real = real_clock.predict(Xte)\n", "mae_real = mean_absolute_error(yte, pred_real)\n", "print(f\"Real-data clock test error (MAE): {mae_real:.1f} years\")\n", "print(f\"Sites used: {int(np.sum(real_clock.coef_!=0))} of {N_SITES_REAL}\")\n", "\n", "plt.figure(figsize=(5.2,5.2))\n", "plt.scatter(yte, pred_real, alpha=0.5, s=18, color=\"#1d3557\")\n", "plt.plot([15,105],[15,105],'k--',lw=1,label=\"perfect prediction\")\n", "plt.xlabel(\"true age (years)\"); plt.ylabel(\"clock's predicted age (years)\")\n", "plt.title(f\"A clock built from REAL blood - off by ~{mae_real:.1f} years\"); plt.legend(); plt.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "id": "8175f643", "metadata": {}, "source": [ "From a slice of one real study, you just built a clock that reads a real stranger's age to within a few years. This is, in miniature, exactly what Hannum and others did in 2013." ] }, { "cell_type": "markdown", "id": "2a3509f9", "metadata": {}, "source": [ "**Step D β€” Now fool the real clock.** The exact same `fake_younger` trick from Part 1 β€” nudge only the sites this clock uses, toward youth, changing no biology." ] }, { "cell_type": "code", "execution_count": null, "id": "5aee3373", "metadata": {}, "outputs": [], "source": [ "X_gamed_real, nudge_real, n_real = fake_younger(real_clock, Xte, FAKE_YEARS)\n", "pred_real_gamed = real_clock.predict(X_gamed_real)\n", "\n", "print(f\"We nudged each of the real clock's {n_real} sites by only {nudge_real:.4f} on the 0-1 scale.\")\n", "print(f\"Mean predicted age BEFORE: {pred_real.mean():.1f} years\")\n", "print(f\"Mean predicted age AFTER : {pred_real_gamed.mean():.1f} years\")\n", "print(f\"Apparent rejuvenation: {(pred_real-pred_real_gamed).mean():.1f} years - on REAL data, with no biology changed.\")\n", "\n", "plt.figure(figsize=(5.4,5.4))\n", "plt.scatter(yte, pred_real, alpha=0.5, s=18, color=\"#1d3557\", label=\"honest reading\")\n", "plt.scatter(yte, pred_real_gamed, alpha=0.5, s=18, color=\"#e76f51\", label=\"after gaming\")\n", "plt.plot([15,105],[15,105],'k--',lw=1)\n", "plt.xlabel(\"true age (years) - UNCHANGED\"); plt.ylabel(\"clock's predicted age (years)\")\n", "plt.title(\"Real clock, real blood - and just as easy to fool\"); plt.legend(); plt.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "id": "56933ba4", "metadata": {}, "source": [ "On real human data the nudge needed is **even tinier** than on our toy β€” because the signal is spread across more sites, so each one has to move only a hair. The warning from the chapter isn't an artifact of a simulation. It's a real property of real clocks." ] }, { "cell_type": "markdown", "id": "04681d4d", "metadata": {}, "source": [ "---\n", "# Both ways, side by side\n", "\n", "You've now built the same clock two ways. Here's why you'd reach for each:\n", "\n", "| | **Part 1 β€” Simulated** | **Part 2 β€” Real (GSE40279)** |\n", "|---|---|---|\n", "| Speed | Instant | ~30 seconds (a download) |\n", "| Reproducible? | Identical for everyone | Depends on the live dataset |\n", "| Realism | A teaching toy | Actual human blood, real ages |\n", "| Best for | *Seeing the machinery* | *Trusting that it's real* |\n", "| The lesson | Predictor β‰  lever | **Same lesson β€” confirmed on real people** |\n", "\n", "The point of doing both: the simulation makes the idea *clear*, and the real data makes it *true*. A scientist needs both habits β€” build the simplest version you can understand, then test it against reality." ] }, { "cell_type": "markdown", "id": "f757d17e", "metadata": {}, "source": [ "## What you just learned (worth reading twice)\n", "\n", "You built a real epigenetic clock β€” twice β€” and rediscovered the chapter's three carry-aways by *doing* them:\n", "\n", "1. **A clock is a machine-learning model** β€” economical, sparse, accurate. *(Step 3 / Step C)*\n", "2. **A measurement that predicts an outcome is not a handle that controls it.** You moved the readout without moving the thing, on both fake and real data. *(Step 5 / Step D)*\n", "3. **Precision doesn't change that** β€” and a bigger model doesn't help. *(Step 4)*\n", "\n", "The promising frontier runs the *other* way: don't optimize the clock β€” **intervene on biology and measure a real outcome.** In a real lab that means editing a single methylation site (a disabled Cas9 fused to a methylation-writing or -erasing enzyme) and watching what the cell actually does β€” the one experiment a clock can never run on its own." ] }, { "cell_type": "markdown", "id": "31150b77", "metadata": {}, "source": [ "## Go further (for the reader who finished hungry)\n", "\n", "- **Pull a different GEO dataset.** Search GEO (ncbi.nlm.nih.gov/geo) for another blood-methylation series with ages and swap in its accession.\n", "- **Stream more sites.** Raise `N_SITES_REAL` toward 450,000 and watch the real clock sharpen (and slow down) β€” you'll *feel* why real pipelines pre-select the most informative sites.\n", "- **Try a neural network** (`MLPRegressor`) and confirm for yourself it's *not* meaningfully better than the elastic net.\n", "- **Change the training target.** Train against a mortality-risk score instead of calendar age β€” the single move that took the field from PhenoAge to GrimAge to DunedinPACE.\n", "- **Build the Chapter 3 companion.** Simulate the senescent-cell tipping point with `dN/dt = P - c*N/(K+N)` and find the threshold where maintenance loses to damage.\n", "\n", "---\n", "*Companion notebook to* **Chapter 2 - Epigenetic Clocks**, *from the \"AI for Drug Discovery for Kids\" aging & longevity volume. Part 1 uses simulated data; Part 2 uses the public Hannum et al. (2013) dataset, GEO accession GSE40279. For teaching.*" ] } ], "metadata": { "colab": { "provenance": [] }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 5 }