Bio Differential Splicing

Version Compatibility

Reference examples tested with: rMATS-turbo 4.3+, SUPPA2 2.4+, leafcutter 0.2.9+, MAJIQ 3.0+, Shiba 0.5+, STAR 2.7.11+, regtools 1.0+, pandas 2.2+, R 4.4+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters
• 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.

Differential Splicing

Detect splicing changes between conditions. Tool choice is a decision about statistical model, annotation dependence, and calibration regime under the specific experimental design — not a preference. Wrong tool for the design produces uncalibrated FDR or systematic effect-size bias.

Statistical Model Taxonomy

The first decision is which regime the design falls into: between-group with replicates, heterogeneous cohort, or single-sample-vs-cohort. Within each regime, tool choice is much smaller (1-2 options).

Comprehensive 2023-2026 benchmarks: Olofsson 2023 Biochem Biophys Res Commun; Tran 2025 WIREs RNA; Kubota 2025 NAR. Methodology evolves — verify benchmarks and tool docs before reporting. Default 2026 recommendation: run two complementary tools (rMATS + leafcutter) and require concordance for high-confidence calls.

Decision Tree by Experimental Design

rMATS-turbo Differential Analysis

Goal: Detect statistically significant differential splicing between two groups from BAMs.

Approach: Run rMATS-turbo without --statoff, then filter by FDR + ΔPSI + per-replicate coverage.

rmats.py \
    --b1 condition1_bams.txt \
    --b2 condition2_bams.txt \
    --gtf annotation.gtf \
    -t paired \
    --readLength 150 \
    --variable-read-length \
    --libType fr-firststrand \
    --nthread 8 \
    --od rmats_output \
    --tmp rmats_tmp \
    --novelSS \
    --cstat 0.05

--cstat 0.05 tests |ΔPSI| > 0.05; raise to 0.10 for stricter discovery. --novelSS enables novel-junction discovery (recommended with STAR 2-pass). For paired designs, add --paired-stats.

import pandas as pd
import numpy as np

se = pd.read_csv('rmats_output/SE.MATS.JC.txt', sep='\t')

def min_per_rep(s):
    return s.str.split(',').apply(lambda x: min(int(v) for v in x))

se['min_inc'] = min_per_rep(se['IJC_SAMPLE_1']).combine(min_per_rep(se['IJC_SAMPLE_2']), min)
se['min_skip'] = min_per_rep(se['SJC_SAMPLE_1']).combine(min_per_rep(se['SJC_SAMPLE_2']), min)

significant = se[
    (se['FDR'] < 0.05) &
    (se['IncLevelDifference'].abs() > 0.10) &
    ((se['min_inc'] + se['min_skip']) >= 10)
].copy()

significant['score'] = -np.log10(significant['FDR']) * significant['IncLevelDifference'].abs()
top = significant.nlargest(50, 'score')

leafcutter Differential Intron Usage

Goal: Detect differential intron-cluster usage annotation-free, capturing novel junctions and complex multi-junction events.

Approach: Extract junctions with regtools, cluster introns by shared splice sites, run cluster-level Dirichlet-multinomial test.

for bam in *.bam; do
    regtools junctions extract -a 8 -m 50 -s XS "$bam" -o "${bam%.bam}.junc"
done
ls *.junc > juncfiles.txt

python leafcutter_cluster_regtools.py \
    -j juncfiles.txt \
    -o leafcutter \
    -m 50 \
    -l 500000

library(leafcutter)

groups <- data.frame(
    sample = c('s1', 's2', 's3', 's4', 's5', 's6'),
    group = c('control', 'control', 'control', 'treatment', 'treatment', 'treatment')
)
write.table(groups, 'groups.txt', sep = '\t', quote = FALSE, row.names = FALSE, col.names = FALSE)

system('leafcutter_ds.R --num_threads 4 --exon_file gencode_exons.txt.gz \
    leafcutter_perind_numers.counts.gz groups.txt -o ds_results')

cluster_sig <- read.table('ds_results_cluster_significance.txt', header = TRUE, sep = '\t')
intron_effects <- read.table('ds_results_effect_sizes.txt', header = TRUE, sep = '\t')

sig_clusters <- subset(cluster_sig, p.adjust < 0.05)

LeafCutter2 (Buen Abad Najar 2025 bioRxiv) extends leafcutter with NMD-aware classification of unproductive splicing — useful when AS-NMD coupling is the question.

MAJIQ V3 Differential Analysis

Goal: Detect differential LSVs with full posterior distributions over ΔPSI; ideal for complex multi-junction events and heterogeneous cohorts.

Approach: Build splice graph -> compute coverage per group -> run deltapsi (replicate-structured) or heterogen (cohort-style).

majiq build annotation.gff3 -c settings.ini -j 8 -o build_output

majiq deltapsi \
    -grp1 build_output/ctrl1.majiq build_output/ctrl2.majiq build_output/ctrl3.majiq \
    -grp2 build_output/trt1.majiq build_output/trt2.majiq build_output/trt3.majiq \
    -n control treatment \
    -o deltapsi_output \
    --minreads 10 --minpos 3 \
    -j 8

majiq heterogen \
    -grp1 build_output/het_ctrl{1..20}.majiq \
    -grp2 build_output/het_trt{1..20}.majiq \
    -n control treatment \
    -o heterogen_output \
    -j 8

