# oricle
Pure-Rust prediction of the bacterial replication origin (**oriC**) from genome
sequence — no database, no network, no alignment. Self-contained, with its own
input/output types and only one dependency (`needletail`, for the CLI's FASTA
parsing).
## Method
Ori-Finder–style, database-free:
1. **Cumulative GC-skew (Z-curve).** Along a bacterial chromosome the skew falls
on the replichore running from terminus to origin and rises on the other, so
its **global minimum locates oriC**.
2. **DnaA-box clustering.** The skew minimum is refined by the nearest cluster of
the DnaA-box 9-mer (`TTATCCACA` consensus + reverse complement, ≤1 mismatch).
3. **dnaA gene proximity** (optional). A supplied `dnaA` gene call refines or —
on genomes with a weak or misleading skew — rescues the prediction.
4. **Displaced-origin resolution.** Only ~43% of oriCs abut `dnaA`. When the box
cluster at the skew minimum sits well away from `dnaA` (*Enterobacteriaceae*,
~42 kb), it wins on skew quality; in the ambiguous 7–20 kb band it wins only
if it also carries an AT-rich **DNA-unwinding element** (*Vibrio*, *Aeromonas*).
The signal is gated so noisy, flat, or plasmid sequence is reported as **no
call** rather than guessed at.
## Library
```rust
use oricle::{detect, GeneHint};
let genome: Vec<u8> = std::fs::read("genome.fasta").map(strip_fasta).unwrap();
// Gene hints are optional; pass an empty slice to skip dnaA refinement.
let hints = [GeneHint { start: 3_882_326, end: 3_883_729, name: "dnaA".into() }];
for o in detect(&genome, &hints) {
println!("oriC {}..{} score {:.2} ({:?})", o.start, o.end, o.score, o.method);
}
# fn strip_fasta(_: Vec<u8>) -> Vec<u8> { Vec::new() }
```
`detect` returns at most one `OriC` per sequence (call it once per replicon).
Each `OriC` carries 1-based inclusive `start`/`end`, a `score` in `0..=1`, a
coarse `conf` class (`Confidence::Pass`/`Review`/`Weak`), the qualitative
`notes` (`Vec<Flag>`), the seeding `skew_min`, the `dnaa_boxes` count, and the
`Method` that produced it.
For non-default topology or thresholds, use `detect_with(seq, genes, &Options)`.
Set `Options::topology = Topology::Linear` for genuinely linear chromosomes
(*Streptomyces*, *Borreliella*), where the origin is internal and the `dnaA`
hint is the reliable anchor.
## CLI
```console
$ oricle genome.fasta
seqname start end score conf skew_min dnaa_boxes method notes
NC_000913.3 3925128 3926636 0.742 pass 3925597 7 SkewAndBoxes displaced
$ oricle --genes genes.tsv --all genome.fasta # dnaA hints + show no-calls
$ oricle --linear streptomyces.fasta # linear chromosome
```
`--genes` reads a tab-separated `seqname start end name` file (1-based); only
`dnaA` rows are used. `--min-sigma` and `--min-antipodal` expose the signal
gates for advanced tuning.
### Confidence and notes
The distribution of errors is **bimodal** — oricle either lands on the origin or
misses by tens of kb — so a single number can hide *why* a call is shaky. Two
columns make that explicit:
- **`conf`** is an action-oriented class derived from `score`, with the cut
points calibrated on the full held-out set against real ≤5 kb accuracy:
| `conf` | score | share of calls | lands ≤5 kb | action |
|---|---|---:|---:|---|
| `pass` | ≥ 0.60 | 67 % | **90.7 %** | trust it |
| `review` | 0.30–0.60 | 21 % | 78.8 % | glance before relying |
| `weak` | < 0.30 | 12 % | 45.8 % | likely wrong — verify or drop |
- **`notes`** is a `;`-separated list (`.` when clean) explaining the class.
Warnings — `multimodal` (rearranged skew, no single origin dip), `weak-skew`
(signal barely cleared the gate), `no-boxes` (skew window only, less precise),
`dnaa-disagree` (skew and `dnaA` conflict), `degenerate-box` (2-mismatch box
match) — are the ones to heed. `displaced` (oriC legitimately far from `dnaA`,
e.g. *E. coli*) and `no-hint` (sequence-only call) are informational.
So `pass` + `displaced` is a confident *E. coli*-style call, while
`weak` + `multimodal;no-boxes` is one to check by hand.
## Accuracy
Benchmarked against [DoriC](https://tubic.org/doric/) curated origins, with
coordinates re-derived on the current RefSeq accessions by exact oriC-sequence
matching (DoriC stores coordinates against older accession versions).
### Held-out accuracy (the number to trust)
Measured on **21,890 RefSeq chromosomes across 4,164 species / 1,363 genera, none
used to tune the thresholds**, with the `dnaA` gene supplied as a hint (as an
annotation pipeline would). This is the honest generalization estimate:
| median error | **0 bp** |
| ≤ 5 kb (genome-weighted) | **82.9 %** |
| ≤ 5 kb (per-species) | **78.3 %** |
| ≤ 5 kb (per-genus) | **74.1 %** |
| ≤ 50 kb | 92.5 % |
Major clinically-relevant genera score high: *Mycobacterium*, *Staphylococcus*,
*Mycobacteroides*, *Listeria* ~100 %; *Bacillus* 99 %, *Klebsiella* 95 %,
*E. coli* 89 %, *Salmonella* 89 %, *Aeromonas* 85 %, *Campylobacter* 81 %,
*Vibrio* (chr1) 79 %. The distribution is bimodal — oricle either lands **on**
the origin (median 0 bp) or misses by tens of kb — so treat a low `score` as a
real warning. Known weak spots (reported low-confidence): flat-skew
cyanobacteria, IS-rearranged *Bordetella* (likely DoriC artifacts), multi-replicon
secondary chromosomes (*Vibrio* chr2 uses RctB, not DnaA), a few clades whose
oriC sits 7–26 kb from *dnaA* without a resolvable DUE (*Neisseria*,
*Haemophilus*), and a subset of genomes whose skew minimum is itself displaced
from oriC toward a decoy box cluster — chiefly some high-GC Actinobacteria
(*Arthrobacter*, *Trueperella*, *Glutamicibacter*) and *Dehalococcoides*, where a
local AT dip beside the wrong cluster can mimic a DUE. Without any gene hint
(pure sequence) the median error is
**2.2 kb**, 57 % within 5 kb.
### In-sample numbers (optimistic — for reference only)
On the **615-replicon panel that the thresholds were tuned on**, oricle reaches
87.8 % ≤ 500 bp and 90.7 % ≤ 5 kb (median 1 bp) versus 0.2 % / 91.1 % (median
684 bp) for a strong annotation-only baseline that simply places oriC at the
annotated `dnaA` gene — strictly closer than that baseline on 558 of 591
replicons. These figures are **in-sample and therefore optimistic**; the held-out
numbers above are the ones to rely on. The gap is instructive: the tuning panel
under-represented *Enterobacteriaceae*, whose oriC genuinely sits ~42 kb from
`dnaA`.
**Diverse panel — 35 chromosomes**, 28–72 % GC and varied phyla (*E. coli*,
*B. subtilis*, *P. aeruginosa*, *Caulobacter*, *Streptomyces*, *M. tuberculosis*,
*Vibrio*, *Helicobacter*, …), scored with no gene hints to isolate the
sequence-only signal:
| 1 kb | 20 % | **49 %** |
| 5 kb | 31 % | **60 %** |
| 20 kb | 37 % | **63 %** |
| false positives (on 40 plasmids / secondary replicons) | 12 | **1** |
*E. coli* K-12 MG1655 is predicted at 3,925,128–3,926,636, overlapping both the
DoriC origin and NCBI's `rep_origin` annotation (3,925,744–3,925,975).
The v0.1.0 → v0.1.1 jump comes mainly from fixing a rotation-variant amplitude
gate that silently dropped every genome deposited starting at oriC (the majority
of RefSeq chromosomes), plus circular-window handling, a scale-free skew
threshold, terminus-antipodality gating, proximity-first box clustering,
linear-topology support, and dnaA-anchored box refinement.
### Comparison with ORCA
[ORCA](https://github.com/ZoyavanMeel/ORCA) is a random-forest oriC predictor
**trained on the full DoriC 12.0 dataset**. Head-to-head:
| 32 DOI-backed **experimental** oriCs (bias-free) | median error | **1 bp** | 11.0 kb |
| — same, midpoint-vs-point (like-for-like) | ≤ 5 kb | **59 %** | 42 % |
| held-out DoriC, *Enterobacteriaceae* (n=20) | ≤ 5 kb | **65 %** | 45 % |
oricle (which never sees DoriC) is more accurate on the bias-free experimental
set and even on ORCA's own DoriC training distribution — most clearly in
precision (median ~0–1 bp vs ~9–11 kb). The remaining oricle weaknesses are the
biologically expected ones: **flat-skew cyanobacteria** (no usable skew anchor,
where a database-trained model like ORCA wins) and a subset of genomes whose skew
minimum is itself displaced from oriC. `dnaA`-displaced clades are handled by
trusting the skew minimum when it is clean and strong (*Enterobacteriaceae*, ~42
kb) or, in the ambiguous 7–20 kb band, when it carries an AT-rich DUE (*Vibrio*,
*Aeromonas*); a residual few (*Neisseria*, *Haemophilus*) whose DUE is not
resolvable at this displacement remain a known weak spot.
## Testing
`cargo test` runs synthetic-profile unit tests (clean, noisy, high-GC,
flat, origin-at-coordinate-0, linear, plasmid-like). A real-genome regression
harness is opt-in:
```console
ORICLE_GENOME_DIR=/path/to/genomes cargo test --test regression -- --nocapture
```
where the directory holds `<accession>.fa` files (and optional
`<accession>.genes.tsv` dnaA hints). Expected origins are in
`tests/data/regression_truth.tsv`.
## License
MIT OR Apache-2.0, at your option.