Bio Methylation Dmr Detection

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

Goal: Map differentially methylated regions to overlapping genes, promoters, and CpG islands for biological interpretation.

Approach: Build a genome annotation set with annotatr, convert DMRs to GRanges, and intersect with genomic features to classify each DMR by functional context.

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

Key Parameters

methylKit tileMethylCounts

bsseq dmrFinder

Related Skills

• methylkit-analysis - Single CpG analysis
• differential-cpg-testing - Per-CpG testing before region-level analysis
• methylation-calling - Generate input files
• pathway-analysis/go-enrichment - Functional annotation of DMR genes
• differential-expression/deseq2-basics - Compare with expression changes