rust-ise 0.2.3

//! `rust-ise build-db` — reproducible-from-source build of the union IS profile DB.
//!
//! Ports the canonical pipeline (`dbbuild/build_db.py` stages 00–01 + the
//! `isdb_build` chain) to pure-Rust orchestration:
//!
//!   ISOSDB.V3.fna (nt) --rustygal meta--> ORFs (fam-tagged)   ┐
//!   IS.faa (ISfinder)  --transposase filter--> ORFs (fam-tagged) ┘ → union_orfs.faa
//!     → mmseqs createdb + cluster (--min-seq-id 0.4 -c 0.8 --cluster-mode 0)
//!     → per-cluster POASTA multiple sequence alignment (IN-PROCESS, pure Rust; replaces clustalo)
//!     → mmseqs convertmsa + msa2profile → mmdb_union/profileDb
//!     → ISfinder-authority relabel → manifest_union.tsv
//!     → FP-control reference DB (module C): fpc/refset (nt P|IS + N|host pool)
//!
//! ## Reproducibility (identical DB for every user)
//!
//! Every determinism input is pinned: sources by commit SHA, rustygal + poasta by
//! crate version, and mmseqs by an exact version requirement (checked at runtime).
//! mmseqs is deterministic given byte-identical, fixed-order input regardless of
//! thread count (verified: `cluster` at 1 vs 16 threads yields identical rep/member
//! maps), so we emit every intermediate in a fixed (sorted / FASTA) order and run
//! mmseqs multi-threaded for speed. The per-profile GTDB thresholds are calibrated
//! ONCE by the maintainer against this reproducible build and shipped embedded in
//! the crate — a user rebuild reproduces the same reps, so the thresholds apply.
//!
//! ## FP-control reference DB (`--host` / `--fpc-neg`)
//!
//! `build-db` also emits `<out>/fpc/refset` — the nt neg-pos reference the scanner's
//! FP-control module ([`crate` `fp_control_module`]) uses to discriminate IS-lineage
//! calls from host-lineage false positives. It is the concat of an IS-lineage POSITIVE
//! pool (prefixed `P|`) and a host-lineage NEGATIVE pool (prefixed `N|`), `mmseqs
//! createdb --dbtype 2`'d (ports `build_fpc_db.sh`). The positive pool is derived
//! automatically from the ISOSDB per-ORF nucleotide CDS assembled during ingest; the
//! negative pool is supplied externally via `--host <host.fna>` (a host-lineage nt
//! FASTA — e.g. non-IS CDS from representative host genomes). Without `--host` the
//! fpc DB is skipped and the scanner silently no-ops the FP-control stage.

use std::io::Write;
use std::path::Path;
use std::process::Command;

/// Exact mmseqs version this DB build is pinned to (reproducibility contract).
pub const MMSEQS_PINNED_VERSION: &str = "18-8cc5c";

// Ingest submodules (raw IS sources → fam-tagged transposase ORFs).
mod isfinder;
mod isosdb;
// Source auto-download (`--fetch-sources`): pinned-SHA + sha256-verified fetch.
pub mod fetch;

/// One transposase ORF entering the union IS profile DB.
#[derive(Clone, Debug)]
pub struct DbOrf {
    /// Stable ORF id. ISOSDB: `ISOSDB<elem>_orf<n>`; ISfinder: `html.2020//<ISname>`.
    pub id: String,
    /// IS family label at ingest time (refined later by the ISfinder-authority relabel).
    pub family: String,
    /// Provenance: `"ISOSDB"` or `"ISfinder"`.
    pub source: &'static str,
    /// Protein residues (uppercase, no stop).
    pub seq: Vec<u8>,
}

/// Read a protein/nucleotide FASTA into `(id, sequence_bytes)` in file order.
/// `id` is the first whitespace-delimited token of the header.
pub fn read_fasta(path: &Path) -> Result<Vec<(String, Vec<u8>)>, String> {
    let data = std::fs::read(path).map_err(|e| format!("read {}: {e}", path.display()))?;
    let mut out: Vec<(String, Vec<u8>)> = Vec::new();
    let mut id: Option<String> = None;
    let mut seq: Vec<u8> = Vec::new();
    for line in data.split(|&b| b == b'\n') {
        if line.first() == Some(&b'>') {
            if let Some(prev) = id.take() {
                out.push((prev, std::mem::take(&mut seq)));
            }
            let hdr = String::from_utf8_lossy(&line[1..]);
            id = Some(hdr.split_whitespace().next().unwrap_or("").to_string());
        } else {
            seq.extend(line.iter().filter(|&&b| b != b'\r' && b != b' '));
        }
    }
    if let Some(prev) = id.take() {
        out.push((prev, seq));
    }
    Ok(out)
}

// ---------------------------------------------------------------- POASTA MSA
// In-process partial-order multiple sequence alignment, replacing the reference
// pipeline's external `clustalo`. Parameters match the `poasta` CLI defaults used
// in the validation experiment (mismatch 4, gap-open 6, gap-extend 2, global),
// which reproduced the clustalo reference profiles to ~100% on the benchmark.

