Skip to content

Technical algorithm

Jump to bottom
Graham edited this page Jul 31, 2026 · 53 revisions

Technical Details

This page provides a detailed technical description of intronIC's algorithm, data flow, and machine learning architecture (current as of v3.1.0).

Quick reference: this page is the authoritative prose description of the scoring stages, the species-adjudicator gate logic, and the parameter defaults. The exact values are defined in the code — the adjudicator parameters and gate logic in src/intronIC/scoring/species_adjudicator.py, also stamped in the bundled model's metadata — and the model-bundle layout is documented in docs/v3_bundle_schema.md.

Pipeline Overview

Input: Genome (FASTA) + Annotation (GFF3/GTF) + Species Name

Note: By default, intronIC uses streaming mode, which processes one chromosome at a time. Streaming and in-memory modes produce bit-identical classifications since v2.4. The choice is purely a memory tradeoff (roughly 50% less RSS on full human at -p 5 for streaming).

  1. Stage 1: Intron Extraction

    • Parse annotation hierarchy (gene → transcript → CDS/exon)
    • Infer intron coordinates from exon gaps
    • Extract intron + flanking exon sequences from genome
    • Filter duplicates, short introns, isoforms (configurable)

  2. Stage 2: PWM Scoring

    • Score 5' splice site, branch point, 3' splice site regions
    • Calculate log-odds ratios: $\log\left(\frac{P(\text{seq}|\text{U12})}{P(\text{seq}|\text{U2})}\right)$
    • Select best branch point position from search window
    • Apply species-specific U2-type background correction (see Species-Specific Background Correction)

  3. Stage 3: Classification

    • A single calibrated RBF SVM ensemble (pmotif_adjudicated, 42 sub-models) scores the 6 raw features (5'_raw, bp_raw, 3'_raw, bp_offset, bp_scan_confidence, support2_raw); there is no per-species z-normalization
    • Per-intron ensemble margin → P_motif = σ(2.796·margin − 1.178): a species-agnostic, Platt-calibrated motif probability

  4. Stage 4: Species adjudication

    • Per-species z_excess (Poisson significance of the strong-call count vs the genome's own U2-type tail) + a per-genome strength gate → motif_category ∈ {DETECTED, INCONCLUSIVE, NOT_DETECTED, UNASSESSABLE}
    • Anchors loss_ceiling_z = 2.60, bearer_floor_z = 5.50; strength gate p_gumbel_p95 ≤ 0.01 (co-fallback cs_p95 ≥ 5.0)

  5. Stage 5: Calling

    • type_id = u12 iff P_motif ≥ 0.5 and motif_category ≠ NOT_DETECTED; only NOT_DETECTED suppresses calls
    • Write adjusted_score = 100·q_eff·P_motif (the 0–100 calling scale) and rel_score = adjusted_score − 90

Output: .score_info.iic, .meta.iic, .bed.iic, .introns.iic, .metrics.iic.json, .tail_model.iic.json, plots

Relationship to the published method (v1)

The stages below describe intronIC v3. For readers familiar with the original method (Moyer et al. 2020), this section consolidates what changed and why; the mechanics are detailed in the linked sections.

v3 retains the PWM motif scoring of v1 (5′SS, branch-point, and 3′SS log-odds) but replaces the per-species z-normalization that v1 used to make scores comparable across genomes. That normalization assumed every genome shared the human reference geometry — a real U12-type mode at a fixed location and a bounded U2-type background — and failed in two ways when that assumption did not hold. The two are the same effect seen from opposite ends:

  • Divergent bearers were under-called. Where a species' U12-type population sits far from the human calibration (e.g. Amborella, Oryza), rescaling placed genuine U12-type introns below the decision boundary, so real introns were missed. v3's species-agnostic P_motif no longer under-calls these.
  • Loss genomes accrued false positives. In genomes lacking the minor spliceosome, per-species rescaling pushes the top of the U2-type distribution into U12-type score territory ("z-inflation"), producing false positives where no U12-type population exists. Leave-clade-out AUC is 0.916 (raw features) versus 0.786 (z-normalized), with essentially the entire difference on this loss-species false-positive class.
v1 (Moyer 2020) v3 (current)
Features 3 motif scores 6 raw features → 9 with interaction terms
Classifier one human-trained SVM 42-model RBF-SVM ensemble, 97 species / 14 clades
Background correction none per-species U2-type subtraction
Cross-species comparability per-species z-normalization species-agnostic P_motif
Per-species logic folded into normalization output-level adjudicator
Final call threshold on the score P_motif ≥ 0.5 and motif_category ≠ NOT_DETECTED

The three motif PWMs are largely shared between v1 and v3, with two motif-layer changes introduced during the v2 series: the branch-point PWM was rebuilt from CoLa-seq empirical branch points (Zeng et al. 2022), and per-species background correction was added. The substantive changes are at the feature, classifier, and adjudication layers described below.

Stage 1: Intron Extraction

Coordinate Inference

Introns are inferred from gaps between consecutive exons (or CDS features) within the same transcript:

Exon 1        Intron 1          Exon 2          Intron 2          Exon 3
[==]------------------------------[==]------------------------------[==]
1-100          101-1600         1601-1700        1701-3200        3201-3300

Priority: CDS features are preferred over exon features when available, as they enable phase calculation.

Mixed CDS/Exon Handling: For transcripts with both CDS and exon features, introns are first generated from CDS features, then exon-only introns (typically in UTR regions) are added if they don't overlap existing CDS-derived introns. All introns are sorted by genomic position and assigned sequential indices, ensuring proper ordering (e.g., 5' UTR introns are numbered before CDS introns in coding direction).

