SimurgArray/src/simurg/pipeline/vcf.py

"""``simurg gtc-to-vcf``: GTC + manifest + FASTA → one bgzipped, indexed SNV VCF per sample (spec CLI-022,
OUT-010 … OUT-042, IN-060 … IN-064; plan §8).

The manifest is mapped to plus-strand sites once (``simurg.harmonize``); each sample is then written by a worker
process with vectorised genotype lookup for single-probe sites and the OUT-031 algorithm for merged sites.
"""

from __future__ import annotations

import hashlib
import math
import multiprocessing as mp
import os
from collections import Counter
from concurrent.futures import ProcessPoolExecutor
from dataclasses import dataclass, field
from pathlib import Path
from typing import Any

import numpy as np
import numpy.typing as npt

from simurg import PRODUCT_NAME, __version__
from simurg.core.atomic import atomic_open, check_outputs_absent
from simurg.core.errors import ExitCode, SimurgError
from simurg.core.log import Reporter
from simurg.core.provenance import RunManifest, sha256_file
from simurg.formats.bgzf import BgzfWriter, TabixBuilder, write_tbi
from simurg.harmonize.combine import genotype as combine_genotype
from simurg.harmonize.combine import gt_text
from simurg.harmonize.sites import AuditRow, ManifestView, Site, group_sites, map_locus

PAR = {  # 1-based inclusive pseudo-autosomal intervals on X (IN-065, human)
    "GRCh38": ((10_001, 2_781_479), (155_701_383, 156_030_895)),
    "GRCh37": ((60_001, 2_699_520), (154_931_044, 155_260_560)),
}
CLASS_AUTO, CLASS_X, CLASS_Y, CLASS_MT, CLASS_PAR = 0, 1, 2, 3, 4


@dataclass
class VcfOptions:
    out: Path
    bpm: Path
    fasta: Path
    gtcs: Path | None = None
    sheet: Path | None = None
    csv: Path | None = None
    bgzip: bool = True
    genome_cache: bool = True
    exclude_loci: Path | None = None
    include_loci: Path | None = None
    unsquash_duplicates: bool = False
    inf1_nocall_informative: bool = False
    drop_ambiguous: bool = False
    filter_gs: float | None = None
    contig_alias: str = "strict"
    ploidy_policy: str = "sex-aware"
    threads: int = 0
    merge: bool = False
    merge_snv_indel: bool = False
    auxiliary_loci: Path | None = None
    info_fields: list[str] = field(default_factory=list)
    extra_format: list[str] = field(default_factory=list)
    egt: Path | None = None
    force: bool = False
    argv: list[str] = field(default_factory=list)


# ------------------------------------------------------------------------------------------------ site table
@dataclass
class SiteTable:
    """Everything a worker needs, precomputed once (picklable)."""

    prefix: list[bytes]  # "CHROM\tPOS\tID\tREF\tALT\t.\t" per site
    single_rows: npt.NDArray[np.int64]  # site indices of single-probe sites
    single_locus: npt.NDArray[np.int64]  # their manifest locus index
    single_gt: npt.NDArray[np.object_]  # (n_single, 2 ploidy, 4 calls) GT text
    single_flip: npt.NDArray[np.bool_]  # BAF must be 1 − BAF
    merged_sites: list[Site]
    merged_rows: npt.NDArray[np.int64]
    site_class: npt.NDArray[np.int8]
    contigs: list[str]
    header: str
    site_beg: npt.NDArray[np.int64]
    site_end: npt.NDArray[np.int64]
    site_contig: list[str]
    site_id: list[str] = field(default_factory=list)
    site_alleles: list[list[str]] = field(default_factory=list)
    single_ab: npt.NDArray[np.int16] | None = None  # (n_single, 2) allele indices of the probe's A and B alleles
    site_info: list[str] = field(default_factory=list)  # INFO text per site ('.' without --info-fields)
    site_info_typed: list[list[tuple[str, str, list[Any]]]] = field(default_factory=list)  # (key, kind, values)
    extra_format: list[str] = field(default_factory=list)
    single_lookup: npt.NDArray[np.int32] | None = None  # normalization-transform index of single-probe loci
    inf1_nc_informative: bool = False
    ploidy_policy: str = "sex-aware"
    filter_gs: float | None = None


INFO_DEFS = {  # OUT-012, one value per contributing probe (BPM order)
    "ASSAY_TYPE": ("Integer", "Assay type per probe: 0 Infinium II, 1 Infinium I red, 2 Infinium I green"),
    "NORM_BIN": ("Integer", "Normalization bin per probe"),
    "GENTRAIN": ("Float", "GenTrain score per probe (cluster file)"),
    "CLUSTER_SEP": ("Float", "Cluster separation score per probe (cluster file)"),
    "ILMN_STRAND": ("String", "Design strand per probe (TOP, BOT, PLUS, MINUS)"),
    "PROBE_IDS": ("String", "Names of the contributing probes"),
}
EXTRA_DEFS = {  # plan §8.3; single-probe records only
    "X": ("Integer", "Raw intensity, X (red) channel"),
    "Y": ("Integer", "Raw intensity, Y (green) channel"),
    "NORMX": ("Float", "Normalized X intensity"),
    "NORMY": ("Float", "Normalized Y intensity"),
    "THETA": ("Float", "Normalized theta = (2/pi) atan2(NORMY, NORMX)"),
    "R": ("Float", "Normalized R = NORMX + NORMY"),
}