use poasta::aligner::config::AffineMinGapCost;
use poasta::aligner::scoring::{AlignmentType, GapAffine};
use poasta::aligner::PoastaAligner;
use poasta::graphs::poa::POAGraph;
use poasta::io::fasta::poa_graph_to_fasta;

/// Multiple-sequence-align `seqs` (`(name, residues)`, **rep first** — it seeds the
/// graph and becomes the match-state reference for `msa2profile --match-mode 1`).
/// Returns `(name, aligned_row)` in MSA column order, gaps as `-`. Deterministic.
pub fn poa_msa(seqs: &[(String, Vec<u8>)]) -> Result<Vec<(String, String)>, String> {
    if seqs.is_empty() {
        return Ok(Vec::new());
    }
    // GapAffine::new(cost_mismatch, cost_gap_extend, cost_gap_open) — poasta CLI defaults.
    let scoring = GapAffine::new(4, 2, 6);
    let aligner = PoastaAligner::new(AffineMinGapCost(scoring), AlignmentType::Global);
    let mut graph: POAGraph<u32> = POAGraph::new();
    for (name, seq) in seqs {
        let weights = vec![1usize; seq.len()];
        if graph.is_empty() {
            graph
                .add_alignment_with_weights(name, seq, None, &weights)
                .map_err(|e| format!("poasta add first seq {name}: {e:?}"))?;
        } else {
            let res = aligner.align::<u32, _>(&graph, seq);
            graph
                .add_alignment_with_weights(name, seq, Some(&res.alignment), &weights)
                .map_err(|e| format!("poasta align {name}: {e:?}"))?;
        }
    }
    let mut buf: Vec<u8> = Vec::new();
    poa_graph_to_fasta(&graph, &mut buf).map_err(|e| format!("poa_graph_to_fasta: {e:?}"))?;
    Ok(parse_aligned_fasta(&buf))
}

/// Parse a tabular aligned-FASTA block (poasta `-O fasta` output) into rows.
fn parse_aligned_fasta(data: &[u8]) -> Vec<(String, String)> {
    let text = String::from_utf8_lossy(data);
    let mut rows: Vec<(String, String)> = Vec::new();
    let mut name: Option<String> = None;
    let mut seq = String::new();
    for line in text.lines() {
        if let Some(hdr) = line.strip_prefix('>') {
            if let Some(n) = name.take() {
                rows.push((n, std::mem::take(&mut seq)));
            }
            name = Some(hdr.split_whitespace().next().unwrap_or("").to_string());
        } else {
            seq.push_str(line.trim());
        }
    }
    if let Some(n) = name.take() {
        rows.push((n, seq));
    }
    rows
}

/// Serialise `msa` rows as one Stockholm block with accession `cid` (the input for
/// `mmseqs convertmsa --identifier-field 0`).
pub fn stockholm_block(cid: &str, msa: &[(String, String)]) -> String {
    let mut out = String::from("# STOCKHOLM 1.0\n");
    out.push_str(&format!("#=GF ID {cid}\n"));
    let w = msa.iter().map(|(n, _)| n.len()).max().unwrap_or(1);
    for (n, s) in msa {
        out.push_str(&format!("{n:<w$}  {s}\n"));
    }
    out.push_str("//\n");
    out
}

// ---------------------------------------------------------------- mmseqs runner
/// Run `mmseqs <args…>`, returning an error string on non-zero exit. `mmseqs` must
/// be on PATH and match [`MMSEQS_PINNED_VERSION`] (see [`check_mmseqs_version`]).
pub fn mmseqs(args: &[&str]) -> Result<(), String> {
    let status = Command::new("mmseqs")
        .args(args)
        .status()
        .map_err(|e| format!("failed to launch mmseqs (is it on PATH?): {e}"))?;
    if !status.success() {
        return Err(format!("mmseqs {} failed ({status})", args.first().unwrap_or(&"")));
    }
    Ok(())
}

/// Verify the on-PATH mmseqs matches the pinned version (reproducibility guard).
/// Returns the detected version string. A mismatch is a hard error unless `strict`
/// is false (then it warns) — different mmseqs versions can cluster differently.
pub fn check_mmseqs_version(strict: bool) -> Result<String, String> {
    let out = Command::new("mmseqs")
        .arg("version")
        .output()
        .map_err(|e| format!("failed to launch mmseqs: {e}"))?;
    let ver = String::from_utf8_lossy(&out.stdout).trim().to_string();
    // `mmseqs version` prints the full commit hash (e.g. `8cc5ce367b56…`), whereas the
    // pin is the release tag `<release>-<hash5>` (e.g. `18-8cc5c`). Accept a match when
    // the pin tag is contained verbatim OR the reported hash starts with the tag's
    // commit-prefix (the part after the last `-`).
    let commit_prefix = MMSEQS_PINNED_VERSION.rsplit('-').next().unwrap_or(MMSEQS_PINNED_VERSION);
    let matches = ver.contains(MMSEQS_PINNED_VERSION) || ver.starts_with(commit_prefix);
    if !matches {
        let msg = format!(
            "mmseqs version '{ver}' != pinned '{MMSEQS_PINNED_VERSION}'. A different \
             version may produce a different clustering/profile DB, so the embedded \
             GTDB thresholds (keyed to the pinned build's reps) may not apply."
        );
        if strict {
            return Err(msg);
        }
        eprintln!("[build-db] WARNING: {msg}");
    }
    Ok(ver)
}

