| name | bio-methylation-methylkit |
| description | Imports Bismark coverage or cytosine-report files into the methylKit object model, then runs the import-to-results spine - filterByCoverage, normalizeCoverage, unite/destrand, calculateDiffMeth, getMethylDiff - for both per-CpG (DMC) and fixed-tile (DMR) differential methylation, plus tileMethylCounts, PCA/correlation/clustering QC, and assocComp/removeComp batch handling. Covers the silent default traps that shape the false-positive rate: overdispersion='none' does no correction while 'MN' forces the F-test (ignoring test='Chisq'), adjust defaults to SLIM not BH, getMethylDiff defaults difference=25/qvalue=0.01, cov.bases=0 admits single-CpG tiles, and pool destroys biological replication. Use when importing bisulfite count tables, filtering/normalizing/uniting methylation samples, running methylKit differential testing, or QC-ing methylomes. For per-site test-choice (count vs continuous) see differential-cpg-testing; for selection-aware region FDR (dmrseq/DSS) see dmr-detection. |
Version Compatibility
Reference examples tested with: methylKit 1.28+, GenomicRanges 1.54+.
Before using code patterns, verify installed versions match. If versions differ:
- R:
packageVersion('<pkg>') then ?function_name to verify parameters
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
The assembly= string (e.g. hg38) is metadata only - methylKit never checks it. The genome build is real elsewhere: coordinates must match the alignment genome, and annotation packages (TxDb.Hsapiens.UCSC.hg38.knownGene, annotatr build_annotations(genome='hg38')) are genome-build-specific. methylKit's overdispersion, test, and adjust defaults have shifted across Bioconductor releases - run ?calculateDiffMeth on the installed build before trusting any default.
methylKit Analysis
"Analyze methylation across my samples" -> Import per-cytosine counts into a methylRawList, then filter, normalize, unite, and test - because each of those steps is a modeling decision that sets which CpGs survive and how many false positives the test emits, not boilerplate.
- R:
methRead(pipeline='bismarkCoverage') -> filterByCoverage() -> normalizeCoverage() -> unite(destrand=) -> calculateDiffMeth(overdispersion='MN') -> getMethylDiff()
Scope: the methylKit OBJECT MODEL and the short-read bisulfite import-to-results workflow, for BOTH per-CpG (DMC) and fixed-tile (DMR) results. Which per-site test to use (count vs continuous, beta vs M) -> differential-cpg-testing. Selection-aware region callers (dmrseq/DSS/metilene) and region FDR -> dmr-detection. Long-read MM/ML modBAM input -> long-read-sequencing/nanopore-methylation (its counts pipe back into this object model).
The Single Most Important Modern Insight -- The Object Model Is the Analysis
Coverage filtering, normalization, destranding, and the overdispersion model are not setup before the "real" test - they ARE the test. Each silently changes which CpGs exist and what the p-value means, and methylKit's defaults are tuned for nothing in particular. Three corollaries every misuse violates:
- The defaults offer no protection.
calculateDiffMeth defaults to overdispersion='none' - a plain logistic LRT that assumes binomial-only variance and over-calls under biological replication. Two healthy replicates differ at a CpG far more than coin-flip sampling predicts; only overdispersion='MN' adds the between-replicate (beta-binomial) layer. The paper discusses overdispersion; the function does not apply it unless told.
- The knobs interact and several are silent.
overdispersion='MN' automatically switches to the F-test, so a passed test='Chisq' is ignored with no warning. adjust defaults to SLIM (methylKit's own q-method), not BH, so counts are not comparable to a DSS/limma BH analysis. tileMethylCounts(cov.bases=0) lets a one-CpG window become a "region." None of these throw an error; the result just quietly changes.
- Coverage is the substrate, not the answer. A single CpG's methylation percentage is a count ratio; below ~10x it is a coin flip. Filtering the low tail (noise) and the high tail (PCR/repeat artifacts) before testing decides the result more than the test does.
Organize the workflow around defending these, not around calling functions in order.
The Import-to-Results Spine
Run these in order. Skipping or reordering them changes the result silently.
1. Import: methRead with the pipeline matching the input
Goal: Load per-cytosine counts into a methylRawList, choosing the parser that matches the Bismark output format.
Approach: pipeline='bismarkCoverage' reads .cov/.cov.gz (chr/start/end/%meth/numC/numT - NO strand, so destranding is limited); pipeline='bismarkCytosineReport' reads the CX/CpG report (carries strand + context, enables proper destranding). treatment is an integer vector (0/1, or 0/1/2 for multi-group). context='CpG' only - never destrand CHG/CHH downstream.
library(methylKit)
file_list <- list('ctrl1.cov.gz', 'ctrl2.cov.gz', 'treat1.cov.gz', 'treat2.cov.gz')
sample_ids <- list('ctrl_1', 'ctrl_2', 'treat_1', 'treat_2')
meth_obj <- methRead(file_list, sample.id=sample_ids, treatment=c(0,0,1,1),
assembly='hg38', context='CpG', pipeline='bismarkCoverage')
2. Filter and normalize BEFORE uniting (and before tiling)
Goal: Drop unreliable and artifactual CpGs per sample, then remove library-size-driven coverage differences so a deeper sample does not look more "confident."
Approach: filterByCoverage(lo.count, hi.perc) removes the noisy low tail and the artifactual high tail; normalizeCoverage scales coverage between samples. Both are per-sample and must precede unite and tileMethylCounts.
meth_filt <- filterByCoverage(meth_obj, lo.count=10, lo.perc=NULL, hi.count=NULL, hi.perc=99.9)
meth_norm <- normalizeCoverage(meth_filt, method='median')
3. Unite: destrand only for CpG, only with strand info
Goal: Build the per-base table of CpGs covered across samples for testing.
Approach: unite keeps CpGs covered in ALL samples; min.per.group=2L relaxes that to >=2 per group (keeps more sites, allows missingness). destrand=TRUE merges the + and - strand counts of a CpG dyad - valid ONLY for symmetric CpG context AND only meaningful when strand is present (cytosine report). On .cov (bismarkCoverage, no strand) destranding is limited; on CHG/CHH it is wrong.
meth_united <- unite(meth_norm, destrand=TRUE)
meth_united <- unite(meth_norm, min.per.group=2L)
4. QC the united object before testing
Goal: Confirm samples cluster by biology, not by batch, before believing any DMC.
Approach: Run correlation, PCA, and clustering on the united (% methylation) object. A control clustering with the treated group, or PC1 tracking sequencing batch, means the contrast is confounded.
getCorrelation(meth_united, plot=TRUE)
PCASamples(meth_united)
clusterSamples(meth_united, dist='correlation', method='ward.D', plot=TRUE)
5. Test: calculateDiffMeth with overdispersion correction
Goal: Test each CpG for a group difference using a model that accounts for between-replicate overdispersion.
Approach: With replicates, set overdispersion='MN', which automatically uses the F-test (the test= argument is then ignored - passing test='Chisq' alongside MN does not produce a chi-square test). Set adjust='BH' if the q-values must be comparable to other tools; the default SLIM is methylKit-specific. getMethylDiff filters by effect size AND q.
diff_meth <- calculateDiffMeth(meth_united, overdispersion='MN', adjust='BH', mc.cores=4)
dmcs <- getMethylDiff(diff_meth, difference=25, qvalue=0.01)
dmcs_hyper <- getMethylDiff(diff_meth, difference=25, qvalue=0.01, type='hyper')
Tile-Based Regions (a fast screen, not selection-corrected inference)
Goal: Aggregate CpGs into fixed windows for a quick region-level scan.
Approach: Tile AFTER filter/normalize, raise cov.bases so a window needs real CpG support, then flow through the same unite -> calculateDiffMeth -> getMethylDiff path. The same getMethylDiff returns DMCs on a per-base object and (window) DMRs on a tiled object.
tiles <- tileMethylCounts(meth_norm, win.size=1000, step.size=1000, cov.bases=3)
tiles_united <- unite(tiles, destrand=FALSE)
diff_tiles <- calculateDiffMeth(tiles_united, overdispersion='MN', adjust='BH', mc.cores=4)
dmrs <- getMethylDiff(diff_tiles, difference=25, qvalue=0.01)
The per-tile q is a per-test SLIM/BH value: it does NOT model correlation between tiles and is NOT corrected for the region-selection step. Fixed windows also split or merge true DMRs at arbitrary boundaries. Treat methylKit tiles as a defensible screen; for rigorous region FDR (a permutation null that survives region selection) go to dmr-detection (dmrseq).
Batch, Multi-Group, and the Single-Factor Limit
methylKit's calculateDiffMeth is a single-factor 2-group test. For known batch, remove the associated principal components before testing; for >2 groups, subset to pairwise contrasts. Complex designs (covariates, multi-factor) exceed what methylKit models - move to dmr-detection (DSS multiFactor / dmrseq covariates) or a continuous limma-on-M path (differential-cpg-testing).
sample_anno <- data.frame(batch=c('a','a','b','b'))
as_comp <- assocComp(meth_united, sample_anno)
meth_corrected <- removeComp(meth_united, comp=1)
meth_AB <- reorganize(meth_united, sample.ids=c('ctrl_1','ctrl_2','treat_1','treat_2'),
treatment=c(0,0,1,1))
pool(meth_united, sample.ids=...) sums replicate counts into one pseudo-sample per group. This DESTROYS biological replication - the test then has no within-group variance estimate and its p-values are meaningless for inference. Use it only for no-replicate exploratory visualization, never for the reported test.
Per-Method Failure Modes
Overdispersion left at the default
Trigger: calculateDiffMeth with replicates and no overdispersion= argument. Mechanism: the default 'none' is a binomial-only logistic LRT that ignores between-replicate variance. Symptom: implausibly many significant CpGs; q-values far smaller than a beta-binomial tool gives on the same data. Fix: overdispersion='MN' (which uses the F-test) whenever replicates exist.
MN plus test='Chisq'
Trigger: passing both overdispersion='MN' and test='Chisq'. Mechanism: MN forces the F-test; test= is silently ignored. Symptom: the reported "chi-square test" was never run. Fix: drop test= when using MN; test= is honored only with overdispersion='none'.
SLIM read as BH
Trigger: comparing methylKit q-value counts to a DSS/limma BH analysis. Mechanism: adjust defaults to SLIM, methylKit's own sliding-linear-model q-method. Symptom: DMC counts disagree with another tool's BH results at the "same" q. Fix: set adjust='BH' for cross-tool comparability and state the method used.
cov.bases=0 tiles
Trigger: tileMethylCounts at the default cov.bases=0. Mechanism: a window with one covered CpG becomes a "region." Symptom: thousands of single-CpG "DMRs," many noisy. Fix: raise cov.bases to >=3.
Tiling or testing the raw object
Trigger: tileMethylCounts(meth_obj, ...) or uniting before filtering/normalizing. Mechanism: low-coverage and library-size artifacts propagate into the tiles and the test. Symptom: artifactual regions; deeper samples look hyper-confident. Fix: filterByCoverage then normalizeCoverage BEFORE tiling/uniting.
Destranding the wrong context or input
Trigger: unite(destrand=TRUE) on .cov input or on CHG/CHH context. Mechanism: .cov carries no strand; non-CpG dyads are not symmetric. Symptom: double-counting or wrong merges. Fix: destrand CpG only, ideally from the cytosine report.
pool() then test
Trigger: pool() followed by calculateDiffMeth. Mechanism: pooling removes within-group variance. Symptom: tiny p-values with no biological meaning. Fix: never pool for the reported test; keep replicates separate.
Quantitative Thresholds
| Threshold | Source | Rationale |
|---|
lo.count = 10 | convention (methylKit tutorial) | below ~10x a single-CpG percentage is a coin flip; not a derived value |
hi.perc = 99.9 | convention | drops the top 0.1% coverage (PCR/repeat artifacts) |
overdispersion = 'MN' with replicates | Akalin 2012 Genome Biol 13:R87 | adds the beta-binomial between-replicate layer; default 'none' over-calls |
adjust = 'BH' (default SLIM) | Akalin 2012 Genome Biol 13:R87 | BH for comparability; SLIM is methylKit-specific |
getMethylDiff difference=25, qvalue=0.01 | methylKit defaults | 25% is tutorial convention, NOT derived; justify per feature/coverage/purity and report it |
tileMethylCounts cov.bases >= 3 | nuance | default 0 admits single-CpG tiles; require real CpG support |
win.size=step.size=1000 (default) | methylKit defaults | step < win gives overlapping (sliding) tiles |
Common Errors
| Error / symptom | Cause | Solution |
|---|
| Implausibly many DMCs | overdispersion='none' under replication | set overdispersion='MN' |
| Reported chi-square never ran | test='Chisq' with MN | MN forces F; drop test= |
| Counts disagree with another tool at same q | default adjust='SLIM' | set adjust='BH', state the method |
| Thousands of single-CpG "regions" | cov.bases=0 | raise cov.bases to >=3 |
| Destrand error / double counts | destrand on .cov or non-CpG | destrand CpG only, from cytosine report |
| Meaningless tiny p-values | pool() before testing | keep replicates; never pool for inference |
| Deeper sample looks more confident | no normalizeCoverage | normalize before unite/test |
References
- Akalin A, Kormaksson M, Li S, Garrett-Bakelman FE, Figueroa ME, Melnick A, Mason CE. 2012. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol 13:R87.
- Krueger F, Andrews SR. 2011. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 27:1571-1572.
- Robinson MD, Kahraman A, Law CW, Lindsay H, Nowicka M, Weber LM, Zhou X. 2014. Statistical methods for detecting differentially methylated loci and regions. Front Genet 5:324.
Related Skills
- methylation-calling - Produces the coverage/cytosine reports read here
- differential-cpg-testing - Per-site statistical model choice (count vs continuous)
- dmr-detection - Selection-aware region callers (dmrseq/DSS) beyond methylKit tiles
- pathway-analysis/go-enrichment - Functional annotation of differentially methylated genes
- long-read-sequencing/nanopore-methylation - Long-read MM/ML calling; pipe counts into this object model
- workflows/methylation-pipeline - End-to-end bisulfite pipeline
