差异甲基化区域DMR检测¶

DMR检测 = 找"甲基化变化热点区域"

思路：
  单个CpG位点的差异（DML）可能是噪音
  多个连续CpG都变化 → 说明这个区域的甲基化真的变了 → DMR

类比：
  DML = 一个路口堵了（可能只是偶然）
  DMR = 一整条路都堵了（这条路肯定出问题了）

DMR的生物学意义：
  启动子区域超甲基化 → 基因沉默（常见于肿瘤抑癌基因）
  启动子区域低甲基化 → 基因激活（常见于癌基因）
  基因体甲基化 → 可能促进转录延伸
  增强子甲基化变化 → 调控远端基因表达

#!/usr/bin/env Rscript
# 使用DSS包检测DMR

# ========== 安装 ==========
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DSS")  # 安装DSS包

library(DSS)  # 加载DSS
library(bsseq)  # 加载bsseq（DSS依赖）

# ========== 准备数据 ==========
# 读取Bismark的CpG报告
# 格式：chr, pos, strand, count_meth, count_unmeth, context, trinucleotide
read_bismark_cov <- function(file, sample_name) {
    df <- read.delim(file, header = FALSE, 
                     col.names = c("chr", "start", "end", "meth_pct", "count_meth", "count_unmeth"))
    df$N <- df$count_meth + df$count_unmeth  # 总覆盖度
    df$X <- df$count_meth  # 甲基化reads数
    return(df[, c("chr", "start", "N", "X")])  # 返回DSS需要的4列
}

# 读取样品数据
ctrl1 <- read_bismark_cov("control_1.bismark.cov", "Control_1")
ctrl2 <- read_bismark_cov("control_2.bismark.cov", "Control_2")
ctrl3 <- read_bismark_cov("control_3.bismark.cov", "Control_3")
treat1 <- read_bismark_cov("treatment_1.bismark.cov", "Treatment_1")
treat2 <- read_bismark_cov("treatment_2.bismark.cov", "Treatment_2")
treat3 <- read_bismark_cov("treatment_3.bismark.cov", "Treatment_3")

# 创建BSseq对象
BSobj <- makeBSseqData(
    list(ctrl1, ctrl2, ctrl3, treat1, treat2, treat3),  # 所有样品数据
    c("C1", "C2", "C3", "T1", "T2", "T3")  # 样品名
)
cat("BSseq对象创建完成\n")
cat("CpG位点数:", nrow(BSobj), "\n")

# ========== DML检测（单位点水平） ==========
dml_test <- DMLtest(
    BSobj,
    group1 = c("C1", "C2", "C3"),  # 对照组样品
    group2 = c("T1", "T2", "T3"),  # 处理组样品
    smoothing = TRUE  # 平滑处理（推荐）
)
cat("DML检测完成\n")

# 提取显著DML
dmls <- callDML(
    dml_test,
    p.threshold = 0.001  # p值阈值
)
cat("显著DML数:", nrow(dmls), "\n")

# ========== DMR检测（区域水平） ==========
dmrs <- callDMR(
    dml_test,
    p.threshold = 0.01,  # 单位点p值阈值
    minlen = 50,  # 最小DMR长度(bp)
    minCG = 3,  # 最少包含3个CpG
    dis.merge = 50,  # 50bp内的相邻DMR合并
    pct.sig = 0.5  # 至少50%的CpG位点显著
)

cat("\n===== DMR检测结果 =====\n")
cat("检测到的DMR数:", nrow(dmrs), "\n")
cat("超甲基化DMR:", sum(dmrs$areaStat > 0), "\n")  # 处理组甲基化升高
cat("低甲基化DMR:", sum(dmrs$areaStat < 0), "\n")  # 处理组甲基化降低

# 保存结果
write.csv(dmrs, "DSS_DMR_results.csv", row.names = FALSE)
cat("结果已保存到 DSS_DMR_results.csv\n")

# ========== 可视化 ==========
# 展示单个DMR区域的甲基化模式
if (nrow(dmrs) > 0) {
    showOneDMR(dml_test, dmrs[1, ])  # 显示第一个DMR
}

#!/usr/bin/env Rscript
# 使用methylKit进行DMR分析

BiocManager::install("methylKit")
library(methylKit)  # 加载methylKit

# ========== 读取数据 ==========
# methylKit可以直接读取Bismark的覆盖度文件
file_list <- list(
    "control_1.bismark.cov",
    "control_2.bismark.cov",
    "treatment_1.bismark.cov",
    "treatment_2.bismark.cov"
)

# 读取并创建methylRawList对象
myobj <- methRead(
    file_list,
    sample.id = list("C1", "C2", "T1", "T2"),  # 样品ID
    assembly = "hg38",  # 参考基因组版本
    treatment = c(0, 0, 1, 1),  # 0=对照组, 1=处理组
    context = "CpG",  # 甲基化上下文
    pipeline = "bismarkCoverage"  # 输入格式
)

# ========== 质量控制 ==========
# 每个样品的甲基化分布
getMethylationStats(myobj[[1]], plot = TRUE, both.strands = FALSE)  # 甲基化分布图
getCoverageStats(myobj[[1]], plot = TRUE, both.strands = FALSE)  # 覆盖度分布图

# ========== 过滤 ==========
# 过滤覆盖度极端的位点
filtered <- filterByCoverage(
    myobj,
    lo.count = 5,  # 最低覆盖度5x
    lo.perc = NULL,
    hi.count = NULL,
    hi.perc = 99.9  # 去掉最高0.1%（可能是PCR偏好）
)

