Bio Single Cell Splicing

Version Compatibility

Reference examples tested with: MARVEL 2.0+, BRIE2 0.2.4+, scQuint 0.1+, SpliZ 0.0.1+, Sierra 1.0+, Psix 0.1+, anndata 0.10+, scanpy 1.10+, pandas 2.2+, scipy 1.13+

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.

Single-Cell Splicing Analysis

The fundamental decision is chemistry, not tool. Most droplet 3' scRNA-seq cannot support transcriptome-wide splicing inference because reverse transcription primes from the poly(A) tail and most reads land in the 3' UTR — far from CDS-region splicing events. Plate-based full-length methods and single-cell long-read sequencing are the chemistries that give per-cell isoform structure across the gene body.

The 10X 3' Problem (Quantified)

Three compounding mechanisms make 10X Chromium 3' (v3.1, GEM-X, v4) hostile to splicing:

1. 3' enrichment: median fragment <1 kb from poly(A); >70% of unique reads fall within 3' UTR.
2. Short R2 (~91 nt): each read straddles at most one junction; usually none, because R2 lands in 3' UTR.
3. PCR concatemers and TSO artifacts: pollute junction detection; UMI collapse is gene-level, not isoform-level.

Quantitative estimate: Only a small fraction of cassette exons sit close enough to the polyA site to be sampled by 3' chemistry (empirical estimates from APA/3'-end atlases — see Tian & Manley 2017 Nat Rev Mol Cell Biol for the 3' UTR isoform landscape). Effective junction read yield from 10X 3' is <0.1 per cell per AS event — vs the 5-10 needed for stable per-cell PSI. Most splicing analyses on 10X 3' data report artifacts.

The 5' kit (10X 5' GEX) does not solve this — it shifts capture from 3' UTR to 5' UTR / TSS-proximal regions. Marginal improvement; not a transcriptome-wide solution. Note that V(D)J recovery requires the 10X Chromium Single Cell Immune Profiling kit (with TCR/BCR-specific enrichment), not 5' GEX alone — postdocs designing immune-repertoire experiments must use the dedicated V(D)J kit.

Decision: Does the Chemistry Support Splicing Analysis?

Tool Selection Matrix

Decision Tree by Goal

MARVEL Plate-Based Workflow

Goal: Run a unified workflow from STAR junctions to cell-type-specific splicing calls.

Approach: Build a wide splice-junction count matrix (rows = junctions keyed by coord.intron, columns = cells), assemble per-event feature tables, then construct MARVEL object with named slots (SpliceJunction, SplicePheno, SpliceFeature, IntronCounts, GeneFeature, Exp, GTF). Quantify PSI per event class, classify modality, test differential splicing.

library(MARVEL); library(Seurat); library(data.table)

seurat_obj <- readRDS('cells.rds')

# Build wide SJ matrix: first column 'coord.intron' (e.g. 'chr1:100007082:100022621'),
# subsequent columns are per-cell sample IDs with junction counts as values.
# This is constructed from STAR SJ.out.tab files (one per cell) merged on intron coord.
sj_files <- list.files('star_pass2/', pattern='SJ.out.tab$', full.names=TRUE)
sj_long <- rbindlist(lapply(sj_files, function(f) {
    d <- fread(f, sep='\t', header=FALSE,
               col.names=c('chr','start','end','strand','motif','annot','unique','multi','overhang'))
    d$coord.intron <- paste(d$chr, d$start, d$end, sep=':')
    d$sample <- gsub('_SJ.out.tab$', '', basename(f))
    d[, .(coord.intron, sample, unique)]
}))
sj <- dcast(sj_long, coord.intron ~ sample, value.var='unique', fill=0)

# SpliceFeature is a NAMED LIST keyed by event class
df.feature.list <- list(
    SE   = read.table('events_SE.txt',   header=TRUE, sep='\t'),
    A5SS = read.table('events_A5SS.txt', header=TRUE, sep='\t'),
    A3SS = read.table('events_A3SS.txt', header=TRUE, sep='\t'),
    MXE  = read.table('events_MXE.txt',  header=TRUE, sep='\t'),
    RI   = read.table('events_RI.txt',   header=TRUE, sep='\t')
)

# SplicePheno: per-cell metadata; sample.id column maps to SpliceJunction column names
df.pheno <- seurat_obj@meta.data
df.pheno$sample.id <- rownames(df.pheno)