// ---------------------------------------------------------------- profile DB build
/// Build a profile DB from `clusters` (each `(cid, rep_first_members)`), using
/// in-process poasta MSAs → one Stockholm file → `mmseqs convertmsa + msa2profile`.
/// `out_prefix` names the resulting `<out_prefix>` profile DB. Clusters are emitted
/// in the given order (caller sorts for determinism).
pub fn build_profiledb_poasta(
    clusters: &[(String, Vec<(String, Vec<u8>)>)],
    out_prefix: &Path,
    work: &Path,
) -> Result<(), String> {
    use rayon::prelude::*;
    std::fs::create_dir_all(work).map_err(|e| format!("mkdir {}: {e}", work.display()))?;
    let sto_path = work.join("union_msas.sto");
    // POA-align every cluster in parallel; `par_iter().collect()` preserves the input
    // (rep-sorted) order, so the concatenated Stockholm is deterministic.
    let blocks: Vec<String> = clusters
        .par_iter()
        .map(|(cid, members)| poa_msa(members).map(|msa| stockholm_block(cid, &msa)))
        .collect::<Result<Vec<_>, String>>()?;
    let mut sto = std::io::BufWriter::new(
        std::fs::File::create(&sto_path).map_err(|e| format!("create {}: {e}", sto_path.display()))?,
    );
    for b in &blocks {
        sto.write_all(b.as_bytes()).map_err(|e| format!("write stockholm: {e}"))?;
    }
    sto.flush().map_err(|e| format!("flush stockholm: {e}"))?;
    drop(sto);

    let msa_db = work.join("msaDb");
    let (sto_s, msa_s, prof_s) = (
        sto_path.to_str().ok_or("non-utf8 path")?,
        msa_db.to_str().ok_or("non-utf8 path")?,
        out_prefix.to_str().ok_or("non-utf8 path")?,
    );
    mmseqs(&["convertmsa", sto_s, msa_s, "--identifier-field", "0", "-v", "0"])?;
    // `--threads 1`: msa2profile computes byte-identical profiles per MSA regardless of
    // thread count, but multi-threaded writes interleave entries in completion order, so
    // the DB's physical byte layout (and .index offsets) varies run-to-run. Single-thread
    // writes entries in index order → byte-identical profileDb for every user (verified:
    // per-cid content identical across 1/16 threads; only physical order differed).
    mmseqs(&["msa2profile", msa_s, prof_s, "--match-mode", "1", "--threads", "1", "-v", "0"])?;
    Ok(())
}

