Orchestrates the end-to-end bulk short-read alternative-splicing pipeline from FASTQ to differential splicing, chaining fastp QC, cohort-consistent STAR 2-pass alignment (one shared junction DB), junction QC, event-level differential splicing (rMATS-turbo + leafcutter, optional MAJIQ V3), parallel isoform-level DTU (Salmon -> tximport dtuScaledTPM -> DRIMSeq/DEXSeq -> stageR), and sashimi visualization. Use when committing the annotation GTF and a shared 2-pass junction database for the whole cohort, keeping the analysis at splice-aware resolution (never collapsing to gene), choosing event-level vs isoform-level DTU and reconciling them, applying the stageR two-stage gene->transcript FDR, or off-ramping to splice-variant / outlier / long-read / single-cell splicing. Hands mechanism to the alternative-splicing component skills; not a re-teach of any single step.
Orchestrates the end-to-end bulk short-read alternative-splicing pipeline from FASTQ to differential splicing, chaining fastp QC, cohort-consistent STAR 2-pass alignment (one shared junction DB), junction QC, event-level differential splicing (rMATS-turbo + leafcutter, optional MAJIQ V3), parallel isoform-level DTU (Salmon -> tximport dtuScaledTPM -> DRIMSeq/DEXSeq -> stageR), and sashimi visualization. Use when committing the annotation GTF and a shared 2-pass junction database for the whole cohort, keeping the analysis at splice-aware resolution (never collapsing to gene), choosing event-level vs isoform-level DTU and reconciling them, applying the stageR two-stage gene->transcript FDR, or off-ramping to splice-variant / outlier / long-read / single-cell splicing. Hands mechanism to the alternative-splicing component skills; not a re-teach of any single step.
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.
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.
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.
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
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.
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
QC/trim — but do NOT over-trim; short reads lose junction-spanning power.
STAR pass 1 (all samples) -> collect every SJ.out.tab.
Build ONE combined junction DB from the concatenated SJ.out.tab (order-trap: skipping this / per-sample DBs makes PSI incomparable).
STAR pass 2 (all samples, same combined DB) -> coordinate-sorted BAMs.
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)
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.
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.
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.)
