はじめに

老化は「避けられない運命」ではなく、データで観測し、AIで定量化できる 生物学的プロセス として扱えるようになってきました。本記事では、この考え方をエンジニア目線で咀嚼し、DNAメチル化データから生物学的年齢を推定する「老化時計（Aging Clock）」をPythonで実装する手順 を具体的に解説します。

参考にした解説記事：「老化」は避けられない運命ではなく、治療できる「病気」だった——AIが挑む人類最古の難問 | AI Robotics Quantum Lab

この記事のゴールは以下のとおりです。

環境構築からデータ取得、モデル学習・評価まで 手を動かして体験 できるようにする
2013年のHorvath時計を踏襲した ElasticNetベースのベースライン実装 を作る
後段の改善（マルチオミクス統合・ディープラーニング化）の 発展の道筋 を示す

対象読者は、PythonやscikitLearnの基本を知っていて、バイオインフォマティクスに触れてみたいエンジニア・データサイエンティストです。

1. 全体像 ― 何を作るのか

作るものは、次のようなシンプルなパイプラインです。

公開データセット（GEOリポジトリ）からDNAメチル化データ（β値）と年齢メタデータを取得
特徴量（CpGサイトのβ値）を前処理
ElasticNet回帰で「年齢」を目的変数として学習
生物学的年齢を予測し、暦年齢との差（Age Acceleration）を算出
評価指標（MAE・相関係数）で精度を確認

[GEO DNAm dataset] → [前処理] → [ElasticNet学習] → [推定年齢] → [Age Acceleration]

2. 環境構築

2.1 推奨環境

OS: macOS / Linux / WSL2
Python: 3.10 以上
メモリ: 16GB 推奨（CpGサイトは数十万次元あるため）

2.2 仮想環境の作成

# プロジェクトディレクトリ作成
mkdir aging-clock && cd aging-clock

# 仮想環境
python3 -m venv .venv
source .venv/bin/activate

# 必須ライブラリ
pip install --upgrade pip
pip install numpy pandas scikit-learn scipy matplotlib seaborn
pip install GEOparse     # GEOデータ取得
pip install pyreadr      # R形式(RData)対応（必要に応じて）
pip install jupyterlab   # ノート実行用

2.3 プロジェクト構成

aging-clock/
├── data/            # 生データ・前処理済みデータ
├── notebooks/       # 実験用ノート
├── src/
│   ├── data.py      # データ取得・前処理
│   ├── model.py     # モデル学習
│   └── evaluate.py  # 評価
└── README.md

3. データ取得 ― GEOからDNAメチル化データを取る

Horvath時計のオリジナル論文では複数のGEOデータセットが用いられています。ここでは導入として 小規模で扱いやすい公開データセット を使う前提で書きます。

⚠️ GEOデータは数GB単位になることもあります。ディスク容量とダウンロード時間に注意してください。

3.1 データ取得スクリプト

# src/data.py
import os
import GEOparse
import pandas as pd

def download_geo(accession: str, dest_dir: str = "./data") -> GEOparse.GSE:
    """GEOからデータセットをダウンロードして読み込む"""
    os.makedirs(dest_dir, exist_ok=True)
    gse = GEOparse.get_GEO(geo=accession, destdir=dest_dir, silent=True)
    return gse


def build_beta_matrix(gse: GEOparse.GSE) -> pd.DataFrame:
    """各サンプルのβ値を行列化（行: CpG、列: サンプル）"""
    frames = []
    for gsm_name, gsm in gse.gsms.items():
        # 列名は "VALUE" のことが多いが、データセットにより異なる
        table = gsm.table.set_index("ID_REF")
        col = "VALUE" if "VALUE" in table.columns else table.columns[0]
        frames.append(table[[col]].rename(columns={col: gsm_name}))
    return pd.concat(frames, axis=1)


def extract_age(gse: GEOparse.GSE, keyword: str = "age") -> pd.Series:
    """メタデータから年齢を抽出（characteristics_ch1を走査）"""
    ages = {}
    for gsm_name, gsm in gse.gsms.items():
        for item in gsm.metadata.get("characteristics_ch1", []):
            if keyword in item.lower():
                value = item.split(":")[-1].strip()
                try:
                    ages[gsm_name] = float(value)
                except ValueError:
                    continue
                break
    return pd.Series(ages, name="age")

