DNA 甲基化 450K/850K 芯片分析¶

一句话概述¶

Illumina 甲基化芯片（450K/EPIC 850K）是检测全基因组 CpG 位点甲基化水平的高通量平台，使用 R 语言的 minfi/ChAMP 包进行分析。

核心知识点表格¶

知识点	说明
450K 芯片	覆盖约 485,000 个 CpG 位点，2011年推出
EPIC/850K 芯片	覆盖约 850,000 个 CpG 位点，2016年推出，增加了增强子区域
EPICv2	2023年推出的最新版本，约 935,000 个 CpG 位点
Beta 值	甲基化水平 0-1，直观但统计性质差
M 值	log2(Beta/(1-Beta))，统计分析更合适
Infinium I/II	两种探针类型，Type II 有技术偏差需要校正
minfi	Bioconductor 最流行的甲基化芯片分析包
ChAMP	整合型分析流水线，封装了多种方法
DMP	Differentially Methylated Position，差异甲基化位点
DMR	Differentially Methylated Region，差异甲基化区域

白话解释原理¶

想象一下： DNA 上的某些位置（CpG）可以被加上一个小标记（甲基化），就像在书本上做标记。被标记的基因通常被"关闭"了。

甲基化芯片做的事情： 1. 用亚硫酸盐处理 DNA → 未甲基化的 C 变成 U（T），甲基化的 C 不变 2. 芯片上有探针分别识别甲基化和未甲基化的 DNA 3. 读取荧光信号 → 计算每个位点的甲基化比例

Beta 值 = 甲基化信号 / (甲基化信号 + 未甲基化信号) - Beta = 0 → 完全未甲基化 - Beta = 1 → 完全甲基化 - Beta = 0.5 → 一半甲基化

各步骤详解¶

第一步：数据读取（R 语言）¶

# 安装必要的包
if (!require("BiocManager")) install.packages("BiocManager")  # Bioconductor管理器
BiocManager::install("minfi")         # 安装minfi
BiocManager::install("ChAMP")        # 安装ChAMP
BiocManager::install("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")  # EPIC注释包

# === 方法1: 使用 minfi ===
library(minfi)  # 加载minfi包

# 读取 idat 文件
# idat 文件是芯片扫描仪输出的原始数据文件（红/绿通道各一个）
baseDir <- "/data/methylation/"       # idat文件目录
targets <- read.metharray.sheet(baseDir)  # 读取样本信息表
rgSet <- read.metharray.exp(targets = targets)  # 读取所有idat文件

# 查看基本信息
cat("样本数:", ncol(rgSet), "\n")     # 打印样本数
cat("探针数:", nrow(rgSet), "\n")     # 打印探针数

# === 方法2: 使用 ChAMP（更简单） ===
library(ChAMP)  # 加载ChAMP包

# ChAMP 一步读取 + QC
myLoad <- champ.load(
    directory = baseDir,              # idat文件目录
    filterDetP = TRUE,                # 过滤检测p值不合格的探针
    filterBeads = TRUE,               # 过滤bead数太少的探针
    filterNoCG = TRUE,                # 过滤非CpG探针
    filterSNPs = TRUE,                # 过滤SNP位置的探针
    filterMultiHit = TRUE,            # 过滤多重比对的探针
    filterXY = TRUE,                  # 过滤性染色体探针
    arraytype = "EPIC"                # 指定芯片类型: "450K" 或 "EPIC"
)
# ChAMP 会自动生成QC报告

第二步：质量控制¶

# === 使用 minfi 做 QC ===
# 1. 检测 p 值过滤（信号可靠性）
detP <- detectionP(rgSet)             # 计算检测p值
failed <- detP > 0.01                 # p > 0.01 的探针不可靠
cat("不合格探针比例:", sum(failed) / length(failed), "\n")  # 打印比例

# 过滤不合格样本（>5%探针失败的样本需移除）
keep_samples <- colMeans(detP) < 0.05  # 保留合格样本
rgSet <- rgSet[, keep_samples]         # 过滤

# 2. 密度图检查（观察 Beta 值分布）
mSet <- preprocessRaw(rgSet)           # 原始预处理
beta_raw <- getBeta(mSet)              # 获取Beta值
densityPlot(beta_raw,                  # 绘制密度图
            sampGroups = targets$Group, # 按组着色
            main = "Raw Beta 值分布")   # 标题
# 应该看到两个峰：0附近（未甲基化）和1附近（甲基化）

# 3. QC 报告
qcReport(rgSet, pdf = "QC_report.pdf")  # 生成PDF报告

# 4. 性别预测（检查样本标记是否正确）
predictedSex <- getSex(mapToGenome(mSet))  # 预测性别
table(predictedSex$predictedSex)            # 查看预测结果

