oricle 0.2.0

//! Real-genome oriC regression harness.
//!
//! This test is **opt-in**: complete bacterial genomes are megabytes each and
//! do not belong in the crate package, so the genomes live outside the repo and
//! the test is a no-op unless `ORICLE_GENOME_DIR` points at a directory holding
//! `<accession>.fa` (or `.fasta` / `.fna`) files. The expected origins are in
//! `tests/data/regression_truth.tsv`, re-derived on current RefSeq accessions
//! from DoriC (see that file's header).
//!
//! ```text
//! ORICLE_GENOME_DIR=/data/genomes cargo test --test regression -- --nocapture
//! ```
//!
//! An optional per-genome `<accession>.genes.tsv` (columns: seqname, start, end,
//! name) supplies a dnaA hint, required for linear chromosomes such as
//! Streptomyces.

use std::fs;
use std::path::{Path, PathBuf};

use oricle::{detect_with, GeneHint, Options, Topology};

struct Truth {
    accession: String,
    start: usize,
    end: usize,
    topology: Topology,
    tol: i64,
    note: String,
}

fn parse_truth() -> Vec<Truth> {
    let path = Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/data/regression_truth.tsv");
    let text = fs::read_to_string(&path).expect("read regression_truth.tsv");
    text.lines()
        .filter(|l| !l.trim_start().starts_with('#') && !l.trim().is_empty())
        .map(|l| {
            let f: Vec<&str> = l.split('\t').collect();
            Truth {
                accession: f[0].to_string(),
                start: f[1].parse().unwrap(),
                end: f[2].parse().unwrap(),
                topology: if f[3] == "linear" {
                    Topology::Linear
                } else {
                    Topology::Circular
                },
                tol: f[4].parse().unwrap(),
                note: f.get(5).unwrap_or(&"").to_string(),
            }
        })
        .collect()
}

fn find_fasta(dir: &Path, acc: &str) -> Option<PathBuf> {
    for ext in ["fa", "fasta", "fna", "fa.gz", "fasta.gz", "fna.gz"] {
        let p = dir.join(format!("{acc}.{ext}"));
        if p.exists() {
            return Some(p);
        }
    }
    None
}

/// Minimal FASTA reader: concatenate sequence lines of the first record.
fn read_fasta(path: &Path) -> Vec<u8> {
    let text = fs::read_to_string(path).expect("read fasta (plain, uncompressed expected)");
    let mut seq = Vec::new();
    for line in text.lines() {
        if line.starts_with('>') {
            if !seq.is_empty() {
                break;
            }
        } else {
            seq.extend(line.trim().bytes());
        }
    }
    seq
}

fn read_genes(dir: &Path, acc: &str) -> Vec<GeneHint> {
    let p = dir.join(format!("{acc}.genes.tsv"));
    let Ok(text) = fs::read_to_string(&p) else {
        return Vec::new();
    };
    text.lines()
        .filter_map(|l| {
            let f: Vec<&str> = l.split('\t').collect();
            if f.len() < 4 {
                return None;
            }
            Some(GeneHint {
                start: f[1].trim().parse().ok()?,
                end: f[2].trim().parse().ok()?,
                name: f[3].trim().to_string(),
            })
        })
        .collect()
}

/// Circular distance from a predicted 1-based inclusive interval to the truth
/// interval (0 if they overlap), handling origin-wrapping intervals.
fn interval_dist(p: (usize, usize), t: (usize, usize), n: i64) -> i64 {
    let arcs = |s: usize, e: usize| -> Vec<(i64, i64)> {
        if s <= e {
            vec![(s as i64, e as i64)]
        } else {
            vec![(s as i64, n), (1, e as i64)]
        }
    };
    let pa = arcs(p.0, p.1);
    let ta = arcs(t.0, t.1);
    for &(ps, pe) in &pa {
        for &(ts, te) in &ta {
            if ps <= te && ts <= pe {
                return 0;
            }
        }
    }
    let pt_arc = |pt: i64, s: i64, e: i64| -> i64 {
        let d = (pt - s).abs().min((pt - e).abs());
        d.min(n - d)
    };
    let mut best = i64::MAX;
    for &(ps, pe) in &pa {
        for pt in [t.0 as i64, t.1 as i64] {
            best = best.min(pt_arc(pt, ps, pe));
        }
    }
    for &(ts, te) in &ta {
        for pt in [p.0 as i64, p.1 as i64] {
            best = best.min(pt_arc(pt, ts, te));
        }
    }
    best
}

#[test]
fn real_genome_regression() {
    let Ok(dir) = std::env::var("ORICLE_GENOME_DIR") else {
        eprintln!(
            "skipping real_genome_regression: set ORICLE_GENOME_DIR to a directory \
             of <accession>.fa genomes to enable it"
        );
        return;
    };
    let dir = PathBuf::from(dir);
    let truths = parse_truth();
    let mut ran = 0;
    let mut failures = Vec::new();

    for t in &truths {
        let Some(fa) = find_fasta(&dir, &t.accession) else {
            eprintln!("  {} — genome not found, skipped", t.accession);
            continue;
        };
        if fa.to_string_lossy().ends_with(".gz") {
            eprintln!(
                "  {} — gzip not supported by the harness reader, skipped",
                t.accession
            );
            continue;
        }
        ran += 1;
        let seq = read_fasta(&fa);
        let genes = read_genes(&dir, &t.accession);
        let opt = Options {
            topology: t.topology,
            ..Default::default()
        };
        let calls = detect_with(&seq, &genes, &opt);
        let n = seq.len() as i64;
        match calls.first() {
            None => failures.push(format!("{} ({}) — NO CALL", t.accession, t.note)),
            Some(o) => {
                let d = interval_dist((o.start, o.end), (t.start, t.end), n);
                let ok = d <= t.tol;
                eprintln!(
                    "  {} — pred {}..{} vs {}..{}  dist={} bp  tol={}  {}",
                    t.accession,
                    o.start,
                    o.end,
                    t.start,
                    t.end,
                    d,
                    t.tol,
                    if ok { "PASS" } else { "FAIL" }
                );
                if !ok {
                    failures.push(format!(
                        "{} ({}) — {} bp off (tol {})",
                        t.accession, t.note, d, t.tol
                    ));
                }
            }
        }
    }

    eprintln!("real_genome_regression: {ran} genomes exercised");
    assert!(
        ran > 0,
        "no genomes found in ORICLE_GENOME_DIR — check the path"
    );
    assert!(
        failures.is_empty(),
        "regressions:\n  {}",
        failures.join("\n  ")
    );
}