M6a Peak Calling

What the source does

• FilesystemNone observedNo reads or writes in the source
• NetworkNone observedNo outbound calls in the source
• ShellNone observedNo commands run in the source
• DependenciesPinnedNothing floats
• SecretsNone foundNo keys in source

Five audit layers

Gates applied: no behavioural pass.

Audit 0.4.1 · source sha c8c1fc3a2b13 · 19 Oct 2026Report an issue or request a re-scan

• PublisherGPTomics GitHub organizationOwner of the repository this skill is read from
• Source folderGPTomics/bioSkills/epitranscriptomics/m6a-peak-callingThe exact folder we audited
• LicenceMITAs the repository states it
• Commit pinned for this verdictc8c1fc3a2b13A new audit runs when the source changes
• Repository statusArchivedLast push 65 days ago
• Sandboxed runNot runSkills are read, not executed; layer L4 is skipped

Its description: Calls m6A peaks from MeRIP-seq / m6A-seq paired IP-vs-input data using exomePeak2 (transcript-aware, GC-bias-corrected Poisson GLM), MeTPeak (HMM over sliding windows), MACS3/MACS2 with --nomodel --broad --keep-dup all…

---
name: bio-epitranscriptomics-m6a-peak-calling
description: Calls m6A peaks from MeRIP-seq / m6A-seq paired IP-vs-input data using exomePeak2 (transcript-aware, GC-bias-corrected Poisson GLM), MeTPeak (HMM over sliding windows), MACS3/MACS2 with --nomodel --broad --keep-dup all (genome-wide broad alternative), and DRACH motif enrichment via HOMER or ggseqlogo as a sanity check (NOT a filter). Covers BED12 vs narrowPeak output, exonic vs intronic peak handling, multi-tool reconciliation (intersection vs union), the m6A-vs-m6Am ambiguity at 5'UTR peaks that antibody methods cannot resolve, and orthogonal validation (miCLIP/GLORI/m6A-SAC-seq/m6Anet). Use when calling peaks from paired IP/input genome BAMs, choosing exomePeak2 (transcript-aware default) vs MACS3 (broad genomic) vs MeTPeak (HMM-smoothed low-coverage), confirming DRACH enrichment as a sanity check on the peak set, reconciling differing peak sets across tools, validating MeRIP peaks against single-base methods, interpreting 5' peaks where m6Am contamination is possible, or recommending a consensus strategy.
tool_type: mixed
primary_tool: exomePeak2
---

## 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 smoot

Excerpt, first 6,000 characters of SKILL.md, verbatim from the audited commit. It is text from someone else’s repository, quoted so you can judge it before your agent reads it.

git clone https://github.com/GPTomics/bioSkills
git -C bioSkills checkout c8c1fc3a2b13

Works in

1. 2026-10-19
   ✓SAFE B 89 / 100

   first audit · source c8c1fc3a2b13

Questions about M6a Peak Calling

• github.com/GPTomics/bioSkills/tree/main/epitranscriptomics/m6a-p — the audited source at commit c8c1fc3a2b13 and every finding aboveObserved · 19 Oct 2026

Spotted something wrong? Report a correction · Methodology