marvel <- CreateMarvelObject(
    SpliceJunction = sj,
    SplicePheno    = df.pheno,
    SpliceFeature  = df.feature.list,
    GeneFeature    = read.table('gene_features.tsv', header=TRUE, sep='\t'),
    Exp            = read.table('tpm.tsv', header=TRUE, sep='\t', row.names=1),
    GTF            = rtracklayer::import('annotation.gtf')
)

marvel <- ComputePSI(marvel, CoverageThreshold=10, EventType='SE')
marvel <- AssignModality(marvel, EventType='SE')
marvel <- CompareValues(
    marvel,
    cell.group.g1 = neurons, cell.group.g2 = glia,
    method = 'wilcox', n.cells = 25, psi.delta = 0.1
)

For 10X droplet data, MARVEL v2+ provides CreateMarvelObject.10x() and AnnotateSJ.10x() constructors. Verify the exact API via ?CreateMarvelObject.10x in installed MARVEL.

MARVEL classifies events into modalities (Song 2017 Mol Cell): included (PSI1), excluded (PSI0), bimodal (mixture at 0/1), middle (peaked ~0.5), multimodal. Bimodality usually reflects mixed cell states or stochastic monoallelic-like bursting. Mid-modality (peaked at 0.5) can be technical (mixed cells in a droplet) — confirm with full-length data.

BRIE2 Bayesian PSI

Goal: Estimate per-cell PSI with informative regulatory-feature prior; test cell-state association via likelihood-ratio testing on covariate effects.

Approach: BRIE2 is a CLI-driven workflow (brie-count for read counting, brie-quant for variational inference + LRT). Prepare a GFF3 of splicing events, count cell-barcoded junction reads, then fit the model with covariate testing.

# 1. Count splicing events per cell
brie-count \
    -a splicing_events.gff3 \
    -S sample_list.tsv \
    -o brie_counts/ \
    -p 16

# 2. Fit BRIE2 with LRT against the cell-type covariate
brie-quant \
    -i brie_counts/brie_count.h5ad \
    -c cell_metadata.tsv \
    -o brie_quant.h5ad \
    --interceptMode gene \
    --LRTindex All \
    --testBase null \
    --MCsize 3 \
    --batchSize 1000000 \
    -p 16

--interceptMode gene fits a gene-specific intercept (recommended); --LRTindex All tests all covariates; --testBase null uses the null model as the LRT reference. Verify exact flag set via brie-quant -h in installed BRIE2.

import scanpy as sc

adata_splice = sc.read_h5ad('brie_quant.h5ad')
# Per-event covariate effects, ELBO values, and LRT statistics live in
# adata_splice.varm and adata_splice.var; column names depend on BRIE2 version.
# Inspect with: print(adata_splice); print(adata_splice.varm.keys())
# Per-event significance is typically derived from LRT delta-ELBO.

BRIE2 (Huang & Sanguinetti 2021 Genome Biol) uses a sequence-derived feature prior (exon length, GC content, splice site strength, motif counts) to regularize PSI estimates in low-coverage cells. The LRT-based covariate test answers "is this event associated with cell state?" without requiring per-cell PSI accuracy. Threshold the delta-ELBO at ~3 (analogous to log-Bayes-factor); confirm against version-specific output keys via the brie-tutorials repo.

SpliZ for Annotation-Free Discovery

Goal: Identify splicing-defined cell populations without an event database.

Approach: Compute per-gene splicing Z-score across cells; test for cell-state association via permutation.

# SpliZ is a Nextflow pipeline (not a standalone CLI). Configure inputs in a .config
# file (dataname, input_file, libraryType, grouping_level_1/2) - either SICILIAN
# output (SICILIAN=true) or BAMs via a samplesheet CSV + metadata + GTF (SICILIAN=false).
nextflow run salzmanlab/spliz -r main -latest -c spliz.config

SpliZ (Olivieri 2022 Nat Methods) is robust to dropout because it pools junction information across the gene; particularly useful for discovering splicing diversity in heterogeneous tumor samples.

Psix for Regulated AS Along Trajectories

Goal: Detect AS that varies coherently with cell state along a developmental trajectory, robust to dropout.

Approach: Score whether observed PSI is smooth on the cell-cell kNN graph from expression-space embedding.

import psix
import scanpy as sc

adata = sc.read_h5ad('cells.h5ad')
sc.pp.neighbors(adata, n_neighbors=30, use_rep='X_pca')

psix_obj = psix.Psix(adata, psi_matrix_path='psi_matrix.tsv')
psix_obj.run_psix()

regulated = psix_obj.psix_results.query('psix_score > 1.5 and pvalue < 0.05')

