oricle 0.2.0

//! Bacterial oriC (replication origin) prediction — DB-free, pure computation.
//!
//! Method (Ori-Finder style): the cumulative GC-skew (Z-curve) of a bacterial
//! chromosome falls along the replichore running from the terminus to the
//! origin and rises along the other, so its **global minimum sits at oriC**.
//! The minimum is refined with clusters of the DnaA-box 9-mer (`TTATCCACA`
//! consensus, plus reverse complement) and, when supplied, proximity to the
//! `dnaA` gene call.
//!
//! No database, no network, no alignment: [`detect`] takes bytes and returns
//! coordinates.
//!
//! ```
//! use oricle::{detect, GeneHint};
//!
//! # fn main() {
//! # let genome: Vec<u8> = Vec::new();
//! let hits = detect(&genome, &[]);          // gene hints are optional
//! for o in &hits {
//!     println!("oriC {}..{} (score {:.2})", o.start, o.end, o.score);
//! }
//! # }
//! ```
//!
//! # Scope and limits
//!
//! * Intended for **complete or near-complete replicons**. The skew curve is
//!   meaningless on short draft contigs, so sequences below
//!   [`Options::min_len`] are rejected outright.
//! * Assumes a **circular** replicon by default. For a linear chromosome
//!   (*Streptomyces*, *Borreliella*) set [`Options::topology`] to
//!   [`Topology::Linear`], which disables the antipodal terminus check.
//! * Some lineages genuinely have a flat or noisy skew (cyanobacteria,
//!   *Deinococcus*). Those are reported as **no call** rather than guessed at;
//!   supply a `dnaA` [`GeneHint`] to enable the gene-anchored fallback.
//! * Plasmids do not carry a DnaA-dependent oriC. The signal gates are tuned to
//!   leave them uncalled.

/// Replicon topology. Bacterial chromosomes are usually circular.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum Topology {
    /// Circular replicon (the default, and the usual bacterial chromosome).
    Circular,
    /// Linear replicon, e.g. *Streptomyces* or *Borreliella* chromosomes.
    Linear,
}

/// How a call was arrived at.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum Method {
    /// Cumulative skew minimum passed the signal gates and DnaA boxes cluster
    /// there. The strongest evidence.
    SkewAndBoxes,
    /// Skew minimum passed the gates but no convincing DnaA-box cluster was
    /// found; the region is a fixed window around the minimum.
    SkewOnly,
    /// The skew signal was too weak or too noisy to trust, and the call is
    /// anchored on the supplied `dnaA` gene instead.
    DnaaAnchored,
}

/// Coarse, action-oriented confidence class for a call, for quick triage and
/// filtering. Derived from [`OriC::score`] against thresholds calibrated on the
/// full held-out set so the class tracks real ≤5 kb accuracy: `Pass` calls are
/// right the large majority of the time, `Weak` calls are wrong more often than
/// not, and `Review` sits in between. The [`OriC::notes`] say *why*.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum Confidence {
    /// Trust it: high ≤5 kb accuracy in calibration.
    Pass,
    /// Borderline — worth a human look before relying on it.
    Review,
    /// Low confidence — likely wrong; verify or discard.
    Weak,
}

impl Confidence {
    /// Lower-case CLI token (`pass` / `review` / `weak`).
    pub fn as_str(self) -> &'static str {
        match self {
            Confidence::Pass => "pass",
            Confidence::Review => "review",
            Confidence::Weak => "weak",
        }
    }
}

/// A qualitative note attached to a call: either a *warning* that should lower a
/// reader's trust, or a contextual *fact* worth surfacing. Zero or more apply to
/// any call; they explain the [`OriC::conf`] class. Each has a stable
/// [`code`](Flag::code) for the CLI `notes` column.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum Flag {
    /// ⚠ The skew has a second trough nearly as deep as the global minimum — a
    /// rearranged genome with no single trustworthy origin dip (*Bordetella*).
    Multimodal,
    /// ⚠ The skew amplitude (or terminus antipodality) barely cleared the gate,
    /// so the skew-based position is shaky.
    WeakSkew,
    /// ⚠ No DnaA-box cluster was found; the region is a plain window around the
    /// skew minimum or the `dnaA` span, so it is less precise.
    NoBoxes,
    /// ⚠ A `dnaA` hint was given but the skew minimum points far elsewhere and
    /// the skew was not clean enough to trust the displacement — the two anchors
    /// genuinely conflict.
    DnaaDisagree,
    /// ⚠ The DnaA boxes were matched only at the relaxed 2-mismatch tolerance
    /// (high-GC drift); usually correct but weaker evidence.
    DegenerateBox,
    /// ℹ The origin sits well away from `dnaA` (E. coli ~42 kb, Vibrio ~10 kb).
    /// Not an error — a confident displaced origin — but worth knowing.
    Displaced,
    /// ℹ No `dnaA` gene hint was supplied; the call is sequence-only (typically
    /// less precise, median ~2 kb rather than 0 bp).
    NoHint,
}

impl Flag {
    /// Stable kebab-case token used in the CLI `notes` column.
    pub fn code(self) -> &'static str {
        match self {
            Flag::Multimodal => "multimodal",
            Flag::WeakSkew => "weak-skew",
            Flag::NoBoxes => "no-boxes",
            Flag::DnaaDisagree => "dnaa-disagree",
            Flag::DegenerateBox => "degenerate-box",
            Flag::Displaced => "displaced",
            Flag::NoHint => "no-hint",
        }
    }

    /// Whether this note is a warning (`true`) or purely informational
    /// (`false`). Warnings are the signals to be cautious about; the two
    /// informational notes (`Displaced`, `NoHint`) describe context.
    pub fn is_warning(self) -> bool {
        !matches!(self, Flag::Displaced | Flag::NoHint)
    }
}

/// Tuning knobs for [`detect_with`]. Start from [`Options::default`].
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Options {
    /// Minimum replicon length to attempt a call. Below this the skew curve is
    /// not meaningful. Default 50 kb.
    pub min_len: usize,
    /// Replicon topology. Default [`Topology::Circular`].
    pub topology: Topology,
    /// Required cumulative-skew amplitude, in units of `sqrt(n)` where `n` is
    /// the number of unambiguous bases.
    ///
    /// A sequence with no strand asymmetry is a random walk whose skew range
    /// grows as `sqrt(n)`, so this threshold is **scale-free**: it means the
    /// same thing for a 0.6 Mb and a 9 Mb replicon. Across a 35-chromosome
    /// panel real bacterial chromosomes score 6–66 on this scale, while
    /// plasmids and shuffled sequence sit near 1–4. Default 6.0.
    pub min_skew_sigma: f64,
    /// For circular replicons, how antipodal the skew maximum (terminus) must
    /// be to the minimum (origin), as a fraction of half the replicon length.
    /// `1.0` means exactly opposite. Real chromosomes with a usable skew score
    /// 0.59–1.00; the worst plasmid noise sits below 0.5. Default 0.55.
    /// Ignored when [`Topology::Linear`].
    pub min_antipodal: f64,
    /// Half-width (bp) of the window around the skew minimum searched for DnaA
    /// boxes. Default 50 kb.
    pub search_half_win: i64,
    /// Maximum gap (bp) between consecutive DnaA boxes in one cluster.
    /// Default 1 kb.
    pub cluster_gap: i64,
    /// Mismatches tolerated when matching the 9-mer DnaA box. Default 1.
    pub box_max_mismatch: u32,
    /// Minimum number of boxes for a cluster to define the reported region.
    /// Isolated single boxes are chance matches (a 9-mer with <=1 mismatch
    /// occurs by chance roughly every 5 kb per strand), so a lone box is not
    /// evidence. Default 3.
    pub min_cluster_boxes: usize,
    /// Half-width (bp) of the reported region when no box cluster anchors it.
    /// Default 250.
    pub region_half_win: i64,
}

