甲基化芯片分析（450K/EPIC芯片数据分析）¶

Beta = M_methylated / (M_methylated + M_unmethylated + alpha)
# alpha通常取100，防止分母为0

M_value = log2(Beta / (1 - Beta))

# 安装必要的R包
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("minfi", "IlluminaHumanMethylation450kanno.ilmn12.hg19",
                         "IlluminaHumanMethylation450kmanifest",
                         "IlluminaHumanMethylationEPICanno.ilm10b4.hg19",
                         "IlluminaHumanMethylationEPICmanifest"))

library(minfi)

# 准备样本信息表（Sample Sheet）
# CSV文件至少包含：Sample_Name, Sentrix_ID, Sentrix_Position, 分组信息
targets <- read.metharray.sheet(baseDir = "path/to/idat/", pattern = "SampleSheet.csv")

# 读取IDAT文件为RGChannelSet对象
rgSet <- read.metharray.exp(targets = targets)

# 查看基本信息
rgSet
pData(rgSet)  # 样本表型信息

BiocManager::install("ChAMP")
library(ChAMP)

# ChAMP一步读取+质控
myLoad <- champ.load(
  directory = "path/to/idat/",
  method = "ChAMP",           # 或 "minfi"
  methValue = "B",             # Beta值
  autoimpute = TRUE,           # 自动填补缺失值
  filterDetP = TRUE,           # 过滤检测P值不显著的探针
  ProbeCutoff = 0,             # 样本中失败探针的比例阈值
  SampleCutoff = 0.1,          # 探针在样本中失败的比例阈值
  detPcut = 0.01,              # 检测P值阈值
  filterBeads = TRUE,          # 过滤bead数不足的探针
  beadCutoff = 0.05,           # bead数阈值
  filterNoCG = TRUE,           # 过滤非CpG探针
  filterSNPs = TRUE,           # 过滤SNP相关探针
  filterMultiHit = TRUE,       # 过滤多重映射探针
  filterXY = TRUE,             # 过滤性染色体探针
  arraytype = "450K"           # 或 "EPIC"
)

# ===== minfi QC =====
library(minfi)

# 1. 检测P值
detP <- detectionP(rgSet)
# 统计每个样本失败探针比例
failed_fraction <- colMeans(detP > 0.01)
print(failed_fraction)
# 去除质量差的样本（>5%探针失败）
keep_samples <- failed_fraction < 0.05
rgSet <- rgSet[, keep_samples]

# 2. QC报告（生成PDF）
qcReport(rgSet, sampNames = pData(rgSet)$Sample_Name,
         sampGroups = pData(rgSet)$Group,
         pdf = "QC_report.pdf")

# 3. 密度图
densityPlot(rgSet, sampGroups = pData(rgSet)$Group, main = "Beta Value Density")

# 4. MDS图（多维缩放）检查批次效应和分组
mdsPlot(rgSet, sampGroups = pData(rgSet)$Group, numPositions = 1000)

# 5. 性别预测
predictedSex <- getSex(mapToGenome(rgSet))
plotSex(predictedSex)

# 6. SNP探针用于样本身份验证
snpBetas <- getSnpBeta(rgSet)
# 可用于检查样本是否搞混

# ===== ChAMP QC =====
# ChAMP在champ.load()时已完成大部分过滤
# 额外QC
champ.QC(
  beta = myLoad$beta,
  pheno = myLoad$pd$Sample_Group,
  mdsPlot = TRUE,
  densityPlot = TRUE,
  dendrogram = TRUE
)

# ===== minfi 预处理 =====

# 方法1：Funnorm（推荐，适合有大的生物学差异时）
grSet_funnorm <- preprocessFunnorm(rgSet)

# 方法2：Quantile归一化
grSet_quantile <- preprocessQuantile(rgSet)

# 方法3：Noob + BMIQ组合
mSet_noob <- preprocessNoob(rgSet)
# 然后可以额外做BMIQ

# 提取Beta值和M值
beta_values <- getBeta(grSet_funnorm)
m_values <- getM(grSet_funnorm)

# 过滤检测P值不合格的探针
detP_filtered <- detP[rownames(beta_values), colnames(beta_values)]
beta_values[detP_filtered > 0.01] <- NA

