Read The epigenetic landscape of transgenerational acclimation to ocean warming to get an idea for potential analysis methods. The authors used different functions in methylKit to obtain DMRs:
"Briefly, the ‘methRead’ function of methylKit reads the mapping results with 10 reads per cytosine as a minimum coverage threshold. High coverage bases (99.9%) were filtered to exclude potential PCR bias and then normalized using ‘filterByCoverage’ and ‘normalizeCoverage’ functions, respectively. Genomic regions were categorized as CpG island, CpG shore, promoter, 5′ untranslated region (UTR), exon, introns, 3′ UTR and repeats. Methylated or unmethylated cytosines in each genomic region were summed for each sample by the ‘regionCounts’ function of methylKit. The P values of methylation differences for each region between two samples were calculated using a chi-squared test in the ‘calculateDiffMeth’ function."
They also generated heatmaps with DMR data, which would be useful in my case as well.
In BEDtools, they used the closest function to pair DMRs and genes:
"The closest gene to a DMR on the same scaffold was identified using ‘closest’ from BEDTools v. 2.2339. This resulted in 1,563 genes from 2,078 CpG DMRs, while 115 DMRs were on scaffolds without annotated genes (Supplementary Table 4)."
Because they paired gene expression data with epigenetic data, they did not do any gene enrichment. I’ll need to refer to Emma’s geoduck paper for those methods.