rust-ise 0.2.2

// 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 / $ISSCAN_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;

// 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) ----
// SIG_EVALUE: the primary significance gate — an ORF's best profile hit must be at least this
// significant to be called an IS at all (family_call). This is the dominant code-side precision
// lever. Empirically calibrated on E. coli MG1655 vs RefSeq (50 mobile_genetic_element truth
// features): loosening to 1e-10 admits ~5 low-significance false positives; tightening past ~1e-40
// starts discarding real IS calls. 1e-30 is the knee — it removes those FPs while retaining every
// true IS element (MG1655: precision 0.72 -> 0.76, recall unchanged at 0.98 / 49-of-50 truth).
// PRECISION/RECALL TRADEOFF: smaller value => fewer FPs but risks dropping divergent IS transposases.
const SIG_EVALUE: f64 = 1e-30;
// SIG_QCOV: minimum query (ORF) coverage of the winning profile hit. A real IS transposase ORF
// usually aligns across most of the profile, so a very partial hit is a weak-evidence FP candidate.
// IMPORTANT / HONEST CAVEAT: on MG1655 query coverage is a WEAK discriminator — false-positive qcov
// (median ~0.88) overlaps true-positive qcov (median ~0.89, min 0.13), so a strict floor (e.g. 0.7)
// would destroy recall (0.98 -> 0.78) for almost no precision gain. This floor is therefore set LOW
// as a conservative safety net that trims only the most extreme fragmentary hits without harming
// recall; the E-value floor above does the real work. Raise only with fresh empirical validation.
const SIG_QCOV: f64 = 0.35;
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 (== constants_isscan.py) ----
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)).
// constants_isscan filters4ssw4trial = (match2, mismatch6, gap_open2, gap_extend2).
const SW_OPEN: i32 = 2;
const SW_EXTEND: i32 = 2;

// ---------------------------------------------------------- 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) -> Vec<(String, Vec<u8>)> {
    let f = BufReader::new(fs::File::open(path).expect("open fasta"));
    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.unwrap();
        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)); }
    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]) -> HashMap<String, Vec<Orf>> {
    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().unwrap();
            let parts: Vec<&str> = rest.split('#').collect();
            let b: i64 = parts[1].trim().parse().unwrap();
            let e: i64 = parts[2].trim().parse().unwrap();
            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);
    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, gpu: bool)
    -> HashMap<String, Vec<Hit>> {
    let profile_db = db_dir.join("mmdb_union").join("profileDb");
    let manifest = db_dir.join("manifest_union.tsv");
    let run = |args: &[&str]| {
        let st = Command::new("mmseqs").args(args)
            .stdout(std::process::Stdio::null()).stderr(std::process::Stdio::null())
            .status().expect("spawn mmseqs");
        assert!(st.success(), "mmseqs {:?} failed", &args[0]);
    };
    let qdb = tmp.join("qdb");
    run(&["createdb", proteome.to_str().unwrap(), qdb.to_str().unwrap(), "-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();
    for line in BufReader::new(fs::File::open(&manifest).unwrap()).lines() {
        let line = line.unwrap();
        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();
    let mut search_args: Vec<&str> = vec!["search", profile_db.to_str().unwrap(),
        qdb.to_str().unwrap(), res.to_str().unwrap(), tp.to_str().unwrap(),
        "-s", "7.5", "-e", "10", "--threads", &nt, "-v", "0"];
    if gpu { search_args.push("--gpu"); search_args.push("1"); }
    run(&search_args);
    run(&["convertalis", profile_db.to_str().unwrap(), qdb.to_str().unwrap(),
        res.to_str().unwrap(), aln.to_str().unwrap(),
        "--format-output", "query,target,evalue,qcov,pident,tstart,tend", "-v", "0"]);
    let mut hits: HashMap<String, Vec<Hit>> = HashMap::new();
    for line in BufReader::new(fs::File::open(&aln).unwrap()).lines() {
        let line = line.unwrap();
        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),
            });
        }
    }
    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();
    // primary significance gate (dominant precision lever): drop weakly-supported ORFs.
    if top.evalue > SIG_EVALUE { return None; }
    // query-coverage floor: reject extreme fragmentary hits (conservative — see SIG_QCOV doc).
    if top.qcov < SIG_QCOV { 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
#[derive(Clone)]
pub struct Tir { pub irid: i64, pub irlen: i64, pub start1: i64, pub end1: i64, pub start2: i64, pub end2: i64,
             pub seq1: Vec<u8>, pub 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().unwrap();
    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().unwrap();
    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
#[derive(Clone)]
pub struct Call {
    pub seqid: String, pub family: String, pub tier: String, pub is_begin: i64, pub is_end: i64, pub is_len: i64,
    pub orf_begin: i64, pub orf_end: i64, pub strand: char, pub evalue: f64, pub qcov: f64, pub pident: f64,
    pub typ: char, pub tir: Option<Tir>, pub ncopy: i64,
    // (C) nt neg-pos FP-control verdict (which lineage the call nt resembles): IS | shared | host | unresolved
    pub 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) {
    let fpc = db_dir.join("fpc");
    if !fpc.join("refset.dbtype").exists() { return; } // 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).unwrap());
        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).unwrap();
                        o.write_all(&sub).unwrap(); o.write_all(b"\n").unwrap();
                        nq += 1;
                    }
                }
            }
        }
    }
    if nq == 0 { return; } // no extractable query nt (0 calls / all too short) -> mmseqs createdb would fail; all calls keep default "IS"
    let run = |args: &[&str]| {
        let st = Command::new("mmseqs").args(args)
            .stdout(std::process::Stdio::null()).stderr(std::process::Stdio::null())
            .status().expect("spawn mmseqs");
        assert!(st.success(), "mmseqs {:?} failed", &args[0]);
    };
    let nt = nthread.to_string();
    let qdb = tmp.join("fpcqdb");
    run(&["createdb", qfa.to_str().unwrap(), qdb.to_str().unwrap(), "--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"));
    run(&["search", qdb.to_str().unwrap(), fpc.join("refset").to_str().unwrap(), res.to_str().unwrap(),
          tp.to_str().unwrap(), "--search-type", "3", "-s", "7", "--max-seqs", "50",
          "-e", "1e-3", "--threads", &nt, "-v", "0"]);
    run(&["convertalis", qdb.to_str().unwrap(), fpc.join("refset").to_str().unwrap(), res.to_str().unwrap(),
          aln.to_str().unwrap(), "--format-output", "query,target,bits", "-v", "0"]);
    for line in BufReader::new(fs::File::open(&aln).unwrap()).lines() {
        let line = line.unwrap(); 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();
    }
}