voila view -p 5000 -j 8 build_output/splicegraph.zarr deltapsi_output/control_treatment.deltapsi.voila -o voila_html

MAJIQ V3 (Aicher, Slaff, Jewell, Barash bioRxiv 2024; public release 2025) uses Zarr storage (splicegraph.zarr); V2's SQLite splicegraph is deprecated. MAJIQ reports posterior probability P(|ΔPSI| > 0.2); thresholds are interpreted differently from FDR. Use HET for n>=10 vs n>=10 cohort designs (clinical, GTEx-style); deltapsi for tightly controlled n=3 vs n=3.

SUPPA2 Differential Analysis

Goal: Quick differential splicing from existing transcript quantifications, useful as a sanity check or pilot.

Approach: Generate per-condition PSI files from Salmon TPM, then run diffSplice with empirical or classical p-values.

suppa.py generateEvents -i annotation.gtf -o events -f ioe -e SE SS MX RI

for ev in SE A5 A3 MX RI; do
    suppa.py psiPerEvent -i events_${ev}_strict.ioe -e ctrl_tpm.tsv -o ctrl_${ev}
    suppa.py psiPerEvent -i events_${ev}_strict.ioe -e trt_tpm.tsv -o trt_${ev}

    suppa.py diffSplice \
        -m empirical \
        -gc \
        -i events_${ev}_strict.ioe \
        -p ctrl_${ev}.psi trt_${ev}.psi \
        -e ctrl_tpm.tsv trt_tpm.tsv \
        -o diff_${ev}
done

For n<=3 designs, switch -m classical (Wilcoxon). Empirical null requires sufficient between-replicate observations to construct.

Shiba for Low-Coverage / Few-Replicate Designs

Goal: Detect differential splicing with explicit junction-imbalance correction — addresses a known false-positive source for rMATS-style methods.

Approach: Shiba is a Snakemake-based pipeline configured via YAML. Install via bioconda, write a config file describing groups + BAMs, then run with snakemake.

conda install -c bioconda shiba

# Edit config.yaml with reference GTF, BAM groups, output dir, thresholds
# Then run the Snakemake workflow:
snakemake -s snakeshiba.smk \
    --configfile config.yaml \
    --cores 8 \
    --use-singularity \
    --singularity-args "--bind $HOME:$HOME"

Shiba (Kubota 2025 NAR) reportedly outperforms rMATS at n=2 vs n=2 by correcting differential mappability between inclusion and skipping junctions; community calibration still emerging. See https://sika-zheng-lab.github.io/Shiba/ for the full config.yaml schema.

Per-Tool Failure Modes

rMATS: Confounder-Blind LRT

Trigger: Sequencing batch, RIN, library prep date, or sex correlates with the comparison of interest.

Mechanism: rMATS' default LRT does not natively accept covariates the way DESeq2 does; the --paired-stats flag handles paired designs but not arbitrary covariates.

Symptom: Many "significant" hits driven by batch rather than biological condition; PCA on PSI matrix shows samples clustering by batch rather than group.

Fix: Either (a) include batch as a stratification (run rMATS within each batch), (b) regress PSI matrix against batch in R, then test residuals, or (c) switch to leafcutter which accepts a confounders argument in differential_splicing.

leafcutter: Cluster Mis-Topology

Trigger: A cluster spans a complex topology (cassette + alternative donor in same cluster).

Mechanism: Cluster-level p-value reports "something in this cluster differs" but doesn't indicate which intron drove the change; downstream analysis needs per-intron effect sizes.

Symptom: Significant cluster, multiple introns with different ΔPSI directions, ambiguous biological interpretation.

Fix: Inspect cluster in leafviz; report per-intron effect sizes from ds_results_effect_sizes.txt; map cluster topology to canonical SE/A5SS/A3SS via flanking exon coordinates manually.

MAJIQ HET: Power vs Type-1 Tradeoff

Trigger: HET module on n=5-10 cohorts (between regimes).

Mechanism: HET assumes between-sample variability dominates; for moderate-replicate designs (n=5-10), HET is conservative and deltapsi is more powerful.

Symptom: HET reports few hits in n=5-10 designs; deltapsi on the same data reports many.

Fix: Use deltapsi for n=3-5; reserve HET for n>=10 with explicit cohort heterogeneity.

SUPPA2 Empirical: Sparse Null at Low Replicate

Trigger: n<=3 vs n<=3 with --method empirical.

Mechanism: Empirical null is ECDF of |ΔPSI| from between-replicate comparisons within each group, binned by transcript expression. Few replicates -> few null observations -> wide confidence on null distribution.

Symptom: Inflated FDR (15-30%); "significant" hits don't replicate or validate.

Fix: Use -m classical (Wilcoxon) for n<=3 vs n<=3; or switch tool entirely (leafcutter, Shiba).

Reconciliation: When Tools Disagree

The two most common short-read tools answer slightly different questions: rMATS classifies on annotated event templates; leafcutter classifies on observed cluster usage. Disagreement is informative.

Operational rule: for high-confidence reporting, require concordant detection in two tools from different algorithmic families (event-based + cluster-based, or LSV + isoform-based). Document both calls and any explainable disagreements.

rMATS Output Columns Reference