Psix (Buen Abad Najar 2022 Genome Res 32:1385) is the principled alternative to imputing PSI: do not impute (it obliterates heterogeneity); test for graph smoothness instead.

Sierra for APA (Not Splicing)

Goal: Detect alternative polyadenylation in 10X 3' data — frequently confounded with AS.

Approach: Peak-call read pile-ups at 3' ends, then DEXSeq-style DTU on 3' UTR isoforms.

library(Sierra)

peak_file <- FindPeaks(
    output.file = 'peaks.txt',
    gtf.file = 'annotation.gtf',
    bam.file = 'possorted_genome_bam.bam'
)

counts <- CountPeaks(
    peak.sites.file = 'peaks.txt',
    gtf.file = 'annotation.gtf',
    bamfile = 'possorted_genome_bam.bam',
    whitelist.file = 'barcodes.tsv'
)

# CountPeaks returns a peak x cell matrix; annotate it and build a peak Seurat
# object before differential-usage testing.
peak.annotations <- AnnotatePeaksFromGTF(
    peak.sites.file = 'peaks.txt',
    gtf.file = 'annotation.gtf',
    output.file = 'peak_annotations.txt'
)

peaks.seurat <- NewPeakSeurat(
    peak.data = counts,
    annot.info = peak.annotations,
    cell.idents = cell_identities
)

apa_results <- DUTest(peaks.seurat, population.1 = ctrl_cells, population.2 = trt_cells)

If only 10X 3' data is available, this is often what is actually wanted. Distinct UTRs change miRNA targeting, RBP binding, and stability — biologically meaningful but not splicing.

Pseudobulk for Statistical Power

Goal: Recover bulk-level statistical power for differential splicing between cell types.

Approach: Sum junction counts across cells of the same cluster, then run leafcutter / rMATS on aggregated counts.

import pandas as pd
import numpy as np

def pseudobulk_junctions(junction_counts, cell_metadata, groupby='cell_type'):
    out = {}
    for group, cells in cell_metadata.groupby(groupby).groups.items():
        mask = junction_counts.columns.isin(cells)
        out[group] = junction_counts.loc[:, mask].sum(axis=1)
    return pd.DataFrame(out)

Use pseudobulk for differential splicing between well-defined cell types; use per-cell methods for within-population heterogeneity (graded splicing along pseudotime, bimodal cell-state mixtures).

Single-Cell Long-Read = Future of Single-Cell Splicing

In 2024-2026, full-length single-cell long-read sequencing has become practical and is the recommended chemistry for splicing-focused single-cell experiments:

• MAS-Iso-seq / PacBio Kinnex: concatenated full-length cDNA arrays, ~16x throughput vs plain Iso-Seq, compatible with 10X 5' libraries (Al'Khafaji 2024 Nat Biotech)
• scISOr-Seq2: hybrid 10X + PacBio for cell typing + isoform structure (Joglekar et al 2024 Nat Neurosci 27:1051-1063, single-cell long-read brain isoform mapping)
• ONT direct cDNA + 10X: lower cost, similar information content
• FLAMES: barcode demultiplexing + isoform quantification + SNV calling for ONT scRNA (Tian 2021 Genome Biol 22:310)

For splicing-specific full-length single-cell analysis, see long-read-splicing skill.

Per-Tool Failure Modes

MARVEL: SpliceJunction Matrix Format

Trigger: Building the SpliceJunction matrix from STAR SJ.out.tab incorrectly (e.g. long-format instead of wide).

Mechanism: MARVEL plate-based CreateMarvelObject(SpliceJunction = ...) expects a wide matrix with first column coord.intron (formatted chr:start:end) and subsequent columns being per-cell sample IDs with integer junction counts. Long-format data.frames or missing coord.intron column cause runtime errors.

Symptom: "no coord.intron column found" errors; or empty PSI tables despite junction reads being present.

Fix: Verify wide-matrix structure; ensure SJ.out.tabs are merged on the chr:start:end key with cells as columns. Use data.table::dcast for the long->wide reshape.

BRIE2: TensorFlow Memory

Trigger: Large cohort (>10k cells) with deep coverage.

Mechanism: Variational inference loads full count matrix; TensorFlow allocates GPU memory aggressively.

Symptom: OOM kills; training stalls.

Fix: Reduce --batchSize from default (500000) to 100000 or 50000; train per-chromosome batch; use CPU mode for very small cohorts. Note flag is camelCase --batchSize, not --batch_size.