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

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

$ 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 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:

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:

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 is a random-forest oriC predictor trained on the full DoriC 12.0 dataset. Head-to-head:

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:

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.