# 标准化覆盖度
normalized <- normalizeCoverage(filtered, method = "median")  # 中位数标准化

# ========== 合并样品 ==========
# 只保留所有样品都覆盖的CpG位点
meth <- unite(normalized, destrand = FALSE)  # 合并数据
cat("共同覆盖的CpG位点:", nrow(meth), "\n")

# ========== 差异分析 ==========
# 计算差异甲基化
diff <- calculateDiffMeth(
    meth,
    overdispersion = "MN",  # 使用beta-二项模型
    adjust = "BH"  # BH多重检验校正
)

# 提取显著差异位点
# 差异>25%且q值<0.01
hyper <- getMethylDiff(diff, difference = 25, qvalue = 0.01, type = "hyper")  # 超甲基化
hypo <- getMethylDiff(diff, difference = 25, qvalue = 0.01, type = "hypo")  # 低甲基化
all_diff <- getMethylDiff(diff, difference = 25, qvalue = 0.01, type = "all")  # 全部

cat("\n===== 差异甲基化结果 =====\n")
cat("超甲基化位点:", nrow(hyper), "\n")
cat("低甲基化位点:", nrow(hypo), "\n")

# ========== 基因注释 ==========
BiocManager::install("genomation")
library(genomation)

# 读取基因注释
gene_bed <- readTranscriptFeatures("hg38_refseq.bed")  # 读取基因注释BED文件

# 注释差异甲基化位点
diffAnn <- annotateWithGeneParts(as(all_diff, "GRanges"), gene_bed)
getTargetAnnotationStats(diffAnn, percentage = TRUE, precedence = TRUE)  # 统计区域分布
plotTargetAnnotation(diffAnn, precedence = TRUE, main = "DMR Genomic Annotation")  # 画饼图

#!/usr/bin/env python3
"""DMR结果可视化"""

import pandas as pd  # 数据处理
import matplotlib.pyplot as plt  # 画图
import numpy as np  # 数值计算

# ========== 读取DMR结果 ==========
dmrs = pd.read_csv("DSS_DMR_results.csv")

# ========== 1. DMR长度分布 ==========
plt.figure(figsize=(10, 6))
plt.hist(dmrs["length"], bins=50, color="steelblue", edgecolor="none", alpha=0.8)
plt.xlabel("DMR Length (bp)", fontsize=14)
plt.ylabel("Count", fontsize=14)
plt.title("DMR Length Distribution", fontsize=16)
plt.tight_layout()
plt.savefig("dmr_length_dist.png", dpi=300)
plt.close()

# ========== 2. 超/低甲基化DMR统计 ==========
hyper = dmrs[dmrs["areaStat"] > 0]  # 超甲基化
hypo = dmrs[dmrs["areaStat"] < 0]  # 低甲基化

fig, ax = plt.subplots(figsize=(6, 6))
ax.bar(["Hyper", "Hypo"], [len(hyper), len(hypo)], color=["red", "blue"], alpha=0.7)
ax.set_ylabel("Number of DMRs", fontsize=14)
ax.set_title("Hyper vs Hypo DMRs", fontsize=16)
plt.tight_layout()
plt.savefig("dmr_hyper_hypo.png", dpi=300)
plt.close()

# ========== 3. 染色体分布 ==========
chr_counts = dmrs["chr"].value_counts()
chr_order = [f"chr{i}" for i in range(1, 23)] + ["chrX", "chrY"]
chr_counts = chr_counts.reindex(chr_order).dropna()

plt.figure(figsize=(14, 6))
plt.bar(range(len(chr_counts)), chr_counts.values, color="coral")
plt.xticks(range(len(chr_counts)), chr_counts.index, rotation=45)
plt.xlabel("Chromosome", fontsize=14)
plt.ylabel("Number of DMRs", fontsize=14)
plt.title("DMR Distribution Across Chromosomes", fontsize=16)
plt.tight_layout()
plt.savefig("dmr_chromosome_dist.png", dpi=300)
plt.close()

print("所有DMR可视化完成")

========================================
差异甲基化区域DMR检测 速查表
========================================

【三大DMR检测工具】
DSS                  → 基于beta-二项分布，推荐
methylKit            → 全功能甲基化分析（检测+注释+可视化）
dmrseq               → 专为WGBS设计

【DSS核心流程】
1. makeBSseqData()   → 创建BSseq对象
2. DMLtest()         → DML检测（单位点水平）
3. callDML()         → 提取显著DML
4. callDMR()         → 检测DMR（区域水平）

【callDMR关键参数】
p.threshold          → 单位点p值阈值（0.01推荐）
minlen               → 最小DMR长度（50bp推荐）
minCG                → 最少CpG数（3-5个推荐）
dis.merge            → 合并距离（50-100bp）
pct.sig              → 显著CpG比例（0.5推荐）

【DMR判断标准】
差异 > 25%           → 一般标准
差异 > 20%           → 宽松标准
差异 > 10%           → 探索性分析
q值 < 0.01           → 显著性阈值

【生物学意义】
启动子超甲基化       → 基因沉默（抑癌基因常见）
启动子低甲基化       → 基因激活
基因体甲基化         → 可能促进转录
CpG岛shore甲基化     → 组织特异性

【面试考点】
Q: DML和DMR的区别？
A: DML=单个CpG位点差异，DMR=连续多个CpG组成的区域差异

Q: 为什么需要DMR而不只看DML？
A: 单位点容易受噪音影响，区域水平更可靠更有生物学意义

Q: 常用的DMR检测方法？
A: DSS（贝叶斯）、methylKit（逻辑回归）、dmrseq（回归+平滑）
========================================