def _site_class(contig: str, pos: int, par_flag: bool, build: str) -> int:
    c = contig[3:] if contig.lower().startswith("chr") else contig
    if c in ("MT", "M"):
        return CLASS_MT
    if c == "Y":
        return CLASS_Y
    if c == "X":
        if par_flag or any(a <= pos <= b for a, b in PAR.get(build or "GRCh38", ())):
            return CLASS_PAR
        return CLASS_X
    return CLASS_AUTO


def ploidy_for(cls: npt.NDArray[np.int8], sex: str, policy: str) -> npt.NDArray[np.int8]:
    """OUT-018 per site; 0 marks 'haploid and always missing' (Y of female samples)."""
    p = np.full(len(cls), 2, dtype=np.int8)
    p[cls == CLASS_MT] = 1
    if policy == "diploid":
        return p
    p[cls == CLASS_Y] = 0 if sex == "F" else 1
    if sex == "M":
        p[cls == CLASS_X] = 1
    return p


def read_locus_list(path: Path) -> list[str]:
    out = []
    for line in path.read_text(encoding="utf-8").splitlines():
        line = line.split("#", 1)[0].strip()
        if line:
            out.append(line)
    return out


def build_sites(
    opts: VcfOptions, rep: Reporter, rm: RunManifest
) -> tuple[SiteTable, list[AuditRow], Counter[str], Any]:
    from simurg.formats.bpm import read_bpm
    from simurg.formats.fasta import ContigMap, Genome

    bpm = read_bpm(opts.bpm)
    csvm = None
    if opts.csv is not None:
        from simurg.formats.csv_manifest import compare_manifests, read_csv_manifest

        csvm = read_csv_manifest(opts.csv)
        problems = compare_manifests(bpm, csvm)
        if problems:
            raise SimurgError("E_BPM_CSV_MISMATCH", "BPM and CSV manifests differ:\n  " + "\n  ".join(problems))
    genome = Genome(opts.fasta, use_cache=opts.genome_cache)
    for n in genome.notes:
        rep.info(n)
    builds = Counter(bpm.genome_builds)
    mbuild = builds.most_common(1)[0][0] if builds else ""
    want = {"37": "GRCh37", "38": "GRCh38"}.get(mbuild.strip(), mbuild.strip())
    if genome.build and want and want != genome.build:
        raise SimurgError(
            "E_BUILD_MISMATCH", f"manifest genome build {want} but FASTA looks like {genome.build} ({opts.fasta})"
        )
    build = genome.build or want
    src = csvm or bpm
    mv = ManifestView(
        bpm.names,
        bpm.snps,
        bpm.chroms,
        bpm.map_infos.tolist(),
        src.ilmn_strands,
        src.source_strands,
        src.ref_strands if any(src.ref_strands) else bpm.ref_strands,
        csvm.source_seqs if csvm else [""] * bpm.num_loci,
        bpm.assay_types.tolist(),
        bpm.intensity_only.tolist(),
    )
    if csvm is None and any(("D" in s) for s in bpm.snps):
        rep.error("W_SKIP_INDEL_CONTEXT", "no CSV manifest: indels are skipped (give --csv with SourceSeq)")
    names = set(bpm.names)
    excl: set[str] = set()
    incl: set[str] | None = None
    for path, kind in ((opts.exclude_loci, "exclude"), (opts.include_loci, "include")):
        if path is None:
            continue
        lst = read_locus_list(path)
        unknown = [x for x in lst if x not in names]
        if unknown:
            rep.warn("W_UNKNOWN_LOCI", f"{len(unknown)} name(s) in --{kind}-loci are not in the manifest")
        if kind == "exclude":
            excl = set(lst)
        else:
            incl = set(lst)
    cmap = ContigMap([e.name for e in genome.entries], opts.contig_alias)
    probes = []
    audit: list[AuditRow] = []
    skips: Counter[str] = Counter()
    for i in range(bpm.num_loci):
        nm = bpm.names[i]
        if nm in excl or (incl is not None and nm not in incl):
            a = AuditRow(nm, bpm.snps[i], "", "", "", "", "", False, False, "W_SKIP_EXCLUDED", build, "")
            p = None
        else:
            p, a = map_locus(i, mv, genome.fetch, cmap.get, build)
            if p is not None and opts.drop_ambiguous and p.confidence == "low":
                a.reason, p = "W_SKIP_AMBIGUOUS", None
        audit.append(a)
        if p is None or a.reason.startswith("W_"):
            skips[a.reason] += 1
        if p is not None:
            probes.append(p)
    for code, count in sorted(skips.items()):
        rep.warn(code, f"{count} locus/loci: {code}")
    if opts.auxiliary_loci is not None:
        from simurg.harmonize.sites import apply_auxiliary, read_auxiliary

        aux = read_auxiliary(str(opts.auxiliary_loci))
        strands = mv.ref_strands
        probes, errors = apply_auxiliary(probes, aux, genome.fetch, cmap.get, lambda i: (strands[i] or "").strip())
        if errors:
            raise SimurgError("E_AUX_INVALID", f"{len(errors)} invalid auxiliary record(s):\n  " + "\n  ".join(errors))
        rep.info(f"{len(aux)} auxiliary locus definition(s) applied")
        names = {r.name for r in aux}
        for a_row in audit:
            if a_row.locus in names:
                a_row.reason = "aux_locus"
    order = {e.name: k for k, e in enumerate(genome.entries)}
    sites = group_sites(probes, order, unsquash=opts.unsquash_duplicates, merge_snv_indel=opts.merge_snv_indel)
    used = {s.contig for s in sites}
    contigs = [e.name for e in genome.entries if e.name in used]
    # static per-site text and single-probe genotype tables
    prefix, cls = [], np.zeros(len(sites), np.int8)
    single_rows, single_locus, single_flip, gts, single_ab = [], [], [], [], []
    merged_sites, merged_rows = [], []
    beg = np.zeros(len(sites), np.int64)
    end = np.zeros(len(sites), np.int64)
    for k, s in enumerate(sites):
        prefix.append(f"{s.contig}\t{s.pos}\t{s.id}\t{s.ref}\t{','.join(s.alts) if s.alts else '.'}\t.\t".encode())
        cls[k] = _site_class(s.contig, s.pos, any(p.par for p in s.probes), build)
        beg[k], end[k] = s.pos - 1, s.pos - 1 + len(s.ref)
        if len(s.probes) == 1:
            from simurg.harmonize.combine import single

            single_rows.append(k)
            single_locus.append(s.probes[0].index)
            al_a, al_b = s.probes[0].alleles
            single_flip.append(al_a != s.ref and al_b == s.ref)
            single_ab.append((s.allele_index(al_a), s.allele_index(al_b)))
            gts.append([[gt_text(single(s, c, pl), pl) for c in range(4)] for pl in (1, 2)])
        else:
            merged_sites.append(s)
            merged_rows.append(k)
    site_info, site_info_typed = _site_info(opts, bpm, sites)
    header = vcf_header(opts, genome, build, contigs, rm)
    table = SiteTable(
        prefix=prefix,
        single_rows=np.array(single_rows, np.int64),
        single_locus=np.array(single_locus, np.int64),
        single_gt=np.array(gts, dtype=object).reshape(len(single_rows), 2, 4),
        single_flip=np.array(single_flip, np.bool_),
        merged_sites=merged_sites,
        merged_rows=np.array(merged_rows, np.int64),
        site_class=cls,
        contigs=contigs,
        header=header,
        site_beg=beg,
        site_end=end,
        site_contig=[s.contig for s in sites],
        site_id=[s.id for s in sites],
        site_alleles=[[s.ref, *s.alts] for s in sites],
        single_ab=np.array(single_ab, np.int16).reshape(-1, 2),
        site_info=site_info,
        site_info_typed=site_info_typed,
        extra_format=list(opts.extra_format),
        single_lookup=bpm.normalization_lookups()[np.array(single_locus, np.int64)] if opts.extra_format else None,
        inf1_nc_informative=opts.inf1_nocall_informative,
        ploidy_policy=opts.ploidy_policy,
        filter_gs=opts.filter_gs,
    )
    return table, audit, skips, genome


