Creates sashimi-style plots showing RNA-seq read coverage and splice junction counts using ggsashimi (general-purpose, condition-grouped overlays), rmats2sashimiplot (rMATS-output-aware), MAJIQ-VOILA (LSV posteriors interactive HTML), leafviz (leafcutter clusters Shiny), Jutils (tool-agnostic heatmaps and sashimi for rMATS/leafcutter/MntJULiP/MAJIQ output), or pyGenomeTracks (multi-track publication figures). Tool choice depends on the upstream differential-splicing tool's output format and the publication vs interactive use case. Use when visualizing specific splicing events, validating differential splicing calls, or producing publication-quality figures.
Creates sashimi-style plots showing RNA-seq read coverage and splice junction counts using ggsashimi (general-purpose, condition-grouped overlays), rmats2sashimiplot (rMATS-output-aware), MAJIQ-VOILA (LSV posteriors interactive HTML), leafviz (leafcutter clusters Shiny), Jutils (tool-agnostic heatmaps and sashimi for rMATS/leafcutter/MntJULiP/MAJIQ output), or pyGenomeTracks (multi-track publication figures). Tool choice depends on the upstream differential-splicing tool's output format and the publication vs interactive use case. Use when visualizing specific splicing events, validating differential splicing calls, or producing publication-quality figures.
Before using code patterns, verify installed versions match. If versions differ:
Python: pip show <package> then help(module.function) to check signatures
CLI: <tool> --version then <tool> --help to confirm flags
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
Sashimi Plot Visualization
Visualize RNA-seq coverage tracks with splice junction arcs labeled by read count. Sashimi plots originated with MISO (Katz 2010 ); modern tools differ in input handling, group aggregation logic, and customization. Tool choice is not interchangeable — some tools work only with specific upstream output formats.
Nat Methods
Tool Selection Matrix
Tool
Best for
Input
Strengths
Fails when
ggsashimi
Publication-quality grouped overlays from any BAM
BAMs + region
--overlay aggregates samples within a group; clean PDFs
No native rMATS/MAJIQ integration; need to extract coords manually
rmats2sashimiplot
One-line plot from rMATS output
rMATS event file + BAMs
No manual coord extraction
rMATS-specific; doesn't handle leafcutter or MAJIQ
MAJIQ-VOILA
Interactive LSV browsing with posterior PSI distributions
--exon_s 1 --intron_s 5 shrinks intron-to-exon visual ratio 5:1 (introns drawn 1/5 their actual length). The --group-info flag (newer versions) allows custom replicate groupings.
MAJIQ-VOILA Interactive HTML
Goal: Browse LSV posterior PSI distributions interactively with splice-graph topology.
Approach: Run voila on MAJIQ output to generate self-contained HTML.
Complete LSV graphs (single source / single target nodes)
Per-junction posterior PSI violin plots
ΔPSI distributions across all conditions
Confidence by junction within an LSV
The only tool that visualizes complex multi-junction LSVs intuitively. For events that don't fit canonical SE/A5SS/A3SS, VOILA is the visualization of choice.
leafviz Shiny App
Goal: Browse leafcutter clusters with intron-level effects, sashimi-like plots, and NMD annotation.
Approach: Prepare leafviz input from leafcutter differential output, then launch Shiny.
(Yang 2021 Bioinformatics) Useful when comparing multiple tools' outputs across publications or doing meta-analysis.
pyGenomeTracks for Multi-Track Figures
Goal: Combine splicing with chromatin or coverage tracks for publication figures.
Approach: Define tracks in an INI file (genes, BAM, BigWig, BED), then run pyGenomeTracks --tracks tracks.ini --region ... -o figure.pdf.
[gene_models]file = annotation.gtf
height = 3title = GENCODE v45
fontsize = 10file_type = gtf
[ctrl_coverage]file = ctrl_merged.bw
title = Control
color = #1f77b4height = 3file_type = bigwig
[trt_coverage]file = trt_merged.bw
title = Treatment
color = #ff7f0eheight = 3file_type = bigwig
[junctions]file = junctions.bedpe
title = Junctions
height = 2file_type = links
links_type = arcs
The junctions.bedpe file must be in BEDPE format (6 columns: chr1 start1 end1 chr2 start2 end2 [+ optional score]). Convert from regtools .bed12 junctions:
Read coverage at each genomic position (depth-normalized in -A mode)
Arc / curve between exons
Junction-spanning reads; arc connects donor to acceptor
Number on arc
Count of junction-spanning reads (raw, not normalized, unless -A set)
Arc thickness
Often proportional to read count (tool-dependent)
Gene model below
Exons (boxes) and introns (lines) from GTF
Multiple parallel tracks
Per-sample (default) or per-group (with -O)
Junction count interpretation: the number on an arc is the absolute count of reads whose CIGAR string contained an N operation matching that intron coordinate. Higher = more usage. Compare counts on inclusion vs skipping arcs to estimate PSI visually.
Color convention: by convention, control = blue (#1f77b4), treatment = orange (#ff7f0e); always document. Use ColorBrewer or matplotlib defaults for >2 groups.
Per-Tool Failure Modes
ggsashimi: Off-Strand Junction Artifacts
Trigger: Stranded RNA-seq library plotted without strand specification.
Mechanism: ggsashimi reads BAM strand from CIGAR + flag; without strand info, antisense junctions appear as artifacts.
Symptom: Implausible junctions in regions with overlapping antisense genes; "noise" arcs at unexpected locations.
Fix: Set library strandedness with -s MATE2_SENSE (dUTP/TruSeq reverse-stranded PE; use -s MATE1_SENSE for forward, -s SENSE/ANTISENSE for single-end); verify orientation with RSeQC infer_experiment.py. Alternatively, pre-filter BAM by strand with samtools view -f 16 / -F 16.
rmats2sashimiplot: Wrong Coordinate Convention
Trigger: Older versions or non-default rMATS output.
Mechanism: rmats2sashimiplot expects 1-based coordinates from rMATS' .MATS.JC.txt; rMATS outputs 0-based half-open in some columns.
Symptom: Plot region shifted by 1 nt; arcs misaligned with gene model.
Fix: Verify rmats2sashimiplot version matches rMATS-turbo output convention; use ggsashimi for cleaner control.
MAJIQ-VOILA: Browser Memory
Trigger: Loading large VOILA HTML in browser (cohort with hundreds of LSVs).
Mechanism: VOILA HTML embeds all LSV data; large cohorts produce >100 MB HTMLs.
Symptom: Browser unresponsive on opening; "page unresponsive" warnings.
Fix: Filter LSVs in MAJIQ before voila step (--changing-pvalue-threshold 0.95 and --changing-between-group-dpsi-threshold 0.2); split into per-gene HTMLs.
leafviz: Annotation Codes Mismatch
Trigger: Using leafviz with annotation_codes from different GENCODE version than leafcutter clusters.
Mechanism: annotation_codes encodes intron-to-event-class mapping per GTF version.
Symptom: Many clusters show as "unannotated" despite being in canonical GTF.
Fix: Generate annotation_codes from the same GENCODE version used in differential analysis.
Customization Reference
Visual goal
ggsashimi flag
Reduce intron whitespace
--shrink
Identical y-axis across groups
--fix-y-scale
Per-group overlay aggregation
-O 3 -A mean_j
Larger figure
--width 12 --height 4
Bigger fonts
--base-size 16
Vector output
-F pdf or -F svg
Custom palette
Edit colors in groups TSV
Filter junction noise
-M 10 (raise to 20+)
Transparency
--alpha 0.25
Restrict to protein-coding
pre-filter the GTF (awk '$0 ~ /protein_coding/'); ggsashimi has no feature-filter flag
Best Practices
Tip
Rationale
Use --shrink for genes with large introns
Keeps exons visible (TTN, brain genes with multi-kb introns)
