P4 (WIP): GenCall-equivalent calling, BAF/LRR, sex, GTC v5, simurg genotype

- genotype/clusters.py: EGT aligned to manifest order (float32 mean R)
- genotype/gencall.py: flair, six-case likelihood, priors/posteriors,
  fuzzy sub-scores and score map, vectorised with the reference's
  float/double evaluation order and libm transcendentals
- genotype/baf_lrr.py, sex.py (rules v1/v2), sample.py (IDAT -> Gtc)
- pipeline/inputs.py shared by ingest and genotype
- pipeline/genotype.py + CLI: GTC per sample, sample_summary.tsv/json
  (OUT-060..063), run manifest
- GTCs from IDATs are byte-identical to the oracle for 36 GSA samples x
  2 normalization versions (NaN sign canonicalised to the x86 pattern)

src/simurg/pipeline/genotype.py

@@ -0,0 +1,344 @@
+"""``simurg genotype``: IDAT pairs + BPM + EGT → GTC per sample + ``sample_summary.{tsv,json}`` (spec CLI-021,
+OUT-001, OUT-060 … OUT-063; phase P4).
+"""
+
+from __future__ import annotations
+
+import datetime as _dt
+import json
+import math
+import multiprocessing as mp
+import os
+from concurrent.futures import ProcessPoolExecutor
+from dataclasses import dataclass, field, replace
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+from simurg import __version__
+from simurg.core.atomic import atomic_write_text, check_outputs_absent
+from simurg.core.errors import ExitCode, SimurgError
+from simurg.core.log import Reporter
+from simurg.core.provenance import RunManifest
+from simurg.formats.bpm import read_bpm
+from simurg.formats.egt import read_egt
+from simurg.formats.gtc import write_gtc
+from simurg.formats.idat import check_pair, read_idat
+from simurg.genotype.baf_lrr import chrom_core
+from simurg.genotype.clusters import align_clusters
+from simurg.genotype.sample import GenotypeParams, ManifestContext, genotype_sample, manifest_context
+from simurg.genotype.sex import SexParams
+from simurg.pipeline.inputs import check_paths, resolve_sample_set, sample_dir, sheet_values
+from simurg.signal.normalization.apply import normalize, polar
+from simurg.signal.raw_xy import LocusMapCache, check_coverage
+
+SUMMARY_SCHEMA = "simurg.sample_summary/1"
+FIXED_COLUMNS = [
+    "sample_id",
+    "sample_name",
+    "sample_folder",
+    "status",
+    "fail_reason",
+    "autosomal_call_rate",
+    "call_rate",
+    "logr_dev",
+    "sex_est",
+    "sex_reported",
+    "sex_check",
+    "gc10",
+    "gc50",
+    "n_calls",
+    "n_nocalls",
+    "n_intensity_only",
+    "p05_x",
+    "p50_x",
+    "p95_x",
+    "p05_y",
+    "p50_y",
+    "p95_y",
+    "het_rate",
+    "contamination",
+]
+QC_COLUMNS = ["functional_status", "control_status", "control_flag_codes", "flags"]  # filled from phase P5
+FLOAT_COLUMNS = {"autosomal_call_rate", "call_rate", "logr_dev", "gc10", "gc50", "het_rate", "contamination"}
+
+
+@dataclass
+class GenotypeOptions:
+    out: Path
+    bpm: Path
+    egt: Path
+    idats: Path | None = None
+    sheet: Path | None = None
+    csv: Path | None = None
+    threads: int = 0
+    force: bool = False
+    layout: str = "nested"
+    on_duplicate_idat: str = "error"
+    cutoff: float = 0.15
+    norm_version: str = "1.2.0"
+    sex_version: int = 2
+    gtc_version: int = 5
+    autocall_date: str | None = None
+    assay_control_probes: tuple[str, ...] = ("TGA_Ctrl_5716",)
+    argv: list[str] = field(default_factory=list)
+
+
+@dataclass
+class Job:
+    sample_id: str
+    sample_name: str
+    grn: str
+    red: str
+    gtc_path: str
+
+
+_W: dict[str, Any] = {}
+
+
+def _init(
+    ctx: ManifestContext,
+    params: GenotypeParams,
+    auto_1_22: np.ndarray[Any, np.dtype[np.bool_]],
+    probe_idx: dict[str, int],
+    started: _dt.datetime,
+) -> None:
+    _W.update(ctx=ctx, params=params, auto=auto_1_22, probes=probe_idx, started=started)
+
+
+def _process(job: Job) -> dict[str, Any]:
+    ctx: ManifestContext = _W["ctx"]
+    params: GenotypeParams = _W["params"]
+    row: dict[str, Any] = {
+        "sample_id": job.sample_id,
+        "sample_name": job.sample_name,
+        "status": "OK",
+        "fail_reason": "",
+        "notes": [],
+    }
+    try:
+        grn, red = read_idat(job.grn), read_idat(job.red)
+        check_pair(grn, red)
+        _, notes = check_coverage(grn, ctx.manifest)
+        row["notes"].extend(notes)
+        g = genotype_sample(grn, red, ctx, params, _W["started"])
+        write_gtc(g, Path(job.gtc_path))
+        auto = _W["auto"] & ~ctx.manifest.intensity_only
+        n_auto = int(auto.sum())
+        row.update(
+            autosomal_call_rate=float(np.float32(int((g.genotypes[auto] > 0).sum())) / np.float32(n_auto))
+            if n_auto
+            else math.nan,
+            call_rate=float(g.call_rate),
+            logr_dev=float(g.logr_dev),
+            sex_est=g.gender_char,
+            gc10=float(g.gc10),
+            gc50=float(g.gc50),
+            n_calls=g.num_calls,
+            n_nocalls=g.num_no_calls,
+            n_intensity_only=g.num_intensity_only,
+            p05_x=int(g.percentiles_x[0]),
+            p50_x=int(g.percentiles_x[1]),
+            p95_x=int(g.percentiles_x[2]),
+            p05_y=int(g.percentiles_y[0]),
+            p50_y=int(g.percentiles_y[1]),
+            p95_y=int(g.percentiles_y[2]),
+        )
+        if _W["probes"]:
+            idx = np.array(list(_W["probes"].values()), dtype=np.int64)
+            nx, ny = normalize(g.raw_x[idx], g.raw_y[idx], ctx.lookup[idx], g.normalization_transforms)
+            _, r = polar(nx, ny)
+            for (name, _), val in zip(_W["probes"].items(), r.tolist(), strict=True):
+                row[f"{name}_norm_r"] = float(val)
+    except SimurgError as e:
+        row.update(status="FAILED", fail_reason=f"{e.code}: {e.message}", error_code=e.code)
+    except Exception as e:
+        row.update(status="FAILED", fail_reason=f"E_INTERNAL: {type(e).__name__}: {e}", error_code="E_INTERNAL")
+    return row
+
+
+def _fmt(col: str, v: Any) -> str:
+    if v is None or v == "":
+        return ""
+    if isinstance(v, float) or col in FLOAT_COLUMNS or col.endswith("_norm_r"):
+        f = float(v)
+        return "NaN" if math.isnan(f) else f"{f:.8f}"
+    return str(v).replace("\t", " ").replace("\n", " ")
+
+
+def _json_value(col: str, v: Any) -> Any:
+    if v == "" or v is None:
+        return None
+    if isinstance(v, float):
+        return None if math.isnan(v) else v
+    return v
+
+
+def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
+    rm = RunManifest(command="genotype", argv=opts.argv)
+    rep.run_id = rm.run_id
+    for k in (
+        "bpm",
+        "egt",
+        "idats",
+        "sheet",
+        "csv",
+        "threads",
+        "cutoff",
+        "norm_version",
+        "sex_version",
+        "gtc_version",
+        "layout",
+        "on_duplicate_idat",
+    ):
+        v = getattr(opts, k)
+        rm.set_option(k, str(v) if v is not None else None, "cli")
+    check_paths(
+        [
+            (opts.bpm, "BPM"),
+            (opts.egt, "cluster file"),
+            (opts.sheet, "sample sheet"),
+            (opts.csv, "CSV manifest"),
+            (opts.idats, "IDAT folder"),
+        ]
+    )
+    if not 0.0 <= opts.cutoff <= 1.0:
+        raise SimurgError("E_CONFIG_KEY", f"--cutoff must be within [0, 1], got {opts.cutoff}")
+
+    sset = resolve_sample_set(opts.idats, opts.sheet, opts.on_duplicate_idat, rep, rm)
+    jobs: list[Job] = []
+    rows: dict[str, dict[str, Any]] = {}
+    for s in sset.samples:
+        folder = ""
+        if s.pair is not None:
+            try:
+                folder = str(s.pair.folder.relative_to(sset.folder))
+            except ValueError:
+                folder = str(s.pair.folder)
+        base = {
+            "sample_id": s.sample_id,
+            "sample_name": s.sample_name,
+            "sample_folder": folder,
+            "sex_reported": getattr(s.sheet_row, "reported_sex", "U") if s.sheet_row else "U",
+            **{f"sheet:{k}": v for k, v in sheet_values(s, sset.sheet).items()},
+        }
+        rows[s.sample_id] = base
+        if s.status != "READY" or s.pair is None or not s.pair.complete:
+            base.update(status="MISSING_INPUT", fail_reason="no complete IDAT pair found")
+            continue
+        stem = sset.stems[s.sample_id]
+        gtc = sample_dir(opts.out, stem, opts.layout) / f"{stem}.gtc"
+        jobs.append(Job(s.sample_id, s.sample_name, str(s.pair.grn), str(s.pair.red), str(gtc)))
+
+    summary_tsv, summary_json = opts.out / "sample_summary.tsv", opts.out / "sample_summary.json"
+    check_outputs_absent([summary_tsv, summary_json] + [Path(j.gtc_path) for j in jobs], force=opts.force)
+
+    rep.info(f"reading manifest {opts.bpm} and cluster file {opts.egt}")
+    manifest = read_bpm(opts.bpm)
+    if opts.csv:
+        from simurg.formats.csv_manifest import compare_manifests, read_csv_manifest
+
+        problems = compare_manifests(manifest, read_csv_manifest(opts.csv, load_sequences=False))
+        if problems:
+            raise SimurgError("E_BPM_CSV_MISMATCH", "BPM and CSV manifests differ:\n  " + "\n  ".join(problems))
+    egt = read_egt(opts.egt)
+    if (
+        egt.manifest_name
+        and manifest.manifest_name
+        and egt.manifest_name.split(".")[0].lower() not in manifest.manifest_name.lower()
+    ):
+        rep.warn(
+            "W_EGT_NAME", f"cluster file was built for {egt.manifest_name!r}, manifest is {manifest.manifest_name!r}"
+        )
+    clusters = align_clusters(egt, manifest)
+    ctx = manifest_context(manifest, clusters, cluster_file_name=Path(opts.egt).name)
+    lite = manifest.lite()
+    ctx = replace(ctx, manifest=lite, cache=LocusMapCache(lite))  # workers get no string columns
+    for role, p in (("bpm", opts.bpm), ("egt", opts.egt), ("csv_manifest", opts.csv), ("sample_sheet", opts.sheet)):
+        if p is not None:
+            rm.add_input(role, p)
+    auto_1_22 = np.array([chrom_core(c) in {str(i) for i in range(1, 23)} for c in manifest.chroms])
+    name_idx = manifest.name_index()
+    probes = {p: name_idx[p] for p in opts.assay_control_probes if p in name_idx}
+    params = GenotypeParams(
+        cutoff=opts.cutoff,
+        norm_version=opts.norm_version,
+        sex=SexParams(version=opts.sex_version),
+        gtc_version=opts.gtc_version,
+        autocall_date=opts.autocall_date,
+    )
+    started = _dt.datetime.now()
+
+    workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(jobs) or 1))
+    rep.info(f"genotyping {len(jobs)} sample(s) with {workers} worker(s)")
+    try:
+        if workers == 1:
+            _init(ctx, params, auto_1_22, probes, started)
+            results = [_process(j) for j in jobs]
+        else:
+            with ProcessPoolExecutor(
+                workers,
+                mp_context=mp.get_context("spawn"),
+                initializer=_init,
+                initargs=(ctx, params, auto_1_22, probes, started),
+            ) as ex:
+                results = list(ex.map(_process, jobs))
+    except KeyboardInterrupt:
+        rm.status = "interrupted"
+        rm.write(opts.out / "run_manifest.json")
+        raise
+    finally:
+        _W.clear()
+
+    for r in results:
+        for code, msg in r.pop("notes", []):
+            rep.warn(code, msg, sample_id=r["sample_id"])
+        if r["status"] == "FAILED":
+            rep.error(
+                r.get("error_code", "E_INTERNAL"), f"{r['sample_id']}: {r['fail_reason']}", sample_id=r["sample_id"]
+            )
+            rm.failed_samples.append({"sample_id": r["sample_id"], "code": r.get("error_code", "")})
+            r.update(
+                autosomal_call_rate=0.0, call_rate=0.0, logr_dev=math.nan, gc10=math.nan, gc50=math.nan, sex_est="U"
+            )  # OUT-061
+            for probe in probes:
+                r[f"{probe}_norm_r"] = math.nan
+        rows[r["sample_id"]].update(r)
+
+    probe_cols = [f"{p}_norm_r" for p in probes]
+    sheet_cols = [f"sheet:{c}" for c in (sset.sheet.passthrough_columns() if sset.sheet else [])]
+    out_rows = []
+    for s in sset.samples:
+        r = rows[s.sample_id]
+        known = {"M", "F"}
+        est, rep_sex = r.get("sex_est", "U"), r.get("sex_reported", "U")
+        r["sex_check"] = ("match" if est == rep_sex else "mismatch") if est in known and rep_sex in known else "unknown"
+        out_rows.append(r)
+    columns = FIXED_COLUMNS + probe_cols + QC_COLUMNS + sheet_cols
+    header = [c.removeprefix("sheet:") for c in columns]
+    lines = ["\t".join(header)] + ["\t".join(_fmt(c, r.get(c, "")) for c in columns) for r in out_rows]
+    atomic_write_text(summary_tsv, "\n".join(lines) + "\n")
+    doc = {
+        "schema": SUMMARY_SCHEMA,
+        "simurg_version": __version__,
+        "columns": header,
+        "samples": [{h: _json_value(c, r.get(c, "")) for c, h in zip(columns, header, strict=True)} for r in out_rows],
+    }
+    atomic_write_text(summary_json, json.dumps(doc, indent=2, ensure_ascii=False) + "\n")
+
+    ok = sum(r.get("status") == "OK" for r in out_rows)
+    rm.summary = {
+        "samples": len(out_rows),
+        "succeeded": ok,
+        "failed_or_missing": len(out_rows) - ok,
+        "warnings": rep.warning_counts(),
+    }
+    rm.add_output(summary_tsv, SUMMARY_SCHEMA)
+    rm.schema_versions["*.gtc"] = f"gtc/{opts.gtc_version}"
+    rm.status = "completed"
+    rm.write(opts.out / "run_manifest.json")
+    rep.info(f"done: {ok}/{len(out_rows)} samples genotyped; summary {summary_tsv}")
+    if ok == len(out_rows):
+        return ExitCode.OK
+    return ExitCode.PARTIAL if ok else ExitCode.ALL_FAILED

