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.

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.

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.