Python API¶
Generated from the docstrings in src/splicescope. The modules are listed in the order
the pipeline uses them; each takes and returns pandas data frames, so any step can be
run on its own or replaced. Private helpers (names starting with _) are omitted.
Input: splicescope.io¶
Reading splice junctions from STAR and reference annotations from GTF.
Coordinate convention (matches STAR SJ.out.tab)
A junction is the intron it spans: start is the first intronic base and
end the last intronic base, both 1-based inclusive. The 5' splice site
(donor) and 3' splice site (acceptor) therefore depend on strand:
+strand: donor =start, acceptor =end-strand: donor =end, acceptor =start
Every junction is keyed by (chrom, start, end, strand).
sample_name_from_path ¶
sample_name_from_path(path: str | Path) -> str
Derive a sample name from an SJ.out.tab filename.
A1_SJ.out.tab, A1.SJ.out.tab, A1SJ.out.tab and A1.tab all give A1.
open_text ¶
open_text(path: str | Path)
Open a plain or gzipped text file. Reference annotations ship as .gz.
donor_acceptor ¶
donor_acceptor(start: int, end: int, strand: str) -> tuple[int, int]
Return (donor_pos, acceptor_pos) for an intron given its strand.
read_star_sj ¶
read_star_sj(path: str | Path, sample: str | None = None) -> pd.DataFrame
Read a STAR SJ.out.tab file into a tidy junction table.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
path
|
str or Path
|
Path to |
required |
sample
|
str
|
Sample label; defaults to the name |
None
|
Returns:
| Type | Description |
|---|---|
DataFrame
|
Columns |
read_many_star_sj ¶
read_many_star_sj(paths: dict[str, str | Path]) -> pd.DataFrame
Read several SJ.out.tab files into one long table.
paths maps sample name -> file path.
read_gtf_junctions ¶
read_gtf_junctions(path: str | Path) -> pd.DataFrame
Derive the set of known introns (junctions) from a GTF annotation.
Exons are grouped per transcript, sorted, and the gaps between consecutive
exons become known introns. The result is a DataFrame with columns
[chrom, start, end, strand, gene_id] (one row per unique junction).
The file may be plain or gzipped; GENCODE and Ensembl both ship .gtf.gz.
read_gmt ¶
read_gmt(path: str | Path) -> dict[str, list[str]]
Read a GMT gene-set file into {term: [genes]}.
GMT lines are tab-separated: term<TAB>description<TAB>gene1<TAB>gene2....
The description column is ignored. Compatible with MSigDB / GO / KEGG exports.
Junction taxonomy: splicescope.annotate¶
Classify observed junctions against a reference annotation.
Each observed junction is labelled by how it relates to the known splice sites:
====================== ==================================================
class meaning
====================== ==================================================
annotated donor+acceptor pair is a known intron
novel_combination both sites known, but not as a pair (e.g. exon skip)
novel_donor acceptor known, 5' site novel
novel_acceptor donor known, 3' site novel
cryptic both sites novel — candidate cryptic splicing
====================== ==================================================
classify_one ¶
classify_one(chrom, start, end, strand, junctions, donors, acceptors) -> str
Classify a single junction. Pure function for easy testing.
resolve_unstranded ¶
resolve_unstranded(observed: DataFrame, known: DataFrame) -> pd.DataFrame
Give strand-undefined junctions the strand their annotation implies.
STAR writes strand code 0 — which this package reads as "." — whenever it cannot
infer the strand from the intron motif, which is routine for non-canonical junctions.
Such a junction can never equal a stranded annotation, so it is classified cryptic
however ordinary it is, loses its gene_id, and — being the only junction at its own
(chrom, position, ".") site — is handed Ψ ≡ 1.0 in every sample. It then occupies a
row in the differential table that can never show a difference.
Where the annotation knows the intron, or either of its splice sites, the strand is
recoverable. Junctions it cannot place are left as ".": the strand genuinely is
unknown, and guessing one would be worse than saying so.
annotate_junctions ¶
annotate_junctions(observed: DataFrame, known: DataFrame) -> pd.DataFrame
Annotate a table of observed junctions.
observed needs columns [chrom, start, end, strand] (extra columns are
preserved). Returns a copy with two new columns: sclass (one of
:data:CLASSES) and gene_id (best-effort assignment by shared splice site).
Strand-undefined junctions are first placed against the annotation where they can
be — see :func:resolve_unstranded.
annotation_summary ¶
annotation_summary(annotated: DataFrame) -> pd.DataFrame
Counts of unique junctions per class (deduplicated across samples).
Splice-site usage: splicescope.quantify¶
Quantify splice-site usage (Ψ) from junction read counts.
We use a well-defined, model-free metric: donor usage. For a junction j with 5' splice site (donor) D, its usage in a sample is
Ψ(j) = count(j) / Σ count(all junctions sharing donor D)
i.e. the fraction of that donor's spliced reads that flow through j. This is the quantity differential-splicing methods build on and needs no event model. Acceptor usage is defined symmetrically and returned alongside.
A minimum coverage threshold avoids dividing tiny counts: sites with fewer than
min_reads total reads in a sample yield NaN (uninformative) rather than a
noisy ratio.
add_unobserved_zeros ¶
add_unobserved_zeros(df: DataFrame) -> pd.DataFrame
Make the zeros an aligner reports by omission explicit.
STAR and its kin list only junctions carrying at least one read, so a junction missing from a sample arrives looking like missing data when it is in fact a measurement: zero reads out of a denominator that is known. For cryptic splicing that is the strongest evidence there is — off in every control, on in every knockdown — and treating it as missing drops exactly the events the tool exists to find.
A zero row is added only where the junction's donor or acceptor site has reads in
that sample, so Ψ has a denominator to be a fraction of; a site nobody sequenced
stays genuinely unobserved. Requires donor and acceptor columns.
compute_psi ¶
compute_psi(annotated: DataFrame, min_reads: int = 10, fill_unobserved: bool = True) -> pd.DataFrame
Add psi_donor and psi_acceptor columns (per sample) to junctions.
annotated must have [chrom, start, end, strand, sample, count]. The
matching donor_total and acceptor_total denominators are kept so a
count-based test can use the reads Ψ was computed from, not just the ratio.
fill_unobserved adds the zero-count rows the aligner left out — see
:func:add_unobserved_zeros. Turn it off only to reproduce the pre-0.9.0
behaviour, which silently discarded junctions absent from a whole group.
psi_matrix ¶
psi_matrix(psi_df: DataFrame, value: str = 'psi_donor') -> pd.DataFrame
Pivot to a junction × sample matrix of Ψ values (junctions as the index).
Events: splicescope.events¶
Event-level splicing: cassette (skipped) exons and their PSI.
Junction usage (see :mod:splicescope.quantify) is model-free but not an event.
Here we reconstruct the canonical cassette-exon event directly from junctions
and compute the standard percent-spliced-in (Ψ).
A cassette exon is defined by three junctions that share endpoints:
inclusion 1 inclusion 2
┌───────────────┐ ┌───────────────┐
[ up-exon ] [ cassette exon ] [ down-exon ]
└──────────────────────────────────────┘
skipping junction
- inclusion 1
(s3, e1)shares the start of the skip and lands on the cassette exon's 5′ boundary, - inclusion 2
(s2, e3)shares the end of the skip and leaves the exon's 3′ boundary, withe1 < s2(the exon is[e1+1, s2-1]), - skipping
(s3, e3)splices straight from up-exon to down-exon.
This coordinate pattern is strand-independent (the strand only has to match), so detection works for both strands. PSI follows the junction-count definition used by rMATS-style methods:
inclusion = (count(inc1) + count(inc2)) / 2
Ψ = inclusion / (inclusion + count(skip))
event_columns ¶
event_columns(types: Sequence[str] = EVENT_TYPES) -> list[str]
The columns a table of these event types has, found or not.
A5SS and A3SS describe a splice site the same way, so asking for both adds
one set of columns, not two.
detect_cassette_events ¶
detect_cassette_events(annotated: DataFrame) -> pd.DataFrame
Find cassette-exon events from a table of (annotated) junctions.
Returns one row per event with the three junction coordinate pairs, the
inferred cassette-exon interval, and a best-effort gene_id.
detect_mxe_events ¶
detect_mxe_events(annotated: DataFrame, max_exon: int = 1000, exclude: set | None = None, max_candidates: int = 50) -> pd.DataFrame
Detect mutually-exclusive-exon (MXE) events.
An MXE is two non-overlapping exons, A (upstream) and B (downstream), that sit
between the same upstream donor s and downstream acceptor e — each
reached by its own pair of inclusion junctions, and normally never included
together. PSI is exon A's share:
Ψ = incl(A) / (incl(A) + incl(B))
where incl(X) = (count(X_5' junction) + count(X_3' junction)) / 2.
max_candidates bounds how many distinct candidate exons one (donor, acceptor)
anchor may hold before it is treated as noise and skipped with a warning; pair
enumeration is quadratic in that number and the candidates are themselves every
junction combination at the anchor.
detect_alt_ss_events ¶
detect_alt_ss_events(annotated: DataFrame, kind: str, exclude: set | None = None) -> pd.DataFrame
Detect alternative 5′ (kind="A5SS") or 3′ (kind="A3SS") splice-site events.
An A5SS event is a splice acceptor served by two or more donors; A3SS is the mirror (one donor, several acceptors). The inclusion isoform is the one whose variable site sits closest to the shared site (the longer exon), following the rMATS convention. One event is emitted per shared site with ≥2 alternatives.
exclude is a set of (chrom, start, end, strand) junctions already explained
by a cassette or MXE event. A site whose junctions are all explained that way is
skipped, since reporting it here would be the same signal told twice. Excluded
junctions still count towards the alternatives at a site they share with a junction
that is not explained: dropping them outright deleted genuine alternative-splice-site
events whose only fault was sharing an anchor with a cassette exon.
detect_events ¶
detect_events(annotated: DataFrame, types: tuple[str, ...] = EVENT_TYPES, max_exon: int = 1000, max_candidates: int = 50) -> pd.DataFrame
Detect all requested event types and return them in one typed table.
Cassette (SE) events are detected first; the junctions they use are then excluded from A5SS/A3SS detection so the same signal is not double-reported.
max_exon bounds how long a candidate MXE exon may be and max_candidates how
many may share one anchor; raise them for loci with unusually long or unusually
numerous mutually exclusive exons.
event_psi ¶
event_psi(annotated: DataFrame, events: DataFrame, min_reads: int = 10) -> pd.DataFrame
Per-sample PSI for every event, dispatched by event_type.
Returns a long table keyed by event_id (feed to
:func:splicescope.diff.differential_splicing with key=["event_id"]).
Differential splicing: splicescope.diff¶
Differential splicing between two groups of samples.
Two tests are available, and the difference between them decides whether the tool can find anything at all on a real dataset.
betabinom (default where read counts are available)
A beta-binomial likelihood-ratio test on the inclusion and total counts Ψ
was computed from. Evidence accumulates with coverage, so a well-covered
event can reach very small p-values even with three replicates per group.
ranksum (legacy, used only when counts are absent)
Two-sided Mann-Whitney U on the per-sample Ψ values. This discards the
counts, and a rank test on n replicates per group cannot return a
p-value below 2 / C(2n, n) — 0.1 for a 3-vs-3 design. After
Benjamini-Hochberg across a genome's worth of junctions nothing can ever be
called significant, whatever the effect size. Keep it only for Ψ tables that
carry no counts.
Both report delta_psi (group B minus group A), pvalue and BH qvalue.
benjamini_hochberg ¶
benjamini_hochberg(pvals: ndarray) -> np.ndarray
Return BH-adjusted q-values for a 1-D array of p-values (NaNs preserved).
differential_splicing ¶
differential_splicing(psi_df: DataFrame, groups: dict[str, str], value: str = 'psi_donor', min_samples: int = 2, key: list[str] | None = None, test: str = 'auto', inc_col: str | None = None, total_col: str | None = None, dispersion: str = 'shared', filter_invariant: bool = False) -> pd.DataFrame
Test each junction (or event) for differential Ψ between two conditions.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
psi_df
|
DataFrame
|
A long Ψ table (e.g. from :func: |
required |
groups
|
dict
|
Mapping |
required |
value
|
str
|
Which Ψ column to test. |
'psi_donor'
|
min_samples
|
int
|
Minimum informative (non-NaN) samples required per group. |
2
|
key
|
list of str
|
Grouping columns identifying a testable unit; defaults to the junction
coordinates |
None
|
test
|
('auto', 'betabinom', 'ranksum')
|
|
"auto"
|
inc_col
|
str
|
Count columns to model; inferred from |
None
|
total_col
|
str
|
Count columns to model; inferred from |
None
|
dispersion
|
('shared', 'per_unit_floor')
|
|
"shared"
|
filter_invariant
|
bool
|
Set aside units whose Ψ is identical in every sample before correcting for
multiple testing. Neither test finds evidence in such a unit
whichever way the samples are labelled, so it cannot be rejected at any
threshold — it is not a test that came out negative but one that could not be
run. Most are structural: a constitutive donor with a single junction has Ψ ≡ 1
because the site has nothing else to splice to, and such sites outnumber the
alternative ones across a genome. Left in, they inflate the Benjamini–Hochberg
denominator and make every real unit's q-value worse: 48 % of the tested units
on a 200-gene simulation, costing a factor of 1.91 in q. This is the degenerate
case of independent filtering (Bourgon, Gentleman & Huber, PNAS 2010) — the
criterion is the spread over all samples, which never looks at the group labels,
and zero spread forces p = 1 whatever the labels are. Off by default because it
changes every q-value in a run; a warning says how many units it set aside.
Reproduce with |
False
|
Returns:
| Type | Description |
|---|---|
DataFrame
|
One row per unit, sorted by q-value. |
significant ¶
significant(diff: DataFrame, q: float = 0.05, min_delta: float = 0.1) -> pd.DataFrame
Filter a differential table to confident hits (qvalue and |ΔΨ|).
Beta-binomial model: splicescope.betabinom¶
Beta-binomial machinery for count-based differential splicing.
Why not a rank test
Comparing per-sample Ψ with Mann-Whitney throws away the read counts that Ψ was
computed from, and a rank test on n replicates per group cannot return a
p-value below 2 / C(2n, n) — 0.1 for the 3-vs-3 designs that dominate
RNA-seq. Nothing survives genome-wide multiple testing, however large the
effect. Modelling the counts removes that floor: evidence grows with coverage,
not only with replicate count.
The model
Inclusion reads k out of n informative reads are beta-binomial with mean
Ψ and precision s = α + β, so that α = Ψs and β = (1 - Ψ)s. The beta
layer absorbs the replicate-to-replicate variability a plain binomial would
mistake for signal. s is shared across events and estimated once by profile
likelihood; per event we then fit Ψ under
H0— one Ψ for every sample, andH1— a separate Ψ per group,
and compare 2(ℓ₁ − ℓ₀) to a chi-square with one degree of freedom.
Everything is vectorised over events: the fits are bisections on a monotone score function, so a few hundred thousand events cost a handful of seconds.
loglik ¶
loglik(k: ndarray, n: ndarray, mask: ndarray, mu: ndarray, s: float) -> np.ndarray
Per-event beta-binomial log-likelihood, up to a constant.
The binomial coefficient is dropped: it does not depend on mu and
cancels in the likelihood ratio.
fit_mu ¶
fit_mu(k: ndarray, n: ndarray, mask: ndarray, s: float) -> np.ndarray
Maximum-likelihood Ψ per event, by vectorised bisection on the score.
The score is decreasing in mu, so a sign change brackets the root. Events
whose score never changes sign are at a boundary (all reads included or all
skipped) and are returned as 0 or 1.
dispersion_has_information ¶
dispersion_has_information(k: ndarray, n: ndarray, mask: ndarray, groups: list[ndarray] | None = None) -> bool
Whether the data carry any information about dispersion.
:func:dispersion_is_estimable asks the same question of the design. It is given
no counts, so it can only weigh informative samples against fitted group means — and
that is not enough. A table of constitutive splice sites passes it with replicates to
spare, because every group is fully covered, and still leaves nothing to measure
dispersion from, because every group sits at Ψ = 1. :func:estimate_precision then
returns its ceiling, which asserts no overdispersion rather than admitting
ignorance, and narrows the test to a binomial one.
Note what this does not report. Data genuinely tighter than binomial also send the
estimate to the ceiling, and that is a measurement rather than the absence of one, so
it is not a failure and this still returns True.
dispersion_is_estimable ¶
dispersion_is_estimable(n: ndarray, mask: ndarray, groups: list[ndarray] | None = None) -> bool
Whether the design carries any information about replicate-to-replicate scatter.
Dispersion is measured by how far replicates fall from their own group's mean, so a
group holding a single informative sample contributes a residual of exactly zero and
says nothing. When no unit anywhere has more informative samples than fitted group
means there is no residual degree of freedom at all, and :func:estimate_precision
falls back to max_precision — which asserts no overdispersion rather than
unknown dispersion, narrowing the test to a plain binomial one.
That failure is quiet and severe: a 1-vs-1 comparison of Ψ 0.300 against 0.360 at
1000 reads returns q = 4.5e-03 on no replication whatsoever. Callers should check
this before trusting a p-value.
estimate_precision ¶
estimate_precision(k: ndarray, n: ndarray, mask: ndarray, groups: list[ndarray] | None = None, min_precision: float = 1.0, max_precision: float = 100000.0) -> float
Shared precision s, from df-corrected Pearson residuals.
Profile likelihood is not used here: fitting one Ψ per event on a handful of
replicates absorbs part of the very variance we are trying to measure, which
biases s upward and makes the test anti-conservative. The moment
estimator below corrects for the fitted means through the residual degrees
of freedom instead.
Under the model E[r²] = 1 + (n - 1) / (s + 1) for Pearson residuals
r, so a df-corrected residual sum determines s in one pass.
A single shared s is a deliberate simplification, and it is not uniformly
conservative. When dispersion is homogeneous it is well calibrated (simulated at
s = 50: estimate 52.3, false-positive rate 0.051 against a nominal 0.05). When
dispersion is heterogeneous the shared value sits between the extremes and the test
is anti-conservative for everything noisier than it: with half the units at s = 200
and half at s = 5 the shared estimate is 10.7, and the loosely dispersed half runs
at a false-positive rate of 0.173 while the tight half runs at 0.000. The overall
rate of 0.087 hides both. :func:estimate_precision_per_unit and
differential_splicing(dispersion="per_unit_floor") trade power for protection
against this; see docs/METHODS.md §5.4.
groups must be the same sample masks the test will use. Residuals are
then taken about each group's own mean: pooling groups that genuinely differ
would charge that difference to dispersion and destroy the power the test is
supposed to have.
A group whose fitted mean sits exactly at 0 or 1 is excluded — from the residual sum
and from the degrees of freedom both. It has every read on one side, so its residuals
are zero by construction and it carries no information about replicate-to-replicate
variability. Most such groups are structural rather than biological: a constitutive
donor with a single junction has Ψ = 1 in every sample because the site has nothing
else to splice to, and across a genome those sites outnumber the alternative ones.
Counting them on one side only drove the estimate from a true 50 to 842 and the
false-positive rate from 0.049 to 0.173 — see validation/invariant_units.py and
docs/METHODS.md §5.5. The exclusion is per group, so a junction switched fully on
in one condition still contributes what the other condition's replicates do, which is
the shape of a real cryptic event.
When the design cannot support an estimate at all this returns max_precision,
which asserts no overdispersion rather than admitting ignorance — check
:func:dispersion_is_estimable first. differential_splicing does.
estimate_precision_per_unit ¶
estimate_precision_per_unit(k: ndarray, n: ndarray, mask: ndarray, groups: list[ndarray] | None = None, min_precision: float = 1.0, max_precision: float = 100000.0) -> np.ndarray
The same moment estimator applied to each unit on its own.
One unit's handful of replicates gives a very noisy estimate, so this is not a replacement for the shared value — it is the input to a floor. A unit that looks more dispersed than the shared estimate says is the case the shared value gets wrong, and taking the smaller of the two can only widen the null.
Units that cannot be estimated return max_precision, which leaves the shared
value in charge of them.
lrt ¶
lrt(k: ndarray, n: ndarray, mask: ndarray, group_a: ndarray, group_b: ndarray, s: float) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]
Likelihood-ratio test of "one Ψ" against "one Ψ per group".
group_a and group_b are boolean masks over the sample axis. Returns
(mu_a, mu_b, statistic, pvalue); the statistic is compared to a
chi-square with one degree of freedom.
Samples in neither group take no part. Both hypotheses must be fitted to the same observations or the ratio is not a likelihood ratio: leaving them in the null while the alternative sees only the two groups charges their entire likelihood to the null and inflates the statistic without bound — with six ungrouped samples alongside a 3-vs-3 comparison, p went from 1.6e-09 to 9.7e-93 on unchanged group data.
min_achievable_rank_pvalue ¶
min_achievable_rank_pvalue(n_a: int, n_b: int) -> float
Smallest two-sided Mann-Whitney p-value possible for these group sizes.
Provided so the limitation this module exists to fix can be shown rather than asserted: with 3 against 3 the floor is 0.1, so no rank test can clear genome-wide correction no matter how large the effect.
This is the exact test's floor, and not what scipy.stats.mannwhitneyu returns
once there are ties: it then uses the normal approximation, which gives 0.081 for an
untied 3 against 3, 0.064 with one group tied, and 0.047 when Psi is 0 in every
control and 1 in every knockdown — the signature this module exists to find, under
the nominal 0.05. What stops the rank test calling it is Benjamini-Hochberg across a
genome's worth of junctions, not the floor. See METHODS section 5.1.
Protein consequence: splicescope.consequence¶
Predict what a cryptic exon does to the protein its host transcript encodes.
Detecting a cryptic exon says nothing about whether it matters. The same event can be silently tolerated, shift the reading frame, or introduce a premature termination codon (PTC) that sends the transcript to nonsense-mediated decay -- which is how TDP-43 cryptic exons deplete proteins such as STMN2 and UNC13A. This module resolves that, one event at a time:
- locate the host transcripts whose intron contains the event,
- inherit the reading frame from the coding sequence upstream of the intron,
- translate the inserted sequence in that frame and look for a stop codon,
- apply the 50-nucleotide rule to decide whether a PTC triggers NMD.
Sequence comes from an indexed genome FASTA read directly through its .fai,
so nothing here needs a compiled sequence library.
Consequence classes
====================== ==================================================
class meaning
====================== ==================================================
ptc_nmd PTC introduced, far enough upstream to trigger NMD
ptc_escape PTC introduced but predicted to escape NMD
frameshift insert length not a multiple of 3, no PTC inside it
in_frame_insertion multiple of 3 and no stop codon: protein extended
utr_insertion event falls outside the coding sequence
non_coding_host host transcript has no annotated CDS
exon_truncation a novel splice site shortens an annotated exon
no_host_transcript no annotated intron contains the event
====================== ==================================================
Two shapes of event are handled. A cassette cryptic exon is given directly as an exon interval. A splice-site shift — a novel donor or acceptor paired with an annotated site on the other side — is given as a junction, and the sequence it adds to or removes from the neighbouring exon is derived from the annotation. Splice-site shifts are the majority of cryptic events reported by junction-level callers, so a layer that only understood cassettes would miss most of them.
GenomeFasta ¶
Random access to an indexed FASTA using only its .fai sidecar.
The index gives, per contig, the byte offset of its first base plus how many bases and how many bytes each line holds, which is enough to seek straight to any coordinate without reading the file.
fetch ¶
fetch(chrom: str, start: int, end: int, strand: str = '+') -> str
Sequence for chrom:start-end, 1-based inclusive, strand-aware.
Transcript
dataclass
¶
Exon and CDS structure of one transcript, in transcription order.
introns
property
¶
introns: list[tuple[int, int]]
Introns in transcription order, 1-based inclusive.
frame_at ¶
frame_at(position: int) -> int
How many bases of the codon spanning position already lie upstream.
cds_length_before ¶
cds_length_before(position: int) -> int
Length of coding sequence transcribed before position.
Consequence
dataclass
¶
Predicted protein-level effect of including one cryptic exon.
load_transcripts ¶
load_transcripts(gtf_path: str | Path, genes: set[str] | None = None) -> dict[str, Transcript]
Read exon and CDS records from a GTF into transcript models.
The file may be plain or gzipped. genes restricts parsing to a set of
gene names or IDs, which matters:
a full GENCODE GTF holds millions of records and a cryptic exon report
normally concerns a few hundred genes.
index_by_gene ¶
index_by_gene(transcripts: dict[str, Transcript]) -> dict[str, list[Transcript]]
Group transcripts by gene, under both the versioned and bare gene id.
host_transcripts ¶
host_transcripts(transcripts: dict[str, Transcript] | list[Transcript], chrom: str, start: int, end: int, strand: str) -> list[Transcript]
Transcripts with an intron fully containing chrom:start-end.
Pass only the candidate transcripts — the host gene's, via
:func:index_by_gene — when the caller knows which gene the event belongs
to; scanning every transcript for every event does not scale.
downstream_sequence ¶
downstream_sequence(tx: Transcript, genome: GenomeFasta, chrom: str, resume: int) -> tuple[str, list[int]]
Mature sequence retained from resume onward, and its exon lengths.
resume is the genomic position at which transcription continues after
the event. Exons are walked in transcription order and clipped to that
point, so the result is the sequence a ribosome would read downstream, with
the length of each remaining exon kept so the last exon-exon junction can be
located.
junction_change ¶
junction_change(tx: Transcript, start: int, end: int) -> tuple[str, int, int] | None
What a novel junction adds to, or removes from, a neighbouring exon.
start and end are the first and last intronic base, 1-based
inclusive. The junction has to share one splice site with an annotated
intron of tx; the other site is the novel one, and the sequence between
the annotated and the novel position is what changes.
Returns ("extension" | "truncation", first_base, last_base), or None
when neither site matches the annotation (an unanchored junction, which the
annotation cannot interpret).
None is also returned when both sites are annotated splice sites of this
transcript but not as a pair: that is an exon skip, not a shifted splice site, and
reading it as one reports the skipped exon and the intron beyond it as a single
contiguous deletion.
junction_hosts ¶
junction_hosts(transcripts: dict[str, Transcript] | list[Transcript], chrom: str, start: int, end: int, strand: str) -> list[tuple[Transcript, str, int, int]]
Transcripts whose annotation can interpret this junction, with the change.
predict_junction_consequence ¶
predict_junction_consequence(tx: Transcript, genome: GenomeFasta, chrom: str, start: int, end: int) -> Consequence | None
Predict the effect of a novel junction, via the exon change it implies.
An extension adds sequence to an exon and is evaluated exactly like a cassette exon. A truncation removes sequence, which cannot introduce a stop inside the removed interval, so only the frame consequence is reported.
find_ptc ¶
find_ptc(sequence: str, frame_offset: int) -> int | None
Offset of the first in-frame stop codon, or None.
frame_offset is how many bases of the codon spanning the exon boundary
already lie upstream, so the first complete codon inside sequence starts
at (3 - frame_offset) % 3.
predict_consequence ¶
predict_consequence(tx: Transcript, genome: GenomeFasta, chrom: str, start: int, end: int, contiguous_downstream: bool = False) -> Consequence
Predict the effect of inserting chrom:start-end into tx.
contiguous_downstream says the sequence is joined to the following exon without
an intervening junction, which is true of an exon extension but not of a cassette
exon. It decides whether that exon can anchor the 50-nucleotide rule.
annotate_junction_consequences ¶
annotate_junction_consequences(events: DataFrame, transcripts: dict[str, Transcript], genome: GenomeFasta, chrom_col: str = 'chrom', start_col: str = 'intron_start', end_col: str = 'intron_end', strand_col: str = 'strand', gene_col: str | None = 'gene_id') -> pd.DataFrame
Add consequence columns to a table of novel junctions.
Coordinates are the first and last intronic base. Each junction is interpreted against its host gene's annotation, and where several transcripts interpret it differently the most disruptive reading is kept.
check_genome ¶
check_genome(events: DataFrame, transcripts: dict[str, Transcript], genome: GenomeFasta, chrom_col: str = 'chrom', end_col: str = 'end') -> None
Refuse a genome the events cannot be read against.
Reading past the end of a contig returns a truncated string rather than an error, so
the wrong assembly does not fail — it finds no stop codon anywhere and calls
everything a frameshift. Measured on the simulator, swapping in a genome 200 bases
long turned six ptc_nmd calls into zero and reported nine events as confidently
as before.
The .fai already carries every contig's name and length, so both halves of the
mistake are cheap to see: a chromosome the FASTA does not have, and a coordinate past
the end of one it does.
Only the chromosomes where a prediction will actually be made are checked — those
with both an event and a transcript to host it. An event on a contig the annotation
says nothing about already has an honest answer (no_host_transcript), and a
genome covering more than the analysis needs is nobody's mistake.
annotate_consequences ¶
annotate_consequences(events: DataFrame, transcripts: dict[str, Transcript], genome: GenomeFasta, chrom_col: str = 'chrom', start_col: str = 'start', end_col: str = 'end', strand_col: str = 'strand', gene_col: str | None = 'gene_id') -> pd.DataFrame
Add consequence columns to a table of cryptic exon candidates.
Where several transcripts host the event the most disruptive prediction is
reported, since that is the one that determines whether protein is lost.
n_host_transcripts records how many were considered.
gene_col names a column holding the event's gene, which restricts the
transcript search to that gene. Without it every transcript is scanned for
every event, which is quadratic and impractical at genome scale.
describe ¶
describe(row) -> str
One sentence saying what a predicted consequence means, for a reader.
The table says ptc_nmd, insert_length=61, distance_to_last_junction=395.
That is the answer, but only to someone who already knows the rule. This spells it
out: what the event adds or removes, whether the reading frame survives, where the
first premature stop falls, and why 50 nucleotides decides the transcript's fate.
Accepts anything with the fields of :class:Consequence — the dataclass itself, or
a row of the table :func:annotate_consequences returns.
Pathway enrichment: splicescope.enrich¶
Gene-set over-representation analysis (ORA) for differentially spliced genes.
After differential splicing, a natural question is which pathways are affected? This module runs the standard hypergeometric (Fisher's exact, one-sided) test used by tools like clusterProfiler: given a set of "hit" genes, a background, and a collection of gene sets (GO terms, KEGG pathways, MSigDB, …), it asks whether each set is over-represented among the hits, and corrects across sets with BH-FDR.
The method is generic — bring any {term: [genes]} mapping (e.g. via
:func:splicescope.io.read_gmt). No gene sets are bundled, because meaningful
enrichment needs real annotations rather than synthetic ones.
normalize_gene_id ¶
normalize_gene_id(gene: str) -> str
Put an identifier in the form both sides of an intersection can agree on.
A GTF gives versioned accessions (ENSG00000141510.16, NM_000546.6); gene-set
files give either symbols or unversioned accessions, so a verbatim comparison
matches nothing at all. The version suffix is dropped and the rest upper-cased.
Only from accessions, though. Stripping .<digits> from everything looks safe for
human symbols and is not safe in general: C. elegans sequence names — the standard
identifier for the majority of worm genes, which have no CGC name — are exactly that
shape, so C42D8.1, C42D8.2 and C42D8.3 would collapse to one gene and take
the background size, the set size, the overlap and the p-value down with them.
selection_propensity ¶
selection_propensity(background: Iterable[str], hits: Iterable[str], weights: Mapping[str, float], bins: int = 10) -> dict[str, float]
Estimated P(hit) for each background gene, given how much opportunity it had.
In splicing data "opportunity" is the number of units tested in that gene: a long,
many-exon gene offers far more chances to contain a significant junction than a
two-exon one, for reasons that have nothing to do with the biology. Genes are binned
by that measure and each bin's observed hit rate becomes the propensity of its
members — the same device goseq uses for gene length in RNA-seq.
over_representation ¶
over_representation(hits: Iterable[str], background: Iterable[str], gene_sets: Mapping[str, Sequence[str]], min_size: int = 2, max_size: int | None = None, weights: Mapping[str, float] | None = None) -> pd.DataFrame
Hypergeometric over-representation of gene_sets among hits.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
hits
|
iterable of str
|
The genes of interest (e.g. differentially spliced). |
required |
background
|
iterable of str
|
All genes that could have been a hit (the universe). |
required |
gene_sets
|
mapping
|
Term -> member genes. |
required |
min_size
|
int
|
Restrict to sets of this size within the background. |
2
|
max_size
|
int
|
Restrict to sets of this size within the background. |
2
|
weights
|
mapping
|
Per-gene opportunity (for splicing, the number of units tested in that
gene). Supplied, the null stops assuming every gene was equally likely to be a
hit: genes are binned by opportunity, each bin's observed hit rate becomes its
members' propensity, and the set's p-value comes from Wallenius' non-central
hypergeometric with odds = (mean of |
None
|
Returns:
| Type | Description |
|---|---|
DataFrame
|
One row per tested set with the 2×2 counts, fold enrichment, p-value and BH
q-value, sorted by q-value, and carries :data: |
Notes
Uses the survival function
P(X ≥ k) = hypergeom.sf(k-1, M, n, N) with M background size, n set
size in background, N number of hits in background, k the overlap, plus
bias_odds when weights were given.
enrich_differential ¶
enrich_differential(diff_table: DataFrame, gene_sets: Mapping[str, Sequence[str]], q: float = 0.05, min_delta: float = 0.1, gene_col: str = 'gene_id', name_col: str = 'gene_name', weight_by_units: bool = True, **kwargs) -> pd.DataFrame
Convenience: ORA of significant genes from a differential table.
Hits are the genes of significant units (q ≤ threshold and |ΔΨ| ≥
min_delta); the background is every gene that was tested.
Gene sets are keyed by symbols (MSigDB *.symbols.gmt, GO, KEGG) about as often
as by accessions, and a GTF supplies both, so whichever of gene_col and
name_col overlaps the sets more is the one used. Identifiers are matched through
:func:normalize_gene_id rather than verbatim.
weight_by_units corrects the bias that makes gene-level ORA on splicing data
misleading: a gene contributes as many chances of being a hit as it has tested
junctions, so large genes are over-represented among the hits for reasons that are
not biological. The number of rows per gene becomes its opportunity measure — see
weights in :func:over_representation. Set it to False for the plain
hypergeometric, which is what every ORA tool does and what this did before.
Junction features: splicescope.cryptic¶
Feature engineering for cryptic-event calling.
Novel junctions are a mix of genuine cryptic splicing and technical noise (mis-alignment, template switching, low-level transcriptional leakage). We turn each junction into a feature vector that a classifier can use to tell them apart:
====================== ==================================================
feature intuition
====================== ==================================================
intron_length cryptic introns have a plausible length distribution
log_max_count real events are reproducibly supported by reads
n_samples_support real events recur across replicates; noise does not
mean_psi_donor real events take a non-trivial share of their donor
canonical_motif genuine splicing tends to use GT/AG (or CT/AC)
dist_known_donor cryptic sites sit near, but not on, annotated sites
dist_known_acceptor same, for the 3' site
is_novel_both both sites unannotated (the strongest cryptic prior)
====================== ==================================================
extract_features ¶
extract_features(psi_annotated: DataFrame, known: DataFrame, novel_only: bool = True) -> pd.DataFrame
Build a per-junction feature table from an annotated+quantified long table.
psi_annotated is the output of :func:quantify.compute_psi run on an
annotated table (so it carries sclass and psi_donor). One row per
unique junction is returned, ready for :mod:splicescope.ml.
Classifier: splicescope.ml¶
Learn to separate true cryptic-splicing junctions from noise.
Novel junctions are abundant and mostly artefactual. Rather than hand-tuning
thresholds, we train a classifier on the engineered features from
:mod:splicescope.cryptic and evaluate it honestly with stratified
cross-validation. The model ships with a model card (metrics + permutation
importances + provenance) so a reviewer can see exactly how it was assessed.
CrypticClassifier
dataclass
¶
A cross-validated classifier for cryptic-junction calling.
evaluate ¶
evaluate(feats: DataFrame) -> dict
Stratified-CV evaluation. Returns ROC-AUC and average precision.
fit ¶
fit(feats: DataFrame) -> CrypticClassifier
Fit on all data and compute permutation importances.
score_table ¶
score_table(feats: DataFrame) -> pd.DataFrame
Return junctions with a cryptic_score column, most-likely first.
Scores tie readily — a forest that is certain gives several junctions exactly 1.0 — and sorting on the score alone left their order to whatever the input order happened to be. The top of this table is what a reader looks at, and a table whose rows move between runs for no reason is one they cannot cite. Ties break on the coordinates.
Plotting: splicescope.plotting¶
Publication-quality plots (matplotlib only, theme-consistent).
plot_annotation_summary ¶
plot_annotation_summary(summary: DataFrame, ax=None)
Bar chart of unique junctions per class.
plot_volcano ¶
plot_volcano(diff, q: float = 0.05, min_delta: float = 0.1, ax=None)
ΔΨ volcano: effect size vs −log10 q-value.
plot_importance ¶
plot_importance(importances: DataFrame, ax=None)
Horizontal bar of permutation importances.
plot_event_summary ¶
plot_event_summary(events: DataFrame, ax=None)
Bar chart of detected events per type (SE / MXE / A5SS / A3SS).
plot_event_volcano ¶
plot_event_volcano(ediff, q: float = 0.05, min_delta: float = 0.1, ax=None)
ΔΨ volcano for events, coloured by event type.
plot_consequence_summary ¶
plot_consequence_summary(consequence: DataFrame, ax=None)
Bar chart of predicted protein consequences, most disruptive class first.
Only classes actually present are drawn, so a cassette-exon run and a splice-site-shift run each get a chart scaled to what they found.
plot_event_psi ¶
plot_event_psi(psi, groups, ax=None, title: str | None = None)
One event's PSI in every sample, split by condition.
The number a differential test reports is a summary; this is the measurement behind it. A reader who cannot see the replicates has to take ΔΨ on trust, and a single outlying sample is exactly what a summary hides.
psi is the rows of :func:splicescope.events.event_psi for one event; groups
maps sample to condition.
plot_enrichment ¶
plot_enrichment(enrich, top: int = 10, ax=None)
Horizontal bar of the most enriched gene sets (−log10 q-value).
savefig ¶
savefig(fig, path: str | Path, dpi: int = 150, close: bool = True)
Write fig to path, creating parent directories, and release it.
pyplot keeps every figure alive until it is closed, so a pipeline that writes a
dozen panels per run leaks all of them; pass close=False to keep the figure.
Simulation: splicescope.simulate¶
Synthetic ground-truth data so the whole toolkit runs anywhere — no downloads.
The generator builds a tiny annotated "genome" and then emits per-sample splice junctions that mimic what STAR reports, including:
- canonical introns (the known annotation),
- true cryptic-exon events — biologically faithful: a cryptic exon inside an
intron produces a
novel_acceptorjunction (sharing the upstream known donor) and anovel_donorjunction (sharing the downstream known acceptor). These recur across replicates, use a canonical GT/AG motif, and are up-regulated in condition B, and - noise novel junctions — sporadic, low-support, often non-canonical.
Every junction carries an is_cryptic_truth label, which is what the ML module
is trained to recover. Fully reproducible via seed.
simulate_dataset ¶
simulate_dataset(n_genes: int = 8, exons_per_gene: int = 5, n_per_group: int = 4, cryptic_fraction: float = 0.5, alt_ss_fraction: float = 0.0, mxe_fraction: float = 0.0, label_noise: float = 0.0, seed: int = 0) -> SimulatedDataset
Generate a reproducible splicing dataset with known ground truth.
label_noise flips is_cryptic_truth for that fraction of junctions,
mimicking imperfect curation. Left at 0 the labels are exact (deterministic);
a small value (e.g. 0.1) makes the downstream ML task realistically hard.
gene_exons ¶
gene_exons(ds: SimulatedDataset) -> dict[str, list[tuple[int, int]]]
Exon blocks per gene, reconstructed from the known introns (120 bp flanks).
simulate_genome ¶
simulate_genome(ds: SimulatedDataset, utr5: int = 30, utr3: int = 60, seed: int = 0) -> tuple[str, dict[str, list[tuple[int, int]]]]
Build a synthetic chromosome consistent with the simulated annotation.
Exonic sequence is drawn from sense codons, so every gene carries a real open
reading frame that ends in a single stop; intron ends carry the canonical
GT/AG dinucleotides; everything else is random. That is enough for the
protein-consequence layer to be exercised end to end — a cryptic exon spliced
into one of these transcripts genuinely shifts the frame and genuinely has (or
has not) an in-frame stop, rather than being asserted to.
Uses its own random stream, so the junction data in ds is untouched.
Returns the chromosome sequence and the CDS blocks of each gene.
write_fasta ¶
write_fasta(sequence: str, path: str | Path, name: str = 'chr1', width: int = 60) -> Path
Write sequence as a FASTA plus the .fai index :class:GenomeFasta needs.
write_dataset ¶
write_dataset(ds: SimulatedDataset, outdir: str | Path, seed: int = 0) -> Path
Write a simulated dataset: GTF, per-sample SJ.out.tab, groups.tsv and a genome.
The GTF carries CDS records and the genome is written with its .fai, so the
protein-consequence layer runs on the built-in data with no downloads.
Everything at once: splicescope.demo¶
The whole pipeline on synthetic data, in one call.
This is what the Streamlit dashboard runs, and it lives here rather than in the app so that it is exercised by the test suite. An interactive demo is the first thing a reader tries, and its analysis had no test at all: a change to any signature it touches would have shown up as a traceback on a public page rather than as a red build.
from splicescope.demo import run_demo
result = run_demo(n_genes=20, cryptic_fraction=0.6)
result.differential.head()
DemoResult
dataclass
¶
Everything the dashboard displays, and why any of it is missing.
run_demo ¶
run_demo(n_genes: int = 20, n_per_group: int = 6, cryptic_fraction: float = 0.6, label_noise: float = 0.12, seed: int = 11, min_reads: int = 5, n_estimators: int = 300, with_consequences: bool = True) -> DemoResult
Simulate a dataset and run annotation, Ψ, the differential test and the classifier.
The classifier is skipped, rather than raising, when the simulated data leaves every
junction in one class: a small number of genes with a low cryptic fraction and no
label noise is a perfectly reasonable thing to ask for, and the junction classes and
differential results are still worth showing. :attr:DemoResult.classifier_note then
says why the rest is absent.
with_consequences also detects events and predicts what each cassette exon does
to the protein — the part that answers so what? about a cryptic junction. It needs
a genome, so the simulated dataset is written to a temporary directory and read back
through the same reader a real run uses.
n_estimators is exposed only so that a test sweeping many parameter combinations
does not have to train a full forest for each; the dashboard leaves it alone.