3.2 実行例

from src.data import download_geo, build_beta_matrix, extract_age

# 例: 適切な血液サンプルのGEO accessionに置き換える
GEO_ID = "GSEXXXXXX"
gse = download_geo(GEO_ID)
beta = build_beta_matrix(gse)      # (CpG, sample)
age = extract_age(gse)             # (sample,)
print(beta.shape, age.shape)

💡 使用するデータセットは研究目的・ライセンスに沿ったものを選択してください。

4. 前処理 ― 欠損・定数列・ロジット変換

Horvath時計では、β値を対数変換（「Horvath transformation」）する工程があります。ここではシンプルな実装に留めます。

# src/data.py に追記
import numpy as np

def preprocess_beta(beta: pd.DataFrame) -> pd.DataFrame:
    # 1) サンプルを行、CpGを列に
    X = beta.T

    # 2) 欠損値を列平均で補完
    X = X.fillna(X.mean(axis=0))

    # 3) 定数列を除外
    X = X.loc[:, X.std(axis=0) > 1e-6]

    # 4) β値を[0.001, 0.999]にクリップしてから logit 変換
    X = X.clip(0.001, 0.999)
    X = np.log(X / (1 - X))
    return X

5. モデル学習 ― ElasticNetでHorvath風の老化時計を作る

Horvath時計のコアは ElasticNet回帰 です。高次元（数十万CpG）かつスパースなので相性が良いアルゴリズムです。

# src/model.py
from dataclasses import dataclass
import numpy as np
import pandas as pd
from sklearn.linear_model import ElasticNetCV
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error


@dataclass
class ClockResult:
    model: ElasticNetCV
    mae: float
    corr: float
    predictions: pd.Series


def train_clock(X: pd.DataFrame, y: pd.Series, random_state: int = 42) -> ClockResult:
    # 共通のサンプルだけを使う
    common = X.index.intersection(y.index)
    X = X.loc[common]
    y = y.loc[common]

    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, random_state=random_state
    )

    model = ElasticNetCV(
        l1_ratio=[0.1, 0.5, 0.7, 0.9, 0.95, 0.99, 1.0],
        n_alphas=50,
        cv=5,
        max_iter=20000,
        n_jobs=-1,
        random_state=random_state,
    )
    model.fit(X_train, y_train)

    pred = pd.Series(model.predict(X_test), index=X_test.index, name="pred_age")
    mae = mean_absolute_error(y_test, pred)
    corr = float(np.corrcoef(y_test.values, pred.values)[0, 1])

    return ClockResult(model=model, mae=mae, corr=corr, predictions=pred)

5.1 実行

from src.data import preprocess_beta
from src.model import train_clock

X = preprocess_beta(beta)
result = train_clock(X, age)

print(f"MAE : {result.mae:.2f} years")
print(f"Corr: {result.corr:.3f}")

血液サンプル主体の小規模データでも、Horvath時計と同程度の MAE 3〜5歳、相関 0.9前後 が目安になります。

6. 生物学的年齢と「年齢加速（Age Acceleration）」

生物学的年齢そのものより、暦年齢との差（residual）=Age Acceleration の方が健康状態と強く相関することが知られています。

# src/evaluate.py
import numpy as np
import pandas as pd


def age_acceleration(pred_age: pd.Series, chrono_age: pd.Series) -> pd.Series:
    """単純な差分に加え、加齢依存のバイアスを回帰で除去"""
    df = pd.concat([pred_age, chrono_age.rename("chrono")], axis=1).dropna()
    # 予測 ~ 暦年齢 の残差
    slope, intercept = np.polyfit(df["chrono"], df["pred_age"], 1)
    residual = df["pred_age"] - (slope * df["chrono"] + intercept)
    return residual.rename("age_acceleration")

