oricle 0.2.0

//! `oricle` CLI: predict bacterial oriC (replication origin) from a FASTA.
//!
//! ```text
//! oricle [OPTIONS] <input.fasta>
//!
//!   --genes <tsv>   Gene hints: seqname, start, end, name (1-based, tab-separated,
//!                   optional header). Only dnaA is used, to refine or rescue a call.
//!   --linear        Treat replicons as linear (Streptomyces, Borreliella).
//!   --all           Also print a row for replicons with no call.
//!   -h, --help
//! ```
//!
//! Prints a TSV: `seqname  start  end  score  conf  skew_min  dnaa_boxes  method
//! notes`.

use std::collections::HashMap;
use std::env;
use std::fs;
use std::process;

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

const USAGE: &str = "\
Usage: oricle [OPTIONS] <input.fasta>

  --genes <tsv>  Gene hints (seqname, start, end, name; 1-based, tab-separated).
                 Only dnaA is used, to refine or rescue a call.
  --linear       Treat replicons as linear (default: circular).
  --all          Print a row for replicons with no call too.
  -h, --help     Show this help.

Output TSV columns: seqname, start, end, score, conf, skew_min, dnaa_boxes,
method, notes

  conf   pass | review | weak  (action-oriented confidence, from score)
  notes  ;-separated warnings/context, or . when clean
         (multimodal, weak-skew, no-boxes, dnaa-disagree, degenerate-box =
          warnings; displaced, no-hint = informational)
";

fn main() {
    let argv: Vec<String> = env::args().skip(1).collect();
    let mut path: Option<String> = None;
    let mut genes_path: Option<String> = None;
    let mut opt = Options::default();
    let mut show_all = false;

    let mut i = 0;
    while i < argv.len() {
        match argv[i].as_str() {
            "-h" | "--help" => {
                print!("{USAGE}");
                return;
            }
            "--linear" => opt.topology = Topology::Linear,
            "--all" => show_all = true,
            "--min-sigma" => {
                i += 1;
                opt.min_skew_sigma = argv
                    .get(i)
                    .and_then(|v| v.parse().ok())
                    .unwrap_or_else(|| fail("--min-sigma needs a number"));
            }
            "--min-antipodal" => {
                i += 1;
                opt.min_antipodal = argv
                    .get(i)
                    .and_then(|v| v.parse().ok())
                    .unwrap_or_else(|| fail("--min-antipodal needs a number"));
            }
            "--genes" => {
                i += 1;
                match argv.get(i) {
                    Some(v) => genes_path = Some(v.clone()),
                    None => fail("--genes needs a file argument"),
                }
            }
            a if a.starts_with('-') => fail(&format!("unknown option: {a}")),
            a => {
                if path.is_some() {
                    fail("expected exactly one input FASTA");
                }
                path = Some(a.to_string());
            }
        }
        i += 1;
    }

    let Some(path) = path else {
        eprint!("{USAGE}");
        process::exit(2);
    };

    let genes = match genes_path {
        Some(p) => load_genes(&p),
        None => HashMap::new(),
    };

    let mut reader = match parse_fastx_file(&path) {
        Ok(r) => r,
        Err(e) => fail(&format!("failed to open {path}: {e}")),
    };

    println!("seqname\tstart\tend\tscore\tconf\tskew_min\tdnaa_boxes\tmethod\tnotes");

    while let Some(rec) = reader.next() {
        let rec = match rec {
            Ok(r) => r,
            Err(e) => fail(&format!("error reading record: {e}")),
        };
        let id = rec.id();
        let name = id.split(|&b| b == b' ' || b == b'\t').next().unwrap_or(id);
        let name = String::from_utf8_lossy(name).into_owned();

        let hints: &[GeneHint] = genes.get(&name).map(Vec::as_slice).unwrap_or(&[]);
        let calls = detect_with(&rec.seq(), hints, &opt);

        if calls.is_empty() && show_all {
            println!("{name}\t.\t.\t.\t.\t.\t.\tnone\t.");
        }
        for o in calls {
            let notes = if o.notes.is_empty() {
                ".".to_string()
            } else {
                o.notes
                    .iter()
                    .map(|f| f.code())
                    .collect::<Vec<_>>()
                    .join(";")
            };
            println!(
                "{}\t{}\t{}\t{:.3}\t{}\t{}\t{}\t{:?}\t{}",
                name,
                o.start,
                o.end,
                o.score,
                o.conf.as_str(),
                o.skew_min,
                o.dnaa_boxes,
                o.method,
                notes
            );
        }
    }
}

/// Parse a gene-hint TSV into per-sequence hints. Blank lines, short lines, and
/// lines whose coordinates do not parse (e.g. a header row) are skipped.
fn load_genes(path: &str) -> HashMap<String, Vec<GeneHint>> {
    let text = match fs::read_to_string(path) {
        Ok(t) => t,
        Err(e) => fail(&format!("failed to read {path}: {e}")),
    };
    let mut out: HashMap<String, Vec<GeneHint>> = HashMap::new();
    for line in text.lines() {
        let f: Vec<&str> = line.split('\t').collect();
        if f.len() < 4 {
            continue;
        }
        let (Ok(start), Ok(end)) = (f[1].trim().parse(), f[2].trim().parse()) else {
            continue;
        };
        out.entry(f[0].trim().to_string())
            .or_default()
            .push(GeneHint {
                start,
                end,
                name: f[3].trim().to_string(),
            });
    }
    out
}

fn fail(msg: &str) -> ! {
    eprintln!("oricle: {msg}");
    process::exit(1);
}