rust-ise 0.1.0

// rust-ise — Rust-native ISEScan-equivalent IS-element scanner.
// Pipeline: rustygal META gene-finding (LINKED LIBRARY, == pyrodigal) + strand-aware edge recovery
//   -> mmseqs2 profile search (system `mmseqs` on PATH; same union DB/params as the Python reference pv2)
//   -> 4-tier family call -> native affine Smith-Waterman-Gotoh TIR -> c/p -> dedup -> TSV/GFF/.sum.
// Coordinates 1-based inclusive. mmseqs2 is the one external engine (call system binary); the IS DB is
// external data resolved via --db / $RUSTISE_DB (not bundled in the crate).
use rayon::prelude::*;
use std::cell::RefCell;
use std::collections::HashMap;
use std::fs;
use std::io::{BufRead, BufReader, Write, BufWriter};
use std::path::{Path, PathBuf};
use std::process::Command;
use std::time::Instant;

// Compile-time perfect-hash map `REP_THRESHOLDS: phf::Map<&str, f32>` (rep -> min pident to accept a
// hit to that profile), GTDB-calibrated, generated by build.rs from data/rep_thresholds.tsv.
include!(concat!(env!("OUT_DIR"), "/rep_thresholds.rs"));

// ---- family-call thresholds (tuned, == v2) ----
const SIG_EVALUE: f64 = 1e-10;
const LIKE_EVALUE: f64 = 1e-30;
const LIKE_QCOV: f64 = 0.7;
const LIKE_SPAN_OVERLAP: f64 = 0.5;
const NOVEL_PIDENT: f64 = 30.0;
// ---- TIR / flank params (== Python reference constants) ----
const MAX_DIST: i64 = 500;
const MIN_DIST_ABS: i64 = 150; // abs(minDist4ter2orf = -150)
const SW_MATCH: i32 = 2;
const SW_MISMATCH: i32 = -6;
// affine gap (Gotoh, parasail convention: a length-g gap costs OPEN + EXTEND*(g-1)).
// Python-reference filters4ssw4trial = (match2, mismatch6, gap_open2, gap_extend2).
const SW_OPEN: i32 = 2;
const SW_EXTEND: i32 = 2;

macro_rules! log { ($t0:expr, $($a:tt)*) => {{
    let _ = &$t0; print!("[{:7.1}s] ", $t0.elapsed().as_secs_f64()); println!($($a)*);
}}}

// ---------------------------------------------------------- sequence utils
fn revcomp(s: &[u8]) -> Vec<u8> {
    s.iter().rev().map(|&c| match c {
        b'A' => b'T', b'T' => b'A', b'G' => b'C', b'C' => b'G',
        b'a' => b't', b't' => b'a', b'g' => b'c', b'c' => b'g',
        b'N' => b'N', b'n' => b'n', _ => b'N',
    }).collect()
}

fn norm_fam(f: &str) -> String {
    f.replace('/', "").replace('_', "").to_uppercase()
}

// standard genetic code (== v2 _CODON), index by 3 uppercase bases
fn translate(s: &[u8]) -> Vec<u8> {
    const TBL: &[u8] = b"FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG";
    let code = |c: u8| -> i32 { match c { b'T'=>0, b'C'=>1, b'A'=>2, b'G'=>3, _=>-1 } };
    let mut out = Vec::with_capacity(s.len()/3);
    let mut i = 0;
    while i + 3 <= s.len() {
        let (a,b,c) = (code(s[i]), code(s[i+1]), code(s[i+2]));
        if a<0 || b<0 || c<0 { out.push(b'X'); }
        else { out.push(TBL[(a*16+b*4+c) as usize]); }
        i += 3;
    }
    out
}

fn read_fasta(path: &str) -> Result<Vec<(String, Vec<u8>)>, Box<dyn std::error::Error>> {
    let f = BufReader::new(fs::File::open(path)
        .map_err(|e| format!("cannot open input FASTA '{}': {}", path, e))?);
    let mut out: Vec<(String, Vec<u8>)> = Vec::new();
    let mut name: Option<String> = None;
    let mut buf: Vec<u8> = Vec::new();
    for line in f.lines() {
        let line = line.map_err(|e| format!("error reading FASTA '{}': {}", path, e))?;
        if let Some(rest) = line.strip_prefix('>') {
            if let Some(n) = name.take() { out.push((n, std::mem::take(&mut buf))); }
            name = Some(rest.split_whitespace().next().unwrap_or("").to_string());
        } else {
            buf.extend(line.trim().bytes().map(|c| c.to_ascii_uppercase()));
        }
    }
    if let Some(n) = name.take() { out.push((n, buf)); }
    Ok(out)
}

// ---------------------------------------------------------- 1. ORFs (rustygal + edge)
struct Orf { begin: i64, end: i64, strand: char, prot: Vec<u8>, edge: bool }

// Parse rustygal -a FASTA bytes: header ">{contig}_{geneidx} # b # e # strand # ID=..." , seq = protein.
// (rv2 gets these bytes from the rustygal LIBRARY `run_meta`, not a file — same format as the binary's -a.)
// Returns contig -> Vec<Orf> (pyrodigal-equivalent genes only; edge added later).
fn parse_rustygal_faa(data: &[u8]) -> Result<HashMap<String, Vec<Orf>>, Box<dyn std::error::Error>> {
    let mut map: HashMap<String, Vec<Orf>> = HashMap::new();
    let mut cur: Option<(String, i64, i64, char)> = None;
    let mut seq: Vec<u8> = Vec::new();
    let flush = |map: &mut HashMap<String, Vec<Orf>>,
                 cur: &Option<(String,i64,i64,char)>, seq: &mut Vec<u8>| {
        if let Some((ctg, b, e, st)) = cur {
            // strip trailing '*' to match v2 g.translate().rstrip("*")
            while seq.last() == Some(&b'*') { seq.pop(); }
            map.entry(ctg.clone()).or_default().push(Orf {
                begin: *b, end: *e, strand: *st, prot: std::mem::take(seq), edge: false });
        } else { seq.clear(); }
    };
    for raw in data.split(|&c| c == b'\n') {
        let line = String::from_utf8_lossy(raw);
        let line = line.trim_end_matches('\r');
        if let Some(rest) = line.strip_prefix('>') {
            flush(&mut map, &cur, &mut seq);
            // rest: "{contig}_{geneidx} # b # e # strand # ID=..."
            let first = rest.split_whitespace().next()
                .ok_or_else(|| format!("malformed rustygal header (empty): '{}'", rest))?;
            let parts: Vec<&str> = rest.split('#').collect();
            if parts.len() < 4 {
                return Err(format!("malformed rustygal header (missing '#' fields): '{}'", rest).into());
            }
            let b: i64 = parts[1].trim().parse()
                .map_err(|e| format!("bad begin coord in rustygal header '{}': {}", rest, e))?;
            let e: i64 = parts[2].trim().parse()
                .map_err(|e| format!("bad end coord in rustygal header '{}': {}", rest, e))?;
            let strand = if parts[3].trim() == "1" { '+' } else { '-' };
            // contig = first token minus the trailing "_{geneidx}"
            let contig = match first.rfind('_') { Some(i) => &first[..i], None => first };
            cur = Some((contig.to_string(), b, e, strand));
        } else {
            seq.extend(line.trim().bytes());
        }
    }
    flush(&mut map, &cur, &mut seq);
    Ok(map)
}