# 去除含NA过多的探针
keep_probes <- rowMeans(is.na(beta_values)) < 0.1
beta_values <- beta_values[keep_probes, ]
m_values <- m_values[keep_probes, ]

# ===== ChAMP归一化 =====
myNorm <- champ.norm(
  beta = myLoad$beta,
  rgSet = myLoad$rgSet,
  mset = myLoad$mset,
  resultsDir = "./CHAMP_Normalization/",
  method = "BMIQ",           # 或 "PBC", "SWAN", "FunctionalNormalization"
  plotBMIQ = TRUE,
  arraytype = "450K",
  cores = 4
)

# ===== ChAMP SVD分析 =====
champ.SVD(
  beta = myNorm,
  pd = myLoad$pd,
  resultsDir = "./CHAMP_SVD/"
)
# SVD热图会显示哪些技术因素与主成分显著相关

# ===== ComBat校正 =====
# 使用ChAMP
myCombat <- champ.runCombat(
  beta = myNorm,
  pd = myLoad$pd,
  batchname = c("Slide"),     # 批次变量名
  variablename = "Sample_Group" # 保护的生物学变量
)

# 使用sva包
library(sva)
# 需要M值做ComBat
m_values_combat <- ComBat(
  dat = m_values,
  batch = pData(grSet_funnorm)$Slide,
  mod = model.matrix(~ Group, data = pData(grSet_funnorm))
)

# ===== minfi + limma =====
library(limma)

# 构建设计矩阵
group <- factor(pData(grSet_funnorm)$Group)
design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)

# 对比矩阵
contrast_matrix <- makeContrasts(
  Disease_vs_Control = Disease - Control,
  levels = design
)

# 用M值做limma分析
fit <- lmFit(m_values, design)
fit2 <- contrasts.fit(fit, contrast_matrix)
fit3 <- eBayes(fit2)

# 提取结果
dmp_results <- topTable(fit3, coef = "Disease_vs_Control", number = Inf, sort.by = "p")

# 计算delta Beta
beta_disease <- rowMeans(beta_values[, group == "Disease"])
beta_control <- rowMeans(beta_values[, group == "Control"])
dmp_results$deltaBeta <- beta_disease[rownames(dmp_results)] - beta_control[rownames(dmp_results)]

# 筛选显著DMP
sig_dmp <- dmp_results[dmp_results$adj.P.Val < 0.05 & abs(dmp_results$deltaBeta) > 0.1, ]
nrow(sig_dmp)

# ===== ChAMP DMP =====
myDMP <- champ.DMP(
  beta = myCombat,
  pheno = myLoad$pd$Sample_Group,
  compare.group = c("Disease", "Control"),
  adjPVal = 0.05,
  adjust.method = "BH",
  arraytype = "450K"
)
# 结果是一个列表，每个比较一个数据框
head(myDMP[[1]])

# ===== 可视化 =====
# 火山图
library(ggplot2)
ggplot(dmp_results, aes(x = deltaBeta, y = -log10(adj.P.Val))) +
  geom_point(aes(color = (adj.P.Val < 0.05 & abs(deltaBeta) > 0.1)), size = 0.5) +
  scale_color_manual(values = c("grey", "red")) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  geom_vline(xintercept = c(-0.1, 0.1), linetype = "dashed") +
  theme_classic() +
  labs(x = "Delta Beta", y = "-log10(adj.P.Val)", title = "DMP Volcano Plot")

# 热图
library(pheatmap)
top50 <- head(sig_dmp, 50)
pheatmap(beta_values[rownames(top50), ],
         annotation_col = data.frame(Group = group, row.names = colnames(beta_values)),
         scale = "row",
         show_rownames = FALSE)

# ===== DMRcate =====
BiocManager::install("DMRcate")
library(DMRcate)

# 构建注释对象
myAnnotation <- cpg.annotate(
  datatype = "array",
  object = m_values,
  what = "M",
  analysis.type = "differential",
  design = design,
  contrasts = TRUE,
  cont.matrix = contrast_matrix,
  coef = "Disease_vs_Control",
  arraytype = "450K",          # 或 "EPIC"
  fdr = 0.05
)

# 提取DMR
dmr_results <- dmrcate(myAnnotation, lambda = 1000, C = 2)

