Bio Epitranscriptomics M6a Differential

Version Compatibility

Reference examples tested with: exomePeak2 1.14+ (Bioconductor 3.18+), QNB 1.1.11 (GitHub lzcyzm/QNB), MeTDiff (GitHub, bundled with MeTPeak), RADAR 0.2.4+ (GitHub scottzijiezhang/RADAR), DESeq2 1.42+, edgeR 4.0+, GenomicFeatures 1.54+, ggplot2 3.5+, GenomicAlignments 1.38+, Rsubread 2.16+.

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: exomePeak2::exomePeak2(bam_ip, bam_input, bam_treated_ip, bam_treated_input, ...) -- integrated peak + differential GLM (modern default)
• 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

Decision Tree by Scenario

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.

library(exomePeak2)
library(GenomicFeatures)
library(BSgenome.Hsapiens.UCSC.hg38)

txdb <- makeTxDbFromGFF('refs/annotation.gtf', format='gtf')

ctrl_ip      <- c('aligned/ctrl_IP1.bam', 'aligned/ctrl_IP2.bam', 'aligned/ctrl_IP3.bam')
ctrl_input   <- c('aligned/ctrl_Input1.bam', 'aligned/ctrl_Input2.bam', 'aligned/ctrl_Input3.bam')
treat_ip     <- c('aligned/treat_IP1.bam', 'aligned/treat_IP2.bam', 'aligned/treat_IP3.bam')
treat_input  <- c('aligned/treat_Input1.bam', 'aligned/treat_Input2.bam', 'aligned/treat_Input3.bam')

result <- exomePeak2(
    bam_ip             = ctrl_ip,
    bam_input          = ctrl_input,
    bam_treated_ip     = treat_ip,
    bam_treated_input  = treat_input,
    txdb               = txdb,
    genome             = BSgenome.Hsapiens.UCSC.hg38,
    paired_end         = TRUE,
    library_type       = 'unstranded',
    peak_calling_mode  = 'exon',
    save_dir           = 'exomepeak2_diff_output',
    experiment_name    = 'ctrl_vs_treat'
)

diff_table <- as.data.frame(result)
nrow(diff_table)
head(diff_table[, c('seqnames', 'start', 'end', 'log2FC', 'pvalue', 'padj')])

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().

library(QNB)
library(Rsubread)
library(rtracklayer)

peaks <- import('exomepeak2_output/m6a_run1/peaks.bed')
peak_saf <- data.frame(
    GeneID = paste0('peak_', seq_along(peaks)),
    Chr    = as.character(seqnames(peaks)),
    Start  = start(peaks),
    End    = end(peaks),
    Strand = as.character(strand(peaks))
)

count_matrix <- function(bam_paths, peak_saf) {
    fc <- featureCounts(
        files       = bam_paths,
        annot.ext   = peak_saf,
        isPairedEnd = TRUE,
        nthreads    = 8,
        allowMultiOverlap = TRUE
    )
    fc$counts
}

ip_ctrl   <- count_matrix(c('aligned/ctrl_IP1.bam', 'aligned/ctrl_IP2.bam', 'aligned/ctrl_IP3.bam'), peak_saf)
ip_treat  <- count_matrix(c('aligned/treat_IP1.bam', 'aligned/treat_IP2.bam', 'aligned/treat_IP3.bam'), peak_saf)
in_ctrl   <- count_matrix(c('aligned/ctrl_Input1.bam', 'aligned/ctrl_Input2.bam', 'aligned/ctrl_Input3.bam'), peak_saf)
in_treat  <- count_matrix(c('aligned/treat_Input1.bam', 'aligned/treat_Input2.bam', 'aligned/treat_Input3.bam'), peak_saf)

qnb_result <- qnbtest(
    control_ip    = ip_ctrl,
    treated_ip    = ip_treat,
    control_input = in_ctrl,
    treated_input = in_treat,
    mode          = 'per-condition'
)

head(qnb_result)
sig <- qnb_result[qnb_result$padj < 0.05 & abs(qnb_result$log2.RR) > 0.5, ]
nrow(sig)

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.

library(RADAR)

radar <- countReads(
    samplenames  = c('ctrl_rep1', 'ctrl_rep2', 'ctrl_rep3', 'treat_rep1', 'treat_rep2', 'treat_rep3'),
    gtf          = 'refs/annotation.gtf',
    bamFolder    = 'aligned_radar/',
    modification = 'm6A',
    strandToKeep = 'opposite',
    threads      = 8
)

radar <- normalizeLibrary(radar)
radar <- adjustExprLevel(radar)

variable(radar) <- data.frame(group = c('ctrl', 'ctrl', 'ctrl', 'treat', 'treat', 'treat'))

radar <- filterBins(radar, minCountsCutOff = 15)
radar <- diffIP(radar)

result <- reportResult(radar, cutoff = 0.1, Beta_cutoff = 0.5)
sig <- result[result$padj < 0.05 & abs(result$logFC) > 0.5, ]
nrow(sig)

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.

library(ggplot2)

diff_table$significance <- with(diff_table,
    ifelse(padj < 0.05 & abs(log2FC) > 0.5, 'differential', 'not_sig'))

ggplot(diff_table, aes(x=log2FC, y=-log10(padj), colour=significance)) +
    geom_point(alpha=0.5, size=0.8) +
    geom_vline(xintercept=c(-0.5, 0.5), linetype='dashed') +
    geom_hline(yintercept=-log10(0.05), linetype='dashed') +
    scale_colour_manual(values=c(differential='red', not_sig='grey60')) +
    labs(x='log2 (treat / ctrl) MeRIP enrichment ratio',
         y='-log10 (FDR)',
         title='Differential m6A peaks: ctrl vs treat',
         caption='Per-window IP/input ratio normalised; not absolute stoichiometry') +
    theme_minimal()

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.