Bio Epitranscriptomics M6a Peak Calling

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

Decision Tree by Scenario

Methodology evolves; before any high-stakes peak-calling analysis, web-search "exomePeak2 Bioconductor release notes" and "MeRIP peak caller benchmark 2024" for current consensus parameters.

exomePeak2 Standard Workflow

Goal: Produce a transcript-aware set of m6A peaks with FDR and IP/input fold-change from paired IP/input genome BAM files, suitable as input to differential analysis, motif scanning, or downstream visualisation.

Approach: Build a TxDb from the matched GTF; pass paired IP/input BAM vectors with bam_ip and bam_input; let exomePeak2 handle GC-bias correction internally; export BED12 + RDS for downstream use; annotate against transcript features.

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

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

result <- exomePeak2(
    bam_ip       = c('aligned/IP_rep1.bam', 'aligned/IP_rep2.bam', 'aligned/IP_rep3.bam'),
    bam_input    = c('aligned/Input_rep1.bam', 'aligned/Input_rep2.bam', 'aligned/Input_rep3.bam'),
    txdb         = txdb,
    genome       = BSgenome.Hsapiens.UCSC.hg38,
    paired_end   = TRUE,
    library_type = 'unstranded',
    save_dir     = 'exomepeak2_output',
    experiment_name = 'm6a_run1'
)

# Inspect peaks: GRanges with metadata
peaks <- result
length(peaks)
head(as.data.frame(peaks))

exomePeak2() writes BED12 + RDS + per-peak fold-change / FDR to save_dir/experiment_name/. The txdb argument is the canonical interface; older tutorials use gff= (path) which is deprecated. BSgenome is required for GC correction; if unavailable, exomePeak2 falls back to a less-accurate GC-uncorrected mode.

MeTPeak HMM-Smoothed Alternative

Goal: Call peaks with HMM smoothing across spatial windows, useful when coverage is low or when peak boundaries are ambiguous from per-window negative-binomial tests alone.

Approach: MeTPeak accepts IP and INPUT BAM vectors plus either a GTF file path (GENE_ANNO_GTF=) or a TxDb object (TXDB=); defaults are WINDOW_WIDTH=50, SLIDING_STEP=50, FRAGMENT_LENGTH=100, MINIMAL_PEAK_LENGTH=FRAGMENT_LENGTH/2, MINIMAL_MAPQ=30. The example below uses GTF input with explicit defaults shown for transparency; adjust only with rationale.

library(MeTPeak)

metpeak(
    IP_BAM              = c('aligned/IP_rep1.bam', 'aligned/IP_rep2.bam', 'aligned/IP_rep3.bam'),
    INPUT_BAM           = c('aligned/Input_rep1.bam', 'aligned/Input_rep2.bam', 'aligned/Input_rep3.bam'),
    GENE_ANNO_GTF       = 'refs/annotation.gtf',
    OUTPUT_DIR          = 'metpeak_output',
    EXPERIMENT_NAME     = 'm6a_run1',
    WINDOW_WIDTH        = 50,
    SLIDING_STEP        = 50,
    FRAGMENT_LENGTH     = 100,
    MINIMAL_PEAK_LENGTH = 50,
    PEAK_CUTOFF_PVALUE  = 1e-5,
    PEAK_CUTOFF_FDR     = 0.05,
    FOLD_ENRICHMENT     = 1
)

MeTPeak accepts either GENE_ANNO_GTF= (a file path) or TXDB= (a TxDb object); the GTF path is the more common usage and is documented in the GitHub README. Output BED12 is at OUTPUT_DIR/EXPERIMENT_NAME/peak.bed. MeTPeak is GitHub-only (compgenomics/MeTPeak) and unversioned; pin a commit SHA in reproducible analyses.

MACS3 Broad-Peak Alternative

Goal: Call broad MeRIP peaks across the genome using the ChIP-seq sliding-window negative-binomial framework; useful for viral genomes where peaks span kilobases, or as a cross-caller second opinion.

Approach: MACS3 callpeak with --nomodel --extsize 150 --keep-dup all --broad. The --nomodel flag disables the ChIP-seq fragment-shift model (designed for DNA); --keep-dup all is REQUIRED because MACS3 default deduplicates and collapses MeRIP signal at high-coverage transcripts.