// strand-aware edge recovery for one contig (== v2 _orf_one edge part).
fn edge_orfs(seq: &[u8], pyr: &[Orf], minaa: usize, edge: i64) -> Vec<Orf> {
    let l = seq.len() as i64;
    let mut out = Vec::new();
    for &strand in &['+', '-'] {
        let s = if strand == '+' { seq.to_vec() } else { revcomp(seq) };
        for frame in 0..3usize {
            let prot = translate(&s[frame.min(s.len())..]);
            // maximal non-stop stretches [a0,a1) (a1 exclusive)
            let mut start = 0usize;
            let mut stretches: Vec<(usize, usize)> = Vec::new();
            for (k, &aa) in prot.iter().enumerate() {
                if aa == b'*' { if k > start { stretches.push((start, k)); } start = k + 1; }
            }
            if start < prot.len() { stretches.push((start, prot.len())); }
            for (a0, a1) in stretches {
                if a1 - a0 < minaa { continue; }
                let ns = frame as i64 + a0 as i64 * 3;
                let ne = frame as i64 + a1 as i64 * 3;
                if !(ns <= edge || ne >= l - edge) { continue; }
                let (cb, ce) = if strand == '+' { (ns + 1, ne) } else { (l - ne + 1, l - ns) };
                let suppressed = pyr.iter().any(|p| p.strand == strand
                    && (ce.min(p.end) - cb.max(p.begin)) as f64 >= 0.5 * (ce - cb) as f64);
                if suppressed { continue; }
                out.push(Orf { begin: cb, end: ce, strand, prot: prot[a0..a1].to_vec(), edge: true });
            }
        }
    }
    out
}

// ---------------------------------------------------------- 2. mmseqs search
#[derive(Clone)]
struct Hit { fam: String, evalue: f64, qcov: f64, pident: f64, tstart: i64, tend: i64 }

fn mmseqs_search(proteome: &Path, db_dir: &Path, tmp: &Path, nthread: usize)
    -> Result<HashMap<String, Vec<Hit>>, Box<dyn std::error::Error>> {
    let profile_db = db_dir.join("mmdb_union").join("profileDb");
    let manifest = db_dir.join("manifest_union.tsv");
    let run = |args: &[&str]| -> Result<(), Box<dyn std::error::Error>> {
        let st = Command::new("mmseqs").args(args)
            .stdout(std::process::Stdio::null()).stderr(std::process::Stdio::null())
            .status().map_err(|e| format!("failed to run mmseqs (is it on PATH?): {}", e))?;
        if !st.success() {
            return Err(format!("mmseqs {:?} failed (exit {:?})",
                args.first().copied().unwrap_or(""), st.code()).into());
        }
        Ok(())
    };
    let qdb = tmp.join("qdb");
    let proteome_s = proteome.to_str().expect("proteome path is valid UTF-8");
    let qdb_s = qdb.to_str().expect("qdb tmp path is valid UTF-8");
    run(&["createdb", proteome_s, qdb_s, "-v", "0"])?;
    // cid -> family and cid -> rep (rep is the stable key for the embedded per-profile threshold)
    let mut cid2fam: HashMap<String, String> = HashMap::new();
    let mut cid2rep: HashMap<String, String> = HashMap::new();
    let mfile = fs::File::open(&manifest)
        .map_err(|e| format!("cannot open DB manifest '{}': {}", manifest.display(), e))?;
    for line in BufReader::new(mfile).lines() {
        let line = line.map_err(|e| format!("error reading DB manifest '{}': {}", manifest.display(), e))?;
        let p: Vec<&str> = line.split('\t').collect();
        if p.len() >= 2 && p[0] != "cid" { cid2fam.insert(p[0].to_string(), p[1].to_string()); }
        if p.len() >= 4 && p[0] != "cid" { cid2rep.insert(p[0].to_string(), p[3].to_string()); }
    }
    let res = tmp.join("pres");
    let aln = tmp.join("pres.tsv");
    let tp = tmp.join("tp");
    let nt = nthread.to_string();
    // TODO (under development): MMseqs2-GPU acceleration. A real GPU search needs a
    // GPU-padded profile DB (`mmseqs makepaddedseqdb`) plus `--gpu 1 --gpu-server 1`
    // on a CUDA GPU. The earlier `--gpu` CLI flag was unvalidated (it only appended
    // `--gpu 1` without the padded DB / server) and has been removed until this path
    // is implemented and verified.
    let profile_db_s = profile_db.to_str().expect("profile DB path is valid UTF-8");
    let res_s = res.to_str().expect("results tmp path is valid UTF-8");
    let aln_s = aln.to_str().expect("alignment tmp path is valid UTF-8");
    let tp_s = tp.to_str().expect("tp tmp path is valid UTF-8");
    let search_args: Vec<&str> = vec!["search", profile_db_s,
        qdb_s, res_s, tp_s,
        "-s", "7.5", "-e", "10", "--threads", &nt, "-v", "0"];
    run(&search_args)?;
    run(&["convertalis", profile_db_s, qdb_s,
        res_s, aln_s,
        "--format-output", "query,target,evalue,qcov,pident,tstart,tend", "-v", "0"])?;
    let mut hits: HashMap<String, Vec<Hit>> = HashMap::new();
    let afile = fs::File::open(&aln)
        .map_err(|e| format!("cannot open mmseqs alignment '{}': {}", aln.display(), e))?;
    for line in BufReader::new(afile).lines() {
        let line = line.map_err(|e| format!("error reading mmseqs alignment '{}': {}", aln.display(), e))?;
        let p: Vec<&str> = line.split('\t').collect();
        if p.len() < 7 { continue; }
        let pident: f64 = p[4].parse().unwrap_or(0.0);
        // per-profile pident gate (GTDB-calibrated, rep-keyed): drop hits to a contaminant /
        // promiscuous profile when their identity is below the profile's threshold. Absent rep =>
        // threshold 0 => accept. This runs before family_call so a gated hit cannot win the ORF.
        if let Some(rep) = cid2rep.get(p[0]) {
            let thr = REP_THRESHOLDS.get(rep.as_str()).copied().unwrap_or(0.0) as f64;
            if pident < thr { continue; }
        }
        if let Some(fam) = cid2fam.get(p[0]) {
            hits.entry(p[1].to_string()).or_default().push(Hit {
                fam: fam.clone(),
                evalue: p[2].parse().unwrap_or(1.0),
                qcov: p[3].parse().unwrap_or(0.0),
                pident,
                tstart: p[5].parse().unwrap_or(0),
                tend: p[6].parse().unwrap_or(0),
            });
        }
    }
    Ok(hits)
}

// ---------------------------------------------------------- 3. family call (4-tier)
// Total order on hits => representative selection is deterministic regardless of the
// (run-to-run varying) mmseqs hit order. evalue asc, then pident/qcov desc, then
// tstart/tend asc, then family — fully canonical so equal-E ties never flip the output.
fn cmp_hit(a: &Hit, b: &Hit) -> std::cmp::Ordering {
    a.evalue.partial_cmp(&b.evalue).unwrap()
        .then(b.pident.partial_cmp(&a.pident).unwrap())
        .then(b.qcov.partial_cmp(&a.qcov).unwrap())
        .then(a.tstart.cmp(&b.tstart))
        .then(a.tend.cmp(&b.tend))
        .then_with(|| a.fam.cmp(&b.fam))
}