def _sha(path: Path | None) -> str:
    return sha256_file(path) if path else ""


def _site_info(
    opts: VcfOptions, bpm: Any, sites: list[Site]
) -> tuple[list[str], list[list[tuple[str, str, list[Any]]]]]:
    """Static INFO per site for --info-fields (OUT-012)."""
    if not opts.info_fields:
        return ["."] * len(sites), [[] for _ in sites]
    gentrain = sep = None
    if {"GENTRAIN", "CLUSTER_SEP"} & set(opts.info_fields):
        from simurg.formats.egt import read_egt

        assert opts.egt is not None
        egt = read_egt(opts.egt)
        eidx = egt.name_index()
        rows = np.array([eidx.get(n, -1) for n in bpm.names], np.int64)
        gentrain = np.where(rows >= 0, egt.gentrain_score[np.maximum(rows, 0)], np.nan)
        sep = np.where(rows >= 0, egt.cluster_sep[np.maximum(rows, 0)], np.nan)
    bins = bpm.norm_bin_ids()
    texts, typed = [], []
    for s in sites:
        idx = [p.index for p in s.probes]
        vals: dict[str, list[Any]] = {
            "ASSAY_TYPE": [int(bpm.assay_types[i]) for i in idx],
            "NORM_BIN": [int(bins[i]) for i in idx],
            "GENTRAIN": [float(gentrain[i]) for i in idx] if gentrain is not None else [],
            "CLUSTER_SEP": [float(sep[i]) for i in idx] if sep is not None else [],
            "ILMN_STRAND": [bpm.ilmn_strands[i] or "." for i in idx],
            "PROBE_IDS": [p.name for p in s.probes],
        }
        parts, tv = [], []
        for key in opts.info_fields:
            v = vals[key]
            kind = INFO_DEFS[key][0]
            text = ",".join(
                "." if (isinstance(x, float) and math.isnan(x)) else (f"{x:.6g}" if isinstance(x, float) else str(x))
                for x in v
            )
            parts.append(f"{key}={text}")
            tv.append((key, kind, v))
        texts.append(";".join(parts))
        typed.append(tv)
    return texts, typed


