Identifies differential m6A methylation between conditions from MeRIP-seq paired IP/input data using exomePeak2 with `bam_ip` + `bam_input` (control arm) and `bam_treated_ip` + `bam_treated_input` (treatment arm) for integrated GC-bias-aware differential calling (Liu 2022 *NAR Genom Bioinform* 4:lqac046), QNB beta-binomial test (Liu 2017 *BMC Bioinformatics* 18:387), MeTDiff HMM-based differential bundled with MeTPeak, RADAR (Zhang 2019 *Genome Biol* 20:294) with its `filterBins -> diffIP -> reportResult` workflow, and the defensible paired-symmetric case of edgeR / DESeq2 on featureCounts-on-peaks matrices when batch / lot covariates need explicit fixed-effect handling (exomePeak2's top-level API does NOT accept arbitrary covariates). Covers paired vs unpaired vs interaction designs, batch confounding (antibody lot, RNA prep, sequencing run) and the per-lot meta-analysis strategy when exomePeak2 is the primary caller, the stoichiometry-vs-expression-vs-IP-efficiency confound that all MeRIP differential metho
Identifies differential m6A methylation between conditions from MeRIP-seq paired IP/input data using exomePeak2 with `bam_ip` + `bam_input` (control arm) and `bam_treated_ip` + `bam_treated_input` (treatment arm) for integrated GC-bias-aware differential calling (Liu 2022 *NAR Genom Bioinform* 4:lqac046), QNB beta-binomial test (Liu 2017 *BMC Bioinformatics* 18:387), MeTDiff HMM-based differential bundled with MeTPeak, RADAR (Zhang 2019 *Genome Biol* 20:294) with its `filterBins -> diffIP -> reportResult` workflow, and the defensible paired-symmetric case of edgeR / DESeq2 on featureCounts-on-peaks matrices when batch / lot covariates need explicit fixed-effect handling (exomePeak2's top-level API does NOT accept arbitrary covariates). Covers paired vs unpaired vs interaction designs, batch confounding (antibody lot, RNA prep, sequencing run) and the per-lot meta-analysis strategy when exomePeak2 is the primary caller, the stoichiometry-vs-expression-vs-IP-efficiency confound that all MeRIP differential methods inherit, normalisation choice (size factor on IP, on input, on per-sample IP/input ratio), the McIntyre 2020 reproducibility caveat, effect-size filtering as a guardrail against under-powered N=2 designs, and orthogonal-validation routes for absolute stoichiometry (GLORI / SAC-seq / m6Anet `mod_ratio`). Use when comparing m6A levels across two or more conditions, choosing between exomePeak2 / QNB / RADAR / MeTDiff for a given design, handling batch confounding when exomePeak2's API is too rigid, normalising against input properly, distinguishing real hyper- / hypo-methylation from expression-level shifts, applying effect-size thresholds, interpreting volcano plots of differential peaks, or planning a follow-up orthogonal stoichiometry validation.
Before using code patterns, verify installed versions match. If versions differ:
R: packageVersion('exomePeak2') then ?exomePeak2 to verify parameters
R: packageVersion('QNB') then ?qnbtest to confirm argument signature
If R throws unused argument or argument is missing, the API moved between Bioconductor minor releases; consult ?exomePeak2 directly. QNB function signature has been stable since 2017 but the GitHub source has occasional changes; pin the commit SHA.
exomePeak2's differential interface is NOT a mode= argument — populate the bam_treated_ip= and bam_treated_input= arguments alongside the standard bam_ip= / bam_input= (control arm) to trigger paired differential calling. peak_calling_mode is a separate argument controlling locus scope (exon | full_transcript | whole_genome). Verify against ?exomePeak2 and Bioconductor 3.20+ release notes. QNB is on GitHub only; install via devtools::install_github('lzcyzm/QNB'). RADAR is on GitHub only; install via devtools::install_github('scottzijiezhang/RADAR'); its differential workflow is countReads -> normalizeLibrary -> adjustExprLevel -> filterBins -> diffIP -> reportResult.
Differential m6A Analysis
"Compare m6A methylation between my conditions" -> Quantify how much each m6A peak's IP/input enrichment shifts between conditions, after normalising for transcript-abundance changes (which all show up in input) and for IP-efficiency drift (which the design matrix and within-run replicates control for). Then apply effect-size filtering to distinguish real biology from the technical noise floor that all MeRIP differential methods inherit (McIntyre 2020 Sci Rep 10:6590: between-study m6A peak overlap is ~45% median; differential calls within that noise envelope routinely fail to replicate). Critically: a higher MeRIP signal in condition A vs B can mean (1) more transcripts of the peak-bearing gene, (2) more methylation per transcript, OR (3) higher IP efficiency in batch A — distinguishing requires careful normalisation or an orthogonal absolute-stoichiometry method.
R: QNB::qnbtest(control_ip, treated_ip, control_input, treated_input) -- beta-binomial for small N (Liu 2017 BMC Bioinformatics 18:387)
R: RADAR::diffIP() then reportResult() -- Poisson-NB on peak windows with TMM normalisation (Zhang 2019 Genome Biol 20:294)
R: MeTDiff::metdiff() -- HMM-based differential paired with MeTPeak
R: edgeR / DESeq2 on featureCounts-on-peaks matrix -- defensible for paired symmetric designs with strong input normalisation; also the route for arbitrary batch / lot covariates exomePeak2's API does not accept
The Single Most Important Modern Insight -- IP fold-change between conditions conflates stoichiometry change with expression change and IP efficiency drift
A higher m6A peak signal in condition A vs B can mean ANY of: (1) more transcripts of the peak-bearing gene (expression up; the input increases proportionally, so the ratio should not change — but residual normalisation noise leaks through), (2) more methylation per transcript (stoichiometry up; the real biology of interest), (3) higher IP efficiency in batch A (technical; antibody lot, IP day, RNA prep). Differential MeRIP WITHOUT per-window input-normalisation OR an orthogonal stoichiometry-aware method (GLORI Liu 2023 Nat Biotechnol 41:355; SAC-seq Hu 2022 Nat Biotechnol 40:1210; MAZTER-seq Garcia-Campos 2019 Cell 178:731; m6Anet per-read modification rate) cannot separate the three. exomePeak2 differential mode, QNB, RADAR, and MeTDiff all implement per-window IP/input ratio modelling, but each has different default normalisations and the choice matters. Equally critical: McIntyre 2020 Sci Rep 10:6590 showed that "differential" m6A peaks from MeRIP-seq routinely do not replicate between independent studies in nominally identical conditions. The empirical noise floor is high; effect-size filtering (|log2FC| >= 0.5 minimum, often >= 1) AND replicate-direction concordance (the change is consistent in direction across replicates) AND minimum N=3 per condition are needed for differential calls to survive replication. For any absolute stoichiometry claim ("this peak is 80% methylated in tumour vs 20% in normal"), require an orthogonal stoichiometry method, not MeRIP alone.
Algorithmic Taxonomy
Tool / mode
Mechanism
Inputs
Output
Strength
Fails when
exomePeak2 differential (Liu 2022)
Transcript-windowed Poisson GLM with GC-bias correction; integrated peak + differential
Sister to QNB from same group; small-N alternative
Less benchmarked than QNB / exomePeak2
edgeR / DESeq2 on peak counts
Generic RNA-seq differential framework applied to featureCounts-on-peaks
peak count matrix + sample sheet
Differential peaks with log2FC + FDR
Familiar; flexible designs; well-tested in RNA-seq
Treats peak counts as RNA counts; loses IP/input pairing structure; defensible only for paired symmetric designs with strong input normalisation
Ratio-of-ratios (heuristic)
Compute per-peak log2 (IP_A / Input_A) - log2 (IP_B / Input_B) per sample, then t-test
per-peak count matrix
per-peak t-test
Transparent; no model assumptions
No multiple-testing correction; ignores overdispersion; not recommended for primary analysis
Decision Tree by Scenario
Scenario
Recommended
Why wrong choices fail
Standard 3-vs-3 paired-design MeRIP differential
exomePeak2 with bam_ip (ctrl) + bam_treated_ip (treat) and matching inputs
QNB usable but designed for smaller N; edgeR/DESeq2 loses IP/input pairing
Very small N (2 vs 2)
QNB (designed for small-sample overdispersion); supplement with exomePeak2 if possible
edgeR / DESeq2 dispersion estimation collapses; exomePeak2 GLM also struggles at N=2
Paired design (patient as blocking factor)
QNB per-pair then aggregate; OR featureCounts-on-peaks -> DESeq2 with ~patient + condition design; exomePeak2 cannot encode patient blocking via its top-level API
Unpaired analysis inflates within-group variance
Interaction design (genotype × treatment)
featureCounts-on-peaks -> DESeq2 / edgeR with interaction term; exomePeak2 top-level API is two-group only
QNB pairwise only; build interaction model from pairwise contrasts manually
Batch confounding (antibody lot, sequencing run, IP day)
Include batch as fixed effect in DESeq2 / edgeR model on featureCounts-on-peaks matrix; OR run exomePeak2 per-batch and meta-analyse
Pooling cross-batch counts without batch term attributes lot-effect to condition
Time-course differential
featureCounts-on-peaks -> DESeq2 / limma with time as numeric covariate; OR pairwise time-point exomePeak2 contrasts
featureCounts-on-peaks -> DESeq2 / edgeR with lot in design; OR run exomePeak2 per-lot then meta-analyse; ideally avoid confounding lot with condition at the experimental-design stage
Lot effect inflates false positives; exomePeak2 top-level API cannot encode lot
Validation of differential calls
Run >=2 differential methods; require concordant direction across replicates AND
log2FC
Visualising differential peaks
Volcano plot with
log2FC
Wanting to test a single gene / locus
Targeted: per-peak boxplot across replicates with condition factor; manual t-test or Wilcoxon at high-coverage peak
Whole-transcriptome differential testing wastes multiple-testing budget for single-locus questions
Methodology evolves; before any high-stakes differential analysis, web-search "exomePeak2 differential mode Bioconductor 3.20" and "MeRIP differential benchmark McIntyre" for current consensus parameters.
exomePeak2 Differential Workflow
Goal: Identify m6A peaks that differ in methylation level between conditions, controlling for transcript-abundance differences (via input normalisation) and GC bias (via internal correction), with an integrated peak-calling + differential pipeline.
Approach: Build TxDb from the matched GTF; pass control IP/input BAM vectors via bam_ip and bam_input AND treatment IP/input BAM vectors via bam_treated_ip and bam_treated_input — populating the treated arms triggers differential mode (there is no separate mode= argument). Output is per-peak log2FC + FDR.
peak_calling_mode accepts 'exon' (transcript-aware, default), 'full_transcript', or 'whole_genome'; the meaning is locus scope, NOT differential-vs-non-differential. For arbitrary covariate adjustment (batch, antibody lot, patient blocking), the exomePeak2 top-level API is insufficient — move counts into DESeq2 / edgeR via the featureCounts-on-peaks route below.
QNB Beta-Binomial for Small-N Designs
Goal: Test differential m6A at pre-called peaks using a quad-negative-binomial model that handles small-N overdispersion better than generic GLM frameworks.
Approach: Count reads in each IP and Input BAM at each peak using featureCounts or summarizeOverlaps; pass the four count matrices (ip1, ip2, input1, input2 — ip/input per group) to qnbtest().
Verify QNB argument names against ?qnbtest for the installed version; older tutorials may show different signatures.
RADAR Reproducibility-Aware Differential
Goal: Test differential m6A using a Poisson-NB framework with TMM normalisation and explicit replicate-variance modeling; useful when replicate variability is a known issue.
Approach: RADAR's documented workflow is countReads -> normalizeLibrary -> adjustExprLevel -> filterBins -> diffIP -> reportResult. CRITICAL: RADAR expects matched BAMs in bamFolder named <sample>.input.bam and <sample>.m6A.bam per replicate; the generic IP / Input naming used elsewhere must be re-conformed or symlinked. variable() is set with a data.frame, NOT a bare factor.
reportResult thresholds (cutoff = p-value cutoff; Beta_cutoff = effect-size cutoff in beta units) are RADAR-specific; convert to the project's standard reporting thresholds downstream. The aligned_radar/ directory should contain BAMs with RADAR's expected naming (<sample>.input.bam, <sample>.m6A.bam).
Volcano Plot of Differential Peaks
Goal: Visualise the differential peak set with effect size on the x-axis and statistical significance on the y-axis; flag peaks passing |log2FC| and FDR thresholds.
Approach: Standard ggplot2 volcano with colour-coded significance and threshold lines.
The caption is intentional: MeRIP differential reports CHANGES IN ENRICHMENT RATIO, NOT changes in absolute stoichiometry. For stoichiometry claims, cross-validate with GLORI / SAC-seq / m6Anet.
Per-Method Failure Modes
Reporting "hyper-methylation" without orthogonal calibration
Trigger: "Peak X shows hyper-methylation in treatment" inferred from MeRIP IP fold-change alone.
Mechanism: MeRIP IP fold-change conflates per-molecule methylation stoichiometry, transcript abundance, and IP efficiency variation between libraries. An IP fold-change increase can reflect any or all of these.
Symptom: Reported m6A "hyper-methylation" tracks RNA-seq expression changes between conditions; reverse-direction effects when properly normalised against input.
Fix: Use "increased / decreased enrichment" terminology for MeRIP-only studies. Reserve "hyper- / hypo-methylated" for studies with absolute quantification orthogonal validation (GLORI, SAC-seq, MAZTER-seq, m6Anet per-read). For high-stakes claims at named loci, run GLORI on a subset of conditions.
Effect-size threshold absence
Trigger: Reporting "1,500 differential m6A peaks" with FDR < 0.05 (uncorrected p-value or naive multiple testing) and no effect-size filter.
Mechanism: With sufficient sequencing depth, MeRIP-seq has high statistical power to detect very small (~1.1-1.2x) IP-ratio changes that lie within antibody / technical noise. McIntyre 2020 Sci Rep 10:6590 showed these changes do not replicate.
Symptom: Differential peak set has many peaks with small effect sizes; replication in an independent study recovers <30% of original calls.
Fix: Apply effect-size filter (|log2FC| >= 0.5 minimum, often >= 1) AND adjusted p-value (FDR < 0.05) AND replicate-direction concordance. Differential peaks should be reported with effect size, not just p-value. Report effect-size distribution alongside peak count.
Underpowered N=2 design
Trigger: Differential m6A study with N=2 IP and N=2 input per condition.
Mechanism: Per McIntyre 2020 and many subsequent benchmarks, MeRIP replicate variance is high; N=2 estimates of dispersion are unreliable; differential calls are unstable.
Symptom: Many "differential" peaks; volcano plot dense; small fraction replicates in held-out replicate.
Fix: Minimum N=3 per condition (per condition per IP/input arm = 12 BAMs for a 2-condition study); N=4-5 preferred for high-stakes claims. Underpowered studies should report effect-size-only filtered subsets (the most extreme peaks) and acknowledge the noise floor explicitly.
Batch confounded with condition (antibody lot, IP day)
Trigger: Control samples processed in batch 1 with antibody lot A; treatment samples processed in batch 2 with antibody lot B.
Mechanism: Anti-m6A antibody lots have batch-to-batch variability in pulldown efficiency and m6A-vs-m6Am cross-reactivity. Pooling cross-batch counts in a differential model attributes batch-effect to condition.
Symptom: "Differential" peaks at high-abundance transcripts; effect sizes track batch rather than condition; reanalysis with batch in the model removes most differential peaks.
Fix: Include antibody_lot / batch / prep_day as a fixed effect in a DESeq2 / edgeR model on featureCounts-on-peaks counts; OR run exomePeak2 separately per lot and meta-analyse the per-lot differential peak sets; OR re-design the experiment to avoid lot-condition confounding. exomePeak2's top-level API does NOT accept arbitrary covariates — DESeq2 / edgeR is the route when covariate handling is required.
exomePeak2 invoked with a fabricated mode= argument
Trigger:exomePeak2(..., mode='differential') OR exomePeak2(..., mode='diff_peak') returns "unused argument" error.
Mechanism: exomePeak2 has NO mode= argument. Differential is triggered by populating bam_treated_ip and bam_treated_input alongside the standard bam_ip and bam_input (control arm). The peak_calling_mode argument is unrelated — it accepts 'exon' | 'full_transcript' | 'whole_genome' and controls locus scope, not differential-vs-non-differential.
Fix: Populate the four BAM-vector arguments (bam_ip, bam_input, bam_treated_ip, bam_treated_input); drop any mode= reference; consult ?exomePeak2 for the authoritative signature in the installed version.
Treating peak count matrices like RNA count matrices in edgeR/DESeq2
Trigger: Compute featureCounts at peaks, build a count matrix, run edgeR / DESeq2 on the matrix as if peaks were genes.
Mechanism: Peak counts reflect both IP enrichment AND transcript abundance. Generic RNA-seq DE on peak counts mixes the two; size-factor normalisation on IP-only counts loses the input-pair information.
Fix: edgeR / DESeq2 on peak counts is defensible ONLY for paired symmetric designs where input is modelled as an offset (per-sample size factor on input counts AND per-sample size factor on IP counts, then differential on the ratio). For most uses, exomePeak2 / QNB / RADAR's purpose-built models are more appropriate.
Counting reads at peaks WITHOUT featureCounts strand-awareness
Trigger:featureCounts(...) invoked without strandSpecific= flag; or with wrong strand setting.
Mechanism: Strand-aware counting matters when the protocol is stranded (most modern MeRIP is unstranded; some are reverse-stranded). Wrong strand counts include antisense reads as if they were sense.
Fix: Verify protocol strandedness from sequencing-core notes or by inspecting featureCounts summary at a few transcripts. Pass strandSpecific=0 (unstranded), 1 (forward), or 2 (reverse) explicitly.
Reconciliation: When Differential Methods Disagree
Pattern
Likely cause
Action
exomePeak2 calls a differential peak; QNB does not
Different dispersion estimates; QNB more conservative at low coverage
Trust intersection; report concordant set as high-confidence
Most "differential" peaks fall at high-abundance housekeeping transcripts
Expression / IP-efficiency confound
Check log2FC vs input log2FC; if correlated, batch effect or expression-driven
RADAR vs exomePeak2 disagree on direction
Different normalisation (TMM vs internal)
Inspect normalisation diagnostic; choose method aligned with experimental design
Differential peaks anti-correlated with RNA-seq DE
Expression conflated with methylation in the differential model
Re-run with stronger input normalisation; consider per-peak ratio normalisation
Single differential peak survives across all methods
High-confidence call
Orthogonally validate (GLORI / SAC-seq) at the named locus
Differential calls scatter randomly across genome
Underpowered; technical noise dominates
Increase N; apply stricter effect-size filter; report null result honestly
Cross-condition peak overlap < 50% before differential
Conditions are biologically very different; OR antibody lot effect
Inspect cross-replicate concordance; check antibody lot metadata
Volcano shows extreme outliers at low-coverage peaks
Per-peak variance dominated by Poisson sampling
Filter peaks by minimum coverage (>=30 reads in IP AND input) before differential
Yes — antibody_lot included as fixed effect in DESeq2 model on featureCounts-on-peaks counts (exomePeak2's top-level API does not accept covariates); per-lot exomePeak2 + meta-analysis when DESeq2 path infeasible; lot-condition confounding assessed
"Does the differential signal track expression changes?"
Cross-checked log2FC vs input log2FC per peak; only peaks with strong IP/input ratio shift reported as differential
"Was orthogonal validation done?"
Top hits orthogonally validated via GLORI / m6Anet / per-locus assay
"Why exomePeak2 over QNB?"
exomePeak2 default for standard 3-vs-3; QNB used for small-N sensitivity analyses; both reported as concordance check
Per-window IP/input ratio normalisation; exomePeak2 internal; alternative TMM (RADAR) compared
"How does this replicate in independent studies?"
Cross-study peak overlap reported; differential subset checked against published m6A-Atlas
"Why not edgeR / DESeq2?"
Generic RNA-seq DE on peak counts loses IP/input pairing structure; used only for sensitivity analysis with paired symmetric design
References
Dominissini D, Moshitch-Moshkovitz S, Schwartz S et al (2012) Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature 485(7397):201-206. doi:10.1038/nature11112
Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR (2012) Comprehensive analysis of mRNA methylation reveals enrichment in 3' UTRs and near stop codons. Cell 149(7):1635-1646. doi:10.1016/j.cell.2012.05.003
Liu L, Zhang SW, Huang Y, Meng J (2017) QNB: differential RNA methylation analysis for count-based small-sample sequencing data with a quad-negative binomial model. BMC Bioinformatics 18(1):387. doi:10.1186/s12859-017-1808-4
Cui X, Meng J, Zhang S, Chen Y, Huang Y (2016) A novel algorithm for calling mRNA m6A peaks by modeling biological variances in MeRIP-seq data. Bioinformatics 32(12):i378-i385. doi:10.1093/bioinformatics/btw281
Liu L, Zhang SW, Gao F et al (2016) DRME: count-based differential RNA methylation analysis at small sample size scenario. Anal Biochem 499:15-23. doi:10.1016/j.ab.2016.01.014
Zhang Z, Zhan Q, Eckert M et al (2019) RADAR: differential analysis of MeRIP-seq data with a random effect model. Genome Biol 20(1):294. doi:10.1186/s13059-019-1915-9
Meng J, Lu Z, Liu H et al (2014) A protocol for RNA methylation differential analysis with MeRIP-Seq data and exomePeak R/Bioconductor package. Methods 69(3):274-281. doi:10.1016/j.ymeth.2014.06.008
Liu J, Zhang Z, Meng J et al (2022) exomePeak2: a peak calling and differential analysis tool for MeRIP-Seq with bias awareness. NAR Genom Bioinform 4(3):lqac046. doi:10.1093/nargab/lqac046
McIntyre ABR, Gokhale NS, Cerchietti L, Jaffrey SR, Horner SM, Mason CE (2020) Limits in the detection of m6A changes using MeRIP/m6A-seq. Sci Rep 10(1):6590. doi:10.1038/s41598-020-63355-3
Liu C, Sun H, Yi Y et al (2023) Absolute quantification of single-base m6A methylation in the mammalian transcriptome using GLORI. Nat Biotechnol 41(3):355-366. doi:10.1038/s41587-022-01487-9
Hu L, Liu S, Peng Y et al (2022) m6A RNA modifications are measured at single-base resolution across the mammalian transcriptome. Nat Biotechnol 40(8):1210-1219. doi:10.1038/s41587-022-01243-z
Garcia-Campos MA, Edelheit S, Toth U et al (2019) Deciphering the m6A code via antibody-independent quantitative profiling. Cell 178(3):731-747.e16. doi:10.1016/j.cell.2019.06.013
Robinson MD, McCarthy DJ, Smyth GK (2010) edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26(1):139-140. doi:10.1093/bioinformatics/btp616
Love MI, Huber W, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15(12):550. doi:10.1186/s13059-014-0550-8
Liao Y, Smyth GK, Shi W (2014) featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30(7):923-930. doi:10.1093/bioinformatics/btt656