// 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').unwrap();
    let exp: i32 = e.parse().unwrap_or(0);
    format!("{}e{}{:02}", m, if exp < 0 { "-" } else { "+" }, exp.abs())
}

pub fn write_tsv(calls: &[Call], path: &Path) {
    let mut o = BufWriter::new(fs::File::create(path).unwrap());
    writeln!(o, "seqID\tfamily\ttier\tisBegin\tisEnd\tisLen\tncopy4is\torfBegin\torfEnd\tstrand\tE-value\tqcov\tpident\ttype\tirLen\tirId\tstart2\tend2\ttir\tfpFlag").unwrap();
    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).unwrap();
    }
}

pub fn write_gff(calls: &[Call], path: &Path) {
    let mut o = BufWriter::new(fs::File::create(path).unwrap());
    writeln!(o, "##gff-version 3").unwrap();
    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, "{}\tisscan\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).unwrap();
        if c.fp_flag != "IS" {
            writeln!(o, "# IS_{} fpFlag={}", id, c.fp_flag).unwrap();
        }
        if let Some(t) = &c.tir {
            writeln!(o, "{}\tisscan\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).unwrap();
            writeln!(o, "{}\tisscan\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).unwrap();
        }
    }
}

pub fn write_sum(calls: &[Call], path: &Path) {
    let mut o = BufWriter::new(fs::File::create(path).unwrap());
    let mut byfam: HashMap<String, Vec<usize>> = HashMap::new();
    for (i, c) in calls.iter().enumerate() {
        byfam.entry(c.family.split('-').next().unwrap().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").unwrap();
    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).unwrap();
        tot += cs.len() as i64; nc += comp; np += part;
    }
    writeln!(o, "TOTAL\t{}\t{}\t{}\t-", tot, nc, np).unwrap();
    // 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").unwrap();
        writeln!(o, "#fpControl\t{}\t{}\t{}\t{}\t{}",
            cnt("IS"), cnt("putative"), cnt("shared"), cnt("host"), cnt("unresolved")).unwrap();
    }
}

// ---------------------------------------------------------- tables (== constants_isscan.py)
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
/// Configuration for [`run`].
pub struct IseConfig {
    /// mmseqs2 thread count.
    pub threads: usize,
    /// Use MMseqs2 GPU acceleration.
    pub gpu: bool,
    /// Precision mode: drop `fpFlag == "host"` (host-lineage FP) calls.
    pub strict: bool,
    /// IS database directory (must contain `mmdb_union/profileDb*` + `manifest_union.tsv`).
    pub db_dir: PathBuf,
}