def vcf_header(opts: VcfOptions, genome: Any, build: str, contigs: list[str], rm: RunManifest) -> str:
    """OUT-011 (without the ##SAMPLE and #CHROM lines, which are per sample)."""

    def lst(p: Path | None) -> str:
        return "none" if p is None else f"{p.name}:{hashlib.sha256(p.read_bytes()).hexdigest()[:16]}"

    lines = [
        "##fileformat=VCFv4.3",
        f"##source={PRODUCT_NAME} {__version__}",
        f"##reference=file://{Path(opts.fasta).resolve()}",
        f"##simurg_genome_build={build}",
        f"##simurg_manifest={opts.bpm.name};sha256={_sha(opts.bpm)}",
    ]
    if opts.csv:
        lines.append(f"##simurg_csv_manifest={opts.csv.name};sha256={_sha(opts.csv)}")
    lines.append("##simurg_cluster_file=__CLUSTER__")
    lines.append(
        "##simurg_options="
        f"unsquash_duplicates={str(opts.unsquash_duplicates).lower()};exclude_loci={lst(opts.exclude_loci)};"
        f"include_loci={lst(opts.include_loci)};inf1_nocall_informative={str(opts.inf1_nocall_informative).lower()};"
        f"auxiliary_loci={lst(opts.auxiliary_loci)};ploidy_policy={opts.ploidy_policy};contig_alias={opts.contig_alias}"
        + (";drop_ambiguous=true" if opts.drop_ambiguous else "")
        + (";merge_snv_indel=true" if opts.merge_snv_indel else "")
        + (";unknown_sex_policy=diploid_x" if opts.ploidy_policy == "sex-aware" else "")
    )
    lines.append(f"##simurg_run={rm.run_id}")
    for c in contigs:
        lines.append(f"##contig=<ID={c},length={genome.contig_length(c)}>")
    lines.append('##FILTER=<ID=PASS,Description="All filters passed">')
    for key in opts.info_fields:
        kind, desc = INFO_DEFS[key]
        lines.append(f'##INFO=<ID={key},Number=.,Type={kind},Description="{desc}">')
    if opts.filter_gs is not None:
        lines.append(f'##FILTER=<ID=LowGS,Description="GenCall score below {opts.filter_gs}">')
    lines += [
        '##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">',
        '##FORMAT=<ID=GS,Number=1,Type=Float,Description="GenCall score. For records merged from several assays, '
        'the minimum GenCall score over the contributing assays">',
        '##FORMAT=<ID=BAF,Number=1,Type=Float,Description="B allele frequency, oriented so that 1 means homozygous '
        'ALT (REF=A, ALT=B relative to the reference genome)">',
        '##FORMAT=<ID=LRR,Number=1,Type=Float,Description="Log R ratio">',
    ]
    for key in opts.extra_format:
        kind, desc = EXTRA_DEFS[key]
        lines.append(f'##FORMAT=<ID={key},Number=1,Type={kind},Description="{desc}">')
    return "\n".join(lines) + "\n"


# ------------------------------------------------------------------------------------------------ per sample
_W: dict[str, Any] = {}


def _init(table: SiteTable) -> None:
    _W["t"] = table


def _f4(v: float) -> str:
    return "." if math.isnan(v) else f"{v:.4f}"


def _g8(v: float) -> str:
    return "." if math.isnan(v) else f"{v:.8g}"