第三步：标准化¶

# === 方法选择指南 ===
# BMIQ — 校正 Type I/II 探针偏差，ChAMP 默认
# Noob (ssNoob) — 背景校正+染料偏差校正，minfi 推荐用于多平台整合
# Funnorm — 功能归一化，适合大规模队列研究
# SeSAMe — 2025年最新研究认为效果最好

# === 使用 minfi 的 Noob + BMIQ ===
# 1. Noob 背景校正
mSetNoob <- preprocessNoob(rgSet)      # Noob预处理

# 2. 获取 Beta 值和 M 值
beta <- getBeta(mSetNoob)              # Beta值 (0-1)
mval <- getM(mSetNoob)                # M值 (log2转换)

# === 使用 ChAMP 标准化（推荐新手） ===
myNorm <- champ.norm(
    beta = myLoad$beta,               # 输入Beta值
    rgSet = rgSet,                     # 原始数据
    mset = mSet,                       # 预处理数据
    method = "BMIQ",                   # 标准化方法（BMIQ推荐）
    arraytype = "EPIC"                 # 芯片类型
)
# 输出 myNorm 是标准化后的 Beta 矩阵

# === 批次效应检测和校正 ===
# 1. SVD 检测批次效应
champ.SVD(
    beta = myNorm,                     # 标准化后Beta值
    pd = myLoad$pd                     # 样本信息
)
# 如果看到 Slide/Array 与主成分显著相关 → 有批次效应

# 2. ComBat 校正
myNorm_combat <- champ.runCombat(
    beta = myNorm,                     # 输入
    pd = myLoad$pd,                    # 样本信息
    batchname = c("Slide")             # 批次变量
)

第四步：差异甲基化分析¶

# === DMP 分析（单个CpG位点） ===
# 使用 ChAMP
myDMP <- champ.DMP(
    beta = myNorm_combat,              # 标准化Beta值
    pheno = myLoad$pd$Group,           # 分组信息
    adjPVal = 0.05,                    # 校正p值阈值
    adjust.method = "BH",             # BH校正
    arraytype = "EPIC"                 # 芯片类型
)
# myDMP 返回每个对比的差异CpG列表

# 使用 limma（更灵活）
library(limma)  # 加载limma

design <- model.matrix(~ 0 + Group, data = myLoad$pd)  # 设计矩阵
colnames(design) <- levels(factor(myLoad$pd$Group))     # 列名
fit <- lmFit(mval, design)                               # 拟合（用M值!）
contrast <- makeContrasts(Treatment - Control, levels = design)  # 对比
fit2 <- contrasts.fit(fit, contrast)                    # 应用对比
fit2 <- eBayes(fit2)                                    # 贝叶斯修正

dmp_results <- topTable(fit2, number = Inf,             # 所有结果
                         adjust.method = "BH")          # BH校正

# 筛选显著 DMP
sig_dmp <- dmp_results[dmp_results$adj.P.Val < 0.05 &   # FDR < 0.05
                        abs(dmp_results$logFC) > 0.5, ]  # |Delta M| > 0.5
cat("显著DMP数量:", nrow(sig_dmp), "\n")

# === DMR 分析（连续区域） ===
# 使用 DMRcate
BiocManager::install("DMRcate")       # 安装DMRcate
library(DMRcate)                       # 加载

# 用 limma 结果做 DMR
cpg_annotated <- cpg.annotate(        # 注释CpG位点
    "array",                           # 数据类型
    mval,                              # M值矩阵
    analysis.type = "differential",    # 差异分析
    design = design,                   # 设计矩阵
    contrasts = TRUE,                  # 有对比
    cont.matrix = contrast,            # 对比矩阵
    coef = "Treatment - Control",      # 对比名称
    arraytype = "EPIC"                 # 芯片类型
)

dmr_results <- dmrcate(cpg_annotated,  # 调用DMR
                        lambda = 1000,  # 平滑窗口1kb
                        C = 2)          # CpG密度参数

# 提取 DMR 坐标
dmr_ranges <- extractRanges(dmr_results, genome = "hg19")  # 提取为GRanges
cat("DMR数量:", length(dmr_ranges), "\n")                  # 打印数量

第五步：可视化和生物学解读¶

# 1. 火山图
library(ggplot2)
ggplot(dmp_results, aes(x = logFC, y = -log10(adj.P.Val))) +
    geom_point(aes(color = adj.P.Val < 0.05 & abs(logFC) > 0.5),
               size = 0.5, alpha = 0.5) +
    scale_color_manual(values = c("grey", "red")) +
    labs(title = "甲基化差异位点火山图",
         x = "Delta M value", y = "-log10(FDR)") +
    theme_minimal()
