name: bio-methylation-dmr-detection
description: Differentially methylated region (DMR) detection using methylKit tiles, bsseq BSmooth, and DMRcate. Use when identifying contiguous genomic regions with methylation differences between experimental conditions or cell types.
tool_type: r
primary_tool: methylKit
measurable_outcome: Execute skill workflow successfully with valid output within 15 minutes.
allowed-tools:
- read_file
- run_shell_command
DMR Detection
methylKit Tile-Based DMRs
library(methylKit)
# Read and process data
meth_obj <- methRead(location = file_list, sample.id = sample_ids, treatment = treatment,
assembly = 'hg38', pipeline = 'bismarkCoverage')
meth_filt <- filterByCoverage(meth_obj, lo.count = 10, hi.perc = 99.9)
# Create tiles (windows)
tiles <- tileMethylCounts(meth_filt, win.size = 1000, step.size = 1000, cov.bases = 3)
tiles_united <- unite(tiles, destrand = TRUE)
# Differential methylation on tiles
diff_tiles <- calculateDiffMeth(tiles_united, overdispersion = 'MN', mc.cores = 4)
# Get significant DMRs
dmrs <- getMethylDiff(diff_tiles, difference = 25, qvalue = 0.01)
dmrs_hyper <- getMethylDiff(diff_tiles, difference = 25, qvalue = 0.01, type = 'hyper')
dmrs_hypo <- getMethylDiff(diff_tiles, difference = 25, qvalue = 0.01, type = 'hypo')
bsseq BSmooth DMRs
library(bsseq)
# Read Bismark cytosine reports
bs <- read.bismark(files = c('sample1.CpG_report.txt.gz', 'sample2.CpG_report.txt.gz'),
sampleNames = c('ctrl', 'treat'),
rmZeroCov = TRUE,
strandCollapse = TRUE)
# Smooth methylation data
bs_smooth <- BSmooth(bs, mc.cores = 4, verbose = TRUE)
# Filter by coverage
bs_cov <- getCoverage(bs_smooth)
keep <- which(rowSums(bs_cov >= 2) == ncol(bs_cov))
bs_filt <- bs_smooth[keep, ]
# Find DMRs with BSmooth
dmrs_bsseq <- dmrFinder(bs_filt, cutoff = c(-0.1, 0.1), stat = 'tstat.corrected')
DMRcate Method
library(DMRcate)
library(minfi)
# From methylation matrix (beta values)
# Rows = CpGs, columns = samples
design <- model.matrix(~ treatment)
# Run DMRcate
myannotation <- cpg.annotate('array', meth_matrix, what = 'Beta', arraytype = 'EPIC',
design = design, coef = 2)
dmr_results <- dmrcate(myannotation, lambda = 1000, C = 2)
dmr_ranges <- extractRanges(dmr_results)
Annotate DMRs with Genes
library(annotatr)
# Build annotations
annots <- build_annotations(genome = 'hg38', annotations = c(
'hg38_basicgenes',
'hg38_genes_promoters',
'hg38_cpg_islands'
))
# Convert DMRs to GRanges
dmr_gr <- as(dmrs, 'GRanges')
# Annotate
dmr_annotated <- annotate_regions(regions = dmr_gr, annotations = annots, ignore.strand = TRUE)
dmr_df <- data.frame(dmr_annotated)
Annotate with genomation
library(genomation)
# Read gene annotations
gene_obj <- readTranscriptFeatures('genes.bed12')
# Annotate DMRs
dmr_gr <- as(dmrs, 'GRanges')
annot_result <- annotateWithGeneParts(dmr_gr, gene_obj)
# Get promoter/exon/intron breakdown
getTargetAnnotationStats(annot_result, percentage = TRUE, precedence = TRUE)
Visualize DMR
library(Gviz)
# Create track for a DMR
chr <- 'chr1'
start <- 1000000
end <- 1010000
# Methylation data track
meth_track <- DataTrack(
range = bs_smooth,
genome = 'hg38',
name = 'Methylation',
type = 'smooth'
)
# Gene annotation track
gene_track <- GeneRegionTrack(TxDb.Hsapiens.UCSC.hg38.knownGene, genome = 'hg38', name = 'Genes')
# Plot
plotTracks(list(meth_track, gene_track), from = start, to = end, chromosome = chr)
Merge Adjacent DMRs
library(GenomicRanges)
dmr_gr <- as(dmrs, 'GRanges')
# Merge DMRs within 500bp
dmr_merged <- reduce(dmr_gr, min.gapwidth = 500)
Export DMRs
# To BED
library(rtracklayer)
export(dmr_gr, 'dmrs.bed', format = 'BED')
# To CSV
dmr_df <- getData(dmrs)
write.csv(dmr_df, 'dmrs.csv', row.names = FALSE)
# To GFF
export(dmr_gr, 'dmrs.gff3', format = 'GFF3')
DMR Comparison Across Methods
| Method |
Package |
Approach |
Best For |
| Tiles |
methylKit |
Fixed windows |
Quick analysis |
| BSmooth |
bsseq |
Smoothing |
WGBS data |
| DMRcate |
DMRcate |
Kernel smoothing |
Array data |
| DSS |
DSS |
Bayesian |
Complex designs |
Key Parameters
methylKit tileMethylCounts
| Parameter |
Default |
Description |
| win.size |
1000 |
Window size (bp) |
| step.size |
1000 |
Step size (bp) |
| cov.bases |
0 |
Min CpGs per tile |
bsseq dmrFinder
| Parameter |
Description |
| cutoff |
Methylation difference threshold |
| stat |
Statistic to use |
| maxGap |
Max gap between CpGs |
Related Skills
- methylkit-analysis - Single CpG analysis
- methylation-calling - Generate input files
- pathway-analysis/go-enrichment - Functional annotation of DMR genes
- differential-expression/deseq2-basics - Compare with expression changes