def write_sample(task: tuple[str, str, str, str, bool]) -> dict[str, Any]:
    """Worker: one GTC → one VCF (+ index). Returns counts for the run summary."""
    from simurg.formats.gtc import read_gtc

    sample_id, sample_name, gtc_path, out_path, bgzip = task
    t: SiteTable = _W["t"]
    res: dict[str, Any] = {"sample_id": sample_id, "records": 0, "reasons": Counter(), "merge_failures": []}
    try:
        g = read_gtc(gtc_path)
        n = len(t.prefix)
        sex = g.gender_char
        ploidy = ploidy_for(t.site_class, sex, t.ploidy_policy)
        gt = np.empty(n, dtype=object)
        gs = np.full(n, np.nan, np.float64)
        baf = np.full(n, np.nan, np.float64)
        lrr = np.full(n, np.nan, np.float64)
        has_baf = len(g.b_allele_freqs) == g.num_snps
        # single-probe sites (vectorised)
        rows, loc = t.single_rows, t.single_locus
        calls = g.genotypes[loc].astype(np.int64)
        pl = ploidy[rows]
        pidx = np.where(pl == 2, 1, 0)
        gt[rows] = t.single_gt[np.arange(len(rows)), pidx, calls]
        gt[rows[pl == 0]] = "."
        gs[rows] = g.genotype_scores[loc].astype(np.float64)
        if has_baf:
            b = g.b_allele_freqs[loc].astype(np.float32)
            b = np.where(t.single_flip, (np.float32(1.0) - b).astype(np.float32), b)
            baf[rows] = b.astype(np.float64)
            lrr[rows] = g.logr_ratios[loc].astype(np.float64)
        extra_cols: list[list[str]] = []
        if t.extra_format:
            from simurg.signal.normalization.apply import normalize, polar

            rx, ry = g.raw_x[loc], g.raw_y[loc]
            nx, ny = normalize(rx, ry, t.single_lookup, g.normalization_transforms)  # type: ignore[arg-type]
            zero = (rx == 0) & (ry == 0)  # no intensity: normalized values are missing (as in the GTC semantics)
            nx = np.where(zero, np.float32(np.nan), nx).astype(np.float32)
            ny = np.where(zero, np.float32(np.nan), ny).astype(np.float32)
            th, rr = polar(nx, ny)
            src: dict[str, npt.NDArray[Any]] = {"X": rx, "Y": ry, "NORMX": nx, "NORMY": ny, "THETA": th, "R": rr}
            for key in t.extra_format:
                col = np.full(n, "", dtype=object)
                vals = src[key]
                col[rows] = [str(int(v)) for v in vals] if key in ("X", "Y") else [_g8(float(v)) for v in vals]
                extra_cols.append(col.tolist())
        het_haploid = int(((calls == 2) & (pl == 1)).sum())
        if het_haploid:
            res["reasons"]["het_in_haploid"] += het_haploid
        # merged sites
        for s, k in zip(t.merged_sites, t.merged_rows.tolist(), strict=True):
            c = tuple(int(g.genotypes[p.index]) for p in s.probes)
            p_k = int(ploidy[k])
            gs[k] = float(min(float(g.genotype_scores[p.index]) for p in s.probes))
            if p_k == 0:
                gt[k] = "."
                continue
            o = combine_genotype(s, c, p_k, t.inf1_nc_informative)
            gt[k] = gt_text(o, p_k)
            if o.reason:
                res["reasons"][o.reason] += 1
                if o.reason.startswith("W_MERGE") and o.reason != "W_MERGE_ALL_NC":
                    res["merge_failures"].append((s.id, c, o.reason))
        merged_mask = np.zeros(n, np.bool_)
        merged_mask[t.merged_rows] = True
        header = t.header.replace("__CLUSTER__", g.cluster_file or "") + (
            f"##SAMPLE=<ID={sample_name},SentrixID={sample_id}>\n"
            f"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{sample_name}\n"
        )
        low = t.filter_gs
        gsl, bafl, lrrl = gs.tolist(), baf.tolist(), lrr.tolist()
        tbx = TabixBuilder()
        out = Path(out_path)
        with atomic_open(out, "wb") as raw:
            if bgzip:
                w = BgzfWriter(raw)  # type: ignore[arg-type]
                w.write_record(header.encode())
            else:
                raw.write(header.encode())  # type: ignore[arg-type]
            for k in range(n):
                filt = "LowGS" if low is not None and not math.isnan(gsl[k]) and gsl[k] < low else "PASS"
                info = t.site_info[k] if t.site_info else "."
                if merged_mask[k]:
                    body = f"{filt}\t{info}\tGT:GS\t{gt[k]}:{_f4(gsl[k])}\n"
                else:
                    extra_keys = "".join(":" + key for key in t.extra_format)
                    extra_vals = "".join(":" + c[k] for c in extra_cols)
                    body = (
                        f"{filt}\t{info}\tGT:GS:BAF:LRR{extra_keys}\t"
                        f"{gt[k]}:{_f4(gsl[k])}:{_g8(bafl[k])}:{_g8(lrrl[k])}{extra_vals}\n"
                    )
                line = t.prefix[k] + body.encode()
                if bgzip:
                    vs, ve = w.write_record(line)
                    tbx.add(t.site_contig[k], int(t.site_beg[k]), int(t.site_end[k]), vs, ve)
                else:
                    raw.write(line)  # type: ignore[arg-type]
            if bgzip:
                w.close()
        if bgzip:
            write_tbi(Path(str(out) + ".tbi"), tbx)
        res["records"] = n
        res["sex"] = sex
    except SimurgError as e:
        res["error"] = f"{e.code}: {e.message}"
    except Exception as e:
        res["error"] = f"E_INTERNAL: {type(e).__name__}: {e}"
    return res