impl Default for Options {
    fn default() -> Self {
        Options {
            min_len: 50_000,
            topology: Topology::Circular,
            min_skew_sigma: 6.0,
            min_antipodal: 0.55,
            search_half_win: 50_000,
            cluster_gap: 1_000,
            box_max_mismatch: 1,
            min_cluster_boxes: 3,
            region_half_win: 250,
        }
    }
}

/// Consensus DnaA box (*E. coli*) and its reverse complement.
const DNAA_BOX: &[u8; 9] = b"TTATCCACA";
const DNAA_BOX_RC: &[u8; 9] = b"TGTGGATAA";

/// Half-width (bp) of the DnaA-box search when anchored on the dnaA gene. oriC
/// abuts dnaA, so the window is tight to keep out distant chance clusters.
const DNAA_BOX_HALF_WIN: i64 = 5_000;
/// Maximum distance (bp) a box cluster may sit from the dnaA gene and still be
/// taken as the oriC. Beyond this the gene span itself is the better estimate.
const DNAA_BOX_MAX_DIST: i64 = 3_000;
/// How far a skew-minimum box cluster must lie from dnaA before it is read as a
/// displaced origin (E. coli, Chlamydia) that overrides the dnaA gene
/// unconditionally.
const DNAA_ORIC_FAR: i64 = 20_000;
/// A box cluster this far from dnaA (but nearer than [`DNAA_ORIC_FAR`]) is read
/// as a displaced origin only if it also carries an AT-rich DUE — the signature
/// that tells a genuine moderately-displaced oriC (Neisseria, Haemophilus,
/// ~7-15 kb from dnaA) from a chance cluster, which is what makes lowering the
/// distance threshold safe.
const DNAA_ORIC_MODERATE: i64 = 7_000;
/// Minimum AT-richness of the DUE window over the genome's baseline AT for a
/// moderately-displaced cluster to count as an origin. A true oriC's DUE runs
/// ~0.05 above baseline; chance clusters sit near 0. See [`due_excess`].
const DUE_MIN_EXCESS: f64 = 0.03;
/// Maximum genome-wide AT fraction for the moderate-displacement DUE test to
/// fire. In an already AT-saturated genome (*Campylobacter*, *Spiroplasma*,
/// *Borreliella*, GC < ~35 %) almost every window clears [`DUE_MIN_EXCESS`] by
/// chance, so the DUE loses all discriminating power and a decoy skew cluster
/// hijacks the (correct) dnaA-adjacent one. The excess signal only earns the
/// moderate branch where the baseline leaves headroom.
const DUE_MAX_GENOME_AT: f64 = 0.62;

/// Score at or above which a call is [`Confidence::Pass`], and below
/// [`CONF_WEAK`] a call is [`Confidence::Weak`] (in between, `Review`).
/// Calibrated on the full 21,890-replicon held-out set against ≤5 kb accuracy:
/// `Pass` (≥0.60) lands within 5 kb 90.7% of the time, `Review` 78.8%, `Weak`
/// (<0.30) only 45.8% — i.e. a `Weak` call is more often wrong than right, the
/// point below which the score stops being trustworthy. See the crate README.
const CONF_PASS: f32 = 0.60;
const CONF_WEAK: f32 = 0.30;

/// AT-richness of the DNA-unwinding element (DUE) beside a DnaA-box cluster, in
/// excess of the genome's baseline AT fraction. Every confirmed DnaA-dependent
/// oriC has an AT-rich DUE flanking its boxes; a chance box cluster does not.
/// Measured, like Ori-Finder's disparity scan, as the maximum AT fraction of a
/// 350 bp window within +-600 bp of `center`, minus `genome_at`. Absolute AT is
/// a weak discriminator (high-AT genomes are AT-rich everywhere); the *excess*
/// over baseline separates true from chance clusters ~2:1.
fn due_excess(seq: &[u8], center: i64, genome_at: f64, n: i64) -> f64 {
    const WIN: i64 = 350;
    let mut best = 0.0f64;
    let mut w0 = center - 600;
    while w0 <= center + 600 {
        let mut at = 0usize;
        for k in 0..WIN {
            match seq[(w0 + k).rem_euclid(n) as usize].to_ascii_uppercase() {
                b'A' | b'T' => at += 1,
                _ => {}
            }
        }
        best = best.max(at as f64 / WIN as f64);
        w0 += 100;
    }
    best - genome_at
}
/// How many more boxes a far skew cluster must hold than the dnaA-adjacent one
/// before it overrides it. A displaced true origin (E. coli oriC has ~9 boxes)
/// beats the handful near dnaA by a wide margin; the margin keeps a marginally
/// larger chance cluster from hijacking an otherwise correct dnaA call.
const DISPLACED_ORIC_MARGIN: usize = 4;
/// Skew-quality bar above which the cumulative-skew minimum is trusted to
/// override the dnaA gene when the two disagree. A clean, strong skew (high
/// amplitude and a near-antipodal terminus) pins oriC physically, so a box
/// cluster there but far from dnaA marks a genuinely displaced origin — the norm
/// across Enterobacteriaceae (E. coli, Klebsiella, Shigella), Bordetella and
/// Neisseria, where oriC sits tens to hundreds of kb from dnaA. These clades
/// dominate RefSeq, so the box-count margin alone (which favours the dnaA decoy
/// when the true cluster is small) systematically mislocated them.
const DISPLACED_SKEW_SIGMA: f64 = 10.0;
const DISPLACED_SKEW_ANTIPODAL: f64 = 0.80;

/// Minimum separation (as a fraction of replicon length) between the global
/// skew minimum and a rival minimum for the rival to count as a *distinct*
/// trough — i.e. evidence of a multimodal, rearranged skew rather than one clean
/// origin dip.
const MULTIMODAL_SEP_FRAC: f64 = 0.1;
/// If a distinct rival trough reaches at least this fraction of the global
/// minimum's depth, the skew is treated as multimodal and the call's confidence
/// is penalised. Heavily rearranged genomes (e.g. IS-element–scrambled
/// *Bordetella pertussis*) have several comparable minima and no single
/// trustworthy origin dip.
const MULTIMODAL_DEPTH_FRAC: f64 = 0.85;

/// A predicted gene position, so the caller can supply the `dnaA` location for
/// proximity refinement. Optional — pass an empty slice to skip.
#[derive(Debug, Clone, PartialEq, Eq)]
pub struct GeneHint {
    /// 1-based inclusive start coordinate.
    pub start: usize,
    /// 1-based inclusive end coordinate.
    pub end: usize,
    /// Gene symbol or product. `dnaA` is recognised by the name containing
    /// "dnaa" (case-insensitive), which also matches the usual product string
    /// "Chromosomal replication initiator protein DnaA".
    pub name: String,
}

/// A predicted replication origin.
#[derive(Debug, Clone, PartialEq)]
pub struct OriC {
    /// 1-based inclusive start coordinate.
    pub start: usize,
    /// 1-based inclusive end coordinate.
    pub end: usize,
    /// Confidence in 0..=1, combining skew strength, terminus antipodality,
    /// DnaA-box cluster size and `dnaA` proximity.
    pub score: f32,
    /// 1-based position of the cumulative skew minimum that seeded this call.
    pub skew_min: usize,
    /// Number of DnaA boxes in the cluster defining the region.
    pub dnaa_boxes: usize,
    /// Which evidence produced the call.
    pub method: Method,
    /// Action-oriented confidence class derived from `score`. See
    /// [`Confidence`].
    pub conf: Confidence,
    /// Qualitative notes explaining the confidence — warnings to heed and
    /// context worth knowing. See [`Flag`]. Empty for a clean, unremarkable
    /// call.
    pub notes: Vec<Flag>,
}