// ---------------------------------------------------------------- FP-control refset
/// Build the FP-control reference DB `<out>/fpc/refset` (mmseqs nt DB, `--dbtype 2`).
///
/// Ports `build_fpc_db.sh` into the build-db flow: the IS-lineage POSITIVE pool
/// (each record prefixed `P|`) and the host-lineage NEGATIVE pool (prefixed `N|`)
/// are concatenated into one FASTA and `mmseqs createdb`'d as a single nt DB. The
/// scanner's `fp_control_module` searches call nt against this DB and reads the
/// `P|`/`N|` prefixes to decide the IS-vs-host lineage verdict.
///
/// * `isosdb_nt`   — auto-derived IS-lineage nt records `(id, nt)`: the ISOSDB per-ORF
///                   CDS assembled during ingest.
/// * `isfinder_nt` — auto-derived ISfinder transposase CDS nt `(html.2020//<ISname>, nt)`,
///                   located in `IS.fna` by 6-frame translation (pure Rust, no external
///                   file) — see [`isfinder::isfinder_nt_pool`].
/// * `pos_extra`   — optional EXTRA IS-lineage nt FASTA (`--fpc-pos-extra`), appended to
///                   the `P|` pool if given (escape hatch; not required — both IS pools
///                   above are auto-derived).
/// * `neg_fna`     — the host-lineage nt FASTA (`--host` / `--fpc-neg`), the `N|` pool.
///
/// Returns `(n_pos, n_neg)` record counts written.
fn build_fpc(
    isosdb_nt: &[(String, Vec<u8>)],
    isfinder_nt: &[(String, Vec<u8>)],
    pos_extra: Option<&Path>,
    neg_fna: &Path,
    out_dir: &Path,
    work: &Path,
) -> Result<(usize, usize), String> {
    let fpc = out_dir.join("fpc");
    std::fs::create_dir_all(&fpc).map_err(|e| format!("mkdir {}: {e}", fpc.display()))?;
    let combined = work.join("fpc_refset.fna");
    let mut w = std::io::BufWriter::new(
        std::fs::File::create(&combined).map_err(|e| format!("create {}: {e}", combined.display()))?,
    );

    let mut write_rec = |prefix: char, id: &str, seq: &[u8]| -> Result<(), String> {
        w.write_all(format!(">{prefix}|{id}\n").as_bytes())
            .and_then(|_| w.write_all(seq))
            .and_then(|_| w.write_all(b"\n"))
            .map_err(|e| format!("write fpc fasta: {e}"))
    };

    // POSITIVE pool (IS-lineage nt), all prefixed P|:
    //   (a) ISOSDB per-ORF CDS (auto-derived during ingest),
    //   (b) ISfinder transposase CDS (auto-derived from IS.faa+IS.fna, pure Rust),
    //   (c) optional --fpc-pos-extra FASTA (escape hatch).
    let mut n_pos = 0usize;
    for (id, seq) in isosdb_nt {
        write_rec('P', id, seq)?;
        n_pos += 1;
    }
    for (id, seq) in isfinder_nt {
        write_rec('P', id, seq)?;
        n_pos += 1;
    }
    if let Some(px) = pos_extra {
        for (id, seq) in read_fasta(px)? {
            write_rec('P', &id, &seq)?;
            n_pos += 1;
        }
    }
    // NEGATIVE pool (host-lineage nt): the external --host FASTA, prefixed N|.
    let neg = read_fasta(neg_fna)?;
    let n_neg = neg.len();
    for (id, seq) in &neg {
        write_rec('N', id, seq)?;
    }
    w.flush().map_err(|e| format!("flush fpc fasta: {e}"))?;
    drop(w);

    let (fa_s, ref_s) = (
        combined.to_str().ok_or("non-utf8 path")?,
        fpc.join("refset").to_str().ok_or("non-utf8 path")?.to_string(),
    );
    mmseqs(&["createdb", fa_s, &ref_s, "--dbtype", "2", "-v", "0"])?;
    Ok((n_pos, n_neg))
}

// ---------------------------------------------------------------- orchestration
use std::collections::HashMap;

/// Inputs + outputs for a `build-db` run (local source paths; a later revision adds
/// pinned-SHA fetching that populates these).
pub struct BuildArgs {
    pub isosdb_fna: std::path::PathBuf,
    pub fam_annot: std::path::PathBuf,
    pub isfinder_faa: std::path::PathBuf,
    pub isfinder_csv: std::path::PathBuf,
    /// ISfinder element DNA (`IS.fna`). Used ONLY to auto-derive the ISfinder half of
    /// the FP-control `P|` pool (transposase CDS located by 6-frame translation). Not
    /// needed for the profile build; only required when building `fpc/refset` (--host).
    pub isfinder_fna: Option<std::path::PathBuf>,
    pub out_dir: std::path::PathBuf,
    pub work: std::path::PathBuf,
    pub strict_version: bool,
    /// Host-lineage NEGATIVE nt pool for the FP-control refset (`--host`/`--fpc-neg`).
    /// `None` → skip building `fpc/refset` (scanner FP-control then no-ops).
    pub fpc_neg: Option<std::path::PathBuf>,
    /// Optional EXTRA IS-lineage POSITIVE nt pool (`--fpc-pos-extra`), appended to the
    /// auto-derived ISOSDB + ISfinder CDS pools. Escape hatch only — NOT required.
    pub fpc_pos_extra: Option<std::path::PathBuf>,
}

/// Write `orfs` as a protein FASTA (one record per ORF, input order preserved).
fn write_union_faa(orfs: &[DbOrf], path: &Path) -> Result<(), String> {
    let mut w = std::io::BufWriter::new(
        std::fs::File::create(path).map_err(|e| format!("create {}: {e}", path.display()))?,
    );
    for o in orfs {
        w.write_all(format!(">{}\n", o.id).as_bytes())
            .and_then(|_| w.write_all(&o.seq))
            .and_then(|_| w.write_all(b"\n"))
            .map_err(|e| format!("write union faa: {e}"))?;
    }
    w.flush().map_err(|e| format!("flush union faa: {e}"))?;
    Ok(())
}

