bio-sashimi-plots

Version Compatibility

Reference examples tested with: ggsashimi 1.1+, rmats2sashimiplot 3.0+, MAJIQ 3.0+, leafcutter 0.2.9+, pyGenomeTracks 3.8+, ggplot2 3.5+, pandas 2.2+

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 Nat Methods); 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.

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	MAJIQ build + psi/deltapsi	Splice-graph topology; LSV-aware; posterior violins	Static figures; non-academic license
leafviz	Cluster-level interactive browsing with NMD annotation	leafcutter differential output	Filter table + sashimi-like plots; NMD-aware	leafcutter-specific
Jutils	Unified output across rMATS, leafcutter, MntJULiP, MAJIQ	Tool-specific differential output	Heatmaps, Venn, sashimi tool-agnostically	Output less polished than ggsashimi
pyGenomeTracks	Multi-track publication figures (RNA-seq + ChIP/ATAC)	BigWig + BED + GTF	Combine RNA with chromatin tracks	Not splicing-specific; configure tracks manually
IGV (interactive)	Quick ad-hoc inspection	BAM + region	Scrollable, instant	Not for publication figures
MISO sashimi	Historical	MISO output	Original sashimi format	MISO unmaintained; no longer recommended

Decision Tree by Goal

Goal	Recommended tool
Validate a specific rMATS hit	rmats2sashimiplot (one-line) or ggsashimi (custom)
Validate a leafcutter cluster	leafviz (interactive) or ggsashimi with cluster coordinates
Validate a MAJIQ LSV (complex topology)	MAJIQ-VOILA (only tool that shows full LSV graph)
Publication-quality two-condition comparison	ggsashimi -O 3 -A mean_j for grouped overlay
Multi-track figure (RNA-seq + H3K4me3 + ATAC)	pyGenomeTracks
Quick ad-hoc browsing during development	IGV sashimi
Tool-agnostic batch heatmap of significant events	Jutils
Interactive cohort-level filtering of leafcutter results	leafviz Shiny

ggsashimi for Publication Overlays

Goal: Generate publication-quality sashimi plot for a region with samples grouped by condition and per-sample tracks aggregated.

Approach: Define samples + groups + colors in a TSV (no header), then call ggsashimi with coordinates, GTF, and visual flags.

import subprocess
import pandas as pd

# ggsashimi input: col1 = sample id, col2 = BAM path, col3 = group (for -O/-C overlay/color)
groups = pd.DataFrame({
    'sample_id': ['ctrl1', 'ctrl2', 'ctrl3', 'trt1', 'trt2', 'trt3'],
    'bam': ['ctrl1.bam', 'ctrl2.bam', 'ctrl3.bam', 'trt1.bam', 'trt2.bam', 'trt3.bam'],
    'group': ['Control', 'Control', 'Control', 'Treatment', 'Treatment', 'Treatment']
})
groups.to_csv('sashimi_groups.tsv', sep='\t', index=False, header=False)

subprocess.run([
    'ggsashimi.py',
    '-b', 'sashimi_groups.tsv',
    '-c', 'chr17:43094000-43125000',
    '-o', 'BRCA1_sashimi',
    '-M', '10',
    '--alpha', '0.25',
    '--height', '3',
    '--width', '10',
    '--shrink',
    '--fix-y-scale',
    '--ann-height', '4',
    '-g', 'gencode_v45.gtf',
    '--base-size', '14',
    '-O', '3',
    '-A', 'mean_j',
    '-F', 'pdf'
], check=True)

Key ggsashimi flags (Garrido-Martin 2018 PLoS Comput Biol):

• --overlay 3 (or -O 3): aggregate multiple samples within a group into a single overlay track with summary statistics — its signature feature
• -A mean_j: junction aggregation method (mean, median, mean_j accounts for sample-wise normalization); use mean_j for biological replicates
• --shrink: rescale long introns (>2x flanking exons) for compact display
• --fix-y-scale: identical y-axis across groups (essential for visual comparison)
• --alpha 0.25: transparency for per-sample coverage in overlay mode
• -M 10: minimum junction reads to display (lower = noisier; 5-10 typical; raise to 20+ for crowded plots)
• --ann-height: gene annotation track height
• -F pdf: output format (pdf, png, svg, eps)

