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

Decision Tree by Goal

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)

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.