macs3 callpeak \
    --treatment aligned/IP_rep1.bam aligned/IP_rep2.bam aligned/IP_rep3.bam \
    --control   aligned/Input_rep1.bam aligned/Input_rep2.bam aligned/Input_rep3.bam \
    --format BAMPE \
    --gsize hs \
    --nomodel \
    --extsize 150 \
    --keep-dup all \
    --broad \
    --broad-cutoff 0.1 \
    --qvalue 0.05 \
    --outdir macs3_output \
    --name m6a_run1

--keep-dup all is non-negotiable for MeRIP; the default --keep-dup 1 (keep one read per position) destroys real signal. --broad and --broad-cutoff 0.1 extend narrow peaks into broader regions, appropriate for MeRIP fragments. Output: narrowPeak + broadPeak + gappedPeak files in macs3_output/.

DRACH Motif Confirmation (Sanity Check, NOT a Filter)

Goal: Confirm that the called peak set is enriched for the DRACH consensus motif relative to genomic background, validating antibody specificity. NEVER post-hoc remove individual non-DRACH peaks.

Approach: Run HOMER findMotifsGenome.pl with -rna mode on peak centres ±50 bp against a length-matched random-shuffled background; expect DRACH-like consensus in the top motifs with E-value < 1e-50 for a well-behaved MeRIP dataset.

findMotifsGenome.pl \
    exomepeak2_output/m6a_run1/peaks.bed \
    hg38 \
    motif_output \
    -rna \
    -size 100 \
    -len 5,6 \
    -p 8

library(Biostrings)
library(BSgenome.Hsapiens.UCSC.hg38)
library(ggseqlogo)
library(rtracklayer)

peaks <- import('exomepeak2_output/m6a_run1/peaks.bed')

peak_centres <- resize(peaks, width=5, fix='center')
genome <- BSgenome.Hsapiens.UCSC.hg38
seqs <- as.character(getSeq(genome, peak_centres))

ggseqlogo(seqs, method='probability') +
    ggplot2::labs(title='Peak-centre 5-mer (DRACH consensus expected)')