// returns (label, tier, best_hit) ; tier in {plain,like,chimera,novel}
fn family_call(orfhits: &[Hit]) -> Option<(String, String, Hit)> {
    if orfhits.is_empty() { return None; }
    let mut perfam: HashMap<String, Hit> = HashMap::new();
    for h in orfhits {
        let nf = norm_fam(&h.fam);
        match perfam.get(&nf) {
            Some(cur) if cmp_hit(cur, h) != std::cmp::Ordering::Greater => {}
            _ => { perfam.insert(nf, h.clone()); }
        }
    }
    let mut ranked: Vec<Hit> = perfam.into_values().collect();
    ranked.sort_by(cmp_hit);
    let top = ranked[0].clone();
    if top.evalue > SIG_EVALUE { return None; }
    if top.pident < NOVEL_PIDENT { return Some(("novel".to_string(), "novel".to_string(), top)); }
    // competitor: best DIFFERENT family, strong + good coverage
    let comp = ranked[1..].iter().find(|x|
        norm_fam(&x.fam) != norm_fam(&top.fam) && x.evalue <= LIKE_EVALUE && x.qcov >= LIKE_QCOV);
    match comp {
        None => Some((top.fam.clone(), "plain".to_string(), top)),
        Some(c) => {
            let ov = top.tend.min(c.tend) - top.tstart.max(c.tstart);
            let span = ((top.tend - top.tstart).min(c.tend - c.tstart)).max(1);
            if ov > 0 && ov as f64 >= LIKE_SPAN_OVERLAP * span as f64 {
                Some((format!("{}-like", top.fam), "like".to_string(), top))
            } else {
                Some((format!("{}-chimera", top.fam), "chimera".to_string(), top))
            }
        }
    }
}

// ---------------------------------------------------------- 4. TIR (native SW) + c/p
struct Tir { irid: i64, irlen: i64, start1: i64, end1: i64, start2: i64, end2: i64,
             seq1: Vec<u8>, seq2: Vec<u8> }

#[derive(Default)]
struct SwBuf {
    ph: Vec<i32>, ch: Vec<i32>,   // H rows (prev / cur)
    pf: Vec<i32>, cf: Vec<i32>,   // F rows (prev / cur) — vertical-gap matrix
    hp: Vec<u8>, ep: Vec<u8>, fp: Vec<u8>,  // traceback pointer matrices
}
thread_local! {
    // Reused across every TIR alignment on this worker thread → no per-call heap allocation.
    // Buffers only grow to the largest flank seen so far.
    static SW_SCRATCH: RefCell<SwBuf> = RefCell::new(SwBuf::default());
}

// Local Smith-Waterman-Gotoh with AFFINE gaps (parasail convention: a length-g gap costs
// SW_OPEN + SW_EXTEND*(g-1)). Two gap matrices: E = horizontal (gap in query / consume ref),
// F = vertical (gap in ref / consume query). When SW_OPEN==SW_EXTEND this is numerically
// identical to the linear case. Returns (score, beg_q, end_q, beg_r, end_r) 0-based inclusive.
fn sw_local(q: &[u8], r: &[u8]) -> Option<(i32, usize, usize, usize, usize)> {
    let (n, m) = (q.len(), r.len());
    if n == 0 || m == 0 { return None; }
    const NEG: i32 = i32::MIN / 4;   // -inf for gap-matrix init (no overflow on -EXTEND)
    SW_SCRATCH.with(|cell| {
        let mut s = cell.borrow_mut();
        let b = &mut *s;             // deref RefMut once so field borrows are provably disjoint
        let stride = m + 1;
        if b.ph.len() < stride {
            b.ph.resize(stride, 0); b.ch.resize(stride, 0);
            b.pf.resize(stride, 0); b.cf.resize(stride, 0);
        }
        let need = (n + 1) * stride;
        if b.hp.len() < need { b.hp.resize(need, 0); b.ep.resize(need, 0); b.fp.resize(need, 0); }
        let (ph, ch, pf, cf) = (&mut b.ph, &mut b.ch, &mut b.pf, &mut b.cf);
        let (hp, ep, fp) = (&mut b.hp, &mut b.ep, &mut b.fp);
        // row 0: H=0, F=-inf
        for x in ph[..stride].iter_mut() { *x = 0; }
        for x in pf[..stride].iter_mut() { *x = NEG; }
        let (mut best, mut bi, mut bj) = (0i32, 0usize, 0usize);
        for i in 1..=n {
            ch[0] = 0; cf[0] = NEG;
            let mut e_prev = NEG;        // E[i][0] = -inf
            let qi = q[i-1];
            let row = i * stride;
            for j in 1..=m {
                let sc = if qi == r[j-1] && matches!(qi, b'A'|b'C'|b'G'|b'T')
                        { SW_MATCH } else { SW_MISMATCH };
                // E[i][j] (horizontal): open from H[i][j-1] or extend from E[i][j-1]
                let e_open = ch[j-1] - SW_OPEN;
                let e_ext  = e_prev - SW_EXTEND;
                let (e_val, e_code) = if e_ext > e_open { (e_ext, 1u8) } else { (e_open, 0u8) };
                // F[i][j] (vertical): open from H[i-1][j] or extend from F[i-1][j]
                let f_open = ph[j] - SW_OPEN;
                let f_ext  = pf[j] - SW_EXTEND;
                let (f_val, f_code) = if f_ext > f_open { (f_ext, 1u8) } else { (f_open, 0u8) };
                // H[i][j]; tie-break order diag > F(up) > E(left) (matches the old linear order)
                let diag = ph[j-1] + sc;
                let mut h = 0i32; let mut hc = 0u8;
                if diag > h { h = diag; hc = 1; }
                if f_val > h { h = f_val; hc = 3; }
                if e_val > h { h = e_val; hc = 2; }
                ch[j] = h; cf[j] = f_val; e_prev = e_val;
                hp[row + j] = hc; ep[row + j] = e_code; fp[row + j] = f_code;
                // tie-break: keep LAST max cell (>=) — closer to parasail's striped traceback
                if h >= best && h > 0 { best = h; bi = i; bj = j; }
            }
            std::mem::swap(ph, ch);
            std::mem::swap(pf, cf);
        }
        if best <= 0 { return None; }
        // 3-state traceback to the origin (where H came from 0). state: 0=H, 1=E, 2=F.
        let (mut i, mut j, mut state) = (bi, bj, 0u8);
        loop {
            match state {
                0 => match hp[i*stride + j] {
                    0 => break,                       // origin
                    1 => { i -= 1; j -= 1; }          // diagonal
                    2 => state = 1,                   // enter horizontal gap (E)
                    _ => state = 2,                   // enter vertical gap (F)
                },
                1 => { let ext = ep[i*stride + j] == 1; j -= 1; if !ext { state = 0; } }
                _ => { let ext = fp[i*stride + j] == 1; i -= 1; if !ext { state = 0; } }
            }
            if i == 0 || j == 0 { break; }
        }
        Some((best, i, bi - 1, j, bj - 1))
    })
}

fn find_tir(contig: &[u8], orf_b: i64, orf_e: i64, family: &str,
            tir_tbl: &HashMap<String,(i64,i64,i64,i64)>) -> Option<Tir> {
    let fam = family.split('-').next().expect("split always yields at least one element");
    let info = lookup(tir_tbl, fam);
    if let Some(t) = info { if t.3 == 0 { return None; } }
    let min_ir = info.map(|t| t.0).unwrap_or(4);
    if min_ir > 1000 { return None; }
    let l = contig.len() as i64;
    let lb = (orf_b - MAX_DIST).max(1);
    let le = (orf_b + MIN_DIST_ABS).min(l);
    let rb = (orf_e - MIN_DIST_ABS).max(1);
    let re = (orf_e + MAX_DIST).min(l);
    if le < lb || re < rb { return None; }
    let left = &contig[(lb-1) as usize..le as usize];
    let right = &contig[(rb-1) as usize..re as usize];
    if (left.len() as i64) < min_ir || (right.len() as i64) < min_ir { return None; }
    let rrc = revcomp(right);
    let (score, qs, qe, ts, te) = sw_local(left, &rrc)?;
    if score < (2 * min_ir) as i32 { return None; }
    let ir_len = qe as i64 - qs as i64 + 1;
    if ir_len < min_ir { return None; }
    let seq1 = left[qs..=qe].to_vec();
    let r_local_end = right.len() - 1 - ts;
    let r_local_beg = right.len() - 1 - te;
    let start1 = lb + qs as i64; let end1 = lb + qe as i64;
    let start2 = rb + r_local_beg as i64; let end2 = rb + r_local_end as i64;
    let seq2 = right[r_local_beg..=r_local_end].to_vec();
    let seq2rc = revcomp(&seq2);
    let ir_id = seq1.iter().zip(seq2rc.iter()).filter(|(a, b)| a == b).count() as i64;
    if ir_len == 0 || (ir_id as f64) / (ir_len as f64) < 0.5 { return None; }
    Some(Tir { irid: ir_id, irlen: ir_len, start1, end1, start2, end2, seq1, seq2 })
}