/// mmseqs createdb + cluster on `faa`, returning `rep_id -> [member_ids]` (rep
/// included). Deterministic given fixed input order (thread-invariant, verified).
fn cluster_union(faa: &Path, work: &Path) -> Result<HashMap<String, Vec<String>>, String> {
    let seq = work.join("seqDb");
    let clu = work.join("cluDb");
    let tmp = work.join("tmp_clu");
    let tsv = work.join("clusters.tsv");
    let (faa_s, seq_s, clu_s, tmp_s, tsv_s) = (
        faa.to_str().ok_or("non-utf8")?,
        seq.to_str().ok_or("non-utf8")?,
        clu.to_str().ok_or("non-utf8")?,
        tmp.to_str().ok_or("non-utf8")?,
        tsv.to_str().ok_or("non-utf8")?,
    );
    mmseqs(&["createdb", faa_s, seq_s, "-v", "0"])?;
    mmseqs(&[
        "cluster", seq_s, clu_s, tmp_s, "--min-seq-id", "0.4", "-c", "0.8",
        "--cluster-mode", "0", "-v", "0",
    ])?;
    mmseqs(&["createtsv", seq_s, seq_s, clu_s, tsv_s, "-v", "0"])?;

    let data = std::fs::read_to_string(&tsv).map_err(|e| format!("read {}: {e}", tsv.display()))?;
    let mut groups: HashMap<String, Vec<String>> = HashMap::new();
    for line in data.lines() {
        let mut it = line.split('\t');
        if let (Some(rep), Some(mem)) = (it.next(), it.next()) {
            groups.entry(rep.to_string()).or_default().push(mem.to_string());
        }
    }
    Ok(groups)
}

// ---------------------------------------------------------------- relabel
/// One manifest row: `(cid, family, source, rep_id, cluster_size)`.
type ManifestRow = (String, String, &'static str, String, usize);

/// Write selected ORFs to a protein FASTA under `>{name}`. `records` yields
/// `(name, seq_bytes)`; emitted in iterator order (deterministic).
fn write_faa<'a>(
    records: impl Iterator<Item = (&'a str, &'a [u8])>,
    path: &Path,
) -> Result<(), String> {
    let mut w = std::io::BufWriter::new(
        std::fs::File::create(path).map_err(|e| format!("create {}: {e}", path.display()))?,
    );
    for (name, seq) in records {
        w.write_all(format!(">{name}\n").as_bytes())
            .and_then(|_| w.write_all(seq))
            .and_then(|_| w.write_all(b"\n"))
            .map_err(|e| format!("write {}: {e}", path.display()))?;
    }
    w.flush().map_err(|e| format!("flush {}: {e}", path.display()))?;
    Ok(())
}

/// `mmseqs easy-search query target out.m8 tmp -s 7.5 -e <evalue>` (fresh tmp).
/// Mirrors `step9_build_full.sh`'s label-propagation search.
fn easy_search(query: &Path, target: &Path, out_m8: &Path, tmp: &Path, evalue: &str) -> Result<(), String> {
    let _ = std::fs::remove_dir_all(tmp);
    std::fs::create_dir_all(tmp).map_err(|e| format!("mkdir {}: {e}", tmp.display()))?;
    let (q, t, o, tmps) = (
        query.to_str().ok_or("non-utf8")?,
        target.to_str().ok_or("non-utf8")?,
        out_m8.to_str().ok_or("non-utf8")?,
        tmp.to_str().ok_or("non-utf8")?,
    );
    mmseqs(&["easy-search", q, t, o, tmps, "-s", "7.5", "-e", evalue, "-v", "0"])
}

/// Parse a BLAST-tab (m8) file → best (lowest-evalue) target per query.
/// m8 columns: query target pident aln mism gap qs qe ts te **evalue** bits (evalue = col 11).
fn best_hits(m8: &Path) -> Result<HashMap<String, (f64, String)>, String> {
    let data = std::fs::read_to_string(m8).map_err(|e| format!("read {}: {e}", m8.display()))?;
    let mut best: HashMap<String, (f64, String)> = HashMap::new();
    for line in data.lines() {
        let p: Vec<&str> = line.split('\t').collect();
        if p.len() < 11 {
            continue;
        }
        let ev: f64 = p[10].parse().unwrap_or(f64::INFINITY);
        match best.get(p[0]) {
            Some((be, _)) if *be <= ev => {}
            _ => {
                best.insert(p[0].to_string(), (ev, p[1].to_string()));
            }
        }
    }
    Ok(best)
}

