171_表观时钟分析¶

# === 从IDAT文件处理甲基化芯片数据 ===
library(minfi)
library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)

# 读取IDAT文件
rgSet <- read.metharray.exp("idat_directory/")

# 质控
qc <- getQC(preprocessRaw(rgSet))
plotQC(qc)

# 归一化 (推荐BMIQ或Noob)
mSet <- preprocessNoob(rgSet)  # 背景校正+染料偏差校正

# 获取Beta值矩阵
beta_matrix <- getBeta(mSet)
# Beta值范围[0,1]: 0=未甲基化, 1=完全甲基化

# 探针过滤
# 去除检测P值>0.01的探针
detP <- detectionP(rgSet)
failed_probes <- rownames(detP)[rowMeans(detP > 0.01) > 0.05]
beta_filtered <- beta_matrix[!rownames(beta_matrix) %in% failed_probes, ]

# 去除SNP探针和cross-reactive探针
# (使用已发表的列表)
beta_clean <- beta_filtered[!rownames(beta_filtered) %in% cross_reactive_probes, ]

cat(sprintf("保留探针数: %d\n", nrow(beta_clean)))

# === 使用methylclock包 ===
# BiocManager::install("methylclock")
library(methylclock)

# 计算多个时钟
clock_results <- DNAmAge(
  beta_clean,
  clocks = c("Horvath", "Hannum", "skinHorvath", "PhenoAge", "Zhang")
)

# === 使用Horvath在线计算器格式 ===
# https://dnamage.clockfoundation.org/
# 需要上传beta值矩阵(行=CpG探针ID, 列=样本)

# === 手动计算Horvath时钟(概念) ===
# Horvath时钟 = 反变换(截距 + Σ(系数_i × 甲基化_i))
# 353个CpG的系数已发表
horvath_coefs <- read.csv("horvath_coefficients.csv")  # cpg_id, coefficient

# 匹配CpG
common_cpgs <- intersect(horvath_coefs$cpg_id, rownames(beta_clean))
beta_horvath <- beta_clean[common_cpgs, ]
coefs <- horvath_coefs$coefficient[match(common_cpgs, horvath_coefs$cpg_id)]

# 线性预测
linear_pred <- as.numeric(t(coefs) %*% beta_horvath) + horvath_intercept

# 反变换(Horvath使用的变换)
horvath_transform_inv <- function(x) {
  ifelse(x < 0, (1 + adult_age) * exp(x) - 1,
         (1 + adult_age) * x + adult_age)
}
dnam_age <- horvath_transform_inv(linear_pred)

# === GrimAge和DunedinPACE ===
# 使用PC-based时钟 (更推荐)
# devtools::install_github("danbelsky/DunedinPACE")
library(DunedinPACE)

# DunedinPACE计算
pace_results <- PACEProjector(beta_clean)
# 输出: 衰老速率(1.0=正常速率; >1=加速衰老; <1=减缓衰老)

# === PC Clocks (Higgins-Chen et al. 2022) ===
# 更鲁棒的主成分版本
# devtools::install_github("MorganLevineLab/PC-Clocks")

library(ggplot2)

# === 计算年龄加速 ===
results <- data.frame(
  SampleID = colnames(beta_clean),
  Chronological_Age = sample_metadata$age,
  Horvath_Age = clock_results$Horvath,
  Hannum_Age = clock_results$Hannum,
  PhenoAge = clock_results$PhenoAge
)

# 方法1: 简单差值
results$AgeAccel_Horvath <- results$Horvath_Age - results$Chronological_Age

# 方法2: 残差法(推荐) - 去除年龄本身的影响
lm_fit <- lm(Horvath_Age ~ Chronological_Age, data = results)
results$AgeAccel_residual <- residuals(lm_fit)
# 正残差=比同龄人"老"; 负残差=比同龄人"年轻"

# === 与疾病/表型关联 ===
# 比较病例vs对照的年龄加速
results$Disease <- sample_metadata$disease_status

ggplot(results, aes(x = Disease, y = AgeAccel_residual, fill = Disease)) +
  geom_boxplot() +
  geom_jitter(width = 0.2, alpha = 0.3) +
  stat_compare_means(method = "wilcox.test") +
  theme_minimal() +
  labs(title = "Epigenetic Age Acceleration by Disease Status",
       y = "Age Acceleration (years)")

# === 生存分析 ===
library(survival)
library(survminer)

results$AgeAccel_group <- ifelse(results$AgeAccel_residual > 0, "Accelerated", "Normal")
fit <- survfit(Surv(OS_time, OS_event) ~ AgeAccel_group, data = results)
ggsurvplot(fit, pval = TRUE, risk.table = TRUE,
           title = "Survival by Epigenetic Age Acceleration")

# Cox回归(连续变量)
cox_fit <- coxph(Surv(OS_time, OS_event) ~ AgeAccel_residual + age + sex, data = results)
summary(cox_fit)

# === DunedinPACE评估抗衰老干预效果 ===
# 干预前后配对比较

pre_pace <- pace_results[metadata$timepoint == "Baseline"]
post_pace <- pace_results[metadata$timepoint == "Post_intervention"]

# 配对t检验
t_test <- t.test(pre_pace, post_pace, paired = TRUE)
cat(sprintf("干预前PACE: %.3f ± %.3f\n", mean(pre_pace), sd(pre_pace)))
cat(sprintf("干预后PACE: %.3f ± %.3f\n", mean(post_pace), sd(post_pace)))
cat(sprintf("变化: %.3f (p = %.4f)\n", mean(post_pace - pre_pace), t_test$p.value))

# DunedinPACE解读:
# 1.0 = 每年老1年(正常)
# 1.2 = 每年老1.2年(加速20%)
# 0.8 = 每年老0.8年(减缓20%)
# 有效干预应使PACE降低

# === 效应量可视化 ===
ggplot(data.frame(pre = pre_pace, post = post_pace), aes(x = pre, y = post)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
  labs(x = "Baseline DunedinPACE", y = "Post-intervention DunedinPACE",
       title = "Intervention Effect on Pace of Aging") +
  theme_minimal()

# === 甲基化时钟分析流程 ===
# 1. IDAT处理
Rscript process_idat.R --idat_dir data/ --output beta_matrix.rds

# 2. 计算时钟
Rscript calculate_clocks.R --beta beta_matrix.rds \
  --clocks "Horvath,Hannum,PhenoAge,GrimAge" --output clock_ages.csv

# 3. 年龄加速分析
Rscript age_acceleration.R --clock_ages clock_ages.csv \
  --metadata clinical.csv --output results/