// O(1): tables are built with normalized keys, so a direct get on norm_fam(fam) suffices.
fn lookup<'a>(tbl: &'a HashMap<String,(i64,i64,i64,i64)>, fam: &str) -> Option<&'a (i64,i64,i64,i64)> {
    tbl.get(&norm_fam(fam))
}

// returns (isBegin, isEnd, type c/p)
fn classify_is(label: &str, b: i64, e: i64, tir: &Option<Tir>,
               tir_tbl: &HashMap<String,(i64,i64,i64,i64)>,
               tpase_tbl: &HashMap<String,(i64,i64,i64)>) -> (i64, i64, char) {
    let fam = label.split('-').next().expect("split always yields at least one element");
    if let Some(t) = tir {
        return (t.start1.min(t.start2), t.end1.max(t.end2), 'c');
    }
    let info = lookup(tir_tbl, fam);
    let has_tir = info.map(|t| t.3).unwrap_or(-1);
    let lenb = tpase_tbl.get(&norm_fam(fam)).copied();
    let orflen = e - b + 1;
    if has_tir == 0 {
        if let Some((lo, hi, _)) = lenb {
            if lo <= orflen && orflen <= hi { return (b, e, 'c'); }
        }
    }
    (b, e, 'p')
}

// ---------------------------------------------------------- 5. assembly
struct Call {
    seqid: String, family: String, tier: String, is_begin: i64, is_end: i64, is_len: i64,
    orf_begin: i64, orf_end: i64, strand: char, evalue: f64, qcov: f64, pident: f64,
    typ: char, tir: Option<Tir>, ncopy: i64,
    // (C) nt neg-pos FP-control verdict (which lineage the call nt resembles): IS | shared | host | unresolved
    fp_flag: String,
}

fn dedup(mut calls: Vec<Call>, min_overlap: f64) -> Vec<Call> {
    let mut by_contig: HashMap<String, Vec<Call>> = HashMap::new();
    for c in calls.drain(..) { by_contig.entry(c.seqid.clone()).or_default().push(c); }
    let mut kept: Vec<Call> = Vec::new();
    for (_, mut cs) in by_contig {
        // total order: best-E first, then canonical (orfBegin,orfEnd,strand) so the greedy
        // keep is deterministic even when two overlapping ORFs tie on E-value.
        cs.sort_by(|a, b| a.evalue.partial_cmp(&b.evalue).unwrap()
            .then(a.orf_begin.cmp(&b.orf_begin))
            .then(a.orf_end.cmp(&b.orf_end))
            .then(a.strand.cmp(&b.strand)));
        let mut acc: Vec<Call> = Vec::new();
        for c in cs {
            let (cb, ce) = (c.orf_begin, c.orf_end);
            let clash = acc.iter().any(|a| {
                let ov = ce.min(a.orf_end) - cb.max(a.orf_begin);
                ov > 0 && ov as f64 >= min_overlap * ((ce - cb).min(a.orf_end - a.orf_begin)) as f64
            });
            if !clash { acc.push(c); }
        }
        kept.extend(acc);
    }
    kept
}

// Merge adjacent co-oriented same-family ORF calls into ONE IS element. Two-ORF IS families
// (IS3/IS1/IS21/IS66 ...) encode their transposase as orfA + orfB on the same strand, joined by a
// programmed -1 frameshift; the two ORFs are immediately consecutive (small gap or slight overlap).
// dedup() only collapses ORFs that OVERLAP by >=min_overlap, so these adjacent pairs survive as two
// calls => the IS is double-counted. ISEScan reports such a pair as a single element. We chain
// consecutive calls of the same normalized family + same strand whose gap <= `gap` bp, then re-derive
// the IS boundary (TIR re-searched over the merged span) and c/p from the merged span; family/E-value/
// identity are taken from the strongest hit in the chain. Singletons pass through unchanged.
fn merge_adjacent(calls: Vec<Call>, gap: i64, contig_map: &HashMap<&str, &Vec<u8>>,
                  tir_tbl: &HashMap<String,(i64,i64,i64,i64)>,
                  tpase_tbl: &HashMap<String,(i64,i64,i64)>) -> Vec<Call> {
    let mut by_contig: HashMap<String, Vec<Call>> = HashMap::new();
    for c in calls { by_contig.entry(c.seqid.clone()).or_default().push(c); }
    let mut out: Vec<Call> = Vec::new();
    for (_, cs) in by_contig {
        let mut cs = cs;
        // total order on genomic position so chaining is deterministic
        cs.sort_by(|a, b| a.orf_begin.cmp(&b.orf_begin)
            .then(a.orf_end.cmp(&b.orf_end))
            .then(a.strand.cmp(&b.strand)));
        let n = cs.len();
        let mut slots: Vec<Option<Call>> = cs.into_iter().map(Some).collect();
        let mut i = 0;
        while i < n {
            let fam = norm_fam(&slots[i].as_ref().unwrap().family);
            let strand = slots[i].as_ref().unwrap().strand;
            let mut chain_end = slots[i].as_ref().unwrap().orf_end;
            let mut j = i + 1;
            while j < n {
                let cj = slots[j].as_ref().unwrap();
                if norm_fam(&cj.family) == fam && cj.strand == strand
                    && cj.orf_begin - chain_end <= gap {
                    chain_end = chain_end.max(cj.orf_end);
                    j += 1;
                } else { break; }
            }
            if j - i == 1 {
                out.push(slots[i].take().unwrap());
            } else {
                let chain: Vec<Call> = (i..j).map(|k| slots[k].take().unwrap()).collect();
                let mb = chain.iter().map(|c| c.orf_begin).min().unwrap();
                let me = chain.iter().map(|c| c.orf_end).max().unwrap();
                // strongest hit = base for family/tier/evalue/qcov/pident (total-order tiebreak)
                let bi = (0..chain.len()).min_by(|&a, &b|
                    chain[a].evalue.partial_cmp(&chain[b].evalue).unwrap()
                        .then(chain[a].orf_begin.cmp(&chain[b].orf_begin))
                        .then(chain[a].orf_end.cmp(&chain[b].orf_end))
                        .then(chain[a].strand.cmp(&chain[b].strand))).unwrap();
                let base = &chain[bi];
                let label = base.family.clone();
                let tir = contig_map.get(base.seqid.as_str())
                    .and_then(|seq| find_tir(seq, mb, me, &label, tir_tbl));
                let (is_b, is_e, typ) = classify_is(&label, mb, me, &tir, tir_tbl, tpase_tbl);
                out.push(Call {
                    seqid: base.seqid.clone(), family: label, tier: base.tier.clone(),
                    is_begin: is_b, is_end: is_e, is_len: is_e - is_b + 1,
                    orf_begin: mb, orf_end: me, strand: base.strand,
                    evalue: base.evalue, qcov: base.qcov, pident: base.pident,
                    typ, tir, ncopy: 1,
                    fp_flag: "IS".into(),
                });
            }
            i = j;
        }
    }
    out
}

