Bio Epitranscriptomics Modification Visualization

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

Decision Tree by Scenario

Methodology evolves; before high-stakes figure generation, web-search "Guitar Bioconductor release notes" and "pyGenomeTracks vs ggcoverage" for current best practice.

Guitar Transcript-Feature Metagene (THE Canonical Plot)

Goal: Render the canonical m6A metagene plot showing peak density along scaled transcript features (5'UTR, CDS, 3'UTR), with the stop-codon-proximal enrichment that is the field's biological QC anchor.

Approach: Load called peaks from BED; pass to GuitarPlot with the matched TxDb; export per-feature density and the PDF.

library(Guitar)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(rtracklayer)

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

# Older Guitar versions: argument is `txdb=`; newer: `txTxdb=`. Verify with ?GuitarPlot.
GuitarPlot(
    txTxdb           = txdb,
    stBedFiles       = list('exomepeak2_output/m6a_run1/peaks.bed'),
    miscOutFilePrefix = 'figures/m6a_metagene',
    enableCI         = FALSE,
    saveToPDFprefix  = 'figures/m6a_metagene'
)

If txTxdb= is rejected, fall back to txdb= and consult ?GuitarPlot. The expected output: stop-codon-proximal enrichment, with peak density rising toward and peaking near the stop codon in the 3'UTR. If the plot does NOT show this pattern, do not proceed to downstream visualisation — investigate IP / antibody / protocol failure in merip-preprocessing.

deepTools Genome-Coordinate Metagene

Goal: Render a genome-coordinate metagene over scaled gene regions using deepTools, useful when transcript-feature scaling is not needed or when the regions of interest are not standard genes (custom BED).

Approach: Build a signal matrix with computeMatrix scale-regions from a bigWig (typically IP-over-Input log2 from merip-preprocessing); render as profile plot + heatmap.

mkdir -p figures

computeMatrix scale-regions \
    --regionsFileName refs/protein_coding.bed \
    --scoreFileName tracks/IP_rep1_over_Input_rep1.bw \
    --regionBodyLength 2000 \
    --upstream 500 \
    --downstream 500 \
    --skipZeros \
    --numberOfProcessors 8 \
    --outFileName figures/m6a_matrix.gz

plotProfile \
    --matrixFile figures/m6a_matrix.gz \
    --outFileName figures/m6a_profile.pdf \
    --plotTitle 'm6A IP over Input metagene' \
    --plotType lines

plotHeatmap \
    --matrixFile figures/m6a_matrix.gz \
    --outFileName figures/m6a_heatmap.pdf \
    --colorMap RdBu_r \
    --plotTitle 'm6A IP over Input heatmap'

scale-regions is the right mode for gene-body metagene (5'-end-to-3'-end scaled to common length); reference-point is for peak-centred plots. For peak-centred, see the next section.

Peak-Centred Heatmap

Goal: Render a peak-centred heatmap of IP-over-Input signal at +/-window around each peak, clustered by condition or by signal pattern, for cross-condition comparison.

Approach: Use deepTools computeMatrix reference-point --referencePoint center for the matrix; ComplexHeatmap for the clustered render.

computeMatrix reference-point \
    --regionsFileName exomepeak2_output/m6a_run1/peaks.bed \
    --scoreFileName tracks/IP_rep1_over_Input_rep1.bw tracks/IP_rep2_over_Input_rep2.bw tracks/IP_rep3_over_Input_rep3.bw \
    --referencePoint center \
    --upstream 500 \
    --downstream 500 \
    --binSize 25 \
    --skipZeros \
    --numberOfProcessors 8 \
    --outFileName figures/peak_centred_matrix.gz \
    --outFileNameMatrix figures/peak_centred_matrix.tab

plotHeatmap \
    --matrixFile figures/peak_centred_matrix.gz \
    --outFileName figures/peak_centred_heatmap.pdf \
    --kmeans 3 \
    --colorMap viridis \
    --plotTitle 'm6A signal centred at peaks'

For multi-condition heatmap with annotations, parse the matrix into R and use ComplexHeatmap:

library(ComplexHeatmap)
library(circlize)

# deepTools --outFileNameMatrix has a 3-line JSON-style header before per-bin numeric columns.
raw <- read.delim('figures/peak_centred_matrix.tab', skip=3, header=FALSE)
mat <- as.matrix(raw[, -(1:6)])

col_fun <- colorRamp2(c(-2, 0, 2), c('blue', 'white', 'red'))

Heatmap(
    mat,
    name              = 'log2 IP / Input',
    col               = col_fun,
    cluster_columns   = FALSE,
    cluster_rows      = TRUE,
    row_km            = 3,
    show_row_names    = FALSE,
    show_column_names = FALSE,
    column_title      = 'Peak-centred (+/- 500 bp)'
)

IP-vs-Input Paired Browser Tracks via pyGenomeTracks

Goal: Render publication-quality genome-browser tracks at a specific locus showing paired IP / Input bigWig tracks, the IP/input log2 ratio, peak calls, and gene annotation; reproducible via INI config.

Approach: Build a tracks.ini config file listing each track type (bigwig, bed, gtf); invoke pyGenomeTracks --tracks tracks.ini --region chr:start-end.

[x-axis]
where = top
fontsize = 12

[IP rep1]
file = tracks/IP_rep1.bw
title = IP rep1
color = #d62728
height = 2
min_value = 0

[Input rep1]
file = tracks/Input_rep1.bw
title = Input rep1
color = #1f77b4
height = 2
min_value = 0

[IP / Input log2]
file = tracks/IP_rep1_over_Input_rep1.bw
title = log2 IP / Input
color = #2ca02c
height = 2

[spacer]