/// Summary of the cumulative GC-skew curve.
#[derive(Debug, Clone, Copy)]
struct SkewProfile {
    /// 0-based index of the global minimum of the cumulative skew.
    min_pos: usize,
    /// 0-based index of the global maximum.
    max_pos: usize,
    /// max - min, the peak-to-trough amplitude. Invariant under rotation of a
    /// circular sequence, unlike `|min|`.
    range: i64,
    /// Count of unambiguous A/C/G/T bases.
    n_valid: usize,
    /// Count of G+C bases, for the genome GC fraction.
    gc_count: usize,
    /// True when a second, well-separated trough reaches nearly the depth of the
    /// global minimum — the signature of a rearranged, multimodal skew with no
    /// single trustworthy origin dip.
    multimodal: bool,
}

/// Detect oriC in one replicon with default [`Options`] (circular topology).
///
/// `genes` supplies gene calls so a `dnaA` hint can refine or rescue the call;
/// it may be empty. Returns at most one [`OriC`], or an empty vector when no
/// candidate clears the signal gates.
pub fn detect(seq: &[u8], genes: &[GeneHint]) -> Vec<OriC> {
    detect_with(seq, genes, &Options::default())
}

/// Detect oriC in one replicon with explicit [`Options`].
pub fn detect_with(seq: &[u8], genes: &[GeneHint], opt: &Options) -> Vec<OriC> {
    match detect_one(seq, genes, opt) {
        Some(o) => vec![o],
        None => Vec::new(),
    }
}

fn detect_one(seq: &[u8], genes: &[GeneHint], opt: &Options) -> Option<OriC> {
    let n = seq.len();
    if n < opt.min_len {
        return None;
    }

    let prof = skew_profile(seq);
    if prof.n_valid < opt.min_len {
        return None;
    }

    let dnaa_span = dnaa_gene(genes);
    let dnaa_center = dnaa_span.map(|(s, e)| (s + e) / 2);

    // --- Signal gates -----------------------------------------------------
    // 1. Amplitude, normalised against the sqrt(n) of an unbiased random walk.
    let sigma = prof.range as f64 / (prof.n_valid as f64).sqrt();
    let strong_amplitude = sigma >= opt.min_skew_sigma;

    // 2. For a circular chromosome the terminus (skew maximum) must sit roughly
    //    opposite the origin. A curve that fails this is dominated by local
    //    composition, not by replication strand bias.
    let antipodal = match opt.topology {
        Topology::Linear => 1.0,
        Topology::Circular => {
            let d = circular_dist(prof.min_pos as i64, prof.max_pos as i64, n as i64);
            d as f64 / (n as f64 / 2.0)
        }
    };
    let skew_usable = strong_amplitude && antipodal >= opt.min_antipodal;

    // --- Locate oriC by reconciling two independent anchors ---------------
    // Two signals point at oriC, each precise in different regimes:
    //
    //  * The `dnaA` gene. In nearly every bacterium oriC abuts dnaA, so an
    //    annotated dnaA plus the DnaA-box cluster beside it pins oriC to within
    //    a few hundred bp. Across a 615-replicon RefSeq panel this reproduces
    //    the curated (DoriC) origin with a *median error of 1 bp*.
    //  * The cumulative-skew minimum. Independent of annotation, accurate to a
    //    few kb, and — crucially — it stays correct for the rare species whose
    //    oriC sits far from dnaA (E. coli ~42 kb, Chlamydia ~440 kb), exactly
    //    where the dnaA anchor fails.
    //
    // Preference order: a DnaA-box cluster hugging the dnaA gene wins; else a
    // box cluster at the skew minimum; else the skew-minimum window; else the
    // dnaA gene span as an annotation-only estimate.
    let high_gc =
        prof.n_valid > 0 && (prof.gc_count as f64 / prof.n_valid as f64) >= BOX_FALLBACK_MIN_GC;
    let dnaa_hit = dnaa_center.and_then(|c| {
        // The 2-mismatch fallback is for the *high-GC, oriC-at-dnaA* case
        // (Actinobacteria: Mycobacterium, Corynebacterium) where the box has
        // drifted. Requiring the skew minimum to also abut dnaA excludes high-GC
        // *displaced*-origin clades (Rhizobiales, Sphingomonadales), where a
        // loose scan near dnaA would fabricate a decoy and hide the real origin.
        let skew_at_dnaa = circular_dist(prof.min_pos as i64, c, n as i64) <= DNAA_ORIC_FAR;
        let allow_loose = high_gc && skew_at_dnaa;
        best_cluster_near(seq, c, DNAA_BOX_HALF_WIN, n as i64, allow_loose, opt).filter(
            |&(s, e, _, _)| {
                circular_dist((s as i64 + e as i64) / 2, c, n as i64) <= DNAA_BOX_MAX_DIST
            },
        )
    });
    // Whether the dnaA-adjacent cluster relied on the relaxed 2-mismatch scan.
    let dnaa_used_loose = dnaa_hit.map(|(_, _, _, loose)| loose).unwrap_or(false);
    let dnaa_cluster = dnaa_hit.map(|(s, e, c, _)| (s, e, c));
    // The dnaA-box search around the skew minimum stays at the strict tolerance:
    // its window is wide (tens of kb), where the looser 2-mismatch match would
    // admit chance clusters far from the true minimum. The 2-mismatch fallback is
    // confined to the tight dnaA-gene window above, where a hit is trustworthy.
    let skew_cluster = if skew_usable {
        let a = prof.min_pos as i64;
        let hits = find_dnaa_boxes_around(seq, a, opt.search_half_win, opt.box_max_mismatch, opt);
        best_cluster(&hits, a, n as i64, opt)
    } else {
        None
    };

    // A box cluster at the skew minimum that sits far from dnaA is the tell-tale
    // of a displaced origin (E. coli, Chlamydia): the boxes, not the gene, mark
    // oriC. Beyond DNAA_ORIC_FAR that is taken at face value; in the ambiguous
    // 7-20 kb band it counts only when the cluster also carries an AT-rich DUE,
    // which distinguishes a genuinely displaced origin (Neisseria, Haemophilus)
    // from a chance cluster. When the skew cluster is near dnaA it is merely
    // corroborating the gene, and the tighter dnaA-anchored estimate is kept.
    let genome_at = if prof.n_valid > 0 {
        1.0 - prof.gc_count as f64 / prof.n_valid as f64
    } else {
        0.0
    };
    let skew_cluster_far = skew_cluster.filter(|&(s, e, _)| {
        dnaa_center.is_none_or(|c| {
            let center = (s as i64 + e as i64) / 2;
            let d = circular_dist(center, c, n as i64);
            d > DNAA_ORIC_FAR
                || (d > DNAA_ORIC_MODERATE
                    && genome_at <= DUE_MAX_GENOME_AT
                    && due_excess(seq, center, genome_at, n as i64) >= DUE_MIN_EXCESS)
        })
    });

    // Decide between the two clusters. A box cluster far from dnaA wins when
    // EITHER the skew that placed it is clean and strong (a physically reliable
    // displaced origin — the Enterobacteriaceae / Bordetella / Neisseria case,
    // where the true oriC cluster can be small yet correct), OR it simply holds
    // many more boxes than the dnaA-adjacent one. Otherwise the tighter
    // dnaA-anchored estimate is kept.
    // Only a *circular* genome earns the skew-quality override: there the
    // antipodal terminus is real evidence the minimum is the origin. On a linear
    // replicon the skew minimum is unreliable (antipodal is forced to 1.0), so
    // the dnaA anchor must keep priority.
    let skew_clean = opt.topology == Topology::Circular
        && sigma >= DISPLACED_SKEW_SIGMA
        && antipodal >= DISPLACED_SKEW_ANTIPODAL;
    // A far skew cluster overrides the dnaA-adjacent one when the skew is clean
    // AND the far cluster is at least as strong, OR when it simply holds many
    // more boxes. The `f.2 >= d.2` guard under `skew_clean` is what keeps a
    // *chance* 1-mismatch cluster (3 boxes) from overriding a genuine, denser
    // dnaA-adjacent cluster recovered at 2 mismatches (Mycobacteroides,
    // Leptospira, high-GC clades) — while still letting E. coli's strong
    // displaced cluster beat its weak dnaA decoy.
    let far_wins = |d: (usize, usize, usize), f: (usize, usize, usize)| {
        (skew_clean && f.2 >= d.2) || f.2 >= d.2 + DISPLACED_ORIC_MARGIN
    };
    let picked = match (dnaa_cluster, skew_cluster_far) {
        (Some(d), Some(f)) if far_wins(d, f) => Some((f, Method::SkewAndBoxes)),
        (Some(d), _) => Some((d, Method::DnaaAnchored)),
        (None, Some(f)) => Some((f, Method::SkewAndBoxes)),
        (None, None) => None,
    };

    let (start0, end0, boxes, method) = if let Some(((s, e, count), m)) = picked {
        (s, e, count, m)
    } else if let Some((gs, ge)) = dnaa_span {
        // oriC abuts dnaA but its boxes are weak/absent: the gene span is the
        // best annotation-only estimate.
        (gs as usize, ge as usize, 0, Method::DnaaAnchored)
    } else if skew_usable {
        match skew_cluster {
            Some((s, e, count)) => (s, e, count, Method::SkewAndBoxes),
            None => {
                let a = prof.min_pos as i64;
                let (s, e) = (a - opt.region_half_win, a + opt.region_half_win);
                (
                    wrap(s, n as i64, opt) as usize,
                    wrap(e, n as i64, opt) as usize,
                    0,
                    Method::SkewOnly,
                )
            }
        }
    } else {
        return None;
    };

    // How well dnaA agrees with the reported region (for scoring only).
    let region_center = (start0 as i64 + end0 as i64) / 2;
    let dnaa_dist = dnaa_center.map(|c| circular_dist(c, region_center, n as i64));
    let score = confidence(sigma, antipodal, boxes, dnaa_dist, prof.multimodal, opt);

    // Qualitative notes: warnings to heed and context worth surfacing. These
    // explain the confidence class; see [`Flag`].
    let mut notes = Vec::new();
    if prof.multimodal {
        notes.push(Flag::Multimodal);
    }
    // A weak skew only matters when the skew actually placed the call.
    if method != Method::DnaaAnchored
        && (sigma < opt.min_skew_sigma * 1.5
            || (opt.topology == Topology::Circular && antipodal < opt.min_antipodal + 0.10))
    {
        notes.push(Flag::WeakSkew);
    }
    if boxes == 0 {
        notes.push(Flag::NoBoxes);
    }
    if dnaa_used_loose && method == Method::DnaaAnchored {
        notes.push(Flag::DegenerateBox);
    }
    match (dnaa_center, dnaa_dist) {
        (None, _) => notes.push(Flag::NoHint),
        (Some(c), d) => {
            let resolved_displaced =
                method == Method::SkewAndBoxes && d.is_some_and(|d| d > DNAA_ORIC_MODERATE);
            if resolved_displaced {
                notes.push(Flag::Displaced);
            } else if circular_dist(prof.min_pos as i64, c, n as i64) > DNAA_ORIC_FAR && !skew_clean
            {
                // The skew minimum and dnaA point far apart and the skew was not
                // clean enough to trust the displacement — an unresolved conflict.
                notes.push(Flag::DnaaDisagree);
            }
        }
    }

    let conf = if score >= CONF_PASS {
        Confidence::Pass
    } else if score < CONF_WEAK {
        Confidence::Weak
    } else {
        Confidence::Review
    };

    Some(OriC {
        start: start0 + 1,
        end: end0 + 1,
        score,
        skew_min: prof.min_pos + 1,
        dnaa_boxes: boxes,
        method,
        conf,
        notes,
    })
}