正の値 → 暦年齢よりも生物学的に老けている（リスクが高い可能性）
負の値 → 暦年齢より若い（健康的な兆候かもしれない）

7. 可視化 ― モデルの挙動を目で確認する

import matplotlib.pyplot as plt

def plot_pred_vs_true(y_true: pd.Series, y_pred: pd.Series):
    fig, ax = plt.subplots(figsize=(5, 5))
    ax.scatter(y_true, y_pred, alpha=0.6)
    lims = [min(y_true.min(), y_pred.min()), max(y_true.max(), y_pred.max())]
    ax.plot(lims, lims, "k--", lw=1)
    ax.set_xlabel("Chronological age")
    ax.set_ylabel("Predicted biological age")
    ax.set_title("Aging clock: prediction vs. ground truth")
    plt.tight_layout()
    plt.show()

散布図が対角線に沿って分布し、外れ値がAge Accelerationの高い（または低い）サンプルになっていればOKです。

8. マルチオミクス統合への発展

本記事のベースラインはDNAメチル化のみを使いますが、実運用では複数のオミクスを組み合わせると精度と解釈性が上がります。

モダリティ	データソース例	役割
エピゲノム (DNAm)	Illumina 450K / EPIC array	老化の長期トレンドを捉える
プロテオーム	SomaLogic、Olink	「34・60・78歳」の転換点を捉える
メタボローム	Nightingale、Metabolon	代謝の瞬間的な状態を反映
マイクロバイオーム	16S / shotgun	腸内環境と炎症の指標

統合のアプローチとしては、次のような段階的な取り組みがわかりやすいです。

各オミクスごとにElasticNetで単独の時計を作る
それぞれの予測値をメタ特徴量として重回帰やGradient Boostingに投入する（スタッキング）
さらに進める場合は、マルチモーダルTransformer やVariational Autoencoderでジョイント表現を学習する

9. ディープラーニング版への差し替え

ベースラインが動いたら、以下のようにモデル部分だけ差し替えられます。

# 例: PyTorchでシンプルなMLP
import torch
from torch import nn

class AgingMLP(nn.Module):
    def __init__(self, in_dim: int, hidden: int = 256):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(in_dim, hidden),
            nn.ReLU(),
            nn.Dropout(0.3),
            nn.Linear(hidden, hidden),
            nn.ReLU(),
            nn.Dropout(0.3),
            nn.Linear(hidden, 1),
        )

    def forward(self, x):
        return self.net(x).squeeze(-1)

学習ループはHuberLossとAdamWあたりから始めると安定しやすいです。過学習しやすい領域なので、必ず独立したテストセット（できれば別コホート）で評価 してください。

10. 再現性と運用のチェックリスト

本番寄りで運用するなら、以下を押さえておくと安心です。

乱数シードを固定（numpy / scikitlearn / torch）
前処理・学習・評価を1コマンドで回せるMakefile or CLIを用意
モデル・前処理パイプラインを joblib などで保存
データセットのバージョン・アクセッション番号をREADMEに明記
推論時にも 学習時と同じ前処理（logit変換・欠損補完の基準値） を適用

11. 倫理・ライセンスに関する注意

公開データセットでも 利用規約・再配布条件 があります。必ず確認してください
個人ゲノム・エピゲノムは極めてセンシティブな情報です。医療用途に転用する場合は、国内の医療機器規制（薬機法・SaMD）や研究倫理審査を前提に設計してください
モデルの結果を 診断・医療アドバイスとして提示しない ことを明確に注記するのが安全です

まとめ

老化は「観測可能なデータ」と「AI」で定量化でき、ElasticNetのような古典的手法でもHorvath時計相当のベースラインが組める
本記事では、GEOからのデータ取得 → 前処理 → ElasticNet学習 → Age Acceleration算出まで、最短ルート を提示した
ベースラインが動いたら、マルチオミクス統合・ディープラーニング化 に進むのが次のステップ
医療応用を視野に入れる場合は、規制と倫理を前提にした設計が必須

「老化を攻略する最初の世代になる」ために、まずは手元のPythonで 自分の老化時計 を動かしてみるところから始めてみてください。