# 转换为GRanges对象
dmr_ranges <- extractRanges(dmr_results, genome = "hg19")
dmr_ranges

# DMR可视化
DMR.plot(
  ranges = dmr_ranges,
  dmr = 1,                     # 第几个DMR
  CpGs = getBeta(grSet_funnorm),
  phen.col = c("blue", "red"),
  what = "Beta",
  arraytype = "450K",
  genome = "hg19"
)

# ===== ChAMP DMR =====
myDMR <- champ.DMR(
  beta = myCombat,
  pheno = myLoad$pd$Sample_Group,
  compare.group = c("Disease", "Control"),
  arraytype = "450K",
  method = "Bumphunter",       # 或 "DMRcate", "ProbeLasso"
  minProbes = 3,               # DMR至少包含的探针数
  adjPvalDmr = 0.05,
  cores = 4,
  resultsDir = "./CHAMP_DMR/"
)

# ===== missMethyl富集分析（推荐） =====
# missMethyl考虑了芯片设计偏倚：大基因上的探针更多，更容易被选中
BiocManager::install("missMethyl")
library(missMethyl)

# GO分析
go_results <- gometh(
  sig.cpg = rownames(sig_dmp),           # 显著DMP
  all.cpg = rownames(m_values),          # 所有检测的CpG
  collection = "GO",
  array.type = "450K",                   # 或 "EPIC"
  plot.bias = TRUE                       # 绘制偏倚校正图
)
# 显示前20个通路
topGO <- topGSA(go_results, number = 20)
print(topGO)

# KEGG分析
kegg_results <- gometh(
  sig.cpg = rownames(sig_dmp),
  all.cpg = rownames(m_values),
  collection = "KEGG",
  array.type = "450K"
)
topKEGG <- topGSA(kegg_results, number = 20)
print(topKEGG)

# ===== ChAMP GSEA =====
myGSEA <- champ.GSEA(
  beta = myCombat,
  DMP = myDMP[[1]],
  DMR = myDMR,
  pheno = myLoad$pd$Sample_Group,
  arraytype = "450K",
  adjPval = 0.05,
  method = "gometh"
)

# ===== 注释信息提取 =====
# 使用芯片注释包
library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
ann450k <- getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)

# 给DMP结果添加注释
sig_dmp_anno <- merge(sig_dmp, ann450k[, c("chr", "pos", "UCSC_RefGene_Name",
                                            "UCSC_RefGene_Group", "Relation_to_Island",
                                            "Regulatory_Feature_Group")],
                       by = "row.names")

# 统计各区域分布
table(sig_dmp_anno$UCSC_RefGene_Group)
table(sig_dmp_anno$Relation_to_Island)

# ===== 细胞组成估计（血液样本） =====
library(FlowSorted.Blood.450k)
# 或 FlowSorted.Blood.EPIC

cellCounts <- estimateCellCounts(
  rgSet,
  compositeCellType = "Blood",
  processMethod = "auto",
  cellTypes = c("CD8T", "CD4T", "NK", "Bcell", "Mono", "Gran")
)
head(cellCounts)

# 将细胞组成纳入差异分析模型（校正混杂）
design_adj <- model.matrix(~ 0 + group + cellCounts)

# ===== 甲基化年龄（Epigenetic Clock） =====
# 使用Horvath时钟或其他时钟
# 推荐使用在线工具：https://dnamage.genetics.ucla.edu/
# 或R包 methylclock
BiocManager::install("methylclock")
library(methylclock)

# 计算多种甲基化年龄
age_results <- DNAmAge(beta_values, clocks = c("Horvath", "Hannum", "PhenoAge", "skinHorvath"))

# ===== 拷贝数变异（CNV）推断 =====
# minfi内置功能
library(conumee)
# 需要参考样本（正常样本）
# cnv_data <- CNV.load(rgSet)
# cnv_fit <- CNV.fit(cnv_data[tumor_samples], cnv_data[normal_samples])
# CNV.genomeplot(cnv_fit)

# === 环境准备 ===
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("minfi", "limma", "DMRcate", "missMethyl",
                         "IlluminaHumanMethylation450kanno.ilmn12.hg19",
                         "IlluminaHumanMethylation450kmanifest",
                         "pheatmap", "ggplot2"))

library(minfi)
library(limma)
library(DMRcate)
library(missMethyl)
library(pheatmap)
library(ggplot2)