[m6A peaks]
file = exomepeak2_output/m6a_run1/peaks.bed
title = m6A peaks (exomePeak2)
color = #ff7f0e
height = 1
labels = false

[genes]
file = refs/annotation.gtf
title = GENCODE genes
height = 2
prefered_name = gene_name
merge_transcripts = true

pyGenomeTracks \
    --tracks tracks.ini \
    --region chr19:54,792,000-54,799,000 \
    --outFileName figures/browser_metti3_locus.pdf \
    --width 14 \
    --plotWidth 12

ggcoverage ggplot2-Native Browser Track

Goal: Render genome-browser tracks in ggplot2 syntax for combining with other ggplot2 layers (annotations, peak highlights, custom theming).

Approach: LoadTrackFile() parses a bigWig / bigBed / BAM input into the dataframe ggcoverage expects; chain ggcoverage() with geom_gene() for transcript annotation; mark.region requires columns start, end, and label.

library(ggcoverage)
library(rtracklayer)

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

track.df <- LoadTrackFile(
    track.file = 'tracks/IP_rep1_over_Input_rep1.bw',
    format     = 'bw',
    region     = 'chr19:54792000-54799000'
)

mark.df <- data.frame(
    start = peaks$start,
    end   = peaks$end,
    label = peaks$name
)

ggcoverage(
    data        = track.df,
    color       = 'auto',
    mark.region = mark.df
) +
    geom_gene(gtf.file = 'refs/annotation.gtf') +
    ggplot2::theme_classic()

DRACH Sequence Logo

Goal: Render a sequence logo of peak-centre 5-mers as a sanity check that the called peak set is enriched for the DRACH consensus motif.

Approach: Resize peak GRanges to fixed 5-nt centred windows; extract genomic sequences; pass to ggseqlogo with method='probability'.

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)',
                  subtitle='Sanity check on antibody specificity — NOT a per-peak filter')

The expected output: a logo showing approximately D-R-A-C-H consensus (D=A/G/U, R=A/G, A=methylated, C, H=A/C/U) with A clearly dominating position 3. If the logo does NOT show DRACH-like enrichment, the IP failed OR the wrong antibody was used.

5'UTR / CDS / 3'UTR Stacked Bar

Goal: Render the Figure 1 "where do my peaks land?" stacked bar showing the fraction of peaks in 5'UTR vs CDS vs 3'UTR vs intron, paired with the metagene to confirm canonical m6A topology.

Approach: Use ChIPseeker annotatePeak() with a matched TxDb; aggregate to per-feature counts; render as ggplot2 stacked bar.

library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(rtracklayer)
library(ggplot2)
library(dplyr)

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
peaks <- import('exomepeak2_output/m6a_run1/peaks.bed')

peak_anno <- annotatePeak(peaks, TxDb=txdb, level='transcript', verbose=FALSE)
anno_df <- as.data.frame(peak_anno@anno)

anno_df$feature <- gsub(' \\(.*\\)', '', anno_df$annotation)

feature_counts <- anno_df %>%
    group_by(feature) %>%
    summarise(n=n()) %>%
    mutate(fraction=n/sum(n)) %>%
    arrange(desc(fraction))

ggplot(feature_counts, aes(x='m6A peaks', y=fraction, fill=feature)) +
    geom_col(width=0.5) +
    scale_y_continuous(labels=scales::percent_format()) +
    labs(x=NULL, y='Fraction of peaks', title='m6A peak distribution across transcript features') +
    theme_minimal()

Per-Method Failure Modes

Metagene without stop-codon enrichment

Trigger: Guitar GuitarPlot of called peaks does NOT show enrichment at and around the stop codon.

Mechanism: Canonical m6A topology (Dominissini 2012 / Meyer 2012) shows stop-codon-proximal enrichment. Absence indicates (1) IP failure, (2) wrong antibody, (3) wrong protocol, (4) sample-RNA degradation, OR (5) the assay captured a different modification (m1A is TSS-centred; m5C distribution differs).

Symptom: Metagene is flat OR peaks at TSS (m1A signature) OR peaks in introns (unusual).

Fix: Do NOT proceed to downstream analysis. Diagnose at the merip-preprocessing layer: plotFingerprint, per-transcript IP/input distribution, antibody-lot QC. Re-do the IP if necessary. The metagene plot is the biological QC anchor; without it, all downstream interpretation is suspect.

DRACH logo not enriched

Trigger: ggseqlogo of peak-centre 5-mers shows no consensus, OR shows a non-DRACH-like motif.

Mechanism: Antibody specificity failure OR wrong protocol. Anti-m6A antibodies have ~70% DRACH-context enrichment on real m6A peaks; a failed IP captures random sequences with no consensus.

Fix: Re-inspect IP enrichment in merip-preprocessing. If IP is clean but DRACH logo is absent, the assay may have captured a different modification — investigate before claiming m6A.

Guitar txdb= vs txTxdb= argument confusion

Trigger: GuitarPlot(stBedFiles=..., txdb=txdb) rejected with "unused argument" error.

Mechanism: Guitar changed the argument name between Bioconductor releases — txdb= (older) vs txTxdb= (newer).

Fix: Try alternative; consult ?GuitarPlot for installed version. Pin Guitar version in reproducible analyses.

deepTools genome-coordinate metagene used where transcript-feature is needed

Trigger: computeMatrix scale-regions over a BED of genes used to show "5'UTR / CDS / 3'UTR enrichment".

Mechanism: deepTools scales by genomic length, NOT by transcript-feature length. A gene with long 5'UTR and short CDS will scale 5'UTR more than CDS; the metagene loses 5'UTR / CDS / 3'UTR semantics.

Fix: Use Guitar (transcript-feature-aware) for 5'UTR / CDS / 3'UTR semantics. Use deepTools for genome-coordinate metagene over arbitrary BED.