-
Notifications
You must be signed in to change notification settings - Fork 1
Technical algorithm
This page provides a detailed technical description of intronIC's algorithm, data flow, and machine learning architecture.
Input: Genome (FASTA) + Annotation (GFF3/GTF) + Species Name
Note: By default, intronIC uses streaming mode, which processes one chromosome at a time for ~85% memory reduction on large genomes.
-
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)
-
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
-
Stage 3: Normalization
- Convert raw log-odds to z-scores via robust scaling
- Fit on reference data or adapt to experimental data
-
Stage 4: Classification
- Feature augmentation (composite and sequence-derived features)
- RBF SVM ensemble with balanced class weights
- Probability calibration via isotonic regression
- Output:
$P(\text{U12-type})$ from 0–100%
Output: .meta.iic, .bed.iic, .introns.iic, .score_info.iic, plots
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.
| 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
|
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 AT-AC PWMs prevent inflated log-ratio scores when scoring AT-AC introns, where the U2 fallback to GT-AG matrices would produce artificially low U2 scores at the dinucleotide positions.
Default scoring windows relative to splice sites:
| Region | Start | End | Length | Description |
|---|---|---|---|---|
| 5' SS | -3 | +9 | 12 bp | Includes last 3 bp of upstream exon |
| Branch point | -55 | -5 (clamped at 3'SS boundary) | up to 50 bp | Search window, clamped to exclude 3'SS region |
| 3' SS | -6 | +4 | 10 bp | Core acceptor only (excludes 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.
For each intron, the terminal dinucleotides (e.g. GT-AG, AT-AC, GC-AG) determine which PWM subtype is used. Both the U12 and U2 PWMs are selected to match the intron's dinucleotide class:
- A GT-AG intron is scored with U12 GT-AG and U2 GT-AG matrices
- An AT-AC intron is scored with U12 AT-AC and U2 AT-AC matrices
- A GC-AG intron is scored with U12 GT-AG (fallback) and U2 GC-AG matrices
When either model lacks a direct match for the dinucleotide class (e.g. GC-AG introns use U12 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.
For each region, the raw score is a log-odds ratio computed by scoring the same sequence with both the U12 and U2 PWMs:
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.
For the branch point region, position selection and scoring are separate steps:
- Position selection: The search window is scanned with a sliding window using the U12 BPS PWM only. The position with the highest raw U12 probability product is selected as the putative branch point.
- Log-ratio scoring: The same sequence at the selected position is then scored with both the U12 and U2 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).
Raw log-odds scores have different ranges and distributions for each region:
- 5'SS scores might range from -50 to +10
- BP scores might range from -20 to +5
- 3'SS scores might range from -5 to +3
Normalization converts these to comparable z-scores for the SVM.
intronIC uses robust z-score normalization via sklearn's RobustScaler:
Where IQR is the interquartile range (
Key properties:
- Robust to outliers: Median and IQR are resistant to extreme values
- Centered: Distribution centered around 0 for each feature
- Comparable scales: All three regions have similar variance after scaling
- Minimal U12 contamination: Rare U12s (~0.5%) don't significantly affect robust statistics
Note: A second normalization step (sklearn StandardScaler) is applied inside the classification pipeline to the full augmented feature vector (base z-scores + composite features + extra features). This ensures all features have zero mean and unit variance before SVM classification.
| Mode | Description | Use Case |
|---|---|---|
human |
Use scales from training (human) data | U12-absent species, cross-species |
adaptive |
Refit scales on experimental data | Species with different GC content |
auto |
Use human if available in model | Default behavior |
Note: The pretrained model includes its scaler, so results are reproducible without additional steps when using --normalizer-mode human or auto.
The default model uses 8 features (8D):
Base features (3D):
-
s5z: 5' splice site z-score -
BPz: Branch point z-score -
s3z: 3' splice site z-score (core acceptor, -6 to +3)
Composite feature (1D):
-
support2: second-largest of max(0, s5z), max(0, BPz), max(0, s3z) — encodes "at least two of the three scoring regions positively support U12." Zero when only one site is positive, directly targeting the one-end-strong false positive pattern.
Extra features (4D):
-
bp_offset: branch point adenosine distance from 3'SS (negative integer). U12-type introns cluster at -10 to -15; U2-type introns at -20 to -35. -
ppt_t_weighted: T-weighted pyrimidine score (T=1.0, C=0.5, purine=0) over a fixed 20 nt window immediately upstream of the 3'SS scoring boundary. Measures PPT presence; higher values indicate U2-like pyrimidine tracts. -
ppt_longest_run: longest uninterrupted C/T run in the same 20 nt window. -
bp_scan_confidence: log2(best/mean) of U12 PWM scores across all positions in the BP search window. Measures BPS motif sharpness; U12-type introns have sharper peaks than U2-type introns in the human reference training data.
Linear coefficient reference (from a parallel 8D linear model, for interpretability):
| Feature | Coefficient | Importance | Interpretation |
|---|---|---|---|
| s5z | +2.97 | 57.0% | Strong donor → U12 |
| bp_offset | +0.64 | 12.3% | Closer BP → U12 |
| BPz | +0.45 | 8.5% | Strong BP motif → U12 |
| support2 | +0.39 | 7.4% | Distributed support → U12 |
| s3z | -0.35 | 6.8% | See note below |
| bp_scan_confidence | +0.22 | 4.2% | Sharper BPS peak → U12 |
| ppt_t_weighted | -0.11 | 2.1% | More pyrimidine → U2 |
| ppt_longest_run | -0.08 | 1.6% | Longer PPT runs → U2 |
Note: the s3z negative coefficient reflects that the core 3'SS AG acceptor does not strongly discriminate U12 from U2 in a linear model. The RBF kernel handles this more appropriately through its nonlinear boundary.
The default model uses an RBF (radial basis function) kernel SVM (sklearn's SVC):
- Kernel: RBF (nonlinear boundary, can learn local patterns)
- Hyperparameters: C=65.88, gamma=0.01 (optimized via grid search)
- Class weights: Balanced to handle ~0.5% U12 prevalence
- Support vectors: ~28-31 per ensemble model
The RBF kernel allows the classifier to learn curved boundaries — for example, requiring both strong donor AND strong BP for confident U12 calls, rather than allowing one to compensate for the other as a linear model would.
The raw SVM outputs decision distances (signed distance from hyperplane). These are converted to probabilities using sklearn's CalibratedClassifierCV. During training, both sigmoid (Platt scaling) and isotonic calibration are evaluated via cross-validation, and the method with lower log-loss is selected.
- Sigmoid (Platt scaling): Fits a logistic function with 2 parameters — more conservative probability estimates, fewer borderline calls
- Isotonic: Non-parametric monotonically increasing step function — more flexible, typically selected when sufficient training data is available
The default pretrained model uses isotonic calibration.
When --n-models > 1, multiple SVMs are trained with different U2 subsamples:
- Each model sees all U12 references but a different 80% of U2 references
- Predictions are averaged across models
- Reduces variance and improves robustness
The default pretrained model (v2.2.0) was trained on:
Training data:
- 472 human U12-type introns from IPA-conserved gold standard with CoLa-seq empirical branch points (Zeng et al. 2022)
- 30,155 human U2-type introns (20,690 original + 9,465 IPA-conserved/hard-negative, deduplicated by sequence hash)
Training process:
- Score all reference introns with PWMs (12 bp BPS motifs derived from CoLa-seq data)
- Normalize using robust scaling (median/IQR)
- Optimize SVM hyperparameters (C, gamma, penalty, class weight multiplier) via iterative grid search with 7-fold cross-validation
- Train 15-model ensemble on all reference data with 80% U2 subsampling per model
- Calibrate probabilities via isotonic regression (selected over sigmoid by log-loss)
Evaluation:
- 7-fold nested cross-validation: mean F1 = 0.9947 ± 0.0065, PR-AUC = 0.9999 ± 0.0002
- Validated across 13 species spanning vertebrates, arthropods, nematodes, plants, fungi, and protists
- False positive assessment in species thought to lack U12-type introns: 0 confident calls in C. elegans, 1 in Ascaris suum
Species with different GC content than human may show shifted score distributions. Options:
-
Adaptive normalization: Refit scaler on experimental data (
--normalizer-mode adaptive) -
Prior adjustment: Adjust base rate expectation (
--species-prior)
For species thought to lack U12-type introns (C. elegans, many fungi):
intronIC -g genome.fa -a annotation.gff -n species \
--normalizer-mode human --species-prior 1e-6The reduced prior shifts probability thresholds to minimize false positives.
The default human-trained model generalizes well to other vertebrates and most eukaryotes with U12-type introns. Performance may degrade for:
- Very distant lineages (plants, protists)
- Lineages with unusual U12-type motifs
- Species with extreme GC bias
Consider providing species-specific reference sequences for best results.
Memory scales with annotation density:
- Loads all intron sequences into memory
- Human genome (~250k introns): ~12-18 GB peak (varies by system and annotation density)
Dramatically reduced memory:
- Writes sequences to temporary SQLite database
- Keeps only scoring motifs in memory
- Human genome: ~2-3 GB peak (~85% reduction)
- Trade-off: Slightly slower I/O
The -p N flag parallelizes PWM scoring:
- Scoring is CPU-bound and embarrassingly parallel
- Linear speedup up to ~8-16 cores
- Diminishing returns beyond that
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/gkl842
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.182477.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.