表观遗传年龄预测 Horvath 时钟¶

# ===== R: 用 Horvath 时钟计算表观遗传年龄 =====

# 方法1：使用 methylclock 包（推荐，简单直接）
# BiocManager::install("methylclock")
library(methylclock)

# 输入数据：甲基化 beta 值矩阵
# 行 = CpG 探针（如 cg00000029），列 = 样本
# beta 值范围 0-1，0 = 未甲基化，1 = 完全甲基化
beta_matrix <- readRDS("beta_values.rds")

# 计算多种时钟的年龄
clock_results <- DNAmAge(beta_matrix,
    clocks = c("Horvath", "Hannum", "PhenoAge", "skinHorvath"))

# 查看结果
head(clock_results)
# 输出包含每个样本的 Horvath 年龄、Hannum 年龄等

# 方法2：手动用 Horvath 的 353 个 CpG 计算
# 下载 Horvath 的系数文件
# 文件包含：CpG 探针名、系数（coefficient）、截距
horvath_coef <- read.csv("horvath_coefficients.csv")

# 提取 353 个 CpG 的 beta 值
cpg_353 <- horvath_coef$CpG[-1]  # 第一行是截距
beta_sub <- beta_matrix[cpg_353, ]

# 年龄转换函数（Horvath 的反函数）
# 模型输出的是转换后的年龄，需要反转换
anti_trafo <- function(x, adult_age = 20) {
    ifelse(x < 0,
        (1 + adult_age) * exp(x) - 1,   # 儿童期（非线性）
        (1 + adult_age) * x + adult_age)  # 成年期（线性）
}

# 计算 DNAm Age
intercept <- horvath_coef$Coefficient[1]  # 截距
coefficients <- horvath_coef$Coefficient[-1]  # 353 个系数
# 线性组合
raw_age <- intercept + as.numeric(coefficients %*% beta_sub)
# 反转换得到最终年龄
dnam_age <- anti_trafo(raw_age)

# ===== 年龄加速分析 =====
# 已知实际年龄
actual_age <- c(25, 30, 45, 60, 70)  # 示例

# 计算年龄加速（残差法更准确）
# 用线性回归去除实际年龄的影响
model <- lm(dnam_age ~ actual_age)
age_acceleration <- residuals(model)
# 正值 = 加速衰老，负值 = 减缓衰老

# 可视化
library(ggplot2)
df <- data.frame(actual = actual_age, dnam = dnam_age,
                 acceleration = age_acceleration)

# 实际年龄 vs DNAm 年龄散点图
ggplot(df, aes(x = actual, y = dnam)) +
    geom_point(size = 3) +               # 散点
    geom_abline(slope = 1, intercept = 0, # 对角线（理想）
                linetype = "dashed", color = "red") +
    labs(x = "实际年龄", y = "DNAm 年龄",
         title = "Horvath 时钟年龄预测") +
    theme_minimal()