/// Cumulative GC-skew: `G` = +1, `C` = -1, other bases ignored. Records where
/// the running sum is globally minimal (oriC) and maximal (terminus).
fn skew_profile(seq: &[u8]) -> SkewProfile {
    let mut skew: i64 = 0;
    let mut min_val: i64 = 0;
    let mut max_val: i64 = 0;
    let mut min_pos: usize = 0;
    let mut max_pos: usize = 0;
    let mut n_valid: usize = 0;
    let mut gc_count: usize = 0;
    for (i, &b) in seq.iter().enumerate() {
        match b.to_ascii_uppercase() {
            b'G' => {
                skew += 1;
                n_valid += 1;
                gc_count += 1;
            }
            b'C' => {
                skew -= 1;
                n_valid += 1;
                gc_count += 1;
            }
            b'A' | b'T' => n_valid += 1,
            _ => {}
        }
        if skew < min_val {
            min_val = skew;
            min_pos = i;
        }
        if skew > max_val {
            max_val = skew;
            max_pos = i;
        }
    }

    // Second pass: the deepest trough well-separated from the global minimum. If
    // it comes close to the global depth the skew is multimodal (rearranged).
    let range = max_val - min_val;
    let n = seq.len() as i64;
    let sep = (seq.len() as f64 * MULTIMODAL_SEP_FRAC) as i64;
    let mut skew2: i64 = 0;
    let mut rival_min: i64 = max_val;
    for (i, &b) in seq.iter().enumerate() {
        match b.to_ascii_uppercase() {
            b'G' => skew2 += 1,
            b'C' => skew2 -= 1,
            _ => {}
        }
        if skew2 < rival_min && circular_dist(i as i64, min_pos as i64, n) > sep {
            rival_min = skew2;
        }
    }
    let multimodal =
        range > 0 && (rival_min - min_val) as f64 <= (1.0 - MULTIMODAL_DEPTH_FRAC) * range as f64;

    SkewProfile {
        min_pos,
        max_pos,
        range,
        n_valid,
        gc_count,
        multimodal,
    }
}

/// Shortest distance between two coordinates, wrapping for circular replicons.
fn circular_dist(a: i64, b: i64, n: i64) -> i64 {
    let d = (a - b).abs();
    d.min(n - d)
}

/// Clamp (linear) or wrap (circular) a coordinate into `0..n`.
fn wrap(x: i64, n: i64, opt: &Options) -> i64 {
    match opt.topology {
        Topology::Circular => x.rem_euclid(n),
        Topology::Linear => x.clamp(0, n - 1),
    }
}

/// Collect DnaA-box positions within `half` bp of `center`.
///
/// The window is materialised so that a circular replicon wraps correctly:
/// when oriC sits at coordinate ~0 (the common convention of depositing a
/// genome starting at the origin) half the window lies at the far end of the
/// sequence, and a boxes-only scan of a clamped slice would miss it.
fn find_dnaa_boxes_around(
    seq: &[u8],
    center: i64,
    half: i64,
    mism: u32,
    opt: &Options,
) -> Vec<usize> {
    let n = seq.len() as i64;
    let (lo, hi) = match opt.topology {
        Topology::Circular => (center - half, center + half),
        Topology::Linear => ((center - half).max(0), (center + half).min(n)),
    };
    let span = (hi - lo).min(n);

    let mut window = Vec::with_capacity(span as usize + 8);
    for k in 0..(span + 8) {
        let idx = match opt.topology {
            Topology::Circular => (lo + k).rem_euclid(n),
            Topology::Linear => {
                if lo + k >= n {
                    break;
                }
                lo + k
            }
        };
        window.push(seq[idx as usize]);
    }

    let mut hits = Vec::new();
    if window.len() < 9 {
        return hits;
    }
    for i in 0..=window.len() - 9 {
        let kmer = &window[i..i + 9];
        if matches_box(kmer, DNAA_BOX, mism) || matches_box(kmer, DNAA_BOX_RC, mism) {
            let abs = match opt.topology {
                Topology::Circular => (lo + i as i64).rem_euclid(n),
                Topology::Linear => lo + i as i64,
            };
            hits.push(abs as usize);
        }
    }
    hits.sort_unstable();
    hits
}

