GPTomics/bio-outlier-splicing-detection

Version Compatibility

Reference examples tested with: FRASER 2.0 (>=1.99.0), OUTRIDER 1.20+, LeafcutterMD via leafcutter 0.2.9+, DROP 1.4+, R 4.4+, BiocManager 1.30+

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

• 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.

Outlier Splicing Detection

For clinical RNA-seq diagnostics in rare disease, the question is not "what differs between groups?" but "what is aberrant in this single patient relative to a panel of unaffected samples?". The statistical framework is single-sample-vs-cohort outlier detection, fundamentally different from two-group differential splicing. Tools in this space are designed for clinical Mendelian diagnostic settings.

Tool Taxonomy

| Tool | Statistic | Test target | Fails when | |------|-----------|-------------|------------| | FRASER 2.0 | Beta-binomial autoencoder on Intron Jaccard Index | Splicing outliers (per-sample, per-junction) | Cohort <20 samples; tissue mismatch | | OUTRIDER | Autoencoder-denoised expression Z-score | Gene-level expression outliers (LoF, monoallelic) | Cohort <20 samples | | LeafcutterMD | Dirichlet-multinomial outlier mode | Annotation-free intron usage | Beta-binomial fits poorly OR few controls | | DROP | Snakemake pipeline | All of above + monoallelic expression | Pipeline complexity for small projects |

Core reference: FRASER 2.0 for splicing outliers, OUTRIDER for expression outliers, DROP to combine. Standard tool in EU rare-disease programs (Solve-RD) and NIH UDN.

Decision Tree by Diagnostic Scenario

| Scenario | Recommended approach | |----------|----------------------| | Single rare-disease patient + panel of n>=50 controls | FRASER 2.0 (Intron Jaccard Index) | | Single patient + small panel (n=20-50) | FRASER 2.0 with auxiliary GTEx controls; tune q carefully | | Patient + cohort <20 | Insufficient for outlier detection; consider differential or recruit more samples | | Outlier expression suspected (loss of function, monoallelic) | OUTRIDER on same cohort | | Annotation-free outlier (cryptic exon, novel junction) | LeafcutterMD | | Integrated diagnostic pipeline (splicing + expression + MAE) | DROP | | TDP-43 ALS post-mortem brain (cryptic exons) | FRASER 2.0; expect UNC13A, STMN2, ATG4B | | SF3B1-mutant cancer sample | FRASER 2.0 with cohort-matched RNA-seq; expect cryptic 3'ss | | Familial dysautonomia (ELP1) | FRASER 2.0 in fibroblast/iPSC; CNS tissue gives strongest signal | | Stargardt deep-intronic ABCA4 | FRASER 2.0 in retina-relevant tissue | | Solid tumor splicing biomarker | Differential splicing (n>=10 vs cohort) — see differential-splicing skill | | RNA validation of SpliceAI hit | FRASER 2.0 + cross-reference with predicted variant location |

When to Use Outlier vs Differential

Outlier regime (this skill):

• Single patient or small case series vs control panel
• Question: "What is aberrant in this patient?"
• Statistical model: single-sample p-value vs cohort distribution

Differential regime (differential-splicing skill):

• Two well-defined groups, n>=3 each
• Question: "What differs between groups?"
• Statistical model: two-group LRT or related

If n>=10 patients with a shared phenotype are available, prefer differential (more power); if single patient or heterogeneous case series, use outlier.

FRASER 2.0 Workflow

Goal: Detect aberrant splicing in patient samples vs cohort using the Intron Jaccard Index.

Approach: Count split reads per junction, compute Intron Jaccard Index per intron, fit a Beta-binomial autoencoder to estimate expected values, then flag outliers by p-value and delta.

library(FRASER); library(BiocParallel)

bam_files <- list.files('bams/', pattern='.bam$', full.names=TRUE)
sample_table <- data.frame(
    sampleID = gsub('.bam', '', basename(bam_files)),
    bamFile = bam_files,
    pairedEnd = TRUE
)