// Target-site duplication (TSD): upon insertion an IS duplicates a short host target, leaving an
// identical DIRECT repeat immediately flanking each end (...[TSD][element][TSD]...). A TSD is thus a
// homology-INDEPENDENT positive signal of transposition, complementary to the TIR. Family TSD lengths
// run ~3-14 bp, so we look for a direct repeat of length >= MIN_TSD whose left copy abuts is_begin and
// right copy abuts is_end, allowing a small shift because the called boundary can be a few bp off.
// Conservative (flush-ish, exact, >= MIN_TSD) to keep the chance-match rate low.
const MIN_TSD: usize = 4;
const MAX_TSD: usize = 14;
fn has_tsd(seq: &[u8], is_begin: i64, is_end: i64) -> bool {
    const W: i64 = 12;      // flank window scanned on each side
    const SHIFT: usize = 2; // boundary-imprecision tolerance
    let n = seq.len() as i64;
    let (lb, le) = ((is_begin - 1 - W).max(0), (is_begin - 1).max(0)); // upstream flank [lb,le)
    let (rb, re) = (is_end.min(n), (is_end + W).min(n));               // downstream flank [rb,re)
    if le <= lb || re <= rb { return false; }
    let left = &seq[lb as usize..le as usize];   // ends at the left boundary
    let right = &seq[rb as usize..re as usize];  // starts at the right boundary
    for l in (MIN_TSD..=MAX_TSD).rev() {
        for ls in 0..=SHIFT {
            if left.len() < l + ls { break; }
            let lc = &left[left.len() - l - ls..left.len() - ls]; // l bp ending ls before boundary
            for rs in 0..=SHIFT {
                if right.len() < l + rs { break; }
                let rc = &right[rs..rs + l];                       // l bp starting rs after boundary
                if lc == rc { return true; }
            }
        }
    }
    false
}

// ---------------------------------------------------------- (C) nt neg-pos FP-control module
// The HTH/DDE fold is shared between IS transposases and host proteins (regulators, nucleases,
// recombinases) at the protein level -> the aa profile search cannot separate them (marA/soxS
// called IS4). But at NUCLEOTIDE level, IS lineages and host lineages diverge, so the call's nt
// discriminates. For each call, compare its nt to a positive IS-nt reference (db_dir/fpc/pos) and
// a negative host-nt reference (db_dir/fpc/neg): pos-lineage match confirms IS; host-lineage match
// with no IS support flags a likely FP. Two homology-INDEPENDENT structural signals (TIR and TSD)
// add positive support that survives when nt homology is silent (divergent/novel IS): a call with
// no nt evidence either way but a TIR/TSD is "putative" (structure suggests IS), not "unresolved";
// and structure protects a host-leaning call from the droppable "host" verdict (never silently drop
// a structurally-supported call -- the neg DB is contamination-fragile, TIR/TSD are IS-specific).
// Validated on MG1655 (marA/soxS/fimB flagged, 49 curated IS retained) and 1500 GTDB strains
// (85% regulator-FP rejection, robust to genus-level holdout => not overfit).
fn fp_control_module(calls: &mut [Call], contig_map: &HashMap<&str, &Vec<u8>>,
                     db_dir: &Path, tmp: &Path, nthread: usize)
    -> Result<(), Box<dyn std::error::Error>> {
    let fpc = db_dir.join("fpc");
    if !fpc.join("refset.dbtype").exists() { return Ok(()); } // no FP-control DB -> skip (module optional)
    // write each call's element nt as >c{index}
    let qfa = tmp.join("fpcq.fna");
    let mut nq = 0usize;
    {
        let mut o = BufWriter::new(fs::File::create(&qfa)
            .map_err(|e| format!("cannot create '{}': {}", qfa.display(), e))?);
        for (i, c) in calls.iter().enumerate() {
            if let Some(seq) = contig_map.get(c.seqid.as_str()) {
                let b = c.is_begin.max(1) as usize; let e = c.is_end as usize;
                if e <= seq.len() && e >= b {
                    let mut sub = seq[b-1..e].to_vec();
                    if c.strand == '-' { sub = revcomp(&sub); }
                    if sub.len() >= 60 {
                        writeln!(o, ">c{}", i)?;
                        o.write_all(&sub)?; o.write_all(b"\n")?;
                        nq += 1;
                    }
                }
            }
        }
    }
    if nq == 0 { return Ok(()); } // no extractable query nt (0 calls / all too short) -> mmseqs createdb would fail; all calls keep default "IS"
    let run = |args: &[&str]| -> Result<(), Box<dyn std::error::Error>> {
        let st = Command::new("mmseqs").args(args)
            .stdout(std::process::Stdio::null()).stderr(std::process::Stdio::null())
            .status().map_err(|e| format!("failed to run mmseqs (is it on PATH?): {}", e))?;
        if !st.success() {
            return Err(format!("mmseqs {:?} failed (exit {:?})",
                args.first().copied().unwrap_or(""), st.code()).into());
        }
        Ok(())
    };
    let nt = nthread.to_string();
    let qdb = tmp.join("fpcqdb");
    let qfa_s = qfa.to_str().expect("fpc query fasta path is valid UTF-8");
    let qdb_s = qdb.to_str().expect("fpc query DB path is valid UTF-8");
    run(&["createdb", qfa_s, qdb_s, "--dbtype", "2", "-v", "0"])?;
    // Best bits per call against the IS(pos) and host(neg) references, in ONE search over the combined
    // refset (targets prefixed "P|" = IS-lineage, "N|" = host-lineage) -> ~40% faster than two searches.
    let mut posb: HashMap<usize,f64> = HashMap::new();
    let mut negb: HashMap<usize,f64> = HashMap::new();
    let (res, aln, tp) = (tmp.join("fpcres"), tmp.join("fpcaln.tsv"), tmp.join("fpctp"));
    let refset = fpc.join("refset");
    let refset_s = refset.to_str().expect("fpc refset path is valid UTF-8");
    let res_s = res.to_str().expect("fpc results path is valid UTF-8");
    let aln_s = aln.to_str().expect("fpc alignment path is valid UTF-8");
    let tp_s = tp.to_str().expect("fpc tp path is valid UTF-8");
    run(&["search", qdb_s, refset_s, res_s,
          tp_s, "--search-type", "3", "-s", "7", "--max-seqs", "50",
          "-e", "1e-3", "--threads", &nt, "-v", "0"])?;
    run(&["convertalis", qdb_s, refset_s, res_s,
          aln_s, "--format-output", "query,target,bits", "-v", "0"])?;
    let afile = fs::File::open(&aln)
        .map_err(|e| format!("cannot open fpc alignment '{}': {}", aln.display(), e))?;
    for line in BufReader::new(afile).lines() {
        let line = line.map_err(|e| format!("error reading fpc alignment '{}': {}", aln.display(), e))?;
        let p: Vec<&str> = line.split('\t').collect();
        if p.len() < 3 { continue; }
        let idx = match p[0].strip_prefix('c').and_then(|s| s.parse::<usize>().ok()) { Some(i)=>i, None=>continue };
        let bits: f64 = p[2].parse().unwrap_or(0.0);
        let tbl = if let Some(_) = p[1].strip_prefix("P|") { &mut posb }
                  else if let Some(_) = p[1].strip_prefix("N|") { &mut negb } else { continue };
        let ent = tbl.entry(idx).or_insert(0.0);
        if bits > *ent { *ent = bits; }
    }
    for (i, c) in calls.iter_mut().enumerate() {
        let pb = *posb.get(&i).unwrap_or(&0.0);
        let nb = *negb.get(&i).unwrap_or(&0.0);
        // homology-independent IS-architecture support (TIR present, or a flanking TSD direct repeat)
        let structural = c.tir.is_some()
            || contig_map.get(c.seqid.as_str()).map(|s| has_tsd(s, c.is_begin, c.is_end)).unwrap_or(false);
        c.fp_flag = fp_verdict(pb, nb, structural).into();
    }
    Ok(())
}

