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

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

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

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

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.