Batch Plotting from rMATS Hits

Goal: Auto-generate sashimi plots for all significant rMATS differential events.

Approach: Parse SE.MATS.JC.txt, expand coordinates to flanking exons + 500nt context, iterate ggsashimi.

import subprocess
import pandas as pd
from pathlib import Path

diff = pd.read_csv('rmats_output/SE.MATS.JC.txt', sep='\t')
sig = diff[(diff['FDR'] < 0.05) & (diff['IncLevelDifference'].abs() > 0.10)]

Path('sashimi_plots').mkdir(exist_ok=True)
for idx, ev in sig.head(25).iterrows():
    region = f'{ev["chr"]}:{ev["upstreamES"] - 500}-{ev["downstreamEE"] + 500}'
    safe_name = f'{ev["geneSymbol"]}_{ev["chr"]}_{ev["upstreamES"]}'
    subprocess.run([
        'ggsashimi.py',
        '-b', 'sashimi_groups.tsv',
        '-c', region,
        '-o', f'sashimi_plots/{safe_name}',
        '-M', '5',
        '--shrink',
        '--fix-y-scale',
        '-O', '3',
        '-A', 'mean_j',
        '-g', 'annotation.gtf',
        '-F', 'pdf'
    ], check=True)

For MXE events, plot from upstreamES of exon 1 to downstreamEE of exon 2 to show both alternative exons in the same figure.

rmats2sashimiplot

Goal: Plot directly from rMATS event coordinates without manual region calculation.

Approach: Pass rMATS event file + BAM lists + event type; rmats2sashimiplot extracts coordinates and produces per-event PDFs.

rmats2sashimiplot \
    --b1 ctrl1.bam,ctrl2.bam,ctrl3.bam \
    --b2 trt1.bam,trt2.bam,trt3.bam \
    -t SE \
    -e rmats_output/SE.MATS.JC.txt \
    --l1 Control \
    --l2 Treatment \
    -o sashimi_rmats \
    --exon_s 1 \
    --intron_s 5 \
    --color '#1f77b4,#ff7f0e' \
    --group-info group_def.txt

--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.

# MAJIQ V3 (June 2025+) uses Zarr-format splicegraph (V2's .sql is deprecated)
voila view -p 5000 -j 8 build/splicegraph.zarr psi_output/sample.psi.voila -o voila_psi_html

voila view -p 5000 -j 8 build/splicegraph.zarr deltapsi_output/group1_group2.deltapsi.voila -o voila_dpsi_html

VOILA shows:

• 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.

prepare_results.R \
    -o leafviz \
    -m groups.txt \
    leafcutter_perind_numers.counts.gz \
    ds_results_cluster_significance.txt \
    ds_results_effect_sizes.txt \
    annotation_codes

library(leafviz)
run_leafviz('leafviz.RData')

Standalone alternative: jackhump/leafviz GitHub repo for the lightweight installable subset. Useful for cohort-level interactive filtering.

Jutils for Tool-Agnostic Output

Goal: Visualize differential splicing output uniformly across rMATS, leafcutter, MntJULiP, and MAJIQ.

Approach: Convert tool output to Jutils' standard format, then plot.

python3 jutils.py convert-results --rmats-dir rmats_output/ --out-dir jutils_out/
python3 jutils.py heatmap --tsv-file jutils_out/rmats.tsv --meta-file meta.tsv --q-value 0.05
python3 jutils.py sashimi --tsv-file jutils_out/rmats.tsv --meta-file meta.tsv \
    --gtf annotation.gtf --coordinate chr1:1000-2000 --bam-list bam_list.tsv
python3 jutils.py venn-diagram --tsv-file-list jutils_out/rmats.tsv,jutils_out/leafcutter.tsv