// The 5-tier fpFlag decision, isolated as a pure fn so it is unit-testable and the semantics are locked.
// pb/nb = best bits vs the IS(pos)/host(neg) nt references; structural = TIR or TSD present.
fn fp_verdict(pb: f64, nb: f64, structural: bool) -> &'static str {
    if pb > 0.0 && pb >= nb {                       // nt matches IS lineage (homology-confirmed)
        "IS"
    } else if nb > pb {                             // nt leans host
        if !structural && pb == 0.0 { "host" }      // host-lineage only, no IS support of any kind -> droppable
        else { "shared" }                           // matches both / structure-vs-nt conflict -> keep
    } else {                                        // pb == 0 && nb == 0: no nt evidence either way
        if structural { "putative" }                // structure (TIR/TSD) suggests IS
        else { "unresolved" }                       // truly undecided
    }
}

// Format like Python's "%.1e" (signed exponent, zero-padded to >=2 digits) so the
// E-value column is byte-identical to pv2's. Rust's {:.1e} drops the sign and padding
// (e.g. 0.0 -> "0.0e0"); Python gives "0.0e+00".
fn py_e(x: f64) -> String {
    let s = format!("{:.1e}", x);
    let (m, e) = s.split_once('e').expect("{:.1e} formatting always contains 'e'");
    let exp: i32 = e.parse().unwrap_or(0);
    format!("{}e{}{:02}", m, if exp < 0 { "-" } else { "+" }, exp.abs())
}

fn write_tsv(calls: &[Call], path: &Path) -> Result<(), Box<dyn std::error::Error>> {
    let mut o = BufWriter::new(fs::File::create(path)
        .map_err(|e| format!("cannot create '{}': {}", path.display(), e))?);
    writeln!(o, "seqID\tfamily\ttier\tisBegin\tisEnd\tisLen\tncopy4is\torfBegin\torfEnd\tstrand\tE-value\tqcov\tpident\ttype\tirLen\tirId\tstart2\tend2\ttir\tfpFlag")?;
    let mut idx: Vec<usize> = (0..calls.len()).collect();
    idx.sort_by(|&a, &b| (calls[a].seqid.as_str(), calls[a].is_begin)
        .cmp(&(calls[b].seqid.as_str(), calls[b].is_begin)));
    for &i in &idx {
        let c = &calls[i];
        let (irlen, irid, s2, e2, tirstr) = match &c.tir {
            Some(t) => (t.irlen, t.irid, t.start2, t.end2,
                        format!("{}:{}", String::from_utf8_lossy(&t.seq1), String::from_utf8_lossy(&t.seq2))),
            None => (0, 0, 0, 0, "-".to_string()),
        };
        writeln!(o, "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{:.2}\t{:.0}\t{}\t{}\t{}\t{}\t{}\t{}\t{}",
            c.seqid, c.family, c.tier, c.is_begin, c.is_end, c.is_len, c.ncopy,
            c.orf_begin, c.orf_end, c.strand, py_e(c.evalue), c.qcov, c.pident, c.typ,
            irlen, irid, s2, e2, tirstr, c.fp_flag)?;
    }
    Ok(())
}

fn write_gff(calls: &[Call], path: &Path) -> Result<(), Box<dyn std::error::Error>> {
    let mut o = BufWriter::new(fs::File::create(path)
        .map_err(|e| format!("cannot create '{}': {}", path.display(), e))?);
    writeln!(o, "##gff-version 3")?;
    let mut idx: Vec<usize> = (0..calls.len()).collect();
    idx.sort_by(|&a, &b| (calls[a].seqid.as_str(), calls[a].is_begin)
        .cmp(&(calls[b].seqid.as_str(), calls[b].is_begin)));
    for (n, &i) in idx.iter().enumerate() {
        let c = &calls[i];
        let id = n + 1;
        writeln!(o, "{}\trust-ise\tinsertion_sequence\t{}\t{}\t.\t{}\t.\tID=IS_{};family={};type={};tier={};ncopy={};evalue={};orf={}..{}",
            c.seqid, c.is_begin, c.is_end, c.strand, id, c.family, c.typ, c.tier, c.ncopy, py_e(c.evalue), c.orf_begin, c.orf_end)?;
        if c.fp_flag != "IS" {
            writeln!(o, "# IS_{} fpFlag={}", id, c.fp_flag)?;
        }
        if let Some(t) = &c.tir {
            writeln!(o, "{}\trust-ise\tterminal_inverted_repeat\t{}\t{}\t.\t{}\t.\tParent=IS_{};irLen={};irId={}",
                c.seqid, t.start1, t.end1, c.strand, id, t.irlen, t.irid)?;
            writeln!(o, "{}\trust-ise\tterminal_inverted_repeat\t{}\t{}\t.\t{}\t.\tParent=IS_{};irLen={};irId={}",
                c.seqid, t.start2, t.end2, c.strand, id, t.irlen, t.irid)?;
        }
    }
    Ok(())
}

fn write_sum(calls: &[Call], path: &Path) -> Result<(), Box<dyn std::error::Error>> {
    let mut o = BufWriter::new(fs::File::create(path)
        .map_err(|e| format!("cannot create '{}': {}", path.display(), e))?);
    let mut byfam: HashMap<String, Vec<usize>> = HashMap::new();
    for (i, c) in calls.iter().enumerate() {
        byfam.entry(c.family.split('-').next()
            .expect("split always yields at least one element").to_string()).or_default().push(i);
    }
    let mut fams: Vec<&String> = byfam.keys().collect();
    // total order: nIS desc, then family name asc — deterministic on count ties.
    fams.sort_by(|a, b| byfam[*b].len().cmp(&byfam[*a].len()).then_with(|| a.cmp(b)));
    writeln!(o, "family\tnIS\tcomplete\tpartial\tmeanIsLen")?;
    let (mut tot, mut nc, mut np) = (0i64, 0i64, 0i64);
    for f in fams {
        let cs = &byfam[f];
        let comp = cs.iter().filter(|&&i| calls[i].typ == 'c').count() as i64;
        let part = cs.len() as i64 - comp;
        let mean: i64 = cs.iter().map(|&i| calls[i].is_len).sum::<i64>() / cs.len() as i64;
        writeln!(o, "{}\t{}\t{}\t{}\t{}", f, cs.len(), comp, part, mean)?;
        tot += cs.len() as i64; nc += comp; np += part;
    }
    writeln!(o, "TOTAL\t{}\t{}\t{}\t-", tot, nc, np)?;
    // FP-control summary: how many calls the nt neg-pos discriminator supports vs flags
    if calls.iter().any(|c| c.fp_flag != "IS") {
        let cnt = |f: &str| calls.iter().filter(|c| c.fp_flag == f).count();
        writeln!(o, "#fpControl\tIS\tputative\tshared\thost\tunresolved")?;
        writeln!(o, "#fpControl\t{}\t{}\t{}\t{}\t{}",
            cnt("IS"), cnt("putative"), cnt("shared"), cnt("host"), cnt("unresolved"))?;
    }
    Ok(())
}