Touching Exons (Annotation Artifacts): Some annotations contain adjacent exon features with no gap between them (zero-length "introns"). These are silently skipped and not included in the intron count (family_size). The intron index sequence remains contiguous (1, 2, 3...) with no gaps for these annotation artifacts.

Filtering Criteria

Filter Default Description
Duplicates Exclude Same coordinates from multiple isoforms
Longest isoform Keep only Can include all with -i
Minimum length 30 bp Adjustable via --min-intron-len
Ambiguous bases Exclude 'N' in scoring regions
Non-canonical Include Exclude with --no-nc

v3 note: -i (--allow-multiple-isoforms) and -d (--include-duplicates) are now honored end-to-end in both streaming and in-memory modes. (A v2.7.x parity bug, where streaming hardcoded longest-isoform-only and dropped -d duplicate sequences before scoring, was fixed in v3, with streaming/in-memory parity holding across all four flag combinations.)

Stage 2: PWM Scoring

Position Weight Matrices

Position weight matrices (PWMs) capture the probability of observing each nucleotide at each position in a motif. intronIC includes PWMs for:

  • U12-type: AT-AC and GT-AG subtypes
  • U2-type: GT-AG, GC-AG, and AT-AC subtypes

Each subtype has PWMs for all three regions (5' splice site, branch point, 3' splice site). The U2-type AT-AC PWMs prevent inflated log-ratio scores when scoring AT-AC introns, where the U2-type fallback to GT-AG matrices would produce artificially low U2-type scores at the dinucleotide positions.

Scoring Regions

Default scoring windows:

Region Relative to Start End Length Description
5' SS Intron 5' end -3 +9 12 bp Includes last 3 bp of upstream exon
Branch point Intron 3' end -55 -5 (clamped at 3'SS boundary) up to 50 bp Search window, clamped to exclude 3'SS region
3' SS Intron 3' end -6 +4 10 bp Core acceptor only (excludes the polypyrimidine tract, PPT)

Non-overlapping regions: The branch point search window is clamped at the 3'SS scoring boundary to prevent feature overlap. For a 100 bp intron, the BP search region covers positions 45–93 and the 3'SS covers 94–100.

PWM Selection

For each intron, the terminal dinucleotides (e.g. GT-AG, AT-AC, GC-AG) determine which PWM subtype is used. Both the U12-type and U2-type PWMs are selected to match the intron's dinucleotide class:

  • A GT-AG intron is scored with U12-type GT-AG and U2-type GT-AG matrices
  • An AT-AC intron is scored with U12-type AT-AC and U2-type AT-AC matrices
  • A GC-AG intron is scored with U12-type GT-AG (fallback) and U2-type GC-AG matrices

When either model lacks a direct match for the dinucleotide class (e.g. GC-AG introns use U12-type GT-AG as a fallback, rare boundary types fall back to GT-AG for both models), the fallback PWM's terminal dinucleotide positions are masked so that the score reflects the surrounding motif context rather than a dinucleotide mismatch penalty.

Log-Odds Ratio Calculation

For each region, the raw score is a log-odds ratio computed by scoring the same sequence with both the U12-type and U2-type PWMs:

$$\text{LLR} = \log_2\left(\frac{\prod_{i} P(b_i | \text{U12 matrix})}{\prod_{i} P(b_i | \text{U2 matrix})}\right)$$

Where:

  • $b_i$ is the nucleotide at position $i$
  • The product is taken over all positions in the scoring window
  • Higher positive values favor U12-type
  • Higher negative values favor U2-type
  • Zero means equally likely under both models

For the 5'SS and 3'SS regions, the scoring window is fixed (see table above) and the log-ratio is computed directly.

Branch Point Selection

For the branch point region, position selection and scoring are separate steps:

  1. Position selection: The search window is scanned with a sliding window using the U12-type BPS PWM only. The position with the highest raw U12-type probability product is selected as the putative branch point.
  2. Log-ratio scoring: The same sequence at the selected position is then scored with both the U12-type and U2-type BPS PWMs, and the log-ratio is computed. This ensures the ratio compares both models at the same location.

The branch point adenosine position (bp_offset) is reported relative to the 3' splice site, derived from the PWM's defined reference position (biological position 0 in the matrix). U12-type branch points cluster tightly at -10 to -15 nt from the 3'SS (median -13), while U2-type branch points are more dispersed (-20 to -35 nt) (Pineda & Bradley 2018).