If DRACH enrichment is NOT observed in the peak set, the IP failed OR the wrong antibody was used OR the assay actually captured a different modification (m6A is centred in coding regions / 3'UTR; m1A is centred at TSS — different metagene signature). Do NOT proceed to differential analysis until DRACH is confirmed on the peak set as a whole.

Multi-Tool Consensus

Goal: Build a high-confidence peak set by intersecting multiple peak callers run on the same data; report both per-tool peak counts and the intersection.

Approach: Use bedtools intersect to combine exomePeak2, MeTPeak, and MACS3 broadPeak outputs; require at least 2-of-3 caller agreement for the high-confidence set.

bedtools intersect \
    -a exomepeak2_output/m6a_run1/peaks.bed \
    -b metpeak_output/m6a_run1/peak.bed macs3_output/m6a_run1_peaks.broadPeak \
    -wa \
    -u \
    -f 0.5 > consensus_at_least_2of3.bed

wc -l \
    exomepeak2_output/m6a_run1/peaks.bed \
    metpeak_output/m6a_run1/peak.bed \
    macs3_output/m6a_run1_peaks.broadPeak \
    consensus_at_least_2of3.bed

Intersection (consensus) is the conservative choice; union inflates the false-positive rate (every caller's idiosyncratic false positives propagate). For published analyses, report both per-tool count AND consensus count; the ratio is a useful caller-agreement metric.

Per-Method Failure Modes

Anti-m6A antibody cross-reacts with m6Am at 5' peaks

Trigger: Peak called at the very 5' end of a transcript (within ~50 nt of TSS), attributed to METTL3-deposited internal m6A.

Mechanism: Anti-m6A antibodies 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). METTL3 / METTL14 do NOT methylate the cap-adjacent position. Linder 2015 Nat Methods 12:767 explicitly noted the cross-reactivity; Mauer 2017 Nature 541:371 deepened the m6Am story.

Symptom: 5' peaks remain in METTL3-KO cells; agent infers METTL3-independent m6A writers when the peaks are really PCIF1 m6Am.

Fix: Flag peaks within ~50 nt of TSS as m6A-or-m6Am ambiguous. For unambiguous internal-m6A studies, exclude TSS-proximal peaks; for m6Am-specific studies, use PCIF1-KO subtraction or m6Am-seq (Sun H 2021); for orthogonal stoichiometric resolution, use GLORI / SAC-seq.

Peaks called at GC-extreme transcripts by MACS3 but not exomePeak2

Trigger: MACS3 broadPeak output has peaks at high-GC or low-GC transcripts that exomePeak2 does NOT call on the same BAMs.

Mechanism: exomePeak2 implements internal GC-bias correction; MACS3 does not. The IP step has GC bias because anti-m6A antibody pull-down efficiency varies with local GC content. MACS3 attributes GC-driven IP enrichment to true methylation; exomePeak2 corrects for it.

Fix: Trust exomePeak2 for transcript-aware analyses. Use MACS3 only as a second opinion or for broad / viral analyses. If a MACS3-only peak is biologically interesting, validate orthogonally.

DRACH filtering applied post-hoc

Trigger: Filtering called MeRIP peaks to retain only those containing a DRACH motif within the peak window.

Mechanism: ~70% of mammalian METTL3-deposited m6A sites sit within DRACH (Linder 2015), but a non-trivial fraction are non-DRACH. MeRIP peaks span ~100-200 nt; the methylated A is somewhere within the window but the peak boundary is not the modification position. Filtering by "no DRACH in window" rejects real peaks where DRACH sits at the edge or where the m6A is non-DRACH. The filter also creates circular validation — the filtered set is enriched for DRACH by construction.

Symptom: Peak count drops 30-50% after DRACH filtering; reviewers note the filter is not a standard convention.

Fix: Report DRACH enrichment as a SANITY CHECK on the peak set (HOMER E-value < 1e-50 expected) — never as a per-peak filter. For single-nucleotide methods (miCLIP, m6Anet, GLORI), the DRACH constraint is built into the calling model and need not be re-applied.

MACS3 default --keep-dup 1 on MeRIP

Trigger: MACS3 invoked without --keep-dup all; default behaviour is to retain one read per position.

Mechanism: MeRIP libraries have NO UMI (standard protocol); positional duplicates in MeRIP are a mix of PCR duplicates and biological re-sampling of high-coverage transcripts. Default --keep-dup 1 collapses the latter and destroys real signal at the most-abundant transcripts.

Symptom: MACS3 reports very few peaks at housekeeping genes (GAPDH, ACTB) despite obvious IP enrichment in IGV; peak count drops 5-20x compared to --keep-dup all.

Fix: Always pass --keep-dup all to MACS3 / MACS2 for MeRIP. This is the same logic as the do-NOT-dedup rule in merip-preprocessing.

MeTPeak called with wrong argument name for annotation

Trigger: metpeak(IP_BAM=..., INPUT_BAM=..., txdb=txdb_object) (lowercase txdb= rather than TXDB=).

Mechanism: MeTPeak's annotation arguments are GENE_ANNO_GTF= (GTF file path) or TXDB= (uppercase, TxDb object). Lowercase txdb= is not recognised. exomePeak2 uses lowercase txdb= instead, which leads to easy confusion between the two skills' APIs.

Fix: Pass either GENE_ANNO_GTF='annotation.gtf' (string path) or TXDB=txdb_object (capitalised). Verify argument casing against the installed MeTPeak source.

Single-tool peak set reported without cross-caller agreement

Trigger: Paper reports "18,234 m6A peaks" from exomePeak2 alone; downstream biology built on this single set.

Mechanism: Cross-tool concordance on MeRIP is empirically ~70% between exomePeak2 and MeTPeak on the same BAMs; single-tool reports include ~30% tool-specific calls. McIntyre 2020 Sci Rep 10:6590 documents this directly.

Fix: Run at least 2 callers (typically exomePeak2 + MeTPeak; or exomePeak2 + MACS3 broad); report the intersection as high-confidence; full per-tool counts as supplementary.