/// True if `kmer` matches `pat` with at most `max_mismatch` mismatches.
fn matches_box(kmer: &[u8], pat: &[u8; 9], max_mismatch: u32) -> bool {
    let mut mism = 0u32;
    for j in 0..9 {
        if kmer[j].to_ascii_uppercase() != pat[j] {
            mism += 1;
            if mism > max_mismatch {
                return false;
            }
        }
    }
    true
}

/// Bonus, in bp, that one extra DnaA box is worth against distance when
/// choosing among clusters. A cluster this much farther from the anchor is
/// preferred only if it holds one more box.
const BOX_WORTH_BP: i64 = 3_000;

/// Looser DnaA-box match tolerance tried only when the strict
/// ([`Options::box_max_mismatch`]) scan finds no qualifying cluster. The DnaA
/// box drifts from the *E. coli* `TTATCCACA` consensus in high-GC and other
/// clades — across a 21k-genome panel the true oriC of *Mycobacterium* (99%),
/// *Mycobacteroides* and *Leptospira* is invisible at 1 mismatch but recovered
/// at 2. Genomes whose oriC is displaced from dnaA (Enterobacteriaceae) have no
/// near-anchor cluster at *any* tolerance, so the fallback leaves them untouched.
const BOX_FALLBACK_MISMATCH: u32 = 2;

/// Genome GC fraction above which the looser [`BOX_FALLBACK_MISMATCH`] tolerance
/// is enabled. The DnaA box drifts from `TTATCCACA` mainly in **high-GC** clades
/// (Mycobacterium, Mycobacteroides ~65% GC), where the true oriC cluster is
/// invisible at 1 mismatch. Restricting the fallback to high-GC genomes leaves
/// ~50% GC lineages (Enterobacteriaceae) untouched — there a 2-mismatch scan near
/// dnaA would otherwise fabricate a decoy cluster and mislocate the displaced
/// origin.
const BOX_FALLBACK_MIN_GC: f64 = 0.60;

/// Best qualifying DnaA-box cluster within `half` bp of `center`, ranked toward
/// `center`. Tries the strict consensus first, then — only for a high-GC genome
/// (`allow_loose`) whose true boxes have drifted — the looser
/// [`BOX_FALLBACK_MISMATCH`] tolerance if the strict scan came up empty.
/// Returns the cluster and, in the `bool`, whether it came from the relaxed
/// 2-mismatch fallback (`true`) rather than the strict scan.
fn best_cluster_near(
    seq: &[u8],
    center: i64,
    half: i64,
    n: i64,
    allow_loose: bool,
    opt: &Options,
) -> Option<(usize, usize, usize, bool)> {
    let strict = find_dnaa_boxes_around(seq, center, half, opt.box_max_mismatch, opt);
    if let Some((s, e, c)) = best_cluster(&strict, center, n, opt) {
        return Some((s, e, c, false));
    }
    if allow_loose && opt.box_max_mismatch < BOX_FALLBACK_MISMATCH {
        let loose = find_dnaa_boxes_around(seq, center, half, BOX_FALLBACK_MISMATCH, opt);
        return best_cluster(&loose, center, n, opt).map(|(s, e, c)| (s, e, c, true));
    }
    None
}

/// Group boxes separated by at most `cluster_gap` and return the best cluster
/// as `(start0, end0, count)`.
///
/// The anchor (skew minimum, or `dnaA` gene) already locates oriC to within a
/// few kb; the DnaA-box cluster only *refines* it. So clusters are ranked
/// primarily by **proximity to the anchor**, with a modest bonus per box: a
/// larger cluster wins only if it is not much farther away. Ranking by raw box
/// count instead lets a chance cluster tens of kb away hijack the call — the
/// bug that put the *Campylobacter* and *Borreliella* origins ~35 kb off.
///
/// Only clusters of at least [`Options::min_cluster_boxes`] are considered.
/// Ranking first and filtering afterwards was a bug: a chance single box sitting
/// right on the anchor has the lowest cost, wins, and is then discarded by the
/// count filter — throwing away a genuine oriC cluster a little farther out (the
/// *Sinorhizobium* / Rhizobiales case, where the origin box cluster sits ~25 kb
/// from the skew minimum).
fn best_cluster(
    hits: &[usize],
    anchor: i64,
    n: i64,
    opt: &Options,
) -> Option<(usize, usize, usize)> {
    if hits.is_empty() {
        return None;
    }
    let mut best: Option<(usize, usize, usize, i64)> = None; // start,end,count,cost
    let mut i = 0;
    while i < hits.len() {
        let mut j = i;
        while j + 1 < hits.len() && (hits[j + 1] as i64 - hits[j] as i64) <= opt.cluster_gap {
            j += 1;
        }
        let (start, end, count) = (hits[i], hits[j] + 8, j - i + 1);
        i = j + 1;
        if count < opt.min_cluster_boxes {
            continue;
        }
        let center = (start as i64 + end as i64) / 2;
        let dist = circular_dist(center, anchor, n);
        let cost = dist - count as i64 * BOX_WORTH_BP;
        let better = match &best {
            None => true,
            Some((_, _, _, bcost)) => cost < *bcost,
        };
        if better {
            best = Some((start, end, count, cost));
        }
    }
    best.map(|(s, e, c, _)| (s, e, c))
}

/// Blend the independent lines of evidence into a 0..=1 confidence.
fn confidence(
    sigma: f64,
    antipodal: f64,
    boxes: usize,
    dnaa_dist: Option<i64>,
    multimodal: bool,
    opt: &Options,
) -> f32 {
    // Skew amplitude: 0 at the gate, saturating well above it.
    let amp_q = norm(sigma, opt.min_skew_sigma, opt.min_skew_sigma * 4.0);
    // Terminus antipodality.
    let anti_q = norm(antipodal, opt.min_antipodal, 1.0);
    // Box cluster size, saturating: 3 boxes ~0.5, 10 boxes ~0.77.
    let box_q = boxes as f64 / (boxes as f64 + 3.0);
    // dnaA proximity, if a hint was given. Absent hint is neutral, not a
    // penalty: most callers have no annotation.
    let dnaa_q = dnaa_dist.map(|d| 1.0 - norm(d as f64, 0.0, opt.search_half_win as f64));

    let (mut sum, mut wsum) = (0.45 * amp_q + 0.20 * anti_q + 0.35 * box_q, 1.0);
    if let Some(q) = dnaa_q {
        sum += 0.30 * q;
        wsum += 0.30;
    }
    let mut score = (sum / wsum) as f32;
    // A multimodal (rearranged) skew has no single trustworthy origin dip;
    // halve the confidence to flag the call as low-quality without suppressing
    // it outright.
    if multimodal {
        score *= 0.5;
    }
    score.clamp(0.0, 1.0)
}

/// Linear ramp from 0 at `lo` to 1 at `hi`, clamped.
fn norm(x: f64, lo: f64, hi: f64) -> f64 {
    if hi <= lo {
        return 1.0;
    }
    ((x - lo) / (hi - lo)).clamp(0.0, 1.0)
}

/// 0-based inclusive span of the first `dnaA` gene call, if present.
fn dnaa_gene(genes: &[GeneHint]) -> Option<(i64, i64)> {
    genes.iter().find(|g| contains_dnaa(&g.name)).map(|g| {
        let s = (g.start as i64 - 1).max(0);
        let e = (g.end as i64 - 1).max(s);
        (s, e)
    })
}

/// Case-insensitive substring test for "dnaa".
fn contains_dnaa(s: &str) -> bool {
    s.to_ascii_lowercase().contains("dnaa")
}

#[cfg(test)]
pub(crate) mod tests_util {
    //! Deterministic synthetic-genome builder for the unit tests. No `rand`
    //! dependency: a tiny xorshift PRNG keeps every genome reproducible.