The scan also computes a BPS confidence score (bp_scan_confidence): log2(best_score / mean_score) across all scan positions. This measures how sharply the best match stands out from background: high values indicate a well-defined motif (U12-like), low values indicate a flat landscape (U2-like).

Species-Specific Background Correction

Before classification, intronIC optionally corrects the U2-type denominator PWMs using the species' own nucleotide composition, so the raw motif log-odds are anchored against the target genome's background.

In intronIC ≤ v2.3 this step was load-bearing: the default model was trained on human reference data, so non-human species with different composition produced systematically shifted score distributions. Background correction reanchored the U2-type log-odds against the target species' own distribution, recovering calibration.

The v3 default is the multispecies pmotif_adjudicated ensemble, trained on 41,333 introns drawn from 97 species across 14 broad clades and scoring the raw motif log-odds directly (no per-species z-normalization). Background correction now plays the role of an inference-time robustness layer for out-of-distribution species: for genomes deep inside the training distribution it is largely redundant with what the multispecies model already handles, but for extreme-composition genomes (e.g., very AT-rich or GC-rich) or clades not represented in training, it remains the cheapest way to keep the U2-type log-odds well-anchored. Because the training corpus was itself scored with background correction on, disabling it at inference would create a train/inference distribution mismatch and is not recommended for the bundled model.

Method

The correction blends empirical per-position nucleotide frequencies from the target species' intron pool into the U2-type PWMs:

$$\text{PWM}_{\text{corrected}} = w \cdot \text{PWM}_{\text{empirical}} + (1 - w) \cdot \text{PWM}_{\text{prior}}$$

The blending weight is computed per dinucleotide subtype (GT-AG, GC-AG, AT-AC):

$$w = \frac{n_{\text{subtype}}}{n_{\text{subtype}} + n_0}$$

Where $n_{\text{subtype}}$ is the number of introns with that terminal dinucleotide and $n_0$ is a shrinkage parameter (default 1000). Rare subtypes (e.g., AT-AC with a few hundred introns) stay close to the prior PWM; common subtypes (e.g., GT-AG with 200k+ introns) shift toward the species' empirical distribution. Only the U2-type denominator PWMs are modified; the U12-type numerator PWMs remain fixed, preserving the signal that defines U12-type introns.

Configuration

Background correction is enabled by default. The key parameter is n0 (Bayesian shrinkage strength); use the same value at classification as was used to score the training corpus (default 1000) for the bundled model.

scoring:
  species_background:
    enabled: true
    n0: 1000          # Shrinkage strength
    min_introns: 500  # Skip correction below this count

Stage 3: Classification — the raw-feature SVM ensemble

v3 has no normalization stage. Earlier versions (v2.3–v2.7) fit a per-species z-scaler before classification. That entire z-normalization layer was removed in v3 after the eval-corpus work showed per-species z-normalization hurts cross-species discrimination: it manufactures false positives in U12-absent species ("z-inflation"). The classifier now operates directly on the background-corrected raw motif log-odds.

Every intron is scored by a single calibrated RBF SVM ensemble (pmotif_adjudicated, model_id = pmotif_adjudicated_v23corpus_C200_g0.001). The per-intron ensemble margin (the SVC decision_function, also called the decision distance; all three names denote the same signed distance from the hyperplane, used throughout as ensemble margin) is mapped through a frozen Platt calibration to P_motif:

$$P_{\text{motif}} = \sigma(2.796 \cdot \text{margin} - 1.178)$$

The constants (2.796, −1.178) are a frozen sigmoid fit of the ensemble margin against the binary U12-type/U2-type training labels; they are preserved across adjudicator gate re-stamps and are not retrained. P_motif is a species-agnostic motif probability: a textbook GTATCC U12-type site scores ~0.99 regardless of how many U12-type siblings its genome has. It is calibrated on the per-intron continuum, with leave-clade-out (each broad clade held out in turn, then scored by a model trained on the rest) expected calibration error (ECE, the mean gap between predicted probability and observed frequency) ≈ 7×10⁻⁴, so a user can threshold it post-hoc at any confidence level (0.75 / 0.90 / 0.99). This is the cross-species-comparable score that per-species z-normalization was originally meant to provide, without the z-inflation that broke on divergent species like Amborella and Oryza.

Feature space

The model scores 6 declared raw features per intron:

Motif log-odds (3):

  • 5'_raw — 5′SS raw motif log-odds (U12/U2), background-corrected
  • bp_raw — branch-point raw motif log-odds, background-corrected
  • 3'_raw — 3′SS raw motif log-odds, background-corrected

Branch-point geometry + support (3):

  • bp_offset — branch-point adenosine distance from the 3′SS (negative integer; U12-type branch points cluster at −10…−15, U2-type more dispersed at −20…−35)
  • bp_scan_confidence — log₂(best/mean) of U12-type BPS PWM scores across the scan window (motif sharpness; U12-type introns peak more sharply)
  • support2_raw — the second-largest of max(0, 5'_raw), max(0, bp_raw), max(0, 3'_raw); encodes "at least two regions positively support a U12-type call". Zero when only one site is positive.