# === 读取数据 ===
targets <- read.metharray.sheet("./idat_dir/", pattern = "SampleSheet.csv")
rgSet <- read.metharray.exp(targets = targets)

# === 质控 ===
detP <- detectionP(rgSet)
keep <- colMeans(detP > 0.01) < 0.05
rgSet <- rgSet[, keep]
detP <- detP[, keep]

# === 归一化 ===
grSet <- preprocessFunnorm(rgSet)
beta <- getBeta(grSet)
mval <- getM(grSet)

# === 过滤 ===
# 检测P值不合格的探针
beta[detP[rownames(beta), colnames(beta)] > 0.01] <- NA
mval[detP[rownames(mval), colnames(mval)] > 0.01] <- NA
keep_probes <- rowMeans(is.na(beta)) == 0
beta <- beta[keep_probes, ]
mval <- mval[keep_probes, ]

# 过滤SNP探针和交叉反应探针
grSet_filtered <- dropLociWithSnps(grSet[keep_probes, ])
beta <- getBeta(grSet_filtered)
mval <- getM(grSet_filtered)

# 过滤性染色体
ann <- getAnnotation(grSet_filtered)
keep_autosome <- !(ann$chr %in% c("chrX", "chrY"))
beta <- beta[keep_autosome, ]
mval <- mval[keep_autosome, ]

# === DMP分析 ===
group <- factor(pData(grSet)$Group)
design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)
contMatrix <- makeContrasts(Disease - Control, levels = design)

fit <- lmFit(mval, design)
fit2 <- contrasts.fit(fit, contMatrix)
fit3 <- eBayes(fit2)
dmp <- topTable(fit3, number = Inf, coef = 1, sort.by = "p")

# 添加delta Beta
dmp$deltaBeta <- rowMeans(beta[rownames(dmp), group == "Disease"]) -
                 rowMeans(beta[rownames(dmp), group == "Control"])
sig_dmp <- dmp[dmp$adj.P.Val < 0.05 & abs(dmp$deltaBeta) > 0.1, ]

# === DMR分析 ===
myAnno <- cpg.annotate("array", mval, what = "M", analysis.type = "differential",
                         design = design, contrasts = TRUE, cont.matrix = contMatrix,
                         coef = 1, arraytype = "450K", fdr = 0.05)
dmr <- dmrcate(myAnno, lambda = 1000, C = 2)
dmr_ranges <- extractRanges(dmr, genome = "hg19")

# === 富集分析 ===
go <- gometh(sig.cpg = rownames(sig_dmp), all.cpg = rownames(mval),
             collection = "GO", array.type = "450K")
kegg <- gometh(sig.cpg = rownames(sig_dmp), all.cpg = rownames(mval),
               collection = "KEGG", array.type = "450K")

library(ChAMP)

# 一步加载+质控+过滤
myLoad <- champ.load(directory = "./idat_dir/", arraytype = "450K")

# 质控可视化
champ.QC(beta = myLoad$beta, pheno = myLoad$pd$Sample_Group)

# 归一化
myNorm <- champ.norm(beta = myLoad$beta, rgSet = myLoad$rgSet,
                      method = "BMIQ", arraytype = "450K", cores = 4)

# 批次效应检测
champ.SVD(beta = myNorm, pd = myLoad$pd)

# 批次校正（如需要）
myCombat <- champ.runCombat(beta = myNorm, pd = myLoad$pd, batchname = c("Slide"))

# DMP分析
myDMP <- champ.DMP(beta = myCombat, pheno = myLoad$pd$Sample_Group)

# DMR分析
myDMR <- champ.DMR(beta = myCombat, pheno = myLoad$pd$Sample_Group,
                    method = "Bumphunter", arraytype = "450K", cores = 4)

# 富集分析
myGSEA <- champ.GSEA(beta = myCombat, DMP = myDMP[[1]], DMR = myDMR,
                      pheno = myLoad$pd$Sample_Group, arraytype = "450K")

# 从GEO下载已有甲基化芯片数据
library(GEOquery)
gset <- getGEO("GSE123456", destdir = "./geo_data/")
# 提取Beta值矩阵
beta_geo <- exprs(gset[[1]])
pheno_geo <- pData(gset[[1]])