/// Relabel cluster reps to lift family accuracy, faithfully porting
/// `step9_label_assign.py` (unknown propagation, E≤1e-3) then
/// `step10_assign_isf.py` (ISfinder authority, E≤1e-5). Mutates `manifest`:
///   1. reps whose family is `unknown` inherit the family of their best hit
///      against the family-labeled **ISOSDB** ORFs (E≤1e-3);
///   2. any rep hitting **ISfinder** at E≤1e-5 takes that ISfinder family and
///      `source="ISfinder"`; every other rep gets `source="ISOSDB"`.
/// Returns `(propagated, isfinder_authoritative, isfinder_changed)` counts.
fn relabel_manifest(
    manifest: &mut [ManifestRow],
    union: &[DbOrf],
    work: &Path,
) -> Result<(usize, usize, usize), String> {
    let seq_by_id: HashMap<&str, &[u8]> =
        union.iter().map(|o| (o.id.as_str(), o.seq.as_slice())).collect();
    let fam_by_id: HashMap<&str, &str> =
        union.iter().map(|o| (o.id.as_str(), o.family.as_str())).collect();

    // ---- Phase 1: unknown-family propagation (E≤1e-3) --------------------
    // query = reps of `unknown` clusters (named by cid); target = family-labeled
    // ISOSDB ORFs (byfam/*.faa is ISOSDB-only — verified: 36449 ids, 0 ISfinder).
    let unk_q = work.join("unknown_reps.faa");
    write_faa(
        manifest.iter().filter(|r| r.1 == "unknown").filter_map(|r| {
            seq_by_id.get(r.3.as_str()).map(|s| (r.0.as_str(), *s))
        }),
        &unk_q,
    )?;
    let lab_t = work.join("labeled_orfs.faa");
    write_faa(
        union
            .iter()
            .filter(|o| o.source == "ISOSDB" && o.family != "unknown")
            .map(|o| (o.id.as_str(), o.seq.as_slice())),
        &lab_t,
    )?;
    let unk_m8 = work.join("unknown_search.m8");
    easy_search(&unk_q, &lab_t, &unk_m8, &work.join("lp_tmp"), "1e-3")?;
    let unk_best = best_hits(&unk_m8)?;
    let mut propagated = 0;
    for row in manifest.iter_mut() {
        if row.1 == "unknown" {
            if let Some((ev, tgt)) = unk_best.get(&row.0) {
                if *ev <= 1e-3 {
                    if let Some(f) = fam_by_id.get(tgt.as_str()) {
                        row.1 = f.to_string();
                        propagated += 1;
                    }
                }
            }
        }
    }

    // ---- Phase 2: ISfinder authority (E≤1e-5) ---------------------------
    // query = ALL reps (by cid); target = ISfinder ORFs. A hit ⇒ ISfinder is
    // authoritative for both family and provenance.
    let all_q = work.join("all_reps.faa");
    write_faa(
        manifest.iter().filter_map(|r| {
            seq_by_id.get(r.3.as_str()).map(|s| (r.0.as_str(), *s))
        }),
        &all_q,
    )?;
    let isf_t = work.join("isfinder_orfs.faa");
    write_faa(
        union
            .iter()
            .filter(|o| o.source == "ISfinder")
            .map(|o| (o.id.as_str(), o.seq.as_slice())),
        &isf_t,
    )?;
    let isf_m8 = work.join("reps_vs_isfinder.m8");
    easy_search(&all_q, &isf_t, &isf_m8, &work.join("isf_tmp"), "1e-5")?;
    let isf_best = best_hits(&isf_m8)?;
    let (mut n_isf, mut changed) = (0usize, 0usize);
    for row in manifest.iter_mut() {
        match isf_best.get(&row.0) {
            Some((ev, tgt)) if *ev <= 1e-5 => {
                if let Some(f) = fam_by_id.get(tgt.as_str()) {
                    if *f != row.1 {
                        changed += 1;
                    }
                    row.1 = f.to_string();
                    row.2 = "ISfinder";
                    n_isf += 1;
                } else {
                    row.2 = "ISOSDB";
                }
            }
            _ => row.2 = "ISOSDB",
        }
    }
    Ok((propagated, n_isf, changed))
}