settings <- FraserDataSet(
    colData = sample_table,
    workingDir = 'fraser_workdir',
    name = 'rare_disease_cohort'
)

settings <- countRNAData(settings, BPPARAM = MulticoreParam(8))
fds <- calculatePSIValues(settings)

fds <- filterExpressionAndVariability(
    fds,
    minDeltaPsi = 0.0,
    minExpressionInOneSample = 20,
    quantile = 0.05,
    quantileMinExpression = 1
)

fitMetrics(fds) <- 'jaccard'
currentType(fds) <- 'jaccard'  # canonical setter for active metric in FRASER 2.0
fds <- FRASER(
    fds,
    q = c(jaccard = 10),
    BPPARAM = MulticoreParam(8)
)

results <- results(
    fds,
    psiType = 'jaccard',
    padjCutoff = 0.05,
    deltaPsiCutoff = 0.1
)

patient_results <- results[results$sampleID == 'PATIENT_001', ]
patient_results <- patient_results[order(patient_results$padjust), ]

FRASER 2.0 changes vs FRASER 1.x:

• Default psiType changed from three metrics (psi5, psi3, theta) to single Intron Jaccard Index
• Default deltaPsiCutoff dropped from 0.3 to 0.1
• Pseudocount and filtering parameter optimization
• Bioconductor package version >=1.99.0 == FRASER 2.0

q = 10 is the autoencoder dimension hyperparameter. Tune via estimateBestQ(fds, type='jaccard') for cohort-specific optimum — too low: confounders not removed; too high: real signal absorbed.

OUTRIDER for Gene-Level Outlier Expression

Goal: Detect genes with aberrantly high or low expression in patient samples.

Approach: Autoencoder denoising of expression matrix; outliers identified by Z-score and adjusted p-value.

library(OUTRIDER); library(BiocParallel)

countTable <- read.table('counts.tsv', header=TRUE, row.names=1)
ods <- OutriderDataSet(countData = countTable)

ods <- filterExpression(ods, minCounts=TRUE, filterGenes=TRUE)
# OUTRIDER's estimateBestQ returns a scalar q (unlike FRASER's, which returns the object)
q_best <- estimateBestQ(ods)
ods <- OUTRIDER(ods, q = q_best, BPPARAM = MulticoreParam(8))

res <- results(ods, padjCutoff = 0.05, zScoreCutoff = 0)
patient_outliers <- res[res$sampleID == 'PATIENT_001', ]

OUTRIDER (Brechtmann 2018 Am J Hum Genet) catches loss-of-function alleles producing transcript collapse, monoallelic effects, and tissue-inappropriate expression — complements splice outlier detection.

LeafcutterMD for Annotation-Free Outlier Intron Usage

Goal: Detect outlier intron usage relative to a control panel without annotation dependence.