    /// xorshift64 — deterministic, seeded, good enough for test noise.
    pub struct Rng(u64);
    impl Rng {
        pub fn new(seed: u64) -> Self {
            Rng(seed ^ 0x9E37_79B9_7F4A_7C15 | 1)
        }
        pub fn next(&mut self) -> u64 {
            let mut x = self.0;
            x ^= x << 13;
            x ^= x >> 7;
            x ^= x << 17;
            self.0 = x;
            x
        }
        /// A random base; `gc` in 0..=100 sets the G/C probability.
        pub fn base(&mut self, gc: u8) -> u8 {
            if (self.next() % 100) < gc as u64 {
                if self.next() & 1 == 0 {
                    b'G'
                } else {
                    b'C'
                }
            } else if self.next() & 1 == 0 {
                b'A'
            } else {
                b'T'
            }
        }
    }

    /// Parameters for a synthetic circular chromosome with a V-shaped skew.
    pub struct Synth {
        /// Length in bp of each replichore arm. Total length ~= 2 * arm.
        pub arm: usize,
        /// Net skew strength: fraction (0..=100) of arm positions that carry the
        /// strand-biasing base (C on the descending arm, G on the ascending).
        pub skew_pct: u8,
        /// Background GC content (0..=100) of the non-biasing positions, so
        /// high-GC genomes can be exercised.
        pub gc_bg: u8,
        /// Number of consensus DnaA boxes to embed at the skew minimum.
        pub boxes: usize,
        /// Rotate the finished sequence left by this many bp, simulating a
        /// genome deposited starting at an arbitrary coordinate.
        pub rot: usize,
        pub seed: u64,
    }

    impl Default for Synth {
        fn default() -> Self {
            Synth {
                arm: 60_000,
                skew_pct: 40,
                gc_bg: 50,
                boxes: 5,
                rot: 0,
                seed: 1,
            }
        }
    }

    impl Synth {
        /// Build the genome and return `(sequence, oriC_position)` where
        /// `oriC_position` is the 0-based index of the embedded box cluster
        /// (== the intended skew minimum) *before* rotation is accounted for.
        pub fn build(&self) -> (Vec<u8>, usize) {
            let mut rng = Rng::new(self.seed);
            let mut seq = Vec::with_capacity(2 * self.arm + 64);
            // Descending arm: biasing base is C.
            for _ in 0..self.arm {
                if (rng.next() % 100) < self.skew_pct as u64 {
                    seq.push(b'C');
                } else {
                    seq.push(rng.base(self.gc_bg));
                }
            }
            let ori = seq.len();
            for _ in 0..self.boxes {
                seq.extend_from_slice(b"TTATCCACA");
                seq.extend_from_slice(b"AACGT"); // spacer
            }
            // Ascending arm: biasing base is G.
            for _ in 0..self.arm {
                if (rng.next() % 100) < self.skew_pct as u64 {
                    seq.push(b'G');
                } else {
                    seq.push(rng.base(self.gc_bg));
                }
            }
            let n = seq.len();
            let rot = self.rot % n;
            seq.rotate_left(rot);
            // oriC index shifts left by rot (mod n).
            let ori = (ori + n - rot) % n;
            (seq, ori)
        }
    }

    /// Convenience: a clean default genome and its oriC index.
    pub fn clean() -> (Vec<u8>, usize) {
        Synth::default().build()
    }
}

#[cfg(test)]
mod tests {
    use super::tests_util::*;
    use super::*;

    /// Circular distance from a 1-based predicted point to a 0-based oriC index.
    fn off(pred_1based: usize, ori_0based: usize, n: usize) -> usize {
        let d = (pred_1based as i64 - 1 - ori_0based as i64).abs();
        d.min(n as i64 - d) as usize
    }

    // --- pure helpers -----------------------------------------------------

    #[test]
    fn dnaa_recognition_variants() {
        assert!(contains_dnaa("dnaA"));
        assert!(contains_dnaa("DnaA"));
        assert!(contains_dnaa(
            "Chromosomal replication initiator protein DnaA"
        ));
        assert!(!contains_dnaa("gyrA"));
        assert!(!contains_dnaa("DNA gyrase subunit A"));
    }

    #[test]
    fn circular_dist_wraps() {
        assert_eq!(circular_dist(10, 90, 100), 20);
        assert_eq!(circular_dist(90, 10, 100), 20);
        assert_eq!(circular_dist(10, 20, 100), 10);
    }

    #[test]
    fn matches_box_tolerates_one_mismatch() {
        assert!(matches_box(b"TTATCCACA", DNAA_BOX, 1));
        assert!(matches_box(b"TTATCCACC", DNAA_BOX, 1)); // 1 mismatch
        assert!(!matches_box(b"TTATCCAGG", DNAA_BOX, 1)); // 2 mismatches
        assert!(!matches_box(b"TTATCCACC", DNAA_BOX, 0)); // strict
    }

    #[test]
    fn best_cluster_prefers_proximity_over_size() {
        // A big cluster far from the anchor must not beat a small one on it —
        // the Campylobacter / Borreliella failure mode.
        let n = 1_000_000;
        let anchor = 500_000;
        let mut hits = vec![anchor as usize, anchor as usize + 30, anchor as usize + 60];
        // a larger, distant chance cluster
        for k in 0..6 {
            hits.push(100_000 + k * 30);
        }
        hits.sort_unstable();
        let opt = Options::default();
        let (s, e, _) = best_cluster(&hits, anchor, n, &opt).unwrap();
        let center = ((s + e) / 2) as i64;
        assert!(
            circular_dist(center, anchor, n) < 5_000,
            "cluster {center} should sit near anchor {anchor}, not at the larger distant group"
        );
    }

    #[test]
    fn best_cluster_skips_subthreshold_near_anchor() {
        // A chance single box right on the anchor must not shadow a genuine
        // cluster a little farther out — the Sinorhizobium / Rhizobiales bug,
        // where the origin box cluster sits ~25 kb from the skew minimum and
        // was discarded because a lone box outscored it, then failed the count
        // filter, leaving no call.
        let n = 4_000_000;
        let anchor = 3_600_000;
        let mut hits = vec![anchor as usize]; // lone chance box on the anchor
        for k in 0..3 {
            hits.push(25_000 + k * 30); // real 3-box oriC cluster ~25 kb away
        }
        hits.sort_unstable();
        let opt = Options::default();
        let (s, e, count) =
            best_cluster(&hits, anchor, n, &opt).expect("must find the real cluster");
        assert_eq!(count, 3);
        let center = ((s + e) / 2) as i64;
        assert!(
            circular_dist(center, 25_000, n) < 1_000,
            "should report the 3-box cluster, not the lone box"
        );
    }

    #[test]
    fn find_boxes_wraps_around_origin() {
        // A box straddling coordinate 0 on a circular replicon must be found
        // when the search window centred at ~0 wraps to the sequence end.
        let mut seq = vec![b'A'; 200_000];
        // place a box at the very end so the window around center=0 wraps to it
        seq[199_995..].copy_from_slice(&b"TTATC"[..]);
        seq[..4].copy_from_slice(&b"CACA"[..]); // ...continues past 0: TTATC|CACA
        let opt = Options::default();
        let hits = find_dnaa_boxes_around(&seq, 0, opt.search_half_win, opt.box_max_mismatch, &opt);
        assert!(
            hits.iter().any(|&h| h >= 199_990 || h <= 5),
            "expected a wrapped box near the origin, got {hits:?}"
        );
    }

    // --- detection behaviour ---------------------------------------------

    #[test]
    fn detects_clean_oric_with_boxes() {
        let (seq, ori) = clean();
        let n = seq.len();
        let out = detect(&seq, &[]);
        assert_eq!(out.len(), 1, "expected exactly one oriC");
        let o = &out[0];
        assert_eq!(o.method, Method::SkewAndBoxes);
        assert!(
            o.dnaa_boxes >= 3,
            "should anchor on the embedded box cluster"
        );
        assert!(
            off(o.start, ori, n) < 2_000,
            "oriC {} should be within 2 kb of the true origin {ori}",
            o.start
        );
        assert!(
            o.score > 0.6,
            "clean signal should score high, got {}",
            o.score
        );
    }