(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 = 3
title = GENCODE v45
fontsize = 10
file_type = gtf

[ctrl_coverage]
file = ctrl_merged.bw
title = Control
color = #1f77b4
height = 3
file_type = bigwig

[trt_coverage]
file = trt_merged.bw
title = Treatment
color = #ff7f0e
height = 3
file_type = bigwig

[junctions]
file = junctions.bedpe
title = Junctions
height = 2
file_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:

# Convert regtools junctions BED12 to BEDPE for pyGenomeTracks.
# regtools BED12 column 11 is blockSizes (anchor_left, anchor_right);
# column 12 is blockStarts (0, intron_length + anchor_left).
# Intron start = chromStart + anchor_left = $2 + a[1]
# Intron end   = chromStart + blockStarts[2] = $2 + b[2]
awk 'BEGIN{OFS="\t"} {split($11,a,","); split($12,b,","); s=$2+a[1]; e=$2+b[2]; print $1, s, s+1, $1, e-1, e, $5}' \
    regtools_junctions.bed > junctions.bedpe

pyGenomeTracks --tracks tracks.ini --region chr17:43094000-43125000 -o figure.pdf

Reading Sashimi Plots (Interpretation Guide)

Visual element	What it represents
Filled coverage track	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)
--fix-y-scale for cross-group comparisons	Otherwise auto-rescaling visually exaggerates differences
Aggregate replicates with -O 3 -A mean_j	Reduces clutter; per-sample variance still shown via alpha
Limit to 3-4 groups per figure	More becomes hard to read
Include 200-500 nt flanking exons	Show full splicing context
For MXE events, plot both alternative exons	Otherwise only half of the event is visible
Check accessibility colors	Use ColorBrewer-safe palettes for color-blind readers
Always include a legend	Sashimi figures without legends are uninformative for non-experts
Specify output format explicitly	PDF for publication; PNG for slides; SVG for editing

Common Errors

Error	Cause	Solution
ggsashimi: 'samtools' not found	samtools not in PATH	Install via conda; which samtools to verify
ggsashimi: empty plot	Region has no reads or wrong chromosome name	Check BAM with samtools view sample.bam chr1:100-200; chrom name match (chr1 vs 1)
rmats2sashimiplot: KeyError 'IJC_SAMPLE_1'	Old rmats2sashimiplot with new rMATS output	Update both to matching versions
voila: out of memory	Large LSV cohort	Filter by deltapsi threshold before voila
pyGenomeTracks: ini parse error	Missing closing bracket or invalid track type	Validate INI syntax; check pyGenomeTracks --listTracks for supported types
leafviz: missing exon file	annotation_codes path wrong	Re-run prepare_results.R with correct paths

Troubleshooting

Issue	Cause	Solution
No junctions shown	Default -M 10 too strict	Lower to -M 3 or -M 5
Plot too crowded	Many samples without aggregation	Use -O 3 to overlay groups
Annotation missing or wrong gene	GTF lacks gene_name attribute or wrong build	Verify GTF version vs BAM reference; pre-filter the GTF to the relevant features
Memory issues on large regions	>100 kb regions with many samples	Plot smaller windows or pre-extract reads with samtools view
Y-axis dominated by one peak	Outlier sample	Use -A mean_j to aggregate; or filter outlier

Related Skills

• differential-splicing - Identify events to plot; sashimi plots are validation
• splicing-quantification - Context for PSI values; sashimi provides visual confirmation
• data-visualization/genome-tracks - Multi-track figure design (pyGenomeTracks, Gviz)
• data-visualization/ggplot2-fundamentals - ggsashimi customization (extends ggplot2)
• data-visualization/color-palettes - Accessible color choices
• data-visualization/volcano-and-ma-plots - Volcano complement to sashimi
• data-visualization/heatmaps-clustering - Heatmap complement to sashimi

References

• Katz et al 2010 Nat Methods - MISO sashimi plot original
• Garrido-Martin et al 2018 PLoS Comput Biol - ggsashimi
• Yang et al 2021 Bioinformatics - Jutils
• Vaquero-Garcia et al 2016 eLife - MAJIQ / VOILA
• Li et al 2018 Nat Genet - leafcutter / leafviz
• Ramirez et al 2018 Nat Commun - pyGenomeTracks