Approach: Run LeafcutterMD (LeafCutter's Dirichlet-multinomial outlier mode for Mendelian disease) against the control panel.

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

leafcutterMD.R \
    --num_threads 4 \
    --output_prefix patient_outlier \
    leafcutter_perind_numers.counts.gz

LeafcutterMD (Jenkinson 2020 Bioinformatics) reports per-sample p-values per intron-cluster; useful when FRASER's Beta-binomial model fits poorly or when novel-junction sensitivity matters.

DROP Pipeline (Integrated Workflow)

Goal: Run FRASER2 + OUTRIDER + monoallelic expression in a unified Snakemake pipeline for clinical diagnostics.

Approach: DROP is distributed via bioconda (not PyPI). Install in a dedicated environment, then configure with patient + control sample sheets; pipeline handles QC, alignment, counting, autoencoding, and reporting.

# Install via bioconda (DROP is not on PyPI)
mamba create -n drop_env -c conda-forge -c bioconda drop --override-channels
conda activate drop_env

drop init my_diagnostic_run
cd my_diagnostic_run

# Edit config.yaml:
#  - sample_table: samples.tsv (patient + controls)
#  - aberrantSplicing: enabled
#  - aberrantExpression: enabled
#  - mae: enabled (monoallelic expression)

snakemake --cores 16 --use-conda

DROP (Yepez 2021 Nat Protocols) is the standard tool in EU rare-disease genome+RNA-seq programs (Solve-RD) and the NIH UDN. v1.4+ uses FRASER 2.0. The MAE module uses a custom z-score test on heterozygous SNPs from RNA-seq (allele-specific expression) — useful for catching dominant-negative or monoallelic LoF that splicing/expression outliers miss. Cohort >=30 samples recommended for confident outlier detection.

Variant + Outlier Integration

Goal: Connect a candidate splice-altering DNA variant to RNA-level confirmation.

Approach: Cross-reference SpliceAI hits with FRASER2 outliers in the same sample.

library(dplyr)

variants <- read.table('spliceai_hits.tsv', header=TRUE, sep='\t')
fraser_hits <- read.table('fraser_results.tsv', header=TRUE, sep='\t')

confirmed <- variants %>%
    filter(delta_max >= 0.2) %>%
    inner_join(
        fraser_hits %>% filter(sampleID == 'PATIENT_001', padjust < 0.05),
        by = c('chrom' = 'seqnames'),
        relationship = 'many-to-many'
    ) %>%
    filter(abs(pos - start) < 1000 | abs(pos - end) < 1000)

A SpliceAI hit + concordant FRASER2 outlier in the patient = strong PS3 functional evidence in the ACMG framework. This integration is the highest-value clinical pipeline step — converts a computational PP3 to functional PS3.

Cohort Size and Power

| Cohort size | Power | Comment | |-------------|-------|---------| | n < 20 | Marginal | High FDR; consider GTEx tissue-matched controls as auxiliary | | n = 20-50 | Acceptable | FRASER autoencoder can fit; tune q carefully | | n >= 50 | Recommended | Standard clinical diagnostic cohort size | | n >= 100 | Optimal | Tissue-matched and batch-matched gives best calibration |

GTEx-derived tissue-matched controls can supplement small in-house cohorts but introduce batch effects; use only when in-house n < 30 and document the pooling strategy.

Tissue Choice for Mendelian RNA-seq

| Tissue | Pros | Cons | Genes captured | |--------|------|------|----------------| | Whole blood (PAXgene) | Easy, standard | Globin contamination; many disease genes silent | ~70-80% of clinical genes | | Fibroblast (skin biopsy) | Reasonable expression | Requires culture; senescence variability | ~75-85% | | Muscle biopsy | Best for muscular dystrophy | Invasive | ~85-90% for muscle disorders | | iPSC-derived neuron / cardiomyocyte | Disease-relevant tissue | Cost, variability | ~95% if differentiation works | | Urine sediment | Non-invasive | Low yield | ~50-60% |

For UDN-style cases: blood first, then fibroblast if blood lacks expression of candidate gene. Critical: a negative blood RNA-seq doesn't rule out a candidate gene that's silent in blood — verify gene expression with GTEx tissue panel before committing to the tissue.

Hyperparameter Tuning

# useOHT=FALSE runs the injection-based q grid so plotEncDimSearch has a curve to show;
# useOHT=TRUE (default) is the fast deterministic OHT estimate but produces no search table to plot.
fds <- estimateBestQ(fds, type='jaccard', useOHT=FALSE, q_param=c(2, 5, 10, 15, 20))
plotEncDimSearch(fds, type='jaccard')

The encoding dimension q should be where the loss curve plateaus. Too low: confounders not removed; too high: real signal absorbed.

For typical 50-100 sample cohorts, q=8-15 is the usual operating range (DROP / FRASER workflow convention; no single primary citation — verify with plotEncDimSearch on the actual cohort). For very small cohorts (n=20-30), q=5-8 is typical.

Per-Tool Failure Modes

FRASER 2.0: Q Hyperparameter Mistuning

Trigger: Default q=10 used without tuning; or wrong q for cohort size.

Mechanism: Q is the autoencoder bottleneck dimension; too small -> confounders leak into outlier signal; too large -> real biological signal absorbed by autoencoder.

Symptom: Either no significant outliers (q too high) or many spurious calls clustering by batch (q too low).

Fix: Run estimateBestQ(fds, type='jaccard', useOHT=TRUE) (Optimal Hard Thresholding default; very fast); use the returned bestQ(fds) value. For exhaustive search, pass useOHT=FALSE, q_param=c(2, 5, 10, 15) and inspect plotEncDimSearch for the plateau.

FRASER 2.0: Tissue Mismatch

Trigger: Patient sample from different tissue than majority of controls.

Mechanism: FRASER autoencoder learns tissue-specific expression patterns; tissue-mismatched patient appears as global outlier.

Symptom: Hundreds of "significant" outliers in patient; not biologically interpretable.

Fix: Strict tissue matching; if controls are mixed-tissue, use only controls from patient's tissue.

OUTRIDER: Few Controls

Trigger: Cohort <20 samples.

Mechanism: Autoencoder needs sufficient samples to learn expression covariance; fails to fit at very small cohort sizes.

Symptom: Convergence warnings; uncalibrated p-values.

Fix: Pool with GTEx auxiliary controls; or use simpler outlier methods (z-score on log-CPM).

LeafcutterMD: Cluster Count Limits

Trigger: Very few clusters in patient sample (low coverage or filtered out).

Mechanism: LeafcutterMD fits a Dirichlet-multinomial (Beta-binomial per-intron) model over cluster counts; few observations -> unstable fit -> unreliable p-values.

Symptom: Inflated or deflated p-values; few significant calls.

Fix: Increase coverage; relax filtering (-m 10 instead of 50); or switch to FRASER2.

DROP: Snakemake Pipeline Failures

Trigger: Missing dependencies or incompatible R/Bioconductor versions.

Mechanism: DROP orchestrates many tools; version mismatches cascade through pipeline.

Symptom: Snakemake step fails partway through; cryptic R errors.

Fix: Use --use-conda flag for environment isolation; pin versions in environment yamls.

Reconciliation: When Outlier Tools Disagree

| Pattern | Likely cause | Action | |---------|--------------|--------| | FRASER2 sig, OUTRIDER not | Splicing change without expression collapse | Standard splicing outlier; report | | OUTRIDER sig, FRASER2 not | Expression LoF without splicing change | Likely promoter / regulatory; not splicing | | Both sig at same gene | LoF allele triggering NMD on splicing-disrupted transcript | Strong combined evidence; expect downstream | | LeafcutterMD sig, FRASER2 not | Novel cryptic event not in annotation | High-priority novel finding; investigate | | All tools null but biology suggests change | Underpowered cohort or wrong tissue | Verify gene expression in tissue; recruit larger cohort |

Disease-Specific Expectations

| Condition | Expected outlier signature | Tissue | |-----------|----------------------------|--------| | ALS / FTD (TDP-43 loss) | Cryptic exons in UNC13A, STMN2, ATG4B | Post-mortem brain ONLY | | SF3B1-mutant MDS / CLL / uveal melanoma | Aberrant 3'ss ~10-30nt upstream of canonical | Bone marrow / tumor tissue | | Spinal muscular atrophy (untreated SMN2) | SMN exon 7 skipping | Fibroblast / iPSC-MN | | Familial dysautonomia (ELP1 c.2204+6T>C) | ELP1 exon 20 skipping (>=99% in CNS, partial elsewhere) | iPSC-neuron > fibroblast > blood | | Deep-intronic CFTR / USH2A / CEP290 | Pseudoexon inclusion | Cognate disease tissue (lung / retina) | | Duchenne muscular dystrophy (DMD) | Out-of-frame exon skipping pattern | Muscle biopsy | | Stargardt (ABCA4) deep-intronic | Pseudoexon in retina | Retinal organoid / iPSC-RPE |

For each, the gene must be expressed in the queried tissue. Verify with GTEx before assuming negative result rules out the gene.

Common Errors

| Error | Cause | Solution | |-------|-------|----------| | FRASER: cohort too small | n<10 | Pool with auxiliary controls; or recruit more patients | | FRASER: countRNAData failed on chromosome X | BAM index missing or corrupted | Re-index BAMs; check samtools idxstats | | estimateBestQ: convergence not reached | Default q range insufficient | Expand q_param=c(2,5,10,15,20); or use useOHT=TRUE for the fast deterministic alternative | | OUTRIDER: encoding-dim search slow | findEncodingDim grids many q values | Use estimateBestQ(ods) for a fast single-q estimate | | DROP: snakemake job failed at FRASER | DROP-FRASER version mismatch | Update DROP to latest; verify FRASER 2.0 compatibility | | LeafcutterMD: insufficient clusters | Cluster filter too strict | Lower -m minimum cluster reads | | Variant integration: chrom format mismatch | VCF uses 1, FRASER uses chr1 (or vice versa) | Standardize with bcftools annotate --rename-chrs |

Common Pitfalls

• Using bulk differential-splicing tools for n=1 vs cohort — rMATS, leafcutter (regular), SUPPA2 are not designed for this. Use FRASER2 / LeafcutterMD.
• Ignoring tissue choice — clinical gene expression varies dramatically across blood / fibroblast / muscle. A negative blood RNA-seq doesn't rule out a candidate gene that's silent in blood.
• Forgetting batch effects — combining in-house and external (GTEx) controls introduces sequencing batch confounding; use ComBat or include batch as covariate.
• Skipping the variant + outlier integration — RNA-only outlier without DNA variant suggests cellular state or technical artifact; DNA-only prediction without RNA confirmation is supporting only (PP3, not PS3).
• Treating all FRASER2 outliers as pathogenic — many are benign tissue-specific variation. Filter against gnomAD splice constraint and disease gene panels.
• Q hyperparameter not tuned — default q=10 works for ~50-100 sample cohorts; tune for outliers.
• Wrong default delta cutoff for FRASER 1.x vs 2.0 — 1.x default 0.3, 2.0 default 0.1; document which version.

Quality Thresholds

| Metric | Recommendation | Source | |--------|----------------|--------| | Cohort size | n>=50 (ideal); n>=20 (minimum) | Solve-RD / UDN convention | | FRASER 2.0 padj | < 0.05 | Standard | | FRASER 2.0 delta Jaccard | >= 0.1 (default in v2.0) | Scheller 2023 AJHG | | OUTRIDER padj | < 0.05 | Brechtmann 2018 AJHG | | OUTRIDER zScore | abs >= 2 | Conservative | | Sequencing depth | >=50M PE reads/sample | Standard for AS analysis | | Tissue match between patient and controls | Required | Critical for FRASER2 calibration | | Batch match | Strongly recommended | Reduces autoencoder confounding |

Related Skills

• splice-variant-prediction - SpliceAI / Pangolin for in-silico prediction; integration target
• differential-splicing - When testing multiple patients vs controls (>=10 vs cohort)
• splicing-qc - Library / depth / tissue prerequisites
• variant-calling/clinical-interpretation - ACMG/AMP framework integration
• workflows/clinical-trial-pipeline - Trial-grade RNA-seq diagnostics

References

• Mertes et al 2021 Nat Commun - FRASER 1.x
• Scheller et al 2023 Am J Hum Genet - FRASER 2.0 (Intron Jaccard Index)
• Brechtmann et al 2018 Am J Hum Genet - OUTRIDER
• Jenkinson et al 2020 Bioinformatics - LeafcutterMD
• Yepez et al 2021 Nat Protocols - DROP pipeline
• Cummings et al 2017 Sci Transl Med - RNA-seq for muscular dystrophy diagnostics
• Kremer et al 2017 Nat Commun - RNA-seq for mitochondrial disease
• Brown et al 2022 Nature - UNC13A cryptic exon (TDP-43 / ALS)
• Klim et al 2019 Nat Neurosci - STMN2 cryptic splicing (ALS)
• Darman et al 2015 Cell Rep - SF3B1 cryptic 3'ss
• Walker et al 2023 Am J Hum Genet - ClinGen SVI splicing recommendations