Visualises RNA-modification data with transcript-feature metagene plots (Guitar GuitarPlot; MetaPlotR; deepTools computeMatrix scale-regions), peak-centred heatmaps (ComplexHeatmap; deepTools plotHeatmap), IP-vs-input paired browser tracks (log2 IP/input bigWig via deepTools bamCompare; pyGenomeTracks; Gviz; IGV/UCSC track hubs), DRACH sequence-logo plots (ggseqlogo; MEME), feature-distribution stacked bars, and volcano/MA plots for differential modification. Establishes stop-codon enrichment in the metagene plot as the biological QC anchor for any MeRIP dataset (Dominissini 2012; Meyer 2012). Use when producing the canonical metagene plot with stop-codon enrichment as a QC anchor, building paired IP/input genome-browser tracks at single-locus resolution, plotting peak-centred heatmaps clustered by condition, summarising peak distribution across transcript features, generating DRACH motif logos as sanity checks, rendering volcano plots of differential m6A, or reproducing the stop-codon enrichment plot.
Visualises RNA-modification data with transcript-feature metagene plots (Guitar GuitarPlot; MetaPlotR; deepTools computeMatrix scale-regions), peak-centred heatmaps (ComplexHeatmap; deepTools plotHeatmap), IP-vs-input paired browser tracks (log2 IP/input bigWig via deepTools bamCompare; pyGenomeTracks; Gviz; IGV/UCSC track hubs), DRACH sequence-logo plots (ggseqlogo; MEME), feature-distribution stacked bars, and volcano/MA plots for differential modification. Establishes stop-codon enrichment in the metagene plot as the biological QC anchor for any MeRIP dataset (Dominissini 2012; Meyer 2012). Use when producing the canonical metagene plot with stop-codon enrichment as a QC anchor, building paired IP/input genome-browser tracks at single-locus resolution, plotting peak-centred heatmaps clustered by condition, summarising peak distribution across transcript features, generating DRACH motif logos as sanity checks, rendering volcano plots of differential m6A, or reproducing the stop-codon enrichment plot.
Before using code patterns, verify installed versions match. If versions differ:
R: packageVersion('Guitar') then to verify parameters
?GuitarPlot
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)
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)
Other tools produce similar but not identical plots
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.
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.
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.
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 = #d62728height = 2min_value = 0[Input rep1]file = tracks/Input_rep1.bw
title = Input rep1
color = #1f77b4height = 2min_value = 0[IP / Input log2]file = tracks/IP_rep1_over_Input_rep1.bw
title = log2 IP / Input
color = #2ca02cheight = 2[spacer][m6A peaks]file = exomepeak2_output/m6a_run1/peaks.bed
title = m6A peaks (exomePeak2)
color = #ff7f0eheight = 1labels = false[genes]file = refs/annotation.gtf
title = GENCODE genes
height = 2prefered_name = gene_name
merge_transcripts = true
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.
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.
IGV-only browser figure in published paper
Trigger: Figure 4 of paper shows an IGV screenshot of one locus; not reproducible from code.
Mechanism: IGV is interactive; the figure cannot be re-generated from a config file. Reviewers cannot reproduce the figure if track files change.
Fix: Use pyGenomeTracks (INI config) or ggcoverage (ggplot2 script) for figures intended for publication. IGV is for ad-hoc inspection only.
Fix: Resize peak ranges to fixed width before extraction: resize(peaks, width=5, fix='center').
pyGenomeTracks INI typo
Trigger: Track section header missing brackets, OR file = instead of file= (whitespace inconsistency).
Fix: Use the documented INI syntax precisely; bracketed section headers; consistent whitespace around =. Validate with pyGenomeTracks --tracks tracks.ini --region 'chr:start-end' --outFileName out.pdf and iterate on errors.
Reconciliation: When Plots Disagree
Pattern
Likely cause
Action
Metagene shows stop-codon enrichment but stacked bar has many 5'UTR peaks
5'UTR peaks are present but density is concentrated at stop codon
Both plots are correct; report together
Guitar metagene and deepTools metagene look very different
Guitar transcript-feature-scaled vs deepTools genome-scaled
Both correct; Guitar for 5'UTR / CDS / 3'UTR semantics; deepTools for genome coordinates
Peak-centred heatmap shows two clusters; condition annotation crosses cluster boundary
Biological signal not condition-aligned; OR clustering driven by per-peak coverage variance
Inspect per-cluster fold-change; consider removing low-coverage peaks before re-clustering
pyGenomeTracks shows IP track higher than Input but bamCompare bigWig shows log2 near zero
Track signal magnitudes are CPM-normalised; absolute counts and log2 ratios are different summaries
Report log2 ratio as the primary; per-track CPM as supplement
DRACH logo enriched but stacked bar shows mostly intronic peaks
Intronic m6A (Louloupi 2018 nascent transcripts)
Genuine biology; report intronic vs exonic separately; consider library prep (poly-A vs ribo-depleted)
Volcano plot symmetric but most differential peaks in 3'UTR
Differential is feature-restricted
Cross-check with per-feature differential testing; biological interpretation
Quantitative Thresholds
Quantity
Threshold
Source / rationale
Guitar metagene -- expected pattern
Stop-codon-proximal enrichment, 3'UTR > CDS > 5'UTR density
Guitar is transcript-feature-scaled (5'UTR / CDS / 3'UTR semantics); deepTools is genome-coordinate; both reported when needed
"Is the DRACH motif enriched?"
ggseqlogo of peak-centre 5-mers; OR HOMER findMotifsGenome.pl E-value reported
"Is the browser figure reproducible?"
Yes — pyGenomeTracks INI config OR ggcoverage R script; not IGV screenshot
"How were peaks annotated to features?"
ChIPseeker annotatePeak with matched TxDb; hierarchy explicit
"Why these specific clusters in the heatmap?"
k-means with k=3 chosen via elbow / silhouette; clusters reflect signal heterogeneity
"Does the metagene differ between conditions?"
Per-condition metagenes plotted alongside; differences quantified at the feature level
"Why is the colour scheme red-blue?"
Standard convention for log2 ratios (red = up, blue = down); colour-blind-safe palette via viridis available
"Was a cross-check with published m6A-Atlas peaks done?"
Common-core overlap reported; cited m6A-Atlas v2
"Are the browser track signal magnitudes comparable across samples?"
bigWig CPM-normalised in merip-preprocessing; documented
References
Dominissini D, Moshitch-Moshkovitz S, Schwartz S et al (2012) Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature 485(7397):201-206. doi:10.1038/nature11112
Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR (2012) Comprehensive analysis of mRNA methylation reveals enrichment in 3' UTRs and near stop codons. Cell 149(7):1635-1646. doi:10.1016/j.cell.2012.05.003
Cui X, Wei Z, Zhang L et al (2016) Guitar: an R/Bioconductor package for gene annotation guided transcriptomic analysis of RNA-related genomic features. Biomed Res Int 2016:8367534. doi:10.1155/2016/8367534
Olarerin-George AO, Jaffrey SR (2017) MetaPlotR: a Perl/R pipeline for plotting metagenes of nucleotide modifications and other transcriptomic sites. Bioinformatics 33(10):1563-1564. doi:10.1093/bioinformatics/btx002
Ramírez F, Ryan DP, Grüning B et al (2016) deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res 44(W1):W160-W165. doi:10.1093/nar/gkw257
Gu Z, Eils R, Schlesner M (2016) Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 32(18):2847-2849. doi:10.1093/bioinformatics/btw313
Lopez-Delisle L, Rabbani L, Wolff J et al (2021) pyGenomeTracks: reproducible plots for multivariate genomic datasets. Bioinformatics 37(3):422-423. doi:10.1093/bioinformatics/btaa692
Song Y, Wang J (2023) ggcoverage: an R package to visualize and annotate genome coverage for various NGS data. BMC Bioinformatics 24(1):309. doi:10.1186/s12859-023-05438-2
Robinson JT, Thorvaldsdóttir H, Winckler W et al (2011) Integrative genomics viewer. Nat Biotechnol 29(1):24-26. doi:10.1038/nbt.1754
Wagih O (2017) ggseqlogo: a versatile R package for drawing sequence logos. Bioinformatics 33(22):3645-3647. doi:10.1093/bioinformatics/btx469
Yu G, Wang LG, He QY (2015) ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics 31(14):2382-2383. doi:10.1093/bioinformatics/btv145
Hahne F, Ivanek R (2016) Visualizing Genomic Data Using Gviz and Bioconductor. Methods Mol Biol 1418:335-351. doi:10.1007/978-1-4939-3578-9_16
Heinz S, Benner C, Spann N et al (2010) Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell 38(4):576-589. doi:10.1016/j.molcel.2010.05.004
Machanick P, Bailey TL (2011) MEME-ChIP: motif analysis of large DNA datasets. Bioinformatics 27(12):1696-1697. doi:10.1093/bioinformatics/btr189
Related Skills
merip-preprocessing - Generates the bigWig tracks (bamCompare log2 IP/input) used here
m6a-peak-calling - Generates the peak BED used for metagene + stacked bar + heatmap
m6a-differential - Differential results visualised via volcano + MA (those plots live in m6a-differential, not duplicated here)
m6anet-analysis - Per-site DRS modification calls; visualisation analogous via metagene
data-visualization/ggplot2-fundamentals - General ggplot2 grammar
data-visualization/multipanel-figures - Combining metagene + heatmap + volcano into figures
data-visualization/heatmaps-clustering - General heatmap clustering patterns
data-visualization/volcano-and-ma-plots - General volcano / MA recipes (modification-specific volcano lives in m6a-differential)
data-visualization/genome-tracks - General genome-track rendering (this skill adds the IP/input pairing specifics)
data-visualization/sequence-logos - General sequence-logo plotting (this skill adds DRACH-specific context)
chip-seq/chipseq-visualization - Closest sibling for browser-track + peak-centred heatmap patterns
pathway-analysis/enrichment-visualization - For visualising downstream pathway results