# ------------------------------------------------------------------------------------------------ command
def run_gtc_to_vcf(opts: VcfOptions, rep: Reporter) -> int:
    from simurg.formats.discovery import discover_gtcs
    from simurg.formats.samplesheet import output_stem, parse_samplesheet
    from simurg.pipeline.inputs import check_paths

    rm = RunManifest(command="gtc-to-vcf", argv=opts.argv)
    rep.run_id = rm.run_id
    check_paths(
        [
            (opts.bpm, "BPM"),
            (opts.fasta, "FASTA"),
            (opts.csv, "CSV manifest"),
            (opts.gtcs, "GTC folder"),
            (opts.sheet, "sample sheet"),
            (opts.exclude_loci, "locus list"),
            (opts.include_loci, "locus list"),
            (opts.auxiliary_loci, "auxiliary loci VCF"),
        ]
    )
    if opts.ploidy_policy not in ("sex-aware", "diploid"):
        raise SimurgError("E_CONFIG_KEY", "--ploidy-policy must be sex-aware or diploid")
    bad = [x for x in opts.info_fields if x not in INFO_DEFS] + [x for x in opts.extra_format if x not in EXTRA_DEFS]
    if bad:
        raise SimurgError(
            "E_CONFIG_KEY",
            f"unknown field(s) {', '.join(bad)}; --info-fields takes {','.join(INFO_DEFS)}; "
            f"--extra-format takes {','.join(EXTRA_DEFS)}",
        )
    if {"GENTRAIN", "CLUSTER_SEP"} & set(opts.info_fields) and opts.egt is None:
        raise SimurgError("E_CONFIG_KEY", "--info-fields GENTRAIN/CLUSTER_SEP need the cluster file (--egt)")
    sheet = parse_samplesheet(opts.sheet) if opts.sheet else None
    root = opts.gtcs or (sheet.root_folder if sheet else None)
    if root is None:
        raise SimurgError("E_NO_INPUT", "give --gtcs or a sample sheet with a RootFolder")
    found = discover_gtcs(root)
    samples: list[tuple[str, str, str]] = []  # (sample_id, sample_name, path)
    if sheet is not None:
        for r in sheet.rows:
            if r.sample_id in found:
                samples.append((r.sample_id, r.sample_name or r.sample_id, str(found[r.sample_id])))
            else:
                rep.warn("W_QC_UNMATCHED_SAMPLES", f"{r.sample_id}: no GTC under {root}", sample_id=r.sample_id)
    else:
        samples = [(k, k, str(v)) for k, v in sorted(found.items())]
    suffix = ".snv.vcf.gz" if opts.bgzip else ".snv.vcf"
    outs = (
        [opts.out / "cohort.snv.bcf"]
        if opts.merge
        else [opts.out / f"{output_stem(name)}{suffix}" for _, name, _ in samples]
    )
    check_outputs_absent([*outs, opts.out / "audit" / "loci.parquet"], force=opts.force)
    if not opts.bgzip:
        rep.warn("W_NO_BGZIP", "writing plain VCF without index; pgx and cn commands need bgzipped, indexed VCFs")
    rep.info(f"mapping manifest {opts.bpm.name} to {opts.fasta.name}")
    table, audit, skips, _genome = build_sites(opts, rep, rm)
    for role, p in (("bpm", opts.bpm), ("csv_manifest", opts.csv), ("fasta", opts.fasta), ("sample_sheet", opts.sheet)):
        if p is not None:
            rm.add_input(role, p)
    opts.out.mkdir(parents=True, exist_ok=True)
    write_audit(opts.out / "audit" / "loci.parquet", audit)
    if opts.merge:
        n = write_cohort_bcf(table, samples, opts.out / "cohort.snv.bcf", opts, rep)
        rm.summary = {"samples": len(samples), "records": n, "skipped_loci_by_code": dict(skips)}
        rm.add_output(opts.out / "cohort.snv.bcf", "bcf/2.2")
        rm.status = "completed"
        rm.write(opts.out / "run_manifest.json")
        rep.info(f"done: cohort BCF with {len(samples)} sample(s), {n} records")
        return ExitCode.OK
    workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(samples)))
    tasks = [(sid, name, path, str(o), opts.bgzip) for (sid, name, path), o in zip(samples, outs, strict=True)]
    rep.info(f"{len(table.prefix)} VCF records; writing {len(tasks)} sample(s) with {workers} worker(s)")
    if workers == 1:
        _init(table)
        results = [write_sample(t) for t in tasks]
    else:
        with ProcessPoolExecutor(
            workers, mp_context=mp.get_context("spawn"), initializer=_init, initargs=(table,)
        ) as ex:
            results = list(ex.map(write_sample, tasks))
    reasons: Counter[str] = Counter()
    failures = []
    for res in results:
        if "error" in res:
            rep.error(res["error"].split(":", 1)[0], f"{res['sample_id']}: {res['error']}", sample_id=res["sample_id"])
            continue
        reasons.update(res["reasons"])
        failures += [(res["sample_id"], *f) for f in res["merge_failures"]]
    for code, n in sorted(reasons.items()):
        if code.startswith("W_"):
            rep.warn(code, f"{n} genotype(s): {code}")
    write_merge_audit(opts.out / "audit" / "merges.parquet", failures)
    ok = sum("error" not in r for r in results)
    rm.summary = {"samples": len(results), "succeeded": ok, "records_per_sample": len(table.prefix),
                  "skipped_loci_by_code": dict(skips), "genotype_notes": dict(reasons)}  # fmt: skip
    for o in outs:
        if o.exists():
            rm.add_output(o, "vcf/4.3")
    rm.status = "completed"
    rm.write(opts.out / "run_manifest.json")
    rep.info(f"done: {ok}/{len(results)} VCF(s) in {opts.out}")
    if ok == len(results):
        return ExitCode.OK
    return ExitCode.PARTIAL if ok else ExitCode.ALL_FAILED