    #[test]
    fn rotation_invariant_when_deposited_at_origin() {
        // v0.1.0 regression: a genome numbered starting at oriC gave |skew_min|
        // ~= 0 and was silently dropped. It must now be found regardless of
        // where the sequence is cut.
        for rot in [0usize, 30_000, 60_000, 90_000, 120_000] {
            let (seq, ori) = Synth {
                rot,
                ..Default::default()
            }
            .build();
            let n = seq.len();
            let out = detect(&seq, &[]);
            assert_eq!(out.len(), 1, "rot={rot}: expected one call");
            assert!(
                off(out[0].start, ori, n) < 2_000,
                "rot={rot}: oriC {} not near origin {ori}",
                out[0].start
            );
        }
    }

    #[test]
    fn high_gc_genome_still_detected() {
        // 70% GC background must not defeat the scale-free amplitude gate.
        let (seq, ori) = Synth {
            gc_bg: 70,
            ..Default::default()
        }
        .build();
        let n = seq.len();
        let out = detect(&seq, &[]);
        assert_eq!(out.len(), 1, "high-GC genome should still be called");
        assert!(off(out[0].start, ori, n) < 3_000);
    }

    #[test]
    fn noisy_skew_still_detected() {
        // A weak (skew_pct=12), noisy V with a real box cluster should still be
        // located near the origin.
        let (seq, ori) = Synth {
            skew_pct: 12,
            boxes: 6,
            seed: 42,
            ..Default::default()
        }
        .build();
        let n = seq.len();
        let out = detect(&seq, &[]);
        assert_eq!(out.len(), 1, "noisy but real signal should be called");
        assert!(
            off(out[0].start, ori, n) < 5_000,
            "oriC {} not within 5 kb of {ori} on noisy genome",
            out[0].start
        );
    }

    #[test]
    fn skew_only_call_without_boxes() {
        // Clean V, no DnaA boxes, no gene: the skew minimum alone is trusted.
        let (seq, _) = Synth {
            boxes: 0,
            ..Default::default()
        }
        .build();
        let out = detect(&seq, &[]);
        assert_eq!(out.len(), 1);
        assert_eq!(out[0].method, Method::SkewOnly);
        assert_eq!(out[0].dnaa_boxes, 0);
    }

    #[test]
    fn flat_skew_declined() {
        // No strand bias at all → a random walk → no trustworthy minimum.
        let mut rng = Rng::new(7);
        let seq: Vec<u8> = (0..200_000).map(|_| rng.base(50)).collect();
        assert!(
            detect(&seq, &[]).is_empty(),
            "a strand-balanced sequence must not yield an oriC"
        );
    }

    #[test]
    fn short_contig_declined() {
        let seq = vec![b'A'; 10_000];
        assert!(detect(&seq, &[]).is_empty());
    }

    #[test]
    fn dnaa_gene_rescues_weak_skew() {
        // Flat skew that fails the gate, but a dnaA gene with a box cluster at
        // it → a dnaA-anchored call.
        let mut rng = Rng::new(3);
        let mut seq: Vec<u8> = (0..300_000).map(|_| rng.base(50)).collect();
        let gene_at = 150_000usize;
        for k in 0..5 {
            let p = gene_at + k * 15;
            seq[p..p + 9].copy_from_slice(DNAA_BOX);
        }
        let genes = [GeneHint {
            start: gene_at + 1,
            end: gene_at + 1200,
            name: "dnaA".into(),
        }];
        let out = detect(&seq, &genes);
        assert_eq!(out.len(), 1, "dnaA gene should rescue the call");
        assert_eq!(out[0].method, Method::DnaaAnchored);
    }

    #[test]
    fn weak_skew_no_gene_declined() {
        // Same flat sequence, no gene, no boxes → declined (precision).
        let mut rng = Rng::new(3);
        let seq: Vec<u8> = (0..300_000).map(|_| rng.base(50)).collect();
        assert!(detect(&seq, &[]).is_empty());
    }

    #[test]
    fn linear_prefers_dnaa_over_wrong_skew_min() {
        // A linear chromosome whose skew minimum is far from the true origin
        // (Streptomyces-like), but a dnaA gene sits at oriC. Under Linear
        // topology the gene must win.
        let (mut seq, _) = Synth {
            arm: 80_000,
            rot: 40_000,
            ..Default::default()
        }
        .build();
        let n = seq.len();
        let gene_at = n / 2;
        for k in 0..5 {
            let p = gene_at + k * 15;
            seq[p..p + 9].copy_from_slice(DNAA_BOX);
        }
        let genes = [GeneHint {
            start: gene_at + 1,
            end: gene_at + 1200,
            name: "dnaA".into(),
        }];
        let opt = Options {
            topology: Topology::Linear,
            ..Default::default()
        };
        let out = detect_with(&seq, &genes, &opt);
        assert_eq!(out.len(), 1);
        assert_eq!(out[0].method, Method::DnaaAnchored);
        let d = (out[0].start as i64 - 1 - gene_at as i64).abs();
        assert!(
            d < 3_000,
            "linear call {} should be near dnaA {gene_at}",
            out[0].start
        );
    }

    #[test]
    fn detect_is_single_per_call() {
        // detect() reports at most one origin; a caller loops over replicons.
        let (chrom, _) = clean();
        let plasmid = vec![b'A'; 60_000]; // >min_len but flat → no call
        assert_eq!(detect(&chrom, &[]).len(), 1);
        assert_eq!(detect(&plasmid, &[]).len(), 0);
    }

    #[test]
    fn confidence_increases_with_evidence() {
        let opt = Options::default();
        let weak = confidence(6.0, 0.55, 0, None, false, &opt);
        let strong = confidence(30.0, 1.0, 8, Some(0), false, &opt);
        assert!(strong > weak);
        assert!((0.0..=1.0).contains(&strong));
        assert!((0.0..=1.0).contains(&weak));
        // A multimodal skew penalises the score.
        assert!(confidence(30.0, 1.0, 8, Some(0), true, &opt) < strong);
    }

    #[test]
    fn dnaa_box_cluster_anchors_the_call() {
        // With a dnaA gene abutting the origin, the call is anchored on the
        // DnaA-box cluster beside it (the common bacterial case, ~median 1 bp
        // against DoriC across 615 replicons).
        let (seq, ori) = clean();
        let n = seq.len();
        let genes = [GeneHint {
            start: ori + 1,
            end: ori + 1200,
            name: "dnaA".into(),
        }];
        let out = detect(&seq, &genes);
        assert_eq!(out.len(), 1);
        assert_eq!(out[0].method, Method::DnaaAnchored);
        assert!(out[0].dnaa_boxes >= 3);
        assert!(off(out[0].start, ori, n) < 2_000);
    }

    #[test]
    fn displaced_origin_overrides_dnaa_decoy() {
        // E. coli case: oriC (a strong box cluster at the skew minimum) sits far
        // from the dnaA gene, which has only a small decoy cluster. The denser,
        // distant cluster must win.
        let (mut seq, ori) = Synth {
            boxes: 9,
            ..Default::default()
        }
        .build();
        let n = seq.len();
        // dnaA gene ~80 kb from the origin, with a 3-box decoy beside it.
        let gene = (ori + 80_000) % n;
        for k in 0..3 {
            let p = (gene + 40 + k * 15) % n;
            seq[p..p + 9].copy_from_slice(DNAA_BOX);
        }
        let genes = [GeneHint {
            start: gene + 1,
            end: gene + 1200,
            name: "dnaA".into(),
        }];
        let out = detect(&seq, &genes);
        assert_eq!(out.len(), 1);
        assert_eq!(out[0].method, Method::SkewAndBoxes);
        assert!(
            off(out[0].start, ori, n) < 2_000,
            "should report the displaced origin at {ori}, got {}",
            out[0].start
        );
        // The displacement is surfaced as an informational note.
        assert!(out[0].notes.contains(&Flag::Displaced));
        assert!(!Flag::Displaced.is_warning());
    }