6 declared vs 9 fitted (caveat). Those are the 6 inputs the predictor feeds and what score_info.iic stores. Inside each sub-model's sklearn pipeline (BothEndsStrongTransformer → StandardScaler → SVC) the transformer deterministically appends 3 engineered interaction terms, so the fitted SVC operates on 9 features (n_features_in_ = 9):

  • corr_5_bp = min(⌊5'⌋₊, ⌊BP⌋₊)
  • corr_bp_3 = min(⌊BP⌋₊, ⌊3'⌋₊)
  • gap_5_bp = max(5' − BP, 0)

These carry no new measurement; they only re-express the 5′/BP/3′ interactions to shape the RBF kernel (there is no per-model feature bagging in this bundle; dropped_feature = None for all 42 sub-models).

RBF SVM

Each sub-model is an RBF (radial basis function) kernel SVC:

  • Kernel: RBF (radial basis function)
  • Hyperparameters: C = 200, γ = 0.001 (jointly grid-searched on the multispecies corpus)
  • Class weights: balanced, to handle the ~0.5% U12-type prevalence

The RBF kernel lets the classifier learn curved boundaries: for example, requiring both a strong donor and a strong branch point for a confident U12-type call, rather than letting one compensate for the other as a linear model would.

Probability calibration

The raw SVC outputs the ensemble margins defined above; two calibrations are derived from them:

  • svm_score (column 3 of score_info.iic) is the ensemble's mean isotonic-calibrated U12-type probability (×100), retained for auditability. Isotonic was selected over sigmoid by log-loss during training.
  • P_motif (the call-driving score) is the Platt calibration of the ensemble margin, not the isotonic probability. The isotonic output saturates at 0/1 and discards the margin-strength information the adjudicator needs.

Ensemble

The bundle is a single RBF SVM ensemble of 42 sub-models / 210 sub-estimators. Each sub-model sees all U12-type references but a different random subsample of the U2-type references (bagging); the final per-intron probability is the mean across sub-models, and ensemble_sigma (their standard deviation) quantifies model agreement.

Training the default model

The v3 default bundle (default_pretrained.model.pkl) is the pmotif_adjudicated raw-feature ensemble, trained on the v3 multispecies corpus.

Training corpus

  • 41,333 introns: 10,003 U12-type positives + 31,330 U2-type and hard-negative records (a hard negative is a U2-type intron whose motif closely mimics a U12-type site, so it sits near the decision boundary and sharpens it), drawn from 97 species across 14 broad evolutionary clades (vertebrates, arthropods, nematodes, monocots, eudicots, fungi, basal animals, protists, and others), with a further 5 species held out for evaluation (790 positives).
  • Per-intron labels are assigned by comparative genomic analysis (orthology across species): each intron's classification is supported by orthology evidence from multiple species, rather than a single reference annotation.
  • The corpus ships as sequences at u{12,2}_reference_multispecies.introns.iic.gz (full intron sequences with 50 bp flanks per side); the legacy human-anchored v2.3 references remain bundled as u{12,2}_reference_human.introns.iic.gz. See Training data and PWMs.
  • PWMs are unchanged from v2.3: the 5′/3′SS PWMs from a ~580-intron human + comparative-genomics-conserved gold standard, and the 12 bp BPS PWMs from CoLa-seq empirical branch-point positions (Zeng et al. 2022). The v3 effort happened entirely at the SVM/adjudicator layer.

Frozen calibration constants

The bundle carries version-pinned calibration constants (a gate change is a re-stamp, not a retrain: the Platt/anchors are preserved):

  • Platt (margin → P_motif): (a, c) = (2.796, −1.178) (the frozen sigmoid fit described in Stage 3).
  • Adjudicator (ADJUDICATOR_PARAMS_VERSION = zexcess_gap_pgumbel_cs_2026-07-03): anchors loss_ceiling_z = 2.60, bearer_floor_z = 5.50; strength gate p_gumbel_threshold = 0.01, cs_point_threshold = 5.0, cs_min_calls = 3 (the minimum number of calls needed to evaluate an upper-tail percentile); min_u2 = 200 (the minimum U2-type introns needed to fit a per-genome EVT tail; below it a genome falls to LOW_N → UNASSESSABLE).
  • EVT tail fit: min_evt_excesses = 20 — the minimum number of U2-type tail exceedances (introns above the 90th-percentile EVT threshold) needed to fit a genome's per-genome Gumbel tail. Below it the tail is treated as unfittable and the genome falls to UNASSESSABLE, distinct from the LOW_N (< min_u2) route to the same category.

The ensemble uses isotonic regression for the svm_score probability (selected over sigmoid by log-loss). Training a fresh ensemble with intronIC train --scoring-mode pmotif_adjudicated emits a valid, runnable bundle, but the frozen Platt/adjudicator constants are ensemble-specific and would need re-fitting for a production-grade custom bundle.

Evaluation

  • Raw motif features beat z-normalization on cross-species discrimination: leave-clade-out real-classifier area under the ROC curve (AUC) 0.916 (raw) vs 0.786 (z), with the entire gain on the z-inflated loss-species false-positive class.
  • Loss-FP control: across 19 snRNA-confirmed U12-absent genomes (every U12-absent clade: green algae, ciliates, Apicomplexa, Dikarya fungi, nematodes, kinetoplastids, microsporidians, diplomonads), P_motif ≥ 0.9 gives a median 3, max 15 motif-strong calls per genome, the irreducible motif floor, with no inflation.
  • Adjudicator: leave-clade-out validation shows 0 cross-errors: held-out genomes never mis-cross the gap; the worst case is an honest INCONCLUSIVE.

Prior architectures remain accessible: legacy z-normalization / mode-separation bundles are reproducible from the pre-zstack-removal git tag (they are rejected at load by the v3 code, which is raw-feature only).

Species adjudicator

All per-species adaptation lives in a single output-level layer that sits on top of the per-intron P_motif scores and makes the per-species U12-type motif population call. It emits two independent quantities, reported separately (never fused into one thresholded score; that fusion was the core mistake of the superseded design):

  1. P_motif — the per-intron, species-agnostic motif probability (Stage 3).
  2. motif_category — a per-species gate over the genome's P_motif calls: one of DETECTED / INCONCLUSIVE / NOT_DETECTED / UNASSESSABLE, reported with the continuous z_excess alongside.

z_excess — the population statistic

The gate's primary driver is z_excess: the Poisson significance, how many standard deviations above expectation under a Poisson count model, of the number of strong calls (P_motif ≥ 0.9) against what the genome's own U2-type background tail predicts. A Gumbel extreme-value null (EVT, extreme-value theory, studies the distribution of the maximum of many independent draws, which converges to a Gumbel distribution; Coles 2001) is fit to the genome's U2-type margin tail and extrapolated into the call region; z_excess asks whether there are significantly more high-margin introns than the U2-type tail would throw up by chance. A loss (a lineage that has lost the minor spliceosome and so carries no U12-type introns) has motif-strong false positives that are the U2-type tail (z_excess ≈ 0); a bearer (a lineage that still carries U12-type introns) has a separate population the U2-type tail cannot explain (large z_excess). Referencing the count to each genome's own U2-type tail is robust at both genome-size extremes (unlike raw call count, which inflates on huge genomes).

The gate

classify_motif_category evaluates, in order:

if z_excess is not finite:                          -> UNASSESSABLE   # no per-genome U2-type null (LOW_N / unfit tail)
if z_excess >= bearer_floor_z (5.50):               -> DETECTED       # count gate: a clear U12-type population
if n_calls >= 3 and (p_gumbel_p95 <= 0.01           # PRIMARY strength gate
                     or cs_p95 >= 5.0):              -> DETECTED       # CO-FALLBACK (null-free)
if 0 < n_calls < 3 and min(bg_fdr) <= 3e-3:         -> INCONCLUSIVE   # LOW-K escape hatch (v3.1); see below
if z_excess <= loss_ceiling_z (2.60):               -> NOT_DETECTED   # calls consistent with the U2-type background
else:                                               -> INCONCLUSIVE   # the gap [2.60, 5.50): abstain, corroborate

The two anchors are empirical, not tuned knobs. loss_ceiling_z = 2.60 is the highest z_excess reached by any snRNA-confirmed U12-loss genome in the calibration panel (Aspergillus coremiiformis, which retains only u4atac/u6atac). bearer_floor_z = 5.50 is a trust threshold rather than a separator: it was widened from an earlier 4.00 because the near-floor zone is genuinely mixed (the loss Monocercomonoides at z_excess ≈ 4.24 sits adjacent to the bearer Blyttiomyces at z_excess ≈ 4.54), so the gap between the two anchors is treated as abstention territory rather than forced into a call.

The low-k escape hatch (v3.1)

The strength gate is hard-gated on n_calls >= 3, deliberately: below that the cs_p95 is the maximum, i.e. a single outlier, which is what the floor exists to refuse. But z_excess is also bounded below the loss ceiling at low call counts by construction — with one call, that call is the call-core, so no call exceeds it and z_excess <= 0; with two, z_excess <= 1. Both sit far under loss_ceiling_z = 2.60.

The consequence was that a genome with one or two strong calls could not reach DETECTED by any path. It fell through to NOT_DETECTED, which zeroes adjusted_score, rel_score and type_id on every row. Such genomes were not adjudicated as losses — they were never adjudicated, and the evidence was discarded, including for genomes that carry a complete minor spliceosome.

bg_fdr is the only statistic in the adjudicator that does not degenerate at one call: cs_p95 needs a distribution of calls, whereas a per-intron q-value against the genome's own fitted U2-type tail needs only one. So when a genome has 0 < n_calls < 3 and its strongest call clears min(bg_fdr) <= 3e-3, the adjudicator now returns INCONCLUSIVE rather than NOT_DETECTED. Its scores and calls survive for downstream corroboration; no population is asserted.

INCONCLUSIVE rather than DETECTED is deliberate. One or two introns is not a population, and the asymmetric loss behind the n_calls >= 3 floor still applies. The hatch's purpose is to stop destroying evidence, not to adjudicate.

The 3e-3 threshold is a preserve-the-evidence trust line, not a bearer/loss separator. It sits above the presumed-loss floor on purpose: because the gate abstains rather than asserts, erring wide is the cheap direction, and the axis that actually separates these genomes — the minor spliceosome machinery — is not visible to intronIC, which is snRNA-blind by design. It also falls in the one empty band in the low-call-count distribution, so any value between roughly 2.7e-3 and 6.0e-3 selects the same genomes.

On a 2,785-genome corpus the hatch moves 8 genomes, all NOT_DETECTED → INCONCLUSIVE, and nothing with three or more calls: four from U12-bearing lineages (three Mucor, Pythium, Seison) and four from loss-prior ones (Leishmania, Babesia, Porospora). That mix is the expected shape for a gate that abstains. Set lowk_gate_enabled = False to restore v3.0 behaviour.

Ordering matters: the strength gate is checked before the loss ceiling, so a genome whose z_excess sits below 2.60 (count-wise loss territory) but which has a few genuinely strong calls is rescued to DETECTED — the divergent-bearer recovery path (e.g. the oomycete Achlya hypogyna: z_excess ≈ 2.13 but cs_p95 ≈ 7.3 → DETECTED, snRNA-corroborated).

Anchor calibration

The two z_excess anchors are frozen from a supervised panel assembled at calibration time: a cross-clade set of minor-spliceosome-confirmed genomes, augmented with 68 freshly-run divergent bearers from the lineages that pin the bearer floor (early-diverging fungi, oomycetes, and basal metazoans). Each panel genome's bearer/loss status is fixed by a minor-snRNA search, never by the motif classifier:

  • A minor snRNA counts as present when its best cmsearch/Infernal hit clears an inclusion threshold of E ≤ 0.01.
  • A genome is scored a bearer when at least three of the four minor snRNAs (U11, U12, U4atac, U6atac) are present. This "3-of-4" rule is deliberately defining-aware: a genome that has lost the minor spliceosome retains at most two of the four, so reaching three necessarily requires one of the defining snRNAs — without forcing the U11+U12 pair specifically, which would wrongly reject lineages (for example some gut fungi) whose U11 is too divergent to detect. loss_ceiling_z = 2.60 is then the highest z_excess reached by any confirmed loss in this panel (Aspergillus coremiiformis, which retains only U4atac/U6atac — 2 of 4).

Because minor-snRNA detectability is itself clade-conditional — U12 and U6atac are robustly detectable across clades, while U11 and U4atac are less reliable in some lineages — snRNA absence is only informative for the snRNAs that are dependable in a given clade. This is a calibration-time and downstream-corroboration consideration; the runtime classifier itself remains motif-only.

The strength gate

z_excess is a count statistic, so it is blind to a divergent bearer with only a few genuinely strong U12-type calls. The strength gate fixes this by evaluating the same per-genome Gumbel U2-type null at the call upper tail: p_gumbel_p95 = P(the genome's own U2-type background produces a max ≥ the call cs_p95), where cs_p95 is the 95th percentile of the un-clipped ensemble margins on the calls. p_gumbel_p95 ≤ 0.01 means "these strong calls are a ~1% outlier against this genome's own U2-type noise" → DETECTED. It is composition-adaptive (judged against each genome's own tail), which is why it replaced the earlier tuned constant cs_p95 ≥ 5.0 as the primary driver (cs_p95 = 5.0 ↔ p_gumbel_p95 ≈ 0.005); cs_p95 ≥ 5.0 is retained as a null-free co-fallback for marginal EVT fits. Recovery is no-regression by construction (the worst loss's p_gumbel_p95 = 0.044, well above 0.01).

What the categories mean (and do not)

motif_category names the motif evidence for a U12-type population, not biological bearer/loss truth:

  • DETECTED — a U12-type motif population is present, by motif. Corroborate downstream — not a proof of functional U12-type introns.
  • NOT_DETECTED — the strong calls are consistent with this genome's U2-type background. This is the only category that suppresses calls (type_id → u2, u12_count → 0). It is not itself a loss call; motif-silent divergent bearers can land here, where motif evidence alone is insufficient to resolve them (minor-snRNA presence or phylogenetic context can).
  • INCONCLUSIVE — in the gap [2.60, 5.50) and not strength-rescued: abstain within support, corroborate.
  • UNASSESSABLE — cannot reference a per-genome U2-type null (see low-N below).

INCONCLUSIVE and UNASSESSABLE do not suppress calls — their strong-motif introns still get type_id = u12; the category is a species-level confidence qualifier, not a per-intron veto.

-q / low-N fallback

The per-genome test presumes the genome has a U2-type background to reference:

  • Full-complement genome / -q (≥ min_u2 = 200 U2-type introns) → the EVT tail fits → normal adjudication.
  • Curated / tiny -q (< 200 U2-type introns) → LOW_N → motif_category = UNASSESSABLE, no bearer/loss call: per-intron P_motif is reported as-is (an absolute strength threshold on a hand-picked set is meaningless). Opt in to self-adjudication against a small genome's own noisier U2-type tail with --adjudicator-min-u2 <N>.
  • A well-powered genome with zero strong calls short-circuits to NOT_DETECTED.

Runtime scope + diagnostics

The adjudicator's runtime scope is motif-only — no phylogeny, no snRNA. Taxonomy and snRNA evidence are used only at calibration time (to source the loss_ceiling/bearer_floor anchors from snRNA-confirmed losses and freshly-run divergent bearers), never by the runtime classifier. Conservatism comes from abstention within support (the gap and the UNASSESSABLE band): a future bearer more extreme than anything sampled lands in the gap and abstains rather than being mislabeled.

Every pmotif_adjudicated run emits a per-species tail-model diagnostic — a sidecar .tail_model.iic.json and a figure .plot.tail_model.iic.png — that makes the decision visible: the genome's U2-type background margins (gray bars), the fitted U2-type tail extrapolated into the call region (gray line), and the U12-type calls that escape it (green; a loss genome's calls sit inside the tail's reach). See Output files.

The contrast is clear when a bearer and a loss genome are placed side by side. In D. melanogaster (DETECTED) the green calls stand well above the extrapolated U2-type tail (z_excess ≈ 9); in C. elegans (NOT_DETECTED) the handful of green calls fall on the fitted tail (z_excess ≈ 0.6) — i.e. they are no more numerous than the genome's own U2-type background predicts at that motif strength, so no U12-type population is called.

D. melanogaster — DETECTED C. elegans — NOT_DETECTED

Calling and score columns

The per-intron call combines the two adjudicator outputs:

  • type_id = u12 iff P_motif ≥ 0.5 (equivalently adjusted_score ≥ 50) AND motif_category ≠ NOT_DETECTED; otherwise u2.
  • adjusted_score = 100·q_eff·P_motif, where q_eff ∈ {0,1} is a binary effective gate that is 0 only when motif_category == NOT_DETECTED (else 1). So in a non-NOT_DETECTED genome adjusted_score = 100·P_motif; a NOT_DETECTED genome has all calls zeroed. adjusted_score is the adjudicated 0–100 calling scale.
  • rel_score = adjusted_score − 90 (range ≈ [−90, +10]): signed and centered on the 90% high-confidence threshold. rel_score > 0 ⇔ adjusted_score > 90 ⇔ P_motif > 0.9 = high-confidence U12-type. rel_score < 0 means "below the 90% threshold", not "anti-U12-type".

The q, P_adj, P_adj_lo and P_adj_hi columns were removed from score_info.iic in v3.1.0 (38 → 35 columns) — they were deterministic functions of P_motif + motif_category, and the _lo/_hi pair was bit-identical to P_adj. Read P_motif + motif_category. type_id moved from column 38 to 35 as a result, so parse score_info.iic by header name rather than by position.

Species-Specific Considerations

GC Content Effects

Species with GC content far from human can show shifted raw-score distributions. Species-specific background correction (on by default) reanchors the U2-type PWMs against the genome's own composition and keeps the raw motif log-odds comparable; --species-prior optionally adjusts the base-rate expectation via Bayes' rule as well.

U12-Absent Lineages

Species lacking U12-type introns are handled at the species-adjudicator level, not by scoring. A genome whose strong P_motif calls are consistent with its own U2-type background gets motif_category = NOT_DETECTED (z_excess ≤ 2.60, and no strength-gate rescue), which is the only category that suppresses calls (type_id → u2). Across the snRNA-confirmed U12-absent panel, P_motif ≥ 0.9 leaves only the irreducible motif floor with no false-positive inflation (see Loss-FP control for the 19-genome breakdown). No special configuration is needed.

Cross-Species Performance

The raw-feature pmotif_adjudicated model generalizes broadly across eukaryotes: P_motif is species-agnostic, so lineages with weaker normalized motif separation (which the old z-normalized model missed, e.g. Amborella, Oryza) are no longer under-called. Where genome-level evidence is genuinely thin, the adjudicator abstains within support (the INCONCLUSIVE gap or UNASSESSABLE) rather than mislabeling. Performance may still be uncertain for:

  • Lineages with unusual U12-type motifs far outside the training distribution
  • Divergent low-count U12-type bearers whose z_excess sits in the [2.60, 5.50) gap and are not strength-rescued (→ INCONCLUSIVE; corroborate with snRNA/phylogeny)
  • Curated / very small -q inputs with too few U2-type introns to reference a per-genome tail (→ UNASSESSABLE)

For species suspected to fall outside the training distribution, inspect the per-species motif_category / z_excess and the .tail_model.iic.png diagnostic, and consider providing custom reference sequences via intronIC train if necessary.

Memory and Performance

--streaming (default) and --in-memory produce bit-identical classifications (covered by tests/integration/test_streaming_equivalence.py); the choice is purely a runtime/memory tradeoff.

Streaming Mode (--streaming, default)

Per-contig pipeline:

  • Writes intron sequences to a temporary SQLite database during extraction
  • Keeps only scoring motifs in memory
  • Each phase (extraction, BG correction, PWM scoring, classification) parallelizes across contigs via multiprocessing.Pool

In-memory Mode (--in-memory)

Loads all intron sequences into memory at extraction time. Used internally by --sequences and --bed input modes (those bypass the per-contig streaming path).

Reference benchmark

Approximate timing on commodity hardware (single workstation, NVMe SSD), -p 5 --streaming. Peak RSS is data-bound (dominated by intron/sequence data, not the model) and roughly version-independent; the wall times are v2.7-era figures — v3's single raw-feature ensemble is typically faster than earlier (v2.6–v2.7) bundles:

Species Genome / annotation Scored introns Wall time Peak RSS
Drosophila melanogaster Release 6 + ISO1_MT (~140 Mb) 47,000 ≲8 min ~0.8 GB
Homo sapiens GRCh38.p13 + NCBI RefSeq GFF (~3.2 Gb) 257,123 ≲40 min ~5.3 GB

On multi-contig genomes --in-memory finishes at essentially the same wall time with roughly 2× the peak memory.

Parallelization

The -p N flag parallelizes PWM scoring and per-contig extraction:

  • Scoring is CPU-bound and parallelizes efficiently
  • Linear speedup up to ~8-16 cores
  • Diminishing returns beyond that
  • BG correction, scoring, and classification all dispatch through Pool workers in streaming mode

References

Original intronIC paper:

Moyer DC, Larue GE, Hershberger CE, Roy SW, Padgett RA. (2020) Comprehensive database and evolutionary dynamics of U12-type introns. Nucleic Acids Research 48(13):7066–7078. doi:10.1093/nar/gkaa464

U12-type intron databases:

Larue GE, Roy SW. (2023) Where the minor things are: a pan-eukaryotic survey suggests neutral processes may explain much of minor intron evolution. Nucleic Acids Research 51(20):10884-10908. doi:10.1093/nar/gkad797

Moyer DC, Larue GE, Hershberger CE, Roy SW, Padgett RA. (2020) Comprehensive database and evolutionary dynamics of U12-type introns. Nucleic Acids Research 48(13):7066-7078. doi:10.1093/nar/gkaa464

Alioto TS. (2007) U12DB: a database of orthologous U12-type spliceosomal introns. Nucleic Acids Research 35:D110-D115. doi:10.1093/nar/gkl796

Branch point mapping and spliceosome profiling:

Burke JE, Longhurst AD, Merkurjev D, Sales-Lee J, Rao B, Moresco JJ, Yates JR III, Li JJ, Madhani HD. (2018) Spliceosome profiling visualizes operations of a dynamic RNP at nucleotide resolution. Cell 173(4):1014-1030.e17. doi:10.1016/j.cell.2018.03.020

Zeng Y, Fair BJ, Zeng H, Krishnamohan A, Hou Y, Hall JM, Ruthenburg AJ, Li YI, Staley JP. (2022) Profiling lariat intermediates reveals genetic determinants of early and late co-transcriptional splicing. Molecular Cell 82(24):4681-4699. doi:10.1016/j.molcel.2022.11.004

Mercer TR, et al. (2015) Genome-wide discovery of human splicing branchpoints. Genome Research 25:290-303. doi:10.1101/gr.182899.114

Pineda JMB, Bradley RK. (2018) Most human introns are recognized via multiple and tissue-specific branchpoints. Genes & Development 32(7-8):577-591. doi:10.1101/gad.312058.118

U12-type intron retention and functional studies:

Niemelä EH, et al. (2014) Global analysis of the nuclear processing of transcripts with unspliced U12-type introns by the exosome. Nucleic Acids Research 42(11):7358-7369. doi:10.1093/nar/gku391

Madan V, et al. (2015) Aberrant splicing of U12-type introns is the hallmark of ZRSR2 mutant myelodysplastic syndrome. Nature Communications 6:6042. doi:10.1038/ncomms7042

Cologne A, et al. (2019) New insights into minor splicing—a transcriptomic analysis of cells derived from TALS patients. RNA 25(9):1130-1149. doi:10.1261/rna.071423.119

SVM probability calibration:

Platt JC. (1999) Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. Advances in Large Margin Classifiers pp. 61-74.

Extreme-value statistics (the Gumbel null behind z_excess / p_gumbel):

Coles S. (2001) An Introduction to Statistical Modeling of Extreme Values. Springer Series in Statistics, Springer. doi:10.1007/978-1-4471-3675-0

Clone this wiki locally