swaruplab/operon/alternative-splicing-splicing-quantification

Version Compatibility

Reference examples tested with: rMATS-turbo 4.3+, SUPPA2 2.4+, leafcutter 0.2.9+, MAJIQ 3.0+, IRFinder-S 2.0+, kallisto 0.50+, Salmon 1.10+, pandas 2.2+, STAR 2.7.11+

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.

Splicing Quantification

Quantify alternative splicing events as PSI (percent spliced in) from RNA-seq. PSI = inclusion read evidence / (inclusion + skipping read evidence), normalized for differential mapping opportunity between isoforms. The choice of quantification unit (event, intron cluster, LSV, transcript) determines which biological questions can be answered and which failure modes apply.

Algorithmic Taxonomy

The agent's first decision is which family the question requires, not which tool. Switching family is the response to a tool failing within a family — switching tool within a family rarely fixes systematic blind spots.

Event Taxonomy (Beyond Standard SE/A5SS/A3SS/MXE/RI)

Tool-agnostic taxonomy reference: Wang 2008 Nature; Vaquero-Garcia 2016 eLife (LSV); Tapial 2017 Genome Res (VastDB).

Tool Selection Matrix

Methodology evolves; verify benchmarks (Olofsson 2023 Brief Bioinform; Kubota 2025 NAR; Tran 2025 WIREs RNA) and tool docs before committing. Default 2026 recommendation: run rMATS-turbo + leafcutter and reconcile; add MAJIQ V3 for complex events / heterogeneous cohorts; switch to Shiba for n=2 vs n=2.

PSI Definition and Effective Length Normalization

For a cassette exon, naive PSI ignores that the inclusion isoform contains more positions where a junction read can map than the skipping isoform. rMATS reports IncFormLen (= 2*(read_length - anchor) for the two flanking junctions, plus exon body bases) and SkipFormLen (= read_length - anchor) and computes:

PSI = (IJC / IncFormLen) / (IJC / IncFormLen + SJC / SkipFormLen)

where IJC = inclusion junction counts, SJC = skipping junction counts. Skipping this normalization biases PSI by ~10-30% depending on read length and exon size. SUPPA2 derives PSI from transcript TPMs, which the upstream Salmon/kallisto already accounts for. leafcutter operates on intron usage proportions within a cluster (a different statistic).

For long-read data, every read carries full isoform identity — effective-length normalization becomes unnecessary because each read counts as one isoform.

Sign Conventions for Alternative Splice Sites

Always record which alternative form ΔPSI > 0 corresponds to in publication-grade reporting. Confusion is the most common reviewer comment for AS papers.

rMATS-turbo Workflow

Goal: Quantify SE/A5SS/A3SS/MXE/RI events from BAMs aligned with STAR 2-pass.

Approach: Group BAMs by condition, run rMATS with --statoff for quantification only, then parse JC.txt files for per-replicate PSI.

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 \
    --statoff

Key flags: --novelSS discovers junctions absent from the GTF (recommended with STAR 2-pass output). --variable-read-length allows mixed read lengths in the cohort. --libType fr-firststrand matches Illumina TruSeq stranded; verify with RSeQC infer_experiment.py. --statoff is for quantification-only runs; omit for differential testing.

import pandas as pd
se_jc = pd.read_csv('rmats_output/SE.MATS.JC.txt', sep='\t')
inc_cols = [c for c in se_jc.columns if c.startswith('IncLevel')]
se_jc['mean_PSI'] = se_jc[inc_cols].mean(axis=1)
per_rep_inc = se_jc['IJC_SAMPLE_1'].str.split(',').apply(lambda x: list(map(int, x)))
per_rep_skip = se_jc['SJC_SAMPLE_1'].str.split(',').apply(lambda x: list(map(int, x)))
min_inc = per_rep_inc.apply(min)
min_skip = per_rep_skip.apply(min)
reliable = se_jc[(min_inc + min_skip) >= 20]

JC vs JCEC files

SE.MATS.JC.txt uses only junction-spanning reads. SE.MATS.JCEC.txt adds reads contained within the alternative exon body as inclusion evidence.

• Prefer JC for clean cassette-exon analysis when reads spanning junctions are sufficient.
• Use JCEC when alternative exons are short (<50nt) and junction-spanning reads are scarce.
• Avoid JCEC when intron retention overlaps the alternative exon body — exon-body reads may come from retained introns, not inclusion isoform.

SUPPA2 Workflow

Goal: Compute event PSI from transcript TPM without alignment; useful when Salmon/kallisto TPMs already exist.

Approach: Generate IOE event definitions from GTF, then aggregate TPMs of transcripts including/excluding each event.

suppa.py generateEvents -i annotation.gtf -o events -f ioe -e SE SS MX RI AF AL
 
for ev in SE A5 A3 MX RI; do
    suppa.py psiPerEvent -i events_${ev}_strict.ioe -e transcript_tpm.tsv -o psi_${ev}
done

SUPPA2 is annotation-bound: events absent from the GTF cannot be quantified. Whether an event is detected depends entirely on which transcripts the upstream Salmon/kallisto index contains. Use GENCODE comprehensive over basic when SUPPA2 detection sensitivity matters.

MAJIQ V3 Workflow

Goal: Quantify LSVs with Bayesian posterior PSI distributions; ideal for complex multi-junction events that don't fit canonical event types.

Approach: Build a splice graph from BAMs + GFF3, compute per-junction coverage with bootstrap, then run majiq psi for posterior PSI per LSV.

majiq build annotation.gff3 -c settings.ini -j 8 -o build_output
majiq psi build_output/sample1.majiq build_output/sample2.majiq -j 4 -o psi_output -n condition_psi
voila view -p 5000 -j 8 build_output/splicegraph.zarr psi_output/condition_psi.psi.voila -o voila_output

MAJIQ V3 (Aicher, Slaff, Jewell, Barash bioRxiv 2024; public release 2025) replaced V2's SQLite splicegraph (splicegraph.sql) with Zarr storage (splicegraph.zarr); the .sql is deprecated. V3 is ~3.2x faster than V2 via xarray/zarr/Dask parallelization. LSV output includes posterior mean PSI plus the full posterior distribution; this enables threshold-based testing (e.g. P(|ΔPSI| > 0.2)) rather than frequentist p-values.