/// Run the IS-element detection pipeline on a genome FASTA, returning the calls.
///
/// `work_dir` is where the intermediate files (`proteome.faa`, mmseqs
/// `isscan_tmp/`) are written; the caller owns it and any cleanup. Requires the
/// `mmseqs` binary on PATH and a valid IS DB at `config.db_dir`. This is the exact
/// pipeline the `rust-ise` binary runs (the binary just parses args and writes the
/// TSV/GFF/.sum from the returned calls).
pub fn run(seqfile: &Path, work_dir: &Path, config: &IseConfig) -> Result<Vec<Call>, String> {
    if !config
        .db_dir
        .join("mmdb_union")
        .join("profileDb.dbtype")
        .exists()
    {
        return Err(format!(
            "'{}' is not a valid isscan DB (missing mmdb_union/profileDb)",
            config.db_dir.display()
        ));
    }
    fs::create_dir_all(work_dir).map_err(|e| format!("create {}: {e}", work_dir.display()))?;
    let tmp = work_dir.join("isscan_tmp");
    fs::create_dir_all(&tmp).map_err(|e| format!("create {}: {e}", tmp.display()))?;

    let seqfile = seqfile.to_str().ok_or("non-UTF8 seqfile path")?;
    let contigs = read_fasta(seqfile);
    let contig_map: HashMap<&str, &Vec<u8>> =
        contigs.iter().map(|(k, v)| (k.as_str(), v)).collect();

    // 1. ORFs: rustygal META gene-finding (rust-pure, no subprocess) + native edge recovery.
    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 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();
    for (cid, mut e) in edges {
        orfs.entry(cid).or_default().append(&mut e);
    }
    // proteome.faa headers {cid}_{b}_{e}_{strand}[_edge], contigs in FASTA order
    // (deterministic byte-order == pv2; mmseqs prefilter is order-sensitive).
    let proteome = work_dir.join("proteome.faa");
    {
        let mut o = BufWriter::new(fs::File::create(&proteome).map_err(|e| e.to_string())?);
        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)
                        .map_err(|e| e.to_string())?;
                    o.write_all(&orf.prot).map_err(|e| e.to_string())?;
                    o.write_all(b"\n").map_err(|e| e.to_string())?;
                }
            }
        }
    }

    // 2. search
    let hits = mmseqs_search(&proteome, &config.db_dir, &tmp, config.threads, config.gpu);

    // 3. family call + TIR + classify (parallel over ORFs)
    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
            };
            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 (see merge_adjacent)
    let mut calls = merge_adjacent(calls, 50, &contig_map, &tir_tbl, &tpase_tbl);

    // 3b. nt neg-pos FP-control: flag host-lineage calls
    fp_control_module(&mut calls, &contig_map, &config.db_dir, &tmp, config.threads);
    if config.strict {
        calls.retain(|c| c.fp_flag != "host"); // precision mode: drop host-lineage FPs
    }
    Ok(calls)
}

// ---------------------------------------------------------- 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");
    }

    // helper: a plausible IS transposase hit, tweak the fields under test
    fn hit(evalue: f64, qcov: f64, pident: f64) -> Hit {
        Hit { fam: "IS3".into(), evalue, qcov, pident, tstart: 1, tend: 300 }
    }

    #[test]
    fn family_call_sig_evalue_floor() {
        // a hit right at the floor is accepted; a hit just weaker than the floor is rejected.
        assert!(family_call(&[hit(SIG_EVALUE, 0.9, 90.0)]).is_some(),
            "E-value == SIG_EVALUE must pass");
        assert!(family_call(&[hit(SIG_EVALUE * 10.0, 0.9, 90.0)]).is_none(),
            "E-value weaker than SIG_EVALUE must be rejected");
        // sanity: the pre-tightening 1e-10 regime is now rejected (that was the FP-admitting setting).
        assert!(family_call(&[hit(1e-10, 0.9, 90.0)]).is_none(),
            "1e-10 hits must no longer qualify after tightening to 1e-30");
    }

    #[test]
    fn family_call_qcov_floor() {
        // strongly significant but below the coverage floor -> rejected as a fragmentary hit.
        assert!(family_call(&[hit(1e-50, SIG_QCOV - 0.01, 90.0)]).is_none(),
            "qcov below SIG_QCOV must be rejected even when E-value is strong");
        // at/above the floor with a strong E-value -> accepted.
        assert!(family_call(&[hit(1e-50, SIG_QCOV, 90.0)]).is_some(),
            "qcov == SIG_QCOV must pass");
        // the floor is deliberately low so mid-coverage true IS calls survive (recall protection).
        assert!(family_call(&[hit(1e-50, 0.55, 90.0)]).is_some(),
            "mid-coverage (0.55) real IS hits must be retained");
    }
}