/// Run the whole union-IS-DB build. Produces `<out>/mmdb_union/profileDb*` and
/// `<out>/manifest_union.tsv`, with families/provenance set by the two-phase
/// ISfinder-authority relabel ([`relabel_manifest`]). When `args.fpc_neg` is set it
/// additionally emits `<out>/fpc/refset*` (the nt neg-pos FP-control reference).
pub fn run(args: &BuildArgs) -> Result<(), String> {
    let ver = check_mmseqs_version(args.strict_version)?;
    eprintln!("[build-db] mmseqs {ver}");
    std::fs::create_dir_all(&args.work).map_err(|e| format!("mkdir work: {e}"))?;
    std::fs::create_dir_all(args.out_dir.join("mmdb_union"))
        .map_err(|e| format!("mkdir out: {e}"))?;

    // 1. Ingest both sources → union ORFs (ISOSDB then ISfinder, deterministic order).
    //    ISOSDB ingest also returns the per-ORF nucleotide CDS (IS-lineage FP-control
    //    positive pool), keyed by the same `ISOSDB<n>_orf<k>` id.
    eprintln!("[build-db] ingesting ISOSDB…");
    let (mut union, isosdb_nt) = isosdb::ingest_isosdb(&args.isosdb_fna, &args.fam_annot)?;
    let n_iso = union.len();
    eprintln!("[build-db] ingesting ISfinder…");
    union.extend(isfinder::ingest_isfinder(&args.isfinder_faa, &args.isfinder_csv)?);
    eprintln!("[build-db] union ORFs: {} ({n_iso} ISOSDB + {} ISfinder)", union.len(), union.len() - n_iso);

    let by_id: HashMap<&str, &DbOrf> = union.iter().map(|o| (o.id.as_str(), o)).collect();
    let union_faa = args.work.join("union_orfs.faa");
    write_union_faa(&union, &union_faa)?;

    // 2. Cluster.
    eprintln!("[build-db] clustering (mmseqs, 0.4/0.8)…");
    let groups = cluster_union(&union_faa, &args.work)?;
    eprintln!("[build-db] clusters: {}", groups.len());

    // 3. Assemble clusters (sorted by rep for determinism; rep first in each), assign
    //    synthetic cids, and record manifest rows.
    let mut reps: Vec<&String> = groups.keys().collect();
    reps.sort();
    let mut clusters: Vec<(String, Vec<(String, Vec<u8>)>)> = Vec::with_capacity(reps.len());
    let mut manifest: Vec<(String, String, &'static str, String, usize)> = Vec::with_capacity(reps.len());
    for (i, rep) in reps.iter().enumerate() {
        let cid = format!("c{i:05}");
        let members = &groups[*rep];
        let mut ordered: Vec<&String> = std::iter::once(*rep)
            .chain(members.iter().filter(|m| *m != *rep))
            .collect();
        // stable: keep rep first, sort the rest
        ordered[1..].sort();
        let seqs: Vec<(String, Vec<u8>)> = ordered
            .iter()
            .filter_map(|m| by_id.get(m.as_str()).map(|o| (o.id.clone(), o.seq.clone())))
            .collect();
        let rep_orf = by_id.get(rep.as_str()).ok_or("rep not in union")?;
        manifest.push((cid.clone(), rep_orf.family.clone(), rep_orf.source, (*rep).clone(), seqs.len()));
        clusters.push((cid, seqs));
    }

    // 4. Relabel: unknown-propagation (E≤1e-3) then ISfinder authority (E≤1e-5).
    eprintln!("[build-db] relabeling (unknown propagation + ISfinder authority)…");
    let (propagated, n_isf, changed) = relabel_manifest(&mut manifest, &union, &args.work)?;
    eprintln!(
        "[build-db] relabel: {propagated} unknowns propagated; \
         {n_isf} ISfinder-authoritative ({changed} families changed)"
    );

    // 5. poasta MSA per cluster → union profile DB.
    eprintln!("[build-db] building profiles (poasta MSA → msa2profile)…");
    let prof = args.out_dir.join("mmdb_union").join("profileDb");
    build_profiledb_poasta(&clusters, &prof, &args.work.join("prof_work"))?;

    // 6. manifest_union.tsv.
    let man_path = args.out_dir.join("manifest_union.tsv");
    let mut mw = std::io::BufWriter::new(
        std::fs::File::create(&man_path).map_err(|e| format!("create manifest: {e}"))?,
    );
    mw.write_all(b"cid\tfamily\tsource\trep\tsize\n").map_err(|e| format!("write: {e}"))?;
    for (cid, fam, src, rep, size) in &manifest {
        mw.write_all(format!("{cid}\t{fam}\t{src}\t{rep}\t{size}\n").as_bytes())
            .map_err(|e| format!("write manifest: {e}"))?;
    }
    mw.flush().map_err(|e| format!("flush manifest: {e}"))?;

    // 7. FP-control reference DB (module C): fpc/refset. Requires a host-lineage
    //    negative pool (--host / --fpc-neg); skipped (with a note) otherwise.
    match &args.fpc_neg {
        Some(neg) => {
            // Auto-derive the ISfinder half of the P| pool from IS.faa + IS.fna (pure
            // Rust 6-frame location), mirroring the auto-derived ISOSDB per-ORF CDS.
            let isfinder_nt = match &args.isfinder_fna {
                Some(fna) => {
                    eprintln!("[build-db] deriving ISfinder nt CDS (6-frame location vs IS.fna)…");
                    let (recs, st) = isfinder::isfinder_nt_pool(&args.isfinder_faa, fna)?;
                    eprintln!(
                        "[build-db] ISfinder nt: {} located ({} primary + {} recovered) \
                         of {} ids (no_nt={} not_located={})",
                        recs.len(), st.primary, st.recovered, st.candidates, st.no_nt, st.not_located
                    );
                    recs
                }
                None => {
                    eprintln!(
                        "[build-db] NOTE: --isfinder-fna not given → ISfinder half of the \
                         P| pool is NOT auto-derived (only ISOSDB CDS + any --fpc-pos-extra)."
                    );
                    Vec::new()
                }
            };
            eprintln!("[build-db] building fpc/refset (nt neg-pos FP-control reference)…");
            let (n_pos, n_neg) = build_fpc(
                &isosdb_nt,
                &isfinder_nt,
                args.fpc_pos_extra.as_deref(),
                neg,
                &args.out_dir,
                &args.work,
            )?;
            eprintln!(
                "[build-db] fpc/refset: {n_pos} IS (P|) + {n_neg} host (N|) refs → {}",
                args.out_dir.join("fpc").join("refset").display()
            );
        }
        None => eprintln!(
            "[build-db] NOTE: no --host given → skipping fpc/refset; the scanner's \
             FP-control (host-lineage discrimination) will no-op for this DB."
        ),
    }

    eprintln!(
        "[build-db] done → {} ({} profiles)",
        args.out_dir.display(),
        clusters.len()
    );
    Ok(())
}

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

    #[test]
    fn poa_msa_aligns_equal_length_rows_rep_first() {
        // Three related short peptides; rep first. All output rows must share the
        // MSA width and preserve input order (rep first).
        let seqs = vec![
            ("rep".to_string(), b"MKTAYIAKQR".to_vec()),
            ("m2".to_string(), b"MKTAYIAKR".to_vec()),
            ("m3".to_string(), b"MKTAYYIAKQR".to_vec()),
        ];
        let msa = poa_msa(&seqs).expect("poa_msa ok");
        assert_eq!(msa.len(), 3, "one row per input sequence");
        assert_eq!(msa[0].0, "rep", "rep stays first");
        let w = msa[0].1.len();
        assert!(w >= 10, "aligned width covers the longest seq");
        for (n, s) in &msa {
            assert_eq!(s.len(), w, "row {n} padded to MSA width");
        }
        // Ungapped residues of each row reconstruct the input sequence.
        let deg: String = msa[0].1.chars().filter(|c| *c != '-').collect();
        assert_eq!(deg, "MKTAYIAKQR");
    }

    #[test]
    fn poa_msa_singleton_is_ungapped() {
        let seqs = vec![("solo".to_string(), b"MKQAR".to_vec())];
        let msa = poa_msa(&seqs).expect("ok");
        assert_eq!(msa, vec![("solo".to_string(), "MKQAR".to_string())]);
    }

    // Local ingest validation against the pinned sources (skipped when absent, so CI-safe).
    // Reference counts (build_db.py §0h): ISfinder 7,121 transposase ORFs, ISOSDB 36,449.
    #[test]
    fn ingest_counts_match_reference() {
        use std::path::PathBuf;
        // Dev-only integration check — gated on env-supplied paths so no developer
        // absolute path is baked into the published source:
        //   RUST_ISE_TEST_ISFINDER_DIR → dir holding IS.faa, IS.csv, IS_fam_annot.txt
        //   RUST_ISE_TEST_ISOSDB_FNA   → path to the unzipped ISOSDB.V3.fna
        // Unset → skip (CI-safe).
        let dir = match std::env::var_os("RUST_ISE_TEST_ISFINDER_DIR") {
            Some(d) => PathBuf::from(d),
            None => { eprintln!("skip: RUST_ISE_TEST_ISFINDER_DIR unset"); return; }
        };
        let iso_fna = match std::env::var_os("RUST_ISE_TEST_ISOSDB_FNA") {
            Some(f) => PathBuf::from(f),
            None => { eprintln!("skip: RUST_ISE_TEST_ISOSDB_FNA unset"); return; }
        };
        let isf_faa = dir.join("IS.faa");
        let isf_csv = dir.join("IS.csv");
        let fam_annot = dir.join("IS_fam_annot.txt");
        if !isf_faa.exists() || !iso_fna.exists() {
            eprintln!("skip: ingest source files not present");
            return;
        }
        let isf = super::isfinder::ingest_isfinder(&isf_faa, &isf_csv).expect("isfinder ingest");
        eprintln!("ISfinder transposase ORFs: {}", isf.len());
        assert_eq!(isf.len(), 7121, "ISfinder transposase count == reference");
        assert!(isf.iter().all(|o| o.source == "ISfinder" && !o.seq.is_empty()));
        let n_unknown = isf.iter().filter(|o| o.family == "unknown").count();
        eprintln!("  ISfinder unknown-family: {n_unknown}");

        let (iso, iso_nt) = super::isosdb::ingest_isosdb(&iso_fna, &fam_annot).expect("isosdb ingest");
        eprintln!("ISOSDB ORFs: {} (nt CDS: {})", iso.len(), iso_nt.len());
        assert!(iso.iter().all(|o| o.source == "ISOSDB" && !o.seq.is_empty()));
        // Every kept protein ORF should have a matching nt CDS record (same ORF set).
        assert_eq!(iso_nt.len(), iso.len(), "nt-CDS pool tracks the kept ORF set");
        assert!(iso_nt.iter().all(|(_, nt)| !nt.is_empty()));
        // Reference produced 36,449; allow a tiny tolerance for rustygal vs pyrodigal edge ORFs.
        assert!(
            (36449i64 - iso.len() as i64).abs() <= 200,
            "ISOSDB ORF count {} within tolerance of reference 36449",
            iso.len()
        );
    }

    #[test]
    fn stockholm_block_has_id_and_terminator() {
        let block = stockholm_block("c001", &[("a".into(), "MK-Q".into()), ("b".into(), "MKAQ".into())]);
        assert!(block.starts_with("# STOCKHOLM 1.0\n#=GF ID c001\n"));
        assert!(block.trim_end().ends_with("//"));
    }
}