    #[test]
    fn high_gc_two_mismatch_boxes_recovered() {
        // High-GC Actinobacteria (Mycobacterium/Mycobacteroides): the oriC DnaA
        // boxes have drifted to 2 mismatches from the consensus, invisible at the
        // strict tolerance. With a dnaA gene at the origin, the GC-gated
        // 2-mismatch fallback recovers them.
        let (mut seq, ori) = Synth {
            gc_bg: 75,
            boxes: 0,
            ..Default::default()
        }
        .build();
        let variant = b"TTAGCCATA"; // 2 mismatches from TTATCCACA (pos 4, 8)
        assert!(!matches_box(variant, DNAA_BOX, 1)); // invisible at 1mm
        for k in 0..5 {
            let p = ori + k * 15;
            seq[p..p + 9].copy_from_slice(variant);
        }
        let genes = [GeneHint {
            start: ori + 1,
            end: ori + 1200,
            name: "dnaA".into(),
        }];
        let out = detect(&seq, &genes);
        assert_eq!(out.len(), 1);
        assert!(
            out[0].dnaa_boxes >= 3,
            "drifted boxes should be recovered at 2 mismatches in a high-GC genome"
        );
        assert!(off(out[0].start, ori, seq.len()) < 2_000);
    }

    #[test]
    fn clean_skew_beats_dnaa_decoy_without_box_margin() {
        // The Enterobacteriaceae fix: the true origin at the skew minimum has a
        // SMALL cluster (3 boxes) that would lose to a 3-box dnaA decoy under the
        // box-count margin. Because the skew is clean and strong, it must win
        // anyway. Regression guard for oriC displaced ~42 kb from dnaA.
        let (mut seq, ori) = Synth {
            boxes: 3,
            ..Default::default()
        }
        .build();
        let n = seq.len();
        let gene = (ori + 42_000) % n;
        for k in 0..3 {
            let p = (gene + 40 + k * 15) % n;
            seq[p..p + 9].copy_from_slice(DNAA_BOX);
        }
        let genes = [GeneHint {
            start: gene + 1,
            end: gene + 1200,
            name: "dnaA".into(),
        }];
        let out = detect(&seq, &genes);
        assert_eq!(out.len(), 1);
        assert_eq!(out[0].method, Method::SkewAndBoxes, "clean skew should win");
        assert!(
            off(out[0].start, ori, n) < 2_000,
            "displaced origin {ori} not recovered, got {}",
            out[0].start
        );
    }

    #[test]
    fn moderate_displacement_recovered_with_due() {
        // Vibrio/Aeromonas case: oriC sits ~10 kb from dnaA — the ambiguous
        // 7-20 kb band, short of DNAA_ORIC_FAR — so the displaced box cluster is
        // taken as the origin only because it carries an AT-rich DUE. GC-rich
        // genome, so the AT-headroom guard leaves the DUE test enabled.
        let (mut seq, ori) = Synth {
            boxes: 5,
            gc_bg: 50,
            ..Default::default()
        }
        .build();
        let n = seq.len();
        // AT-rich DUE immediately upstream of the origin boxes.
        for k in 0..400 {
            seq[(ori + n - 400 + k) % n] = if k & 1 == 0 { b'A' } else { b'T' };
        }
        // dnaA gene 10 kb away with a 3-box decoy beside it.
        let gene = (ori + 10_000) % n;
        for k in 0..3 {
            let p = (gene + 40 + k * 15) % n;
            seq[p..p + 9].copy_from_slice(DNAA_BOX);
        }
        let genes = [GeneHint {
            start: gene + 1,
            end: gene + 1200,
            name: "dnaA".into(),
        }];
        let out = detect(&seq, &genes);
        assert_eq!(out.len(), 1);
        assert_eq!(
            out[0].method,
            Method::SkewAndBoxes,
            "the DUE-bearing displaced cluster should override the dnaA decoy"
        );
        assert!(
            off(out[0].start, ori, n) < 2_000,
            "moderately-displaced origin {ori} not recovered, got {}",
            out[0].start
        );
    }

    #[test]
    fn at_rich_genome_suppresses_moderate_due_override() {
        // Campylobacter protection: in an AT-saturated genome (GC ~24 %) a DUE
        // window clears DUE_MIN_EXCESS almost anywhere, so the moderate-band DUE
        // test loses all discriminating power. The AT-headroom guard
        // (DUE_MAX_GENOME_AT) disables it, so a chance skew cluster 10 kb from
        // dnaA no longer hijacks the (correct) dnaA-adjacent call. Identical
        // geometry to the test above; only the genome AT fraction differs.
        let (mut seq, ori) = Synth {
            boxes: 5,
            skew_pct: 20,
            gc_bg: 5,
            ..Default::default()
        }
        .build();
        let n = seq.len();
        for k in 0..400 {
            seq[(ori + n - 400 + k) % n] = if k & 1 == 0 { b'A' } else { b'T' };
        }
        let gene = (ori + 10_000) % n;
        for k in 0..3 {
            let p = (gene + 40 + k * 15) % n;
            seq[p..p + 9].copy_from_slice(DNAA_BOX);
        }
        let genes = [GeneHint {
            start: gene + 1,
            end: gene + 1200,
            name: "dnaA".into(),
        }];
        let out = detect(&seq, &genes);
        assert_eq!(out.len(), 1);
        assert_eq!(
            out[0].method,
            Method::DnaaAnchored,
            "the moderate-band DUE override must be suppressed in AT-saturated genomes"
        );
    }

    #[test]
    fn clean_call_is_pass_with_no_warnings() {
        // A clean genome with dnaA at the origin is a high-confidence call: Pass,
        // and none of its notes are warnings.
        let (seq, ori) = clean();
        let genes = [GeneHint {
            start: ori + 1,
            end: ori + 1200,
            name: "dnaA".into(),
        }];
        let out = detect(&seq, &genes);
        assert_eq!(out[0].conf, Confidence::Pass);
        assert!(out[0].score >= CONF_PASS);
        assert!(
            out[0].notes.iter().all(|f| !f.is_warning()),
            "clean call should carry no warning notes, got {:?}",
            out[0].notes
        );
    }

    #[test]
    fn no_hint_is_flagged() {
        // Called with no dnaA gene, the sequence-only mode is surfaced.
        let (seq, _) = clean();
        let out = detect(&seq, &[]);
        assert!(out[0].notes.contains(&Flag::NoHint));
        assert!(!Flag::NoHint.is_warning());
    }

    #[test]
    fn missing_box_cluster_is_flagged() {
        // Strong skew but no DnaA-box cluster: a SkewOnly call, flagged no-boxes
        // (lower precision) and, with no gene, no-hint.
        let (seq, _) = Synth {
            boxes: 0,
            ..Default::default()
        }
        .build();
        let out = detect(&seq, &[]);
        assert_eq!(out.len(), 1);
        assert_eq!(out[0].method, Method::SkewOnly);
        assert!(out[0].notes.contains(&Flag::NoBoxes));
        assert!(Flag::NoBoxes.is_warning());
    }

    #[test]
    fn confidence_tiers_track_score() {
        // The tier boundaries are exactly CONF_PASS / CONF_WEAK.
        let tier = |s: f32| {
            if s >= CONF_PASS {
                Confidence::Pass
            } else if s < CONF_WEAK {
                Confidence::Weak
            } else {
                Confidence::Review
            }
        };
        assert_eq!(tier(0.95), Confidence::Pass);
        assert_eq!(tier(CONF_PASS), Confidence::Pass);
        assert_eq!(tier(0.45), Confidence::Review);
        assert_eq!(tier(CONF_WEAK), Confidence::Review);
        assert_eq!(tier(0.10), Confidence::Weak);
    }
}