bio-epitranscriptomics-m6a-peak-calling

What this skill does

## Version Compatibility

Reference examples tested with: exomePeak2 1.14+ (Bioconductor 3.18+), MeTPeak (GitHub commit SHA-pinned; no Bioconductor release), MACS3 3.0+, MACS2 2.2.9+, samtools 1.19+, GenomicFeatures 1.54+, BSgenome.Hsapiens.UCSC.hg38 1.4+, HOMER 4.11+, ggseqlogo 0.2+, rtracklayer 1.62+, GenomicRanges 1.54+, ChIPseeker 1.38+.

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('exomePeak2')` then `?exomePeak2` to verify parameters
- CLI: `macs3 callpeak --help`, `findMotifsGenome.pl` to confirm flags

If R throws `unused argument` or `argument is missing`, the exomePeak2 / MeTPeak API may have moved between releases; consult `?exomePeak2` and the installed package's NAMESPACE.

exomePeak2 differential is triggered by populating `bam_treated_ip` / `bam_treated_input` alongside the control-arm `bam_ip` / `bam_input` (see m6a-differential); there is NO `mode=` argument. MeTPeak is GitHub-only and unversioned — pin a commit SHA in any reproducible analysis; defaults are `WINDOW_WIDTH=50, SLIDING_STEP=50, FRAGMENT_LENGTH=100`. MACS3 supersedes MACS2 in development activity but both are widely used; default `--keep-dup` is 1 in BOTH and MUST be overridden to `all` for MeRIP. HOMER `findMotifsGenome.pl` is mature; the `-rna` flag is the RNA-mode entry.

# m6A Peak Calling

**"Find m6A sites in my MeRIP data"** -> Compare paired IP and input read distributions per transcript window, call windows with statistically significant IP enrichment as peaks, annotate against transcript features (5'UTR / CDS / 3'UTR / stop-codon), and confirm DRACH motif enrichment on the peak set relative to background as a sanity check that the antibody-IP worked. CRITICAL: DRACH is a sanity check on the WHOLE peak set, NOT a filter on individual peaks (filtering by DRACH drops 5-10% of real m6A sites and creates circular validation). Equally critical: peaks within ~50 nt of TSS may be PCIF1 m6Am, not METTL3 m6A — anti-m6A antibodies cross-react.

- R: `exomePeak2::exomePeak2()` -- transcript-aware GC-corrected GLM (field default)
- R: `MeTPeak::metpeak()` -- HMM-smoothed sliding-window alternative
- CLI: `macs3 callpeak --nomodel --keep-dup all --broad` -- broad genomic alternative
- CLI: `findMotifsGenome.pl peaks.bed -rna` (HOMER) -- DRACH enrichment sanity check
- R: `ggseqlogo::ggseqlogo()` -- peak-centre 5-mer sequence logo

## The Single Most Important Modern Insight -- MeRIP cannot distinguish m6A from m6Am near the 5' end, and DRACH is a sanity check not a filter

Anti-m6A antibodies (Synaptic Systems 202-003, Abcam ab151230, NEB EpiMark E1610, Cell Signaling 56593, Active Motif 61755) cross-react with m6Am — the cap-adjacent N6,2'-O-dimethyladenosine at the +1 nucleotide of capped mRNAs, installed by PCIF1 / CAPAM (Akichika 2019 *Science* 363:eaav0080; Boulias 2019 *Mol Cell* 75:631). METTL3 / METTL14 do NOT methylate the cap-adjacent position; PCIF1 does. Peaks within the first ~50 nt of a transcript are ambiguous between m6A and m6Am. METTL3-KO will REMOVE internal m6A peaks but LEAVE the cap m6Am peaks intact, which has caused multiple papers to mis-attribute METTL3-independent peaks to non-canonical writers when they are really PCIF1 m6Am. Linder 2015 *Nat Methods* 12:767 (miCLIP, the antibody-based single-base method) explicitly notes the cross-reactivity. Separately: while ~70% of mammalian METTL3-deposited m6A sites sit within the DRACH consensus (D=A/G/U, R=A/G, A=methylated, C=C, H=A/C/U), a non-trivial fraction are non-DRACH. Post-hoc filtering peaks by DRACH content drops real m6A peaks; DRACH should be reported as enrichment-relative-to-background (HOMER / MEME / ggseqlogo) on the PEAK SET as a whole, NOT as a per-peak filter. For unambiguous internal-m6A studies, restrict analysis to peaks past the first ~50 nt AND validate at high-stakes sites with an orthogonal method (miCLIP for single-base antibody validation; GLORI Liu 2023 *Nat Biotechnol* 41:355 for absolute stoichiometry; m6A-SAC-seq Hu 2022 *Nat Biotechnol* 40:1210; m6Anet for ONT direct-RNA confirmation).

## Algorithmic Taxonomy

| Tool / mode | Mechanism | Inputs | Output | Strength | Fails when |
|-------------|-----------|--------|--------|----------|------------|
| exomePeak2 (Liu 2022 *NAR Genom Bioinform* 4:lqac046) | Transcript-windowed Poisson GLM with on-the-fly GC-bias correction; supersedes exomePeak v1 | paired IP/input BAM + TxDb / GTF | BED12 + RDS + per-peak fold-change / FDR | Transcript-aware; GC-aware; integrates motif annotation; modern field default | Slow on very large datasets; argument signatures shift between Bioconductor minor releases — verify against `?exomePeak2` |
| MeTPeak (Cui 2016 *Bioinformatics* 32:i378) | HMM over sliding windows with Beta-binomial emission per window | paired IP/input BAM + GTF (or TxDb) | BED12 | HMM smooths spatial dependency; better at low coverage | GitHub-only; unversioned; default window/step is 50/50 not the small values some tutorials suggest |
| MACS3 / MACS2 broad mode (Zhang 2008 *Genome Biol* 9:R137) | Sliding-window negative-binomial test; `--broad` extends; `--nomodel` disables ChIP fragment-shift model | paired IP/input BAM | narrowPeak / broadPeak | Battle-tested ChIP-seq lineage; very fast | Not transcript-aware; misses GC-confounded peaks; default `--keep-dup 1` collapses MeRIP signal at high-coverage transcripts |
| MeRIPtools (R wrapper) | Bundles exomePeak / MeTPeak / motif / annotation steps | FASTQ -> peaks pipeline | full report | Reproducible end-to-end | Less flexibility than calling tools separately |
| MoAIMS | Mixture model alternative | paired IP/input BAM | peaks | Smaller user base; less benchmarked | Niche use |
| m6Aboost (R) | Boost peak-calling sensitivity by leveraging DRACH motif as a prior | paired IP/input BAM + motif file | refined peak set | Improves sensitivity in low-coverage regions | Builds DRACH into the calling — DON'T use as evidence for DRACH enrichment downstream (circular) |
| m6ACali (recent ML peak filter; verify current citation against the project repo) | ML-based false-peak filter trained on exomePeak2 + MACS2 outputs across many cell lines | called peak set + IP/input BAM | refined peak set | Modern QC layer; cuts antibody artifact peaks | Trained on specific antibody clones; verify it generalises to the antibody used |

## Decision Tree by Scenario

| Scenario | Recommended | Why wrong choices fail |
|----------|-------------|------------------------|
| Standard mammalian MeRIP, 3+ replicates per arm | exomePeak2 + DRACH confirmation; reconcile with MeTPeak as a second opinion | MACS3 misses GC-confounded peaks; MeTPeak alone gives less GC awareness |
| Viral / kilobase-broad peaks | MACS3 `--broad --broad-cutoff 0.1 --keep-dup all` | exomePeak2 splits broad enrichment into many small per-window peaks |
| Low-coverage / scarce-sample MeRIP | MeTPeak HMM smooths; exomePeak2 as cross-check | MACS3 default loses sensitivity at low coverage |
| Need single-nucleotide resolution | NOT MeRIP -- switch to miCLIP, m6Anet, GLORI, SAC-seq | MeRIP windows are ~100-200 nt; cannot resolve to base |
| 5'UTR / cap-proximal peaks | Run normally BUT flag 5' peaks (within ~50 nt of TSS) as m6A-or-m6Am ambiguous; validate with PCIF1-KO if available | Antibody cross-reacts with m6Am; cannot assign without orthogonal data |
| Cross-tool reconciliation | Call with exomePeak2 + MeTPeak + MACS3 broad; intersect (NOT union) for high-confidence; report each separately | Union inflates false-positive rate; single-tool reports under-call ~30% of consensus peaks |
| Validation of high-stakes peak set | Cross-check against published m6A-Atlas / REPIC databases; orthogonal validation at top hits (miCLIP / GLORI / m6A-SAC-seq) | Single-method-single-study peaks have ~50% inter-study overlap (McIntyre 2020) |
| Wanting absolute stoichiometry | NOT MeRIP -- use GLORI (Liu 2023 *Nat Biotechnol*), SAC-seq (Hu 2022), eTAM-seq (Xiao 2023) | MeRIP IP fold-cha