--- name: dna-methylation-differential-analysis description: Use when you have merged methylation call data across multiple biological replicates (samples per group ≥2) with base-pair-level coverage information, and you need to identify loci where methylation levels differ significantly between treatment groups. license: CC-BY-4.0 metadata: edam_operation: http://edamontology.org/operation_3233 edam_topics: - http://edamontology.org/topic_3295 - http://edamontology.org/topic_0654 tools: - R - genomation - GenomicFeatures - methylKit - Bismark - MethylDackel derived_from: - doi: 10.1186/gb-2012-13-10-r87 title: methylkit evidence_spans: - packageVersion('methylKit') - We can annotate our differentially methylated regions/bases based on gene annotation using genomation package - This annotation operation will tell us what percentage of our differentially methylated regions are on promoters/introns/exons/intergenic region. In this case we read annotation from a BED file, claims: [] provenance: collection: https://w3id.org/holobiomicslab/asb-skill/collection/epigenomics/v1 assembled_by: scripts/collect_metabolomics_collection.py sources: - build: coll_methylkit doi: 10.1186/gb-2012-13-10-r87 title: methylkit dedup_kept_from: coll_methylkit schema_version: 0.2.0 --- # dna-methylation-differential-analysis ## Summary Identify and quantify differentially methylated bases or regions between sample groups using statistical tests that account for methylation heterogeneity and overdispersion. This skill applies Fisher's exact test or logistic regression to bisulfite sequencing methylation calls, filtered by q-value and percent methylation difference thresholds, to distinguish hyper-methylated from hypo-methylated loci. ## When to use Apply this skill when you have merged methylation call data across multiple biological replicates (samples per group ≥2) with base-pair-level coverage information, and you need to identify loci where methylation levels differ significantly between treatment groups. This is the core comparative analysis step in bisulfite sequencing workflows after quality filtering and sample merging. ## When NOT to use - Input data lacks technical replicates or biological replication structure—statistical tests require ≥2 samples per group to estimate variance. - Bases have not been pre-filtered for minimum coverage (typically 10X default in methRead())—low-coverage bases produce unreliable methylation percentages. - You seek to analyze regional methylation (e.g., DMRs across promoters) rather than base-resolution differential calls—use regional or tiling window methods instead. ## Inputs - methylBase object (merged sample data from unite() with base-level coverage and methylation percentages) - methylRawList objects (optional, for overdispersion estimation; generated by methRead() from bisulfite alignment outputs) ## Outputs - methylDiff object containing all bases with calculated test statistics and q-values - hyper-methylated bases subset (type='hyper' from getMethylDiff()) - hypo-methylated bases subset (type='hypo' from getMethylDiff()) - numerical counts and percentages of significant bases by direction ## How to apply Execute the calculateDiffMeth() function on a methylBase object (created by merging samples with unite()) to compute differential methylation statistics. The function automatically selects Fisher's exact test for small sample sizes or logistic regression for larger cohorts. For studies showing overdispersion (variance exceeding binomial expectations typical in methylation data), apply the overdispersion='MN' parameter, which calculates a scaling factor φ = X²/(N-P) to adjust variance as φ·n_i·π̂_i·(1-π̂_i) and switches to an F-test, producing more stringent (higher) q-values. Extract differentially methylated bases using getMethylDiff() with dual filtering: q-value threshold (typically < 0.01) and percent methylation difference cutoff (typically > 25%). Separate hyper-methylated (high methylation in treatment) from hypo-methylated (low methylation in treatment) bases using the type parameter to assess directional changes. ## Related tools - **methylKit** (Core R package providing calculateDiffMeth(), getMethylDiff(), methRead(), unite(), and filterByCoverage() functions for differential methylation analysis and statistical testing) — https://github.com/al2na/methylKit - **Bismark** (Bisulfite sequencing aligner and methylation caller that produces the CpG call files (input to methRead()) used to construct methylRawList objects) — https://github.com/FelixKrueger/Bismark - **MethylDackel** (Alternative methylation extractor from BAM/CRAM alignments; produces bedGraph-format methylation calls compatible with methylKit workflows) — https://github.com/dpryan79/MethylDackel - **R** (Runtime environment and statistical framework for methylKit differential analysis functions) ## Examples ``` calculateDiffMeth(methylBase_obj, overdispersion='MN', test='Chisq'); diff_meth <- getMethylDiff(methylDiff_obj, difference=25, qvalue=0.01); hyper <- getMethylDiff(methylDiff_obj, difference=25, qvalue=0.01, type='hyper'); hypo <- getMethylDiff(methylDiff_obj, difference=25, qvalue=0.01, type='hypo') ``` ## Evaluation signals - Hyper-methylated and hypo-methylated base counts match vignette-reported values for the same example dataset and filtering thresholds (q < 0.01, Δmeth > 25%). - Q-value distribution from overdispersion-corrected run (overdispersion='MN') shows higher mean and median q-values than uncorrected run, confirming more stringent correction. - methylDiff object contains non-null test statistics (stat), p-values, and q-values for all rows; no NaN or Inf values for bases meeting coverage criteria. - Extracted hyper and hypo subsets are mutually exclusive and collectively account for all bases passing the q-value and methylation difference thresholds. - Sample size (number of replicates per group) correctly triggers selection of statistical test: Fisher's exact for n<30 per group, logistic regression for n≥30. ## Limitations - Statistical tests assume independence of bases and do not account for spatial autocorrelation across the genome—neighboring bases are treated as independent observations. - Overdispersion correction ('MN') requires sufficient within-group replicates to estimate variance reliably; with very small sample sizes, the scaling factor may be unstable. - Methylation percentages from bisulfite sequencing reflect sample heterogeneity and sequencing errors; binary methylation calls (fully methylated vs. unmethylated) are not produced, limiting interpretation of single-molecule states. - The skill requires pre-filtering of bases by coverage (default ≥10 reads) to ensure reliable methylation calls; bases below this threshold are discarded before statistical testing and may lead to loss of weakly covered regions. - Q-value filtering at 0.01 and 25% methylation difference cutoffs are data-dependent thresholds; very stringent cutoffs may miss subtle but biologically relevant changes, while lenient cutoffs increase false discovery rate. ## Evidence - [intro] The calculateDiffMeth() function is the main function to calculate differential methylation. Depending on the sample size per each set it will either use Fisher's exact or logistic regression: "The calculateDiffMeth() function is the main function to calculate differential methylation. Depending on the sample size per each set it will either use Fisher's exact or logistic regression" - [intro] After q-value calculation, we can select the differentially methylated regions/bases based on q-value and percent methylation difference cutoffs: "After q-value calculation, we can select the differentially methylated regions/bases based on q-value and percent methylation difference cutoffs" - [intro] Using getMethylDiff() with q-value < 0.01 and 25% methylation difference cutoffs yields separate hyper-methylated and hypo-methylated base objects, extractable by specifying type parameter: "Using getMethylDiff() with q-value < 0.01 and 25% methylation difference cutoffs on example methylKit data yields separate hyper-methylated and hypo-methylated base objects" - [intro] The calculateDiffMeth() function with overdispersion='MN' parameter calculates a scaling parameter φ to adjust variance, making statistical tests more stringent and automatically switching to F-test: "The calculateDiffMeth() function with overdispersion='MN' parameter calculates a scaling parameter φ = X²/(N-P) to adjust variance as φ·n_i·π̂_i·(1-π̂_i), which makes statistical tests more stringent" - [intro] Usually, the methylation status of a base determined by a high-throughput bisulfite sequencing will not be a binary score, but it will be a percentage: "Usually, the methylation status of a base determined by a high-throughput bisulfite sequencing will not be a binary score, but it will be a percentage" - [intro] By default methRead requires a minimum coverage of 10 reads per base to ensure good quality of the data and a high confidence methylation percentage: "By default methRead requires a minimum coverage of 10 reads per base to ensure good quality of the data and a high confidence methylation percentage" - [intro] In order to do further analysis, we will need to get the bases covered in all samples. The following function will merge all samples to one object for base-pair locations that are covered in all: "In order to do further analysis, we will need to get the bases covered in all samples. The following function will merge all samples to one object for base-pair locations" - [readme] methylKit is an R package for DNA methylation analysis and annotation from high-throughput bisulfite sequencing: "methylKit is an R package for DNA methylation analysis and annotation from high-throughput bisulfite sequencing"