ggsave("DMP_volcano.pdf", width = 8, height = 6)

# 2. 热图（top DMPs）
library(pheatmap)
top_dmps <- head(rownames(sig_dmp), 50)  # 前50个显著DMP
pheatmap(
    beta[top_dmps, ],                    # Beta值子集
    annotation_col = data.frame(         # 样本注释
        Group = myLoad$pd$Group,
        row.names = colnames(beta)
    ),
    scale = "row",                        # 按行标准化
    show_rownames = FALSE,                # 不显示行名（太多）
    main = "Top 50 DMPs"                  # 标题
)

# 3. 基因组特征分布
# 用 ChAMP 做基因组注释分析
champ.GSEA(
    beta = myNorm_combat,                 # Beta值
    DMP = myDMP[[1]],                     # 差异位点
    DMR = dmr_ranges,                     # 差异区域
    arraytype = "EPIC"                    # 芯片类型
)

常见报错与解决¶

报错	原因	解决方案
`cannot open file .idat`	idat 文件路径或文件名格式错误	确认文件名格式为 SampleID_Grn.idat/SampleID_Red.idat
`different array types`	450K 和 EPIC 混用	使用 combineArrays() 或 convertArray() 统一平台
`BMIQ: beta values out of range`	Beta 值含 0 或 1	用 beta[beta == 0] <- 0.001; beta[beta == 1] <- 0.999
`Error in ChAMP: no pd file`	缺少样本信息表	创建 CSV 格式的 pd 文件，至少含 Sample_Name 和 Sentrix_ID
`Too many DMPs (>100k)`	样本量大或阈值太松	收紧 Delta Beta 阈值（如 0.1 或 0.2）
`Memory error`	EPIC 芯片数据量大	使用 read.metharray.exp(extended = FALSE) 减少内存

速查表¶

┌──────────────────────────────────────────────────────────┐
│              DNA 甲基化芯片分析速查                        │
├──────────────────────────────────────────────────────────┤
│ 芯片类型:                                                 │
│   450K: ~485k CpGs    EPIC: ~850k CpGs                   │
│   EPICv2: ~935k CpGs                                     │
│                                                          │
│ 分析流程:                                                 │
│   idat → QC → 标准化 → 批次校正 → DMP/DMR → 富集          │
│                                                          │
│ 核心R包:                                                  │
│   minfi — 底层分析（预处理+标准化+QC）                    │
│   ChAMP — 一站式流水线（推荐新手）                        │
│   DMRcate — DMR检测                                      │
│   limma — 差异分析                                       │
│                                                          │
│ 标准化方法:                                               │
│   BMIQ — 校正Type I/II偏差（ChAMP默认）                  │
│   Noob/ssNoob — 背景校正（多平台整合推荐）                │
│   SeSAMe — 2025年研究认为最优                            │
│                                                          │
│ 差异分析标准:                                              │
│   DMP: FDR < 0.05 + |ΔBeta| > 0.1 (或0.2)               │
│   DMR: ≥3个CpG + 1kb内                                  │
│   统计用M值, 展示用Beta值                                 │
│                                                          │
│ 批次效应: SVD检测 → ComBat校正                           │
└──────────────────────────────────────────────────────────┘

面试高频问题¶

Beta 值和 M 值有什么区别？分析时用哪个？ 答：Beta 值 = methylated/(methylated + unmethylated)，范围 0-1，生物学意义直观但在极端值（0或1附近）时方差不均匀（异方差性）。M 值 = log2(Beta/(1-Beta))，范围 -Inf 到 +Inf，满足正态分布假设，更适合统计检验。最佳实践：统计分析用 M 值，结果展示和解读用 Beta 值。
450K 和 EPIC 芯片数据能合并分析吗？ 答：可以，但需要注意：(1) 只取两个平台共有的探针（约 450K 的 90%+都在 EPIC 中）；(2) 使用 minfi 的 convertArray() 将 EPIC 转为虚拟 450K 格式；(3) 使用 ssNoob（single-sample Noob）标准化，它适合增量处理单个样本。(4) 合并后仍需要 ComBat 校正平台间的系统偏差。
什么是 Type I 和 Type II 探针？为什么需要校正？ 答：Type I 探针用两个珠子分别检测甲基化和未甲基化状态（一个颜色通道），覆盖高 CpG 密度区域；Type II 探针用一个珠子的两个颜色通道，覆盖低 CpG 密度区域。Type II 探针的动态范围更窄，Beta 值分布被压缩。BMIQ 等方法通过将 Type II 分布映射到 Type I 分布来校正这个偏差。