07 DNA 甲基化分析入门

DNA 甲基化（主要是 CpG 位点的 5-甲基胞嘧啶）是最稳定的表观修饰之一。和 ChIP/ATAC 不同，甲基化分析不做 peak calling，而是直接量化每个 CpG 位点的甲基化率（0~100%）。

实验类型

技术	覆盖度	成本	适用
WGBS	全基因组 ~28M CpG	高	全景图、发现新 DMR
RRBS	富集 CpG 岛附近 ~2M CpG	中	启动子甲基化
450K / EPIC 芯片	固定位点 450K~850K	低	大样本量、TCGA 数据

BioF3 这一章聚焦 WGBS/RRBS 的 bisulfite-seq 分析。芯片数据用 minfi 包处理，思路不同。

分析流程

FASTQ → Bismark (bisulfite-aware alignment) → methylation extraction
     → per-CpG methylation table → DMR calling (methylKit / DSS / dmrseq)

Bismark 比对

# 建立 bisulfite 索引（一次性）
bismark_genome_preparation --bowtie2 reference/

# 比对
bismark --genome reference/ -1 sample_R1.fq.gz -2 sample_R2.fq.gz

# 去重
deduplicate_bismark sample_pe.bam

# 提取甲基化信息
bismark_methylation_extractor --paired-end --comprehensive --cytosine_report \
  --genome_folder reference/ sample_pe.deduplicated.bam

输出的 CpG_context 文件每行一个 CpG 位点，包含染色体、位置、甲基化 reads 数、非甲基化 reads 数。

R 里做差异甲基化

library(methylKit)

# 读入 Bismark 的 CpG 报告
file_list <- list("sample1.CpG_report.txt", "sample2.CpG_report.txt",
                  "sample3.CpG_report.txt", "sample4.CpG_report.txt")
obj <- methRead(file_list,
                sample.id = list("ctrl1","ctrl2","trt1","trt2"),
                assembly = "hg38",
                treatment = c(0, 0, 1, 1),
                context = "CpG",
                mincov = 10)

# 合并所有样本的 CpG 位点
meth <- unite(obj, destrand = FALSE)

# 差异甲基化位点
diff <- calculateDiffMeth(meth)
diff_25 <- getMethylDiff(diff, difference = 25, qvalue = 0.01)

# 差异甲基化区域（DMR）
# 用 tileMethylCounts 把基因组分成 1kb 窗口再做差异
tiles <- tileMethylCounts(obj, win.size = 1000, step.size = 1000)
meth_tiles <- unite(tiles)
diff_tiles <- calculateDiffMeth(meth_tiles)

关键概念

• 甲基化率：某个 CpG 位点被甲基化的 reads 占总 reads 的比例
• DMC（Differentially Methylated Cytosine）：单个 CpG 位点的差异
• DMR（Differentially Methylated Region）：连续多个 CpG 位点一起变化的区域
• 覆盖度过滤：覆盖度 < 10x 的位点噪声太大，通常过滤掉

和表达的关系

启动子区域的高甲基化通常和基因沉默相关；基因体内的甲基化和活跃转录正相关。把 DMR 和 RNA-seq 的差异基因做交叉，能找到"甲基化变化驱动表达变化"的候选基因。

但要注意：甲基化和表达的关系不是简单一对一。

区域	甲基化 ↑	一般预期
启动子（CpG island）	表达 ↓	经典抑制
启动子（非 CpG island）	表达关系不明确	不要硬解读
基因 body	表达 ↑（活跃转录基因）	反直觉
增强子	表达 ↓	类似启动子
重复元件	沉默 retrotransposon	调控功能

实务上：先把 DMR 注释到 promoter / gene body / enhancer / intergenic 各自分类，再分别和 RNA-seq 对比。整体混在一起看会得到无意义的"弱相关"。

常见坑

坑 1：没做 bisulfite 转化率 QC

WGBS 实验的 bisulfite 转化率应该 > 99%。未甲基化的 lambda phage spike-in 是金标准，转化率低意味着大量未甲基化 C 被当成 5mC，所有结果系统偏移。Bismark 输出会自动算这个值，看一眼。

坑 2：覆盖度不足直接做差异

WGBS 单 CpG 位点的可信差异需要每位点至少 10× 覆盖度。一份 30× WGBS 看上去深，但 CpG 只占基因组 1%，单位点覆盖度会被稀释。先 methRead(..., mincov = 10) 过滤。

坑 3：用 DMC 当 DMR 报告

单个 CpG（DMC）噪声大，真正可靠的是连续多个 CpG 一起变化的 DMR。methylKit tileMethylCounts 或 dmrseq 找区域，而不是把 DMC 直接报告。

坑 4：跨平台数据强行合并

WGBS / RRBS / 450K / EPIC 的位点集合不重合（450K 只覆盖 ~2% CpG）。跨平台 meta 分析需要先取交集 + 再 normalize，直接 cbind 是错的。

坑 5：differential = 25% 的阈值死搬

methylKit 默认 difference = 25 意为甲基化率差 ≥ 25%。对小效应（药物处理、衰老）来说 25% 太严，差不到 5% 但显著的位点经常有意义。从分布看实际差异范围再选阈值，不要照搬教程。

下一步

接着深入：

• 08 表观组与转录组的整合 — DMR 最有价值的解读是和差异表达基因的关系
• 多组学整合 module01 — 甲基化是整合分析里最常用的 omics 之一

横向延伸：

• Bismark 完整文档 — bisulfite 比对的细节
• methylKit 官方教程

参考资源

• Bismark 手册
• methylKit 教程
• DSS 差异甲基化
• dmrseq
• minfi（芯片数据）