src/simurg/pipeline/ingest.py

@@ -36,9 +36,9 @@ from simurg.core.schemas import (
     write_parquet,
 )
 from simurg.formats.bpm import Manifest, read_bpm
-from simurg.formats.discovery import ResolvedSample, discover_idats, resolve_input_folder, resolve_samples
+from simurg.formats.discovery import ResolvedSample
 from simurg.formats.idat import check_pair, parse_idat_filename, read_idat
-from simurg.formats.samplesheet import SampleSheet, output_stem, parse_samplesheet
+from simurg.pipeline.inputs import check_paths, resolve_sample_set, sample_dir
 from simurg.signal.raw_xy import LocusMapCache, assemble_raw_xy, check_coverage
 
 SUMMARY_COLUMNS = [
@@ -93,7 +93,7 @@ def _init_worker(manifest: Manifest, manifest_sha256: str) -> None:
 
 
 def _sample_paths(out: Path, stem: str, layout: str) -> tuple[Path, Path]:
-    d = out / "samples" / stem if layout == "nested" else out
+    d = sample_dir(out, stem, layout)
     return d / f"{stem}.intensity.parquet", d / f"{stem}.controls.parquet"
 
 
@@ -202,61 +202,18 @@ def run_ingest(opts: IngestOptions, rep: Reporter) -> int:
         v = getattr(opts, k)
         rm.set_option(k, str(v) if v is not None else None, "cli")
 
-    for p, what in (
-        (opts.bpm, "BPM"),
-        (opts.sheet, "sample sheet"),
-        (opts.csv, "CSV manifest"),
-        (opts.idats, "IDAT folder"),
-    ):
-        if p is not None and not Path(p).exists():
-            raise SimurgError("E_FILE_NOT_FOUND", f"{what} not found: {p}", path=str(p))
-
-    # --- inputs (IN-001 … IN-022) --------------------------------------------------------------------
-    sheet: SampleSheet | None = parse_samplesheet(opts.sheet) if opts.sheet else None
-    if sheet:
-        rm.sheet_header = sheet.header
-        for code, msg in sheet.warnings:
-            rep.warn(code, msg)
-    folder, notes = resolve_input_folder(opts.idats, sheet)
-    for code, msg in notes:
-        rep.warn(code, msg) if code.startswith("W_") else rep.info(msg, code)
-    disc = discover_idats(folder, on_duplicate=opts.on_duplicate_idat)
-    for code, msg in disc.warnings:
-        rep.warn(code, msg)
-    samples, notes = resolve_samples(disc, sheet)
-    for code, msg in notes:
-        rep.info(msg, code)
-    if not samples:
-        raise SimurgError("E_NO_INPUT", f"no IDAT pairs found under {folder}")
-
-    # --- names and output paths (IN-010, IN-012) -----------------------------------------------------
-    use_idat_names = sheet is None or not sheet.has_sample_name
-    stems: dict[str, str] = {}
+    check_paths(
+        [(opts.bpm, "BPM"), (opts.sheet, "sample sheet"), (opts.csv, "CSV manifest"), (opts.idats, "IDAT folder")]
+    )
+    sset = resolve_sample_set(opts.idats, opts.sheet, opts.on_duplicate_idat, rep, rm)
+    samples = sset.samples
     jobs: list[SampleJob] = []
     rows: list[dict[str, Any]] = []
     for s in samples:
-        if use_idat_names and s.pair is not None and s.pair.grn is not None:
-            # sheet Sample_Name > IDAT sample name (field 406) > sample_id
-            try:
-                idat_name = read_idat(s.pair.grn, load_arrays=False).sample_id.strip()
-            except SimurgError:
-                idat_name = ""
-            if idat_name:
-                s.sample_name = idat_name
-        stem = output_stem(s.sample_name)
-        if stem in stems and stems[stem] != s.sample_id:
-            if not use_idat_names:  # sheet names that sanitise to the same stem
-                raise SimurgError(
-                    "E_SHEET_NAME_COLLISION",
-                    f"samples {stems[stem]} and {s.sample_id} both map to output name {stem!r}",
-                )
-            rep.info(f"{s.sample_id}: IDAT sample name {s.sample_name!r} is not unique; using the sample ID")
-            s.sample_name = s.sample_id
-            stem = output_stem(s.sample_id)
-        stems[stem] = s.sample_id
         if s.status != "READY" or s.pair is None or not s.pair.complete:
             rows.append(_missing_row(s))
             continue
+        stem = sset.stems[s.sample_id]
         ip, cp = _sample_paths(opts.out, stem, opts.layout)
         jobs.append(
             SampleJob(s.sample_id, s.sample_name, stem, str(s.pair.grn), str(s.pair.red), str(ip), str(cp), opts.resume)

src/simurg/pipeline/inputs.py

@@ -0,0 +1,87 @@
+"""Shared sample-set resolution for IDAT-based commands (spec CLI-006, IN-001 … IN-022, IN-010/IN-012)."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from pathlib import Path
+from typing import Any
+
+from simurg.core.errors import SimurgError
+from simurg.core.log import Reporter
+from simurg.core.provenance import RunManifest
+from simurg.formats.discovery import ResolvedSample, discover_idats, resolve_input_folder, resolve_samples
+from simurg.formats.idat import read_idat
+from simurg.formats.samplesheet import SampleSheet, output_stem, parse_samplesheet
+
+
+@dataclass
+class SampleSet:
+    folder: Path
+    sheet: SampleSheet | None
+    samples: list[ResolvedSample]
+    stems: dict[str, str]  # sample_id -> output stem
+
+
+def check_paths(pairs: list[tuple[Path | None, str]]) -> None:
+    for p, what in pairs:
+        if p is not None and not Path(p).exists():
+            raise SimurgError("E_FILE_NOT_FOUND", f"{what} not found: {p}", path=str(p))
+
+
+def resolve_sample_set(
+    idats: Path | None, sheet_path: Path | None, on_duplicate_idat: str, rep: Reporter, rm: RunManifest
+) -> SampleSet:
+    sheet = parse_samplesheet(sheet_path) if sheet_path else None
+    if sheet:
+        rm.sheet_header = sheet.header
+        for code, msg in sheet.warnings:
+            rep.warn(code, msg)
+    folder, notes = resolve_input_folder(idats, sheet)
+    for code, msg in notes:
+        rep.warn(code, msg) if code.startswith("W_") else rep.info(msg, code)
+    disc = discover_idats(folder, on_duplicate=on_duplicate_idat)
+    for code, msg in disc.warnings:
+        rep.warn(code, msg)
+    samples, notes = resolve_samples(disc, sheet)
+    for code, msg in notes:
+        rep.info(msg, code)
+    if not samples:
+        raise SimurgError("E_NO_INPUT", f"no IDAT pairs found under {folder}")
+
+    use_idat_names = sheet is None or not sheet.has_sample_name
+    stems: dict[str, str] = {}
+    taken: dict[str, str] = {}
+    for s in samples:
+        if use_idat_names and s.pair is not None and s.pair.grn is not None:
+            # IN-010: sheet Sample_Name > IDAT sample name (field 406) > sample_id
+            try:
+                idat_name = read_idat(s.pair.grn, load_arrays=False).sample_id.strip()
+            except SimurgError:
+                idat_name = ""
+            if idat_name:
+                s.sample_name = idat_name
+        stem = output_stem(s.sample_name)
+        if stem in taken and taken[stem] != s.sample_id:
+            if not use_idat_names:  # IN-012: sheet names that sanitise to the same stem
+                raise SimurgError(
+                    "E_SHEET_NAME_COLLISION",
+                    f"samples {taken[stem]} and {s.sample_id} both map to output name {stem!r}",
+                )
+            rep.info(f"{s.sample_id}: IDAT sample name {s.sample_name!r} is not unique; using the sample ID")
+            s.sample_name = s.sample_id
+            stem = output_stem(s.sample_id)
+        taken[stem] = s.sample_id
+        stems[s.sample_id] = stem
+    return SampleSet(folder=folder, sheet=sheet, samples=samples, stems=stems)
+
+
+def sample_dir(out: Path, stem: str, layout: str) -> Path:
+    return out / "samples" / stem if layout == "nested" else out
+
+
+def sheet_values(s: ResolvedSample, sheet: SampleSheet | None) -> dict[str, Any]:
+    """Pass-through sheet columns for summaries (IN-013)."""
+    if sheet is None or s.sheet_row is None:
+        return {}
+    row = s.sheet_row
+    return {c: row.values.get(c, "") for c in sheet.passthrough_columns()}  # type: ignore[attr-defined]