GPTomics/bio-workflows-splicing-pipeline

Version Compatibility

Reference examples tested with: STAR 2.7.11+, fastp 0.23+, rMATS-turbo 4.3+, leafcutter 0.2.9+, Salmon 1.10+, tximport 1.30+, DRIMSeq 1.30+, DEXSeq 1.48+, stageR 1.24+, IsoformSwitchAnalyzeR 2.2+, RSeQC 5.0+, ggsashimi 1.1+ (numpy 1.26+, 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
• 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.

Note: rMATS-turbo --readLength must match the trimmed read length (add --variable-read-length if trimming produced a range); IsoformSwitchAnalyzeR's importRdata argument is genuinely spelled isoformExonAnnoation (a package typo, not an error here). Confirm in-tool before quoting.

Alternative Splicing Analysis Pipeline

"Find differential alternative splicing between my two conditions" -> Chain QC/trim, cohort-consistent 2-pass alignment, junction QC, event-level and (parallel) isoform-level differential testing, and sashimi visualization.

• CLI + R: fastp -> STAR 2-pass (shared SJ DB) -> junction QC -> rMATS-turbo + leafcutter -> [Salmon -> tximport dtuScaledTPM -> DRIMSeq/DEXSeq -> stageR] -> ggsashimi

This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step. Every step below cross-references the component skill that teaches its mechanism.

The governing principle

Splicing quantification is comparative by construction, so the trustworthy result is decided at four seams, not inside any one caller.

1. Splice-aware alignment must be cohort-consistent — one shared junction database for the whole experiment. STAR 2-pass where pass 2 uses the COMBINED SJ.out.tab from all samples is a made-once commitment, the splicing analogue of a shared reference. Junctions discovered per-sample and applied unevenly make PSI values non-comparable across samples, so a "differential" event is really a coverage artifact. Never run 2-pass independently per sample for a cohort.
2. The analysis lives at splice-aware resolution and does NOT collapse to gene. Feeding the standard gene-level tximport import (txOut=FALSE) averages the isoform-switch signal away. The DTU branch uses the OPPOSITE import from workflows/rnaseq-to-de: countsFromAbundance="dtuScaledTPM" + txOut=TRUE.
3. Event-level and isoform-level tests answer different questions and are reconciled, not merged. rMATS/leafcutter test exon-inclusion (ΔPSI); DTU tests within-gene isoform-proportion shifts. A gene can be significant for one and not the other; treating their p-values as interchangeable is a category error.
4. Transcript-level DTU FDR is only honest through the stageR two-stage test. Report per-transcript q-values only after a gene-level screen has passed; raw transcript-level FDR is inflated. stageR is the seam that converts gene-screened-then-transcript-confirmed tests into calibrated FDR.

Pipeline map

FASTQ (paired)
  | [1] QC & trim ---------------> fastp                (read-qc/fastp-workflow)
  v
  | [2] STAR 2-pass -------------> pass1 (all) -> COMBINED SJ.out.tab -> pass2 (all)  (read-alignment/star-alignment)
  v     ^-- commitment: annotation GTF + ONE shared junction DB + readLength
  | [3] Junction QC -------------> saturation plateau, entropy   (alternative-splicing/splicing-qc)
  v
  +---------------------------+-----------------------------------+
  | EVENT level               | ISOFORM level (parallel, optional)|
  v                           v
[4a] rMATS-turbo + leafcutter [4b] Salmon -> tximport dtuScaledTPM(txOut=TRUE)
     (alternative-splicing/       -> DRIMSeq/DEXSeq -> stageR
      differential-splicing)      (alternative-splicing/isoform-switching)
  |   ^-- reconcile, don't merge      ^-- gene screen BEFORE transcript q
  v
[5] Sashimi on top events -----> ggsashimi   (alternative-splicing/sashimi-plots)

Made-once commitments

Decided before the first alignment; every downstream PSI/DTU value inherits them.

| Commitment | Choice | Consequence inherited downstream | |------------|--------|----------------------------------| | Annotation GTF | ONE GTF used by rMATS, leafcutter, and IsoformSwitchAnalyzeR | Different GTFs make their events irreconcilable; fixes the event universe | | 2-pass junction DB | Combined SJ.out.tab from ALL samples fed into pass 2 | Per-sample junctions -> non-comparable PSI, false "differential" events from coverage | | Measurement level | Splice-aware (event PSI and/or transcript DTU); NEVER gene collapse | Gene-level tximport (txOut=FALSE) averages away the switch signal | | --readLength | Matches the trimmed read length (or --variable-read-length) | A wrong value miscomputes inclusion-junction lengths and biases PSI |

The canonical order and why

1. QC/trim — but do NOT over-trim; short reads lose junction-spanning power.
2. STAR pass 1 (all samples) -> collect every SJ.out.tab.
3. Build ONE combined junction DB from the concatenated SJ.out.tab (order-trap: skipping this / per-sample DBs makes PSI incomparable).
4. STAR pass 2 (all samples, same combined DB) -> coordinate-sorted BAMs.
5. Junction QC — saturation curves must plateau (else deeper sequencing), entropy sane.
6. Event-level differential splicing — rMATS-turbo plus leafcutter for concordance.
7. (Parallel) isoform-level DTU — Salmon transcript quant -> tximport dtuScaledTPM+txOut=TRUE -> DRIMSeq/DEXSeq -> stageR (order-trap: reporting transcript q-values without the stageR gene screen inflates FDR).
8. Sashimi on the top reconciled events.

Order-traps that silently produce wrong results: per-sample (not cohort) 2-pass junctions; collapsing to gene before splicing analysis; transcript FDR without stageR; mixing ΔPSI significance with DTU proportion significance as if interchangeable.

Choosing event-level vs isoform-level (and the caller)

Pipeline-level selection only; mechanism lives in the component skills.

| Fork | Lean toward | Hand off to | |------|-------------|-------------| | Event-level vs isoform-level | rMATS/leafcutter (which exon?) vs IsoformSwitchAnalyzeR/DRIMSeq+DEXSeq (which isoform, + NMD/ORF/domain consequences) | alternative-splicing/differential-splicing, alternative-splicing/isoform-switching | | rMATS vs leafcutter vs MAJIQ | rMATS (known event types, replicate-based) + leafcutter (annotation-free intron clusters) for concordance; MAJIQ V3 HET for complex/heterogeneous cohorts | alternative-splicing/differential-splicing | | DTU import scale | Salmon -> tximport dtuScaledTPM + txOut=TRUE (the ONLY correct DTU import) | rna-quantification/tximport-workflow | | Transcript FDR | stageR two-stage: gene-level screen then transcript confirmation | alternative-splicing/isoform-switching |

Primary path: STAR 2-pass + rMATS-turbo (+ leafcutter)

Goal: produce cohort-comparable PSI and a filtered differential-event table.

Approach: align all samples through one shared junction DB, QC junction saturation, then run rMATS-turbo (and leafcutter for concordance). Full runnable script: examples/splicing_pipeline.sh.

# Pass 1 (all samples) -> collect junctions
STAR --runThreadN 8 --genomeDir star_index/ --readFilesIn ${s}_R1.fq.gz ${s}_R2.fq.gz \
    --readFilesCommand zcat --outSAMtype BAM Unsorted --outFileNamePrefix ${s}_pass1_
cat *_pass1_SJ.out.tab > combined_SJ.out.tab          # ONE shared DB for the whole cohort
# Pass 2 (all samples, same combined DB) -> comparable coordinates
STAR --runThreadN 8 --genomeDir star_index/ --readFilesIn ${s}_R1.fq.gz ${s}_R2.fq.gz \
    --readFilesCommand zcat --sjdbFileChrStartEnd combined_SJ.out.tab \
    --outSAMtype BAM SortedByCoordinate --outFileNamePrefix ${s}_

# Differential splicing. --readLength MUST match the trimmed reads, and --variable-read-length is
# required because the fastp step above trims to a RANGE, not a single length; without it rMATS
# miscomputes inclusion-junction lengths and biases PSI.
rmats.py --b1 cond1_bams.txt --b2 cond2_bams.txt --gtf annotation.gtf \
    -t paired --readLength 150 --variable-read-length --nthread 8 --od rmats_output --tmp rmats_tmp

Filter events on |IncLevelDifference| > 0.1, FDR < 0.05, and >=10 supporting junction reads averaged per replicate (sum IJC_SAMPLE_1, SJC_SAMPLE_1, IJC_SAMPLE_2, SJC_SAMPLE_2 — each a comma-separated per-replicate list — and divide by the replicate count); rank by -log10(FDR) * |IncLevelDifference| (clamp FDR with max(FDR, 1e-300)). Resolve every column by header name: MXE carries two extra coordinate columns, so fixed positions silently read the wrong field. The read floor is not cosmetic — PSI is a ratio, so a 2-read event can reach |dPSI| = 0.9 and pass FDR on noise alone.

Parallel path: isoform-level DTU (Salmon -> stageR)

Goal: detect within-gene isoform-proportion switches with honest transcript-level FDR.

Approach: quantify transcripts with Salmon, import at DTU scale, test proportions with DRIMSeq/DEXSeq, and gate transcript q-values through stageR.

library(tximport)
# DTU import is the OPPOSITE of gene-level DGE: keep transcripts, dtuScaledTPM.
# dtuScaledTPM scales by median tx length AMONG a gene's isoforms, so tx2gene is required even with txOut.
txi <- tximport(files, type = 'salmon', txOut = TRUE, countsFromAbundance = 'dtuScaledTPM', tx2gene = tx2gene)

# Canonical two-stage FDR route (Love, Soneson & Patro 2018): DRIMSeq/DEXSeq proportion test, then
# stageR gene-level SCREEN -> transcript-level CONFIRM. This is the "stageR seam" the principle names;
# mechanism lives in alternative-splicing/isoform-switching.

# Alternative route -- IsoformSwitchAnalyzeR adds NMD/ORF/protein-domain consequences. Its gene-level
# q-value comes from DEXSeq's perGeneQValue (min-p aggregation), NOT the stageR package:
library(IsoformSwitchAnalyzeR)
sList <- importRdata(isoformCountMatrix = txi$counts, isoformRepExpression = txi$abundance,
                     designMatrix = design, isoformExonAnnoation = 'annotation.gtf',
                     isoformNtFasta = 'transcripts.fa')      # 'isoformExonAnnoation' is the real (typo'd) arg
sList <- isoformSwitchTestDEXSeq(sList, reduceToSwitchingGenes = TRUE)   # gene q via DEXSeq perGeneQValue

When NOT to use this pipeline (regime off-ramps)

This pipeline targets bulk short-read differential splicing between two groups. For other regimes, use the dedicated skill.

| Question | Use instead | |----------|-------------| | "Does this DNA variant alter splicing?" | alternative-splicing/splice-variant-prediction (SpliceAI, Pangolin, MMSplice) | | "What is aberrant in this single rare-disease patient?" | alternative-splicing/outlier-splicing-detection (FRASER 2.0, OUTRIDER, DROP) | | "Full-isoform analysis from PacBio Iso-Seq / ONT" | alternative-splicing/long-read-splicing (FLAIR, IsoQuant, Bambu, SQANTI3) | | "Single-cell splicing analysis" | alternative-splicing/single-cell-splicing (chemistry-first; MARVEL, BRIE2) | | "Heterogeneous cohort, n>=10 vs n>=10" | This pipeline + MAJIQ V3 HET (alternative-splicing/differential-splicing) | | "Microexon-focused (3-27 nt)" | This pipeline with VAST-TOOLS or MicroExonator (alternative-splicing/splicing-quantification) |

QC checkpoints between steps

| After | Gate | Interpretation | |-------|------|----------------| | QC/trim | Q30 >80%, adapter <5%, reads not over-trimmed | Aggressive trimming below ~75 nt weakens junction-spanning evidence | | Alignment | Uniquely-mapped >80%; junction-saturation curves plateau | Still-rising curves = under-sequenced for splicing; deeper reads needed | | Differential | |ΔPSI| >0.1 / >0.2, FDR <0.05, >=10 junction reads | Low read support = PSI is noise; require the read floor per event | | DTU | stageR gene-level screen passed | Transcript q-values are only valid after the gene screen (alternative-splicing/isoform-switching) |

Common Errors

| Symptom | Cause | Fix | |---------|-------|-----| | "Differential" events that are really coverage differences | Per-sample 2-pass junctions, not a shared DB | Concatenate all SJ.out.tab, feed the combined DB to pass 2 for every sample | | Isoform-switch signal disappears | Gene-level tximport (txOut=FALSE) before splicing | Keep transcripts: txOut=TRUE + dtuScaledTPM on the DTU branch | | Too many "significant" transcripts | Transcript q-values reported without stageR | Run the stageR two-stage gene->transcript test | | Biased PSI across samples | --readLength != trimmed length | Set --readLength to the real length or add --variable-read-length | | Event- and isoform-level calls "disagree" | Treating ΔPSI and DTU proportion tests as the same question | Reconcile as complementary; they test different things |

Pipeline map (hand-offs)

• read-qc/fastp-workflow - QC/trim without over-trimming junction-spanning reads
• read-alignment/star-alignment - STAR 2-pass cohort-style, the shared junction DB
• alternative-splicing/splicing-qc - junction saturation/entropy, depth thresholds
• alternative-splicing/splicing-quantification - PSI computation, event taxonomy, sign conventions
• alternative-splicing/differential-splicing - rMATS/leafcutter/MAJIQ selection and reconciliation
• rna-quantification/alignment-free-quant - Salmon transcript quant for the DTU branch
• rna-quantification/tximport-workflow - dtuScaledTPM + txOut import
• alternative-splicing/isoform-switching - DTU + stageR + NMD/ORF/domain consequences
• alternative-splicing/sashimi-plots - ggsashimi/leafviz visualization

The complete runnable script is in this skill's examples/ (splicing_pipeline.sh).

Related Skills

• read-qc/fastp-workflow - QC/trim options
• read-alignment/star-alignment - STAR 2-pass cohort-style configuration
• alternative-splicing/splicing-quantification - PSI computation, event taxonomy, sign conventions
• alternative-splicing/differential-splicing - Tool selection, MAJIQ V3, leafcutter, reconciliation
• alternative-splicing/isoform-switching - DTU + NMD/ORF/domain consequences (IsoformSwitchAnalyzeR v2, stageR)
• alternative-splicing/sashimi-plots - ggsashimi, MAJIQ-VOILA, leafviz visualization
• alternative-splicing/splicing-qc - STAR 2-pass cohort-style, library prep, depth thresholds
• alternative-splicing/single-cell-splicing - 10X chemistry decision; plate-based and long-read SC
• alternative-splicing/splice-variant-prediction - SpliceAI / Pangolin / MMSplice variant interpretation
• alternative-splicing/outlier-splicing-detection - FRASER 2.0 / DROP rare-disease workflow
• alternative-splicing/long-read-splicing - PacBio HiFi / ONT full-isoform analysis
• rna-quantification/alignment-free-quant - Salmon TPM for SUPPA2 and DTU pipelines
• rna-quantification/tximport-workflow - dtuScaledTPM + txOut DTU import

References

• Shen S, Park JW, Lu ZX, et al (2014) rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. PNAS 111:E5593-E5601. DOI 10.1073/pnas.1419161111.
• Li YI, Knowles DA, Humphrey J, et al (2018) Annotation-free quantification of RNA splicing using LeafCutter. Nature Genetics 50:151-158. DOI 10.1038/s41588-017-0004-9.
• Vitting-Seerup K, Sandelin A (2019) IsoformSwitchAnalyzeR: analysis of changes in genome-wide patterns of alternative splicing and its functional consequences. Bioinformatics 35:4469-4471. DOI 10.1093/bioinformatics/btz247.
• Van den Berge K, Soneson C, Robinson MD, Clement L (2017) stageR: a general stage-wise method for controlling the gene-level false discovery rate in differential expression and differential transcript usage. Genome Biology 18:151. DOI 10.1186/s13059-017-1277-0.
• Love MI, Soneson C, Patro R (2018) Swimming downstream: statistical analysis of differential transcript usage following Salmon quantification. F1000Research 7:952. DOI 10.12688/f1000research.15398.3. (the dtuScaledTPM -> DRIMSeq/DEXSeq -> stageR two-stage DTU workflow.)