bio-epitranscriptomics-modification-visualization

What this skill does

## Version Compatibility

Reference examples tested with: Guitar 2.18+ (Bioconductor), MetaPlotR (GitHub, unversioned), deepTools 3.5+, ggcoverage 1.4+, pyGenomeTracks 3.8+, Gviz 1.46+, ComplexHeatmap 2.18+, ggseqlogo 0.2+, ggplot2 3.5+, rtracklayer 1.62+, GenomicFeatures 1.54+, BSgenome.Hsapiens.UCSC.hg38 1.4+, TxDb.Hsapiens.UCSC.hg38.knownGene 3.18+, ChIPseeker 1.38+.

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('Guitar')` then `?GuitarPlot` to verify parameters
- CLI: `<tool> --help`; `deepTools <tool> --help`

If R throws `unused argument` or `argument is missing`, the API moved between Bioconductor minor releases. Guitar's older API used `txdb=` while newer uses `txTxdb=`; verify with `?GuitarPlot`. deepTools `bamCompare --operation log2` is the modern syntax (older `--ratio log2` is being phased out).

# RNA Modification Visualisation

**"Make the canonical m6A metagene plot for my paper"** -> Render the 5'UTR / CDS / 3'UTR transcript-feature metagene with stop-codon enrichment for visual confirmation of the canonical m6A topology (THE smoke test that the antibody-IP captured real m6A), the peak-feature-distribution stacked bar ("where do my peaks land?" Figure 1), peak-centred heatmaps clustered by condition, IP-over-Input paired browser tracks at specific loci, the DRACH sequence logo as a sanity check on antibody specificity, and the volcano / MA plots for differential modification. CRITICAL: the stop-codon enrichment plot is a biological-QC anchor — a MeRIP dataset that does NOT show enrichment at and around the stop codon indicates IP failure, wrong antibody, wrong protocol, or the assay captured a different modification (e.g., m1A is centred at TSS, not stop). Build this plot FIRST as a smoke test BEFORE any downstream visualisation.

- R: `Guitar::GuitarPlot(peakBedFiles, txTxdb=...)` -- canonical transcript-feature metagene (THE field-standard plot)
- CLI: `deeptools computeMatrix scale-regions ...` + `plotProfile` -- genome-coordinate metagene
- R: `ComplexHeatmap::Heatmap()` on peak-centred signal matrix -- peak-centred heatmap clustered by condition
- CLI: `pyGenomeTracks --tracks tracks.ini --region chr:start-end` -- multi-track browser plot
- R: `ggcoverage::ggcoverage(data=track, mark.region=peaks)` + `geom_gene()` -- ggplot2-based browser tracks
- R: `ggseqlogo::ggseqlogo(seqs)` -- DRACH sequence logo from peak-centre 5-mers
- R: `ChIPseeker::annotatePeak()` + `ggplot2` stacked bar -- 5'UTR / CDS / 3'UTR feature distribution

## The Single Most Important Modern Insight -- Stop-codon enrichment in the metagene plot is the biological QC anchor — a MeRIP dataset without it is suspect

Dominissini 2012 *Nature* 485:201 and Meyer 2012 *Cell* 149:1635 — concurrent papers from different labs using different protocols (MeRIP-seq and m6A-seq, respectively) — both independently showed m6A enrichment at and around the stop codon (3'UTR-proximal end of the CDS). This is the most robust biological signal in MeRIP-seq, reproduced across cell types, conditions, and decades. A metagene plot from a MeRIP library that does NOT show stop-codon enrichment indicates: (1) IP failure (antibody did not bind), (2) wrong antibody, (3) wrong protocol (e.g., the assay actually captured a different modification — m1A is centred at TSS, not stop), or (4) sample-RNA degradation. Build the Guitar metagene plot FIRST after peak calling, BEFORE any downstream visualisation, as a smoke test. The corollary: the 5'UTR / CDS / 3'UTR feature-distribution stacked bar (the "where do m6A peaks land?" Figure 1) should always be PAIRED with the metagene plot, because the stacked bar can look right (peaks land in 3'UTR / stop area) while the metagene is off (peak DENSITY not concentrated at the codon itself), or vice versa. The two plots are complementary biology-QC anchors, not redundant.

## Algorithmic Taxonomy

| Tool / plot | Mechanism | Inputs | Output | Strength | Fails when |
|-------------|-----------|--------|--------|----------|------------|
| Guitar GuitarPlot (Cui 2016 *Biomed Res Int* 2016:8367534) | Per-peak distance computation relative to TSS / start / stop / TES; feature-scaled rendering | BED + TxDb | PDF + per-feature density | THE field-standard m6A metagene; transcript-feature-aware | Older / newer API uses `txdb=` vs `txTxdb=` argument |
| MetaPlotR (Olarerin-George 2017 *Bioinformatics* 33:1563) | Alternative metagene approach with different segment-rescaling | BED + GTF | PDF | Alternative segment rescaling philosophy | GitHub-only; less actively maintained |
| deepTools computeMatrix + plotProfile (Ramírez 2016 *NAR* 44:W160) | Genome-coordinate signal aggregation over scaled gene regions | bigWig + BED | PDF + matrix | Fast; flexible region/anchor choice; well-tested ChIP-seq lineage | Not transcript-feature-aware; misses 5'UTR / CDS / 3'UTR distinction |
| deepTools plotHeatmap | Same matrix as plotProfile; rendered as heatmap with row clustering | bigWig matrix | PDF | Standard peak-centred heatmap | Single-condition view (use ComplexHeatmap for multi-condition clustering) |
| ComplexHeatmap (Gu 2016 *Bioinformatics* 32:2847) | General-purpose heatmap with multi-dimensional clustering | numeric matrix | PDF / interactive | Multi-condition cluster + annotation; publication-quality | Requires pre-computed signal matrix |
| pyGenomeTracks (Lopez-Delisle 2021 *Bioinformatics* 37:422) | Config-file (INI) driven multi-track browser plot | bigWig / bed / GTF + INI config | PDF / PNG | Reproducible browser figures via config; multi-track stacking | INI config syntax error-prone for new users |
| ggcoverage (Song & Wang 2023 *BMC Bioinformatics* 24:309) | ggplot2-native track plotting | bigWig / BAM + GTF | ggplot2 object | ggplot2 syntax; combinable with annotation layers | Newer tool; smaller user base than Gviz |
| Gviz (Hahne & Ivanek 2016 *Methods Mol Biol* 1418:335) | R/Bioconductor general genome track | bigWig / BAM / GRanges + TxDb | PDF / R plot | Most flexible R-native; long Bioconductor history | Heavier than ggcoverage; steeper learning curve |
| IGV / IGV.js (Robinson 2011 *Nat Biotechnol* 29:24) | Interactive browser via Java / JS | bigWig / BAM / bed | interactive | Standard for ad-hoc inspection | Not reproducible for figure generation |
| UCSC Track Hubs | UCSC genome browser display | bigWig + hub.txt + genomes.txt | URL | Public-display standard | Setup overhead for short projects |
| ggseqlogo (Wagih 2017 *Bioinformatics* 33:3645) | Sequence-logo rendering in ggplot2 | character vector of equal-length sequences | ggplot2 object | Native ggplot2; method='bits' or 'probability' | Requires equal-length sequences |
| MEME-ChIP (Machanick & Bailey 2011 *Bioinformatics* 27:1696) | Motif discovery + logo | BED + genome | HTML + PWM | Comprehensive motif suite; gold standard | Heavier than ggseqlogo for simple visualisation |
| HOMER findMotifsGenome.pl (Heinz 2010 *Mol Cell* 38:576) | Motif discovery via cumulative hypergeometric on shuffled background | BED + genome | HTML + motif files | Battle-tested; RNA mode via `-rna` | Less flexible output formatting than MEME |

## Decision Tree by Scenario

| Scenario | Recommended | Why wrong choices fail |
|----------|-------------|------------------------|
| Canonical Figure 1 metagene plot for m6A paper | Guitar GuitarPlot with TxDb -- transcript-feature scaled | deepTools computeMatrix is genome-coordinate only; misses 5'UTR / CDS / 3'UTR distinction |
| QC: does my MeRIP show stop-codon enrichment? | Guitar GuitarPlot FIRST -- the canonical biological-QC anchor | Skipping this and rushing to peak counting is the most common QC failure |
| Browser-track Figure for a specific locus | pyGenomeTracks (config-file driven; reproducible) OR ggcoverage (ggplot2-native) | IGV is interactive but not reproducible for figures |
| Multi-condition peak-centred heatmap with clustering | ComplexHeatmap; row-cluster k-means; annotate by condition |