// ---------------------------------------------------------- tables (== Python reference constants)
fn tir_table() -> HashMap<String,(i64,i64,i64,i64)> {
    [("IS1",(8,67,14,1)),("IS110",(2,31,14,-1)),("IS1182",(8,44,10,1)),("IS1380",(7,39,10,1)),
     ("IS1595",(10,43,15,1)),("IS1634",(11,32,12,1)),("IS200/IS605",(10000,0,10000,0)),
     ("IS21",(8,76,10,1)),("IS256",(8,48,15,1)),("IS3",(7,54,10,-1)),("IS30",(11,50,12,1)),
     ("IS4",(8,67,12,1)),("IS481",(5,52,10,1)),("IS5",(7,45,14,1)),("IS6",(12,36,14,1)),
     ("IS607",(12,46,12,-1)),("IS630",(3,92,11,1)),("IS66",(11,144,11,1)),("IS701",(12,38,12,1)),
     ("IS91",(11,21,11,-1)),("IS982",(11,35,11,1)),("ISAS1",(12,34,12,1)),("ISAZO13",(18,48,18,1)),
     ("ISH3",(11,31,15,1)),("ISKRA4",(15,40,18,1)),("ISL3",(6,50,11,1)),("ISNCY",(4,52,13,-1)),
     ("new",(10,50,20,-1))]
    .iter().map(|(k,v)| (norm_fam(k), *v)).collect()   // normalized keys -> O(1) lookup
}
fn tpase_table() -> HashMap<String,(i64,i64,i64)> {
    [("IS1",(666,1119,252)),("IS110",(603,1380,156)),("IS1182",(822,1731,570)),("IS1380",(1158,1554,1158)),
     ("IS1595",(576,1158,426)),("IS1634",(1314,1875,1314)),("IS200/IS605",(366,1482,147)),("IS21",(882,1758,231)),
     ("IS256",(990,1389,990)),("IS3",(441,1581,120)),("IS30",(540,1419,189)),("IS4",(570,1629,219)),
     ("IS481",(447,1794,447)),("IS5",(360,1908,75)),("IS6",(528,1062,246)),("IS607",(768,1653,453)),
     ("IS630",(510,1194,318)),("IS66",(354,1695,165)),("IS701",(921,1410,921)),("IS91",(648,1548,648)),
     ("IS982",(627,981,429)),("ISAS1",(594,1329,189)),("ISAZO13",(1203,2094,513)),("ISH3",(573,1206,549)),
     ("ISKRA4",(1047,1719,114)),("ISL3",(414,1716,408)),("ISNCY",(573,1815,123)),("new",(300,2100,50))]
    .iter().map(|(k,v)| (norm_fam(k), *v)).collect()   // normalized keys -> O(1) lookup
}

// ---------------------------------------------------------- main
const USAGE: &str = "\
rust-ise \u{2014} ISEScan-equivalent insertion-sequence (IS) scanner (Rust-native)

USAGE:
    rust-ise -i <contigs.fasta> -o <out_dir> [OPTIONS]

REQUIRED:
    -i, --seqfile <FILE>   input contigs FASTA
    -o, --output  <DIR>    output directory

OPTIONS:
    -t, --threads <N>      worker threads (default: 8)
        --db <DIR>         IS database directory (else $RUSTISE_DB)
        --keep             keep intermediate files
        --strict           drop host-lineage false-positive (fpFlag=host) calls
    -h, --help             print this help and exit

The IS database (external MMseqs2 profile DB) must be provided via --db or the
RUSTISE_DB environment variable, and `mmseqs` must be on PATH. See README.
";