def write_audit(path: Path, audit: list[AuditRow]) -> None:
    import pyarrow as pa

    from simurg.core.schemas import write_parquet

    cols = list(AuditRow.__dataclass_fields__)
    table = pa.table({c: [getattr(a, c) for a in audit] for c in cols})
    path.parent.mkdir(parents=True, exist_ok=True)
    write_parquet(table, path)


def write_merge_audit(path: Path, rows: list[tuple[str, str, tuple[int, ...], str]]) -> None:
    import pyarrow as pa

    from simurg.core.schemas import write_parquet

    table = pa.table(
        {
            "sample_id": [r[0] for r in rows],
            "site_id": [r[1] for r in rows],
            "calls": ["/".join(str(c) for c in r[2]) for r in rows],
            "reason": [r[3] for r in rows],
        }
    )
    path.parent.mkdir(parents=True, exist_ok=True)
    write_parquet(table, path)


# ------------------------------------------------------------------------------------------------ cohort BCF
def _f32_payload(v: npt.NDArray[np.float32]) -> bytes:
    u = np.asarray(v, dtype="<f4").view("<u4").copy()
    u[np.isnan(np.asarray(v, dtype=np.float32))] = 0x7F800001
    return u.tobytes()


def write_cohort_bcf(
    t: SiteTable, samples: list[tuple[str, str, str]], out: Path, opts: VcfOptions, rep: Reporter
) -> int:
    """OUT-040: one multi-sample BCF (+ .csi), same records and FORMAT rules as the per-sample VCFs."""
    from simurg.formats.bcf import (
        FLOAT,
        INT8,
        INT32,
        BcfWriter,
        CsiBuilder,
        _type_byte,
        encode_record,
        header_dictionaries,
        typed_floats,
        typed_ints,
        typed_string,
    )
    from simurg.formats.bpm import read_bpm
    from simurg.pipeline.qc import default_tmp, load_cohort
    from simurg.signal.normalization.apply import polar

    bpm = read_bpm(opts.bpm)
    workers = max(1, opts.threads or (os.cpu_count() or 1))
    store, infos = load_cohort(
        [(sid, path) for sid, _, path in samples],
        bpm.num_loci,
        default_tmp(out.parent),
        baf=True,
        scores=True,
        lrr=True,
        workers=workers,
        lookup=bpm.normalization_lookups() if t.extra_format else None,
    )
    try:
        for s in infos:
            if s.error:
                raise SimurgError("E_GTC_UNSUITABLE", f"{s.sample_id}: {s.error}")
        names = [name for _, name, _ in samples]
        cluster = read_cluster_name(samples[0][2])
        header = t.header.replace("__CLUSTER__", cluster)
        header += "".join(f"##SAMPLE=<ID={name},SentrixID={sid}>\n" for sid, name, _ in samples)
        header += "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t" + "\t".join(names) + "\n"
        strings, contigs = header_dictionaries(header)
        k_gt, k_gs, k_baf, k_lrr = strings["GT"], strings["GS"], strings["BAF"], strings["LRR"]
        f_pass, f_low = strings["PASS"], strings.get("LowGS", -1)
        n = len(names)
        sexes = [s.sex_est for s in infos]
        ploidy = np.stack([ploidy_for(t.site_class, sx, t.ploidy_policy) for sx in sexes], axis=1)  # sites × samples
        single_pos = np.full(len(t.prefix), -1, np.int64)
        single_pos[t.single_rows] = np.arange(len(t.single_rows))
        merged_pos = {int(k): j for j, k in enumerate(t.merged_rows.tolist())}
        csi = CsiBuilder(len(contigs))
        fdesc = _type_byte(1, FLOAT)
        written = 0
        with atomic_open(out, "wb") as raw:
            w = BcfWriter(raw, header)
            chunk = 20000
            for start in range(0, len(t.prefix), chunk):
                rows = range(start, min(len(t.prefix), start + chunk))
                need = {int(t.single_locus[single_pos[k]]) for k in rows if single_pos[k] >= 0}
                for k in rows:
                    if k in merged_pos:
                        need.update(p.index for p in t.merged_sites[merged_pos[k]].probes)
                loci = sorted(need)
                col = {loc: j for j, loc in enumerate(loci)}
                sel = np.array(loci, np.int64)
                gt = np.asarray(store.gt[:, sel]) if len(sel) else np.zeros((n, 0), np.uint8)
                sc = np.asarray(store.score[:, sel]) if len(sel) else np.zeros((n, 0), np.float32)  # type: ignore[index]
                bf = np.asarray(store.baf[:, sel]) if len(sel) else np.zeros((n, 0), np.float32)  # type: ignore[index]
                lr = np.asarray(store.lrr[:, sel]) if len(sel) else np.zeros((n, 0), np.float32)  # type: ignore[index]
                if t.extra_format and len(sel):
                    ex_rx, ex_ry = np.asarray(store.rx[:, sel]), np.asarray(store.ry[:, sel])  # type: ignore[index]
                    ex_nx, ex_ny = np.asarray(store.nx[:, sel]), np.asarray(store.ny[:, sel])  # type: ignore[index]
                for k in rows:
                    pl = ploidy[k]
                    if single_pos[k] >= 0:
                        j = int(single_pos[k])
                        c = col[int(t.single_locus[j])]
                        calls = gt[:, c].astype(np.int64)
                        ia, ib = (int(x) for x in t.single_ab[j])  # type: ignore[index]
                        a1 = np.where(calls == 3, ib, ia)
                        a2 = np.where(calls == 1, ia, ib)
                        lo, hi = np.minimum(a1, a2), np.maximum(a1, a2)
                        width = max(1, int(pl.max()))
                        cells = np.full((n, width), -127, np.int8)
                        miss = (calls == 0) | (pl == 0) | ((pl == 1) & (calls == 2))
                        hap = pl <= 1
                        cells[:, 0] = np.where(miss, 0, ((np.where(hap & (calls == 3), ib, lo) + 1) << 1)).astype(
                            np.int8
                        )
                        if width == 2:
                            second = np.where(miss, 0, (hi + 1) << 1).astype(np.int8)
                            cells[:, 1] = np.where(hap, np.int8(-127), second)
                        gs = sc[:, c]
                        b = bf[:, c]
                        b = np.where(t.single_flip[j], (np.float32(1.0) - b).astype(np.float32), b)
                        fmt = [(k_gt, _type_byte(width, INT8), cells.tobytes()), (k_gs, fdesc, _f32_payload(gs)),
                               (k_baf, fdesc, _f32_payload(b)), (k_lrr, fdesc, _f32_payload(lr[:, c]))]  # fmt: skip
                        if t.extra_format:
                            zero = (ex_rx[:, c] == 0) & (ex_ry[:, c] == 0)
                            enx = np.where(zero, np.float32(np.nan), ex_nx[:, c]).astype(np.float32)
                            eny = np.where(zero, np.float32(np.nan), ex_ny[:, c]).astype(np.float32)
                            th, rr = polar(enx, eny)
                            xsrc = {
                                "X": ex_rx[:, c],
                                "Y": ex_ry[:, c],
                                "NORMX": enx,
                                "NORMY": eny,
                                "THETA": th,
                                "R": rr,
                            }
                            for key in t.extra_format:
                                if key in ("X", "Y"):
                                    payload = xsrc[key].astype("<i4").tobytes()
                                    fmt.append((strings[key], _type_byte(1, INT32), payload))
                                else:
                                    fmt.append((strings[key], fdesc, _f32_payload(xsrc[key].astype(np.float32))))
                        gs_min = gs
                    else:
                        site = t.merged_sites[merged_pos[k]]
                        cs = [col[p.index] for p in site.probes]
                        g_all = gt[:, cs].astype(np.int64)
                        gs_min = sc[:, cs].min(axis=1)
                        width = max(1, int(pl.max()))
                        cells = np.full((n, width), -127, np.int8)
                        for si in range(n):
                            p_k = int(pl[si])
                            if p_k == 0:
                                cells[si, 0] = 0
                                continue
                            o = combine_genotype(site, tuple(int(x) for x in g_all[si]), p_k, t.inf1_nc_informative)
                            if o.alleles is None:
                                cells[si, :p_k] = 0
                            else:
                                for q, a in enumerate(o.alleles):
                                    cells[si, q] = (a + 1) << 1
                        fmt = [(k_gt, _type_byte(width, INT8), cells.tobytes()), (k_gs, fdesc, _f32_payload(gs_min))]
                    low = opts.filter_gs is not None and bool(np.nanmin(gs_min) < opts.filter_gs) if n else False
                    info = []
                    for key, kind, vals in t.site_info_typed[k] if t.site_info_typed else []:
                        if kind == "Integer":
                            enc = typed_ints([int(v) for v in vals])
                        elif kind == "Float":
                            enc = typed_floats([float(v) for v in vals])
                        else:
                            enc = typed_string(",".join(str(v) for v in vals))
                        info.append((strings[key], enc))
                    rec = encode_record(contigs[t.site_contig[k]], int(t.site_beg[k]),
                                        int(t.site_end[k] - t.site_beg[k]), t.site_id[k], t.site_alleles[k],
                                        f_low if low and f_low >= 0 else f_pass, n, fmt, info)  # fmt: skip
                    vs, ve = w.write(rec)
                    csi.add(contigs[t.site_contig[k]], int(t.site_beg[k]), int(t.site_end[k]), vs, ve)
                    written += 1
            w.close()
        Path(str(out) + ".csi").write_bytes(csi.encode())
    finally:
        store.close()
    return written


def read_cluster_name(gtc_path: str) -> str:
    from simurg.formats.gtc import read_gtc

    return read_gtc(gtc_path, check_write_complete=False).cluster_file or ""