// Fetch a required option value at args[i], or print usage and exit if missing.
fn req(args: &[String], i: usize, flag: &str) -> String {
    args.get(i).cloned().unwrap_or_else(|| {
        eprintln!("error: {} requires a value\n", flag);
        eprint!("{}", USAGE);
        std::process::exit(2);
    })
}

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let t0 = Instant::now();
    let args: Vec<String> = std::env::args().collect();
    let mut seqfile = String::new(); let mut output = String::new();
    let mut threads = 8usize; let mut keep = false; let mut strict = false;
    let mut db_dir: Option<PathBuf> = None;   // resolve from --db, then $RUSTISE_DB
    let mut i = 1;
    while i < args.len() {
        match args[i].as_str() {
            "-h" | "--help" => { print!("{}", USAGE); return Ok(()); }
            "-i" | "--seqfile" => { i+=1; seqfile = req(&args, i, "-i/--seqfile"); }
            "-o" | "--output"  => { i+=1; output  = req(&args, i, "-o/--output"); }
            "-t" | "--threads" => {
                i+=1;
                threads = req(&args, i, "-t/--threads").parse().unwrap_or_else(|_| {
                    eprintln!("error: -t/--threads expects an integer\n");
                    eprint!("{}", USAGE);
                    std::process::exit(2);
                });
            }
            "--keep" => { keep = true; }
            "--strict" => { strict = true; } // drop fpFlag=host calls (host-lineage FPs) from output
            "--db" => { i+=1; db_dir = Some(PathBuf::from(req(&args, i, "--db"))); }
            other => { eprintln!("error: unknown argument '{}'\n", other); eprint!("{}", USAGE); std::process::exit(2); }
        }
        i += 1;
    }
    if seqfile.is_empty() || output.is_empty() {
        eprintln!("error: -i/--seqfile and -o/--output are required\n");
        eprint!("{}", USAGE);
        std::process::exit(2);
    }
    // DB is external data (not bundled in the crate). Resolve --db, else $RUSTISE_DB, else error.
    let db_dir = db_dir
        .or_else(|| std::env::var_os("RUSTISE_DB").map(PathBuf::from))
        .unwrap_or_else(|| {
            eprintln!("error: IS database not set. Pass --db <dir> or set RUSTISE_DB (see README).");
            std::process::exit(2);
        });
    if !db_dir.join("mmdb_union").join("profileDb.dbtype").exists() {
        eprintln!("error: '{}' is not a valid rust-ise DB (missing mmdb_union/profileDb). See README.",
            db_dir.display());
        std::process::exit(2);
    }
    rayon::ThreadPoolBuilder::new().num_threads(threads).build_global().ok();
    fs::create_dir_all(&output)
        .map_err(|e| format!("cannot create output dir '{}': {}", output, e))?;
    let tmp = Path::new(&output).join("rust-ise_tmp");
    fs::create_dir_all(&tmp)
        .map_err(|e| format!("cannot create tmp dir '{}': {}", tmp.display(), e))?;

    let contigs = read_fasta(&seqfile)?;
    let contig_map: HashMap<&str, &Vec<u8>> = contigs.iter().map(|(k,v)| (k.as_str(), v)).collect();
    log!(t0, "contigs: {}", contigs.len());

    // 1. ORFs: rustygal LIBRARY meta gene-finding (rust-pure, no subprocess) + native edge recovery.
    // Call run_meta per contig (parallel), concatenate the -a protein FASTA bytes in FASTA order
    // (== the binary's -a output), then parse exactly as before.
    let t1 = Instant::now();
    let meta = rustygal::meta_api::meta_bins();
    let per_contig: Vec<Vec<u8>> = contigs.par_iter().enumerate()
        .map(|(i, (hdr, dna))| rustygal::meta_api::run_meta(i as i32 + 1, hdr, dna, &meta).trans_faa)
        .collect();
    let faa: Vec<u8> = per_contig.concat();
    let mut orfs = parse_rustygal_faa(&faa)?;
    let n_pyr: usize = orfs.values().map(|v| v.len()).sum();
    // edge recovery in parallel over contigs
    let edges: Vec<(String, Vec<Orf>)> = contigs.par_iter().map(|(cid, seq)| {
        let empty = Vec::new();
        let pyr = orfs.get(cid).unwrap_or(&empty);
        (cid.clone(), edge_orfs(seq, pyr, 30, 3))
    }).collect();
    let mut n_edge = 0usize;
    for (cid, mut e) in edges { n_edge += e.len(); orfs.entry(cid).or_default().append(&mut e); }
    // write proteome.faa with headers {cid}_{b}_{e}_{strand}[_edge]
    let proteome = Path::new(&output).join("proteome.faa");
    {
        let mut o = BufWriter::new(fs::File::create(&proteome)
            .map_err(|e| format!("cannot create '{}': {}", proteome.display(), e))?);
        // Iterate contigs in FASTA order (NOT HashMap order) so the proteome byte-order is
        // deterministic and identical to pv2's — mmseqs prefilter is order-sensitive at the
        // --max-seqs margin, so a fixed input order removes the run-to-run / cross-tool jitter.
        for (cid, _) in &contigs {
            if let Some(list) = orfs.get(cid) {
                for orf in list {
                    let suffix = if orf.edge { "_edge" } else { "" };
                    writeln!(o, ">{}_{}_{}_{}{}", cid, orf.begin, orf.end, orf.strand, suffix)?;
                    o.write_all(&orf.prot)?;
                    o.write_all(b"\n")?;
                }
            }
        }
    }
    log!(t0, "ORFs: {} rustygal + {} edge  ({:.1}s)", n_pyr, n_edge, t1.elapsed().as_secs_f64());

    // 2. search
    let t2 = Instant::now();
    let hits = mmseqs_search(&proteome, &db_dir, &tmp, threads)?;
    log!(t0, "search: {} ORFs hit  ({:.1}s)", hits.len(), t2.elapsed().as_secs_f64());

    // 3. family call + TIR + classify (parallel over ORFs)
    let t3 = Instant::now();
    let tir_tbl = tir_table(); let tpase_tbl = tpase_table();
    let hit_vec: Vec<(&String, &Vec<Hit>)> = hits.iter().collect();
    let raw_calls: Vec<Call> = hit_vec.par_iter().filter_map(|(orf, oh)| {
        let (label, tier, best) = family_call(oh)?;
        if tier == "novel" { return None; }
        // parse orf id: contig_begin_end_strand[_edge]
        let core: &str = if orf.ends_with("_edge") { &orf[..orf.len()-5] } else { orf };
        // rsplit into 3 from right: strand, end, begin
        let mut it = core.rsplitn(4, '_');
        let strand = it.next()?; let e: i64 = it.next()?.parse().ok()?;
        let b: i64 = it.next()?.parse().ok()?; let contig = it.next()?;
        let strand_c = strand.chars().next().unwrap_or('+');
        let tir = contig_map.get(contig)
            .and_then(|seq| find_tir(seq, b, e, &label, &tir_tbl));
        let (is_b, is_e, typ) = classify_is(&label, b, e, &tir, &tir_tbl, &tpase_tbl);
        Some(Call {
            seqid: contig.to_string(), family: label, tier, is_begin: is_b, is_end: is_e,
            is_len: is_e - is_b + 1, orf_begin: b, orf_end: e, strand: strand_c,
            evalue: best.evalue, qcov: best.qcov, pident: best.pident, typ, tir, ncopy: 1,
            fp_flag: "IS".into(),
        })
    }).collect();
    let calls = dedup(raw_calls, 0.5);
    // collapse two-ORF (orfA+orfB) IS into one element so they are not double-counted (see merge_adjacent)
    let mut calls = merge_adjacent(calls, 50, &contig_map, &tir_tbl, &tpase_tbl);
    log!(t0, "IS calls: {} after dedup+merge  ({:.1}s for family-call+TIR+dedup)", calls.len(), t3.elapsed().as_secs_f64());

    // 3b. (C) nt neg-pos FP-control: flag host-lineage calls (marA/soxS-type FPs) using nt discrimination
    let t_c = Instant::now();
    fp_control_module(&mut calls, &contig_map, &db_dir, &tmp, threads)?;
    let cnt = |f: &str| calls.iter().filter(|c| c.fp_flag == f).count();
    log!(t0, "fp-control: {} host, {} shared, {} putative, {} unresolved (nt neg-pos + TIR/TSD)  ({:.1}s)",
        cnt("host"), cnt("shared"), cnt("putative"), cnt("unresolved"), t_c.elapsed().as_secs_f64());
    if strict {
        let before = calls.len();
        calls.retain(|c| c.fp_flag != "host"); // precision mode: drop host-lineage FPs
        log!(t0, "strict: dropped {} host-flagged calls ({} -> {})", before - calls.len(), before, calls.len());
    }

    // 4. outputs
    let stem = Path::new(&output);
    write_tsv(&calls, &stem.join("rust-ise.tsv"))?;
    write_gff(&calls, &stem.join("rust-ise.gff"))?;
    write_sum(&calls, &stem.join("rust-ise.sum"))?;
    if !keep { let _ = fs::remove_dir_all(&tmp); }
    log!(t0, "DONE  total {:.1}s -> {}/rust-ise.tsv  ({} IS)", t0.elapsed().as_secs_f64(), output, calls.len());
    Ok(())
}

// ---------------------------------------------------------- unit tests (golden: lock the logic)
#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn fp_verdict_five_tiers() {
        // IS: nt matches IS lineage (pb>0 & pb>=nb), incl. the tie pb==nb>0
        assert_eq!(fp_verdict(100.0, 0.0, false), "IS");
        assert_eq!(fp_verdict(100.0, 100.0, false), "IS");
        // shared: nt leans host but pb>0 (matches both, host closer)
        assert_eq!(fp_verdict(50.0, 100.0, false), "shared");
        // host: host-lineage only, pb==0, no structure -> the only droppable tier
        assert_eq!(fp_verdict(0.0, 100.0, false), "host");
        // structure PROTECTS a host-leaning call from "host" -> "shared" (never silently dropped)
        assert_eq!(fp_verdict(0.0, 100.0, true), "shared");
        // no nt evidence either way: structure -> putative, else unresolved
        assert_eq!(fp_verdict(0.0, 0.0, true), "putative");
        assert_eq!(fp_verdict(0.0, 0.0, false), "unresolved");
    }

    #[test]
    fn tsd_detected_and_rejected() {
        // ...[GGGGGG][ACGTAC]<TTTTTTTTTT>[ACGTAC][CCCCCC]... : 6bp TSD flanking element at 1-based [13,22]
        let s = b"GGGGGGACGTACTTTTTTTTTTACGTACCCCCC";
        assert!(has_tsd(s, 13, 22), "flanking 6bp direct repeat should be found");
        // no flanking direct repeat -> false
        let s2 = b"GGGGGGGGGGGGTTTTTTTTTTCCCCCCCCCCCC";
        assert!(!has_tsd(s2, 13, 22), "distinct flanks have no TSD");
        // 3bp repeat is below MIN_TSD(4) -> not counted
        let s3 = b"GGGGGGGGGACGTTTTTTTTTTACGCCCCCCCCC";
        assert!(!has_tsd(s3, 13, 22), "3bp repeat must not count as TSD");
    }

    #[test]
    fn revcomp_translate_norm() {
        assert_eq!(revcomp(b"ACGT"), b"ACGT");           // palindrome
        assert_eq!(revcomp(b"AAAA"), b"TTTT");
        assert_eq!(revcomp(b"ATGC"), b"GCAT");
        assert_eq!(translate(b"ATGAAA"), b"MK");         // ATG=Met, AAA=Lys
        assert_eq!(norm_fam("IS200/IS605"), "IS200IS605");
        assert_eq!(norm_fam("is_3"), "IS3");
    }

    #[test]
    fn py_e_matches_python_format() {
        assert_eq!(py_e(0.0), "0.0e+00");
        assert_eq!(py_e(1e-5), "1.0e-05");
    }
}