P4: unit tests for calling/BAF/LRR/sex/CLI; vendor validation script

- tests/unit/test_genotype.py: fuzzy memberships incl. NaN rules, calls
  at cluster centres, intensity threshold, BAF/LRR interpolation and
  NaN pattern, LogR-dev locus set, sex rules v1/v2, TOP-strand base
  calls, invariant dates, end-to-end `simurg genotype` (summary schema,
  pass-through columns, thread/layout determinism, overwrite refusal)
- tests/oracle/test_p4_genotype.py: 36 GSA samples x 2 versions
  byte-identical to the oracle
- scripts/validate_p4_large.py + cloud task p4: compare every vendor
  GTC field by field and as whole files; reports which sex rule matches

scripts/cloud/tasks/p4.sh

@@ -0,0 +1,13 @@
+#!/bin/bash
+# Phase P4: genotype every Omni2.5 IDAT pair and compare the GTC with the vendor's (fields + whole file).
+# Args: <data dir with bpm/egt and data/>  <results dir>
+set -uo pipefail
+D=$1; OUT=$2
+SRC=s3://1000genomes/release/20130502/supporting/hd_genotype_chip/broad_intensities
+mkdir -p $D/egt
+aws s3 cp --no-sign-request "$SRC/HumanOmni2.5-4v1-Multi_H.egt" $D/egt/ --only-show-errors
+cp $D/HumanOmni2.5-4v1-Multi_B.egt $D/egt/
+/usr/bin/time -v .venv/bin/python scripts/validate_p4_large.py --root $D/data --bpm $D/HumanOmni2.5-4v1-Multi_B.bpm \
+  --egt-dir $D/egt --out $OUT/summary.json --norm-version 1.1.2 --sex-version 1 --workers "$(nproc)" \
+  2> $OUT/time.txt || echo "P4 validation reported differences"
+tail -25 $OUT/time.txt

scripts/validate_p4_large.py

@@ -0,0 +1,248 @@
+#!/usr/bin/env python3
+"""Phase P4 validation against vendor GTCs: genotype every IDAT pair and compare the produced GTC with the vendor's,
+field by field and as whole files (after aligning free-text metadata that cannot match).
+
+Usage: validate_p4_large.py --root DIR --bpm F --egt-dir DIR --out summary.json
+                            [--norm-version 1.1.2] [--sex-version 1] [--workers N] [--limit N]
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import multiprocessing as mp
+import os
+import re
+import sys
+import time
+import traceback
+from concurrent.futures import ProcessPoolExecutor
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+NEG_NAN = np.frombuffer(np.uint32(0xFFC00000).tobytes(), dtype=np.float32)[0]
+ALIGN = (
+    "cluster_file",
+    "snp_manifest",
+    "autocall_date",
+    "autocall_version",
+    "imaging_date",
+    "sample_name",
+    "sample_plate",
+    "sample_well",
+    "scanner_data",
+)
+_W: dict[str, Any] = {}
+
+
+def _init(bpm: str, egt_dir: str, norm: str, sex_version: int) -> None:
+    _W.update(bpm=bpm, egt_dir=egt_dir, norm=norm, sex_version=sex_version, ctx={})
+
+
+def _ctx(cluster_field: str) -> Any:
+    from simurg.formats.bpm import read_bpm
+    from simurg.formats.egt import read_egt
+    from simurg.genotype.clusters import align_clusters
+    from simurg.genotype.sample import manifest_context
+
+    name = re.sub(r"\(MD5Checksum=.*\)$", "", cluster_field).strip()
+    if name not in _W["ctx"]:
+        if "manifest" not in _W:
+            _W["manifest"] = read_bpm(_W["bpm"])
+        path = Path(_W["egt_dir"]) / name
+        if not path.exists():
+            raise FileNotFoundError(f"cluster file {name!r} named in the GTC is not in {_W['egt_dir']}")
+        m = _W["manifest"]
+        _W["ctx"][name] = manifest_context(m, align_clusters(read_egt(path), m), cluster_file_name=name)
+    return _W["ctx"][name]
+
+
+def _canon(a: np.ndarray) -> np.ndarray:  # type: ignore[type-arg]
+    a = np.asarray(a).copy()
+    if a.dtype.kind == "f":
+        a[np.isnan(a)] = NEG_NAN
+    return a
+
+
+def check(args: tuple[str, str, str]) -> dict[str, Any]:
+    from simurg.formats.gtc import encode_gtc, read_gtc
+    from simurg.formats.idat import read_idat
+    from simurg.genotype.sample import GenotypeParams, genotype_sample
+    from simurg.genotype.sex import SexParams
+
+    gtc, grn, red = args
+    res: dict[str, Any] = {"gtc": gtc}
+    try:
+        ref = read_gtc(gtc, check_write_complete=False)
+        ctx = _ctx(ref.cluster_file)
+        params = GenotypeParams(
+            norm_version=_W["norm"], sex=SexParams(version=_W["sex_version"]), gtc_version=ref.version
+        )
+        g = genotype_sample(read_idat(grn), read_idat(red), ctx, params)
+        res["cluster_file"] = ctx.cluster_file_name
+        diffs: dict[str, Any] = {}
+        for k in ("raw_x", "raw_y", "genotypes", "genotype_scores", "base_calls", "controls_x", "controls_y"):
+            a, b = _canon(getattr(g, k)), _canon(getattr(ref, k))
+            if a.shape != b.shape:
+                diffs[k] = f"shape {a.shape} vs {b.shape}"
+            elif a.tobytes() != b.tobytes():
+                n = int((a.reshape(len(a), -1) != b.reshape(len(b), -1)).any(axis=1).sum())
+                diffs[k] = n
+        if g.normalization_transforms.tobytes() != ref.normalization_transforms.tobytes():
+            diffs["normalization_transforms"] = True
+        if ref.version >= 4:
+            for k in ("b_allele_freqs", "logr_ratios"):
+                if _canon(getattr(g, k)).tobytes() != _canon(getattr(ref, k)).tobytes():
+                    diffs[k] = int(
+                        (_canon(getattr(g, k)).view(np.uint32) != _canon(getattr(ref, k)).view(np.uint32)).sum()
+                    )
+        for k in (
+            "call_rate",
+            "gender",
+            "logr_dev",
+            "gc10",
+            "gc50",
+            "num_calls",
+            "num_no_calls",
+            "num_intensity_only",
+            "ploidy",
+            "ploidy_type",
+            "dx",
+            "slide_identifier",
+        ):
+            a, b = getattr(g, k), getattr(ref, k)
+            same = a == b or (isinstance(a, float) and np.isnan(a) and np.isnan(b))
+            if not same:
+                diffs[k] = [repr(a), repr(b)]
+        if "genotype_scores" in diffs:
+            d = np.abs(g.genotype_scores.astype(np.float64) - ref.genotype_scores.astype(np.float64))
+            res["max_score_abs_diff"] = float(np.nanmax(d))
+        res["field_diffs"] = diffs
+        # also evaluate the other sex-estimate generation on the same calls (cheap), to identify the vendor's rule
+        from simurg.genotype.sex import estimate_sex
+        from simurg.signal.normalization.apply import normalize, polar
+
+        nx, ny = normalize(g.raw_x, g.raw_y, ctx.lookup, g.normalization_transforms)
+        _, r_all = polar(nx, ny)
+        res["sex_by_version"] = {
+            v: estimate_sex(
+                g.genotypes, r_all, ctx.clusters.gentrain, ctx.clusters.present, ctx.sex_classes, SexParams(version=v)
+            )
+            for v in (1, 2)
+        }
+        res["sex_ref"] = ref.gender_char
+        for k in ALIGN:
+            setattr(g, k, getattr(ref, k))
+        for k in ("b_allele_freqs", "logr_ratios", "genotype_scores"):
+            setattr(g, k, _canon(getattr(g, k)))
+            setattr(ref, k, _canon(getattr(ref, k)))
+        res["byte_identical"] = encode_gtc(g) == encode_gtc(ref)
+        res["ref_layout_canonical"] = encode_gtc(ref) == Path(gtc).read_bytes() or None
+        res["metadata"] = {k: str(getattr(ref, k))[:80] for k in ("autocall_version", "cluster_file", "imaging_date")}
+    except Exception as e:
+        res["error"] = f"{type(e).__name__}: {e}"
+        res["trace"] = traceback.format_exc(limit=4)
+    return res
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--root", required=True, type=Path)
+    ap.add_argument("--bpm", required=True, type=Path)
+    ap.add_argument("--egt-dir", required=True, type=Path)
+    ap.add_argument("--out", required=True, type=Path)
+    ap.add_argument("--norm-version", default="1.1.2")
+    ap.add_argument("--sex-version", type=int, default=1)
+    ap.add_argument("--workers", type=int, default=os.cpu_count() or 1)
+    ap.add_argument("--limit", type=int, default=0)
+    a = ap.parse_args()
+
+    from simurg.formats.discovery import walk_files
+
+    t0 = time.time()
+    files = walk_files(a.root, (".gtc", ".idat"))
+    key_re = re.compile(r"(\d+_r\d{2}c\d{2})", re.IGNORECASE)
+    grn, red, gtcs = {}, {}, []
+    for p in files:
+        low = p.name.lower()
+        m = key_re.search(p.name)
+        if not m:
+            continue
+        k = m.group(1).lower()
+        if low.endswith("_grn.idat"):
+            grn[k] = str(p)
+        elif low.endswith("_red.idat"):
+            red[k] = str(p)
+        elif low.endswith(".gtc"):
+            gtcs.append((k, str(p)))
+    triples = [(g, grn[k], red[k]) for k, g in gtcs if k in grn and k in red]
+    if a.limit:
+        triples = triples[: a.limit]
+    with ProcessPoolExecutor(
+        a.workers,
+        mp_context=mp.get_context("spawn"),
+        initializer=_init,
+        initargs=(str(a.bpm), str(a.egt_dir), a.norm_version, a.sex_version),
+    ) as ex:
+        results = list(ex.map(check, triples, chunksize=2))
+
+    field_counts: dict[str, int] = {}
+    for r in results:
+        for k in r.get("field_diffs", {}):
+            field_counts[k] = field_counts.get(k, 0) + 1
+    summary: dict[str, Any] = {
+        "n_gtc": len(gtcs),
+        "n_compared": len(results),
+        "errors": [r for r in results if "error" in r][:20],
+        "byte_identical": sum(bool(r.get("byte_identical")) for r in results),
+        "samples_with_no_field_diff": sum(not r.get("field_diffs") and "error" not in r for r in results),
+        "samples_with_field_diff": field_counts,
+        "genotype_diff_loci_total": sum(
+            int(r.get("field_diffs", {}).get("genotypes", 0) or 0)
+            if isinstance(r.get("field_diffs", {}).get("genotypes", 0), int)
+            else 0
+            for r in results
+        ),
+        "score_diff_loci_total": sum(
+            int(r.get("field_diffs", {}).get("genotype_scores", 0) or 0)
+            if isinstance(r.get("field_diffs", {}).get("genotype_scores", 0), int)
+            else 0
+            for r in results
+        ),
+        "max_score_abs_diff": max((r.get("max_score_abs_diff", 0.0) for r in results), default=0.0),
+        "cluster_files": sorted({r.get("cluster_file", "") for r in results}),
+        "sex_agreement_by_version": {
+            str(v): sum(r.get("sex_by_version", {}).get(v) == r.get("sex_ref") for r in results) for v in (1, 2)
+        },
+        "ref_sex_counts": {s: sum(r.get("sex_ref") == s for r in results) for s in ("M", "F", "U")},
+        "examples": [r for r in results if r.get("field_diffs")][:15],
+        "params": {"norm_version": a.norm_version, "sex_version": a.sex_version},
+        "seconds": round(time.time() - t0, 1),
+    }
+    summary["passed"] = summary["byte_identical"] == len(results) > 0
+    a.out.parent.mkdir(parents=True, exist_ok=True)
+    a.out.write_text(json.dumps(summary, indent=2, default=str) + "\n")
+    print(
+        json.dumps(
+            {
+                k: summary[k]
+                for k in (
+                    "n_compared",
+                    "byte_identical",
+                    "samples_with_no_field_diff",
+                    "samples_with_field_diff",
+                    "genotype_diff_loci_total",
+                    "score_diff_loci_total",
+                    "passed",
+                    "seconds",
+                )
+            }
+        )
+    )
+    return 0 if summary["passed"] else 1
+
+
+if __name__ == "__main__":
+    sys.exit(main())

tests/unit/test_genotype.py

@@ -0,0 +1,235 @@
+"""Phase P4: calling, BAF/LRR, sex, base calls and `simurg genotype` on synthetic data (plan §7.5–§7.8)."""
+
+from __future__ import annotations
+
+import json
+from pathlib import Path
+
+import numpy as np
+import pytest
+from typer.testing import CliRunner
+
+from simurg.formats.bpm import read_bpm
+from simurg.formats.egt import read_egt
+from simurg.formats.gtc import read_gtc
+from simurg.genotype.baf_lrr import NEG_NAN, autosomal_mask, baf_lrr, logr_dev
+from simurg.genotype.clusters import ClusterArrays, align_clusters
+from simurg.genotype.gencall import call_genotypes, call_rate, smf, zmf
+from simurg.genotype.sample import base_calls, invariant_datetime, top_alleles
+from simurg.genotype.sex import SexParams, chromosome_classes, estimate_sex
+from simurg.synth.files import SynthIdat, SynthLocus, encode_bpm, encode_egt, encode_idat
+
+F32 = np.float32
+
+
+def clusters(n: int, *, gentrain: float = 0.8, threshold: float = np.nan) -> ClusterArrays:
+    tm = np.tile(np.array([0.05, 0.5, 0.95], np.float32), (n, 1))
+    rm = np.ones((n, 3), np.float32)
+    return ClusterArrays(
+        present=np.ones(n, bool),
+        n=np.full((n, 3), 100, np.int32),
+        theta_mean=tm,
+        theta_dev=np.full((n, 3), 0.03, np.float32),
+        r_mean=rm,
+        r_dev=np.full((n, 3), 0.1, np.float32),
+        intensity_threshold=np.full(n, threshold, np.float32),
+        gentrain=np.full(n, gentrain, np.float32),
+        mean_r=np.ones(n, np.float32),
+    )
+
+
+def test_fuzzy_memberships() -> None:
+    x = np.array([-1, 0, 0.25, 0.5, 0.75, 1, 2, np.nan], np.float32)
+    s = smf(x, F32(0), F32(1))
+    assert s.tolist()[:7] == [0.0, 0.0, 0.125, 0.5, 0.875, 1.0, 1.0]
+    assert s[7] == 1.0  # the reference returns 1 for NaN (all comparisons false)
+    z = zmf(x, F32(0), F32(1))
+    assert z.tolist()[:7] == [1.0, 1.0, 0.875, 0.5, 0.125, 0.0, 0.0] and z[7] == 0.0
+    assert smf(np.array([0.4, 0.6], np.float32), F32(1), F32(0)).tolist() == [0.0, 1.0]  # a >= b rule
+
+
+def test_calls_at_cluster_centres() -> None:
+    n = 6
+    t = np.array([0.05, 0.5, 0.95, 0.5, 0.5, 0.5], np.float32)
+    r = np.ones(n, np.float32)
+    io = np.array([0, 0, 0, 1, 0, 0], bool)
+    raw_zero = np.array([0, 0, 0, 0, 1, 0], bool)
+    c = clusters(n)
+    c.intensity_threshold[5] = 2.0  # R below the cluster file's intensity threshold -> no-call
+    calls = call_genotypes(t, r, raw_zero, c, io)
+    assert calls.genotypes.tolist() == [1, 2, 3, 0, 0, 0]
+    assert (calls.scores[:3] > 0.5).all() and calls.scores[4] == 0 and calls.scores[5] == 0
+    assert calls.num_calls == 3 and calls.num_intensity_only == 1
+    assert call_rate(calls, n) == F32(3 / 5)
+
+
+def test_low_gentrain_and_outliers_score_low() -> None:
+    n = 2
+    t = np.array([0.5, 0.3], np.float32)  # second point between clusters
+    calls_good = call_genotypes(t, np.ones(n, np.float32), np.zeros(n, bool), clusters(n), np.zeros(n, bool))
+    calls_bad = call_genotypes(
+        t, np.ones(n, np.float32), np.zeros(n, bool), clusters(n, gentrain=0.05), np.zeros(n, bool)
+    )
+    assert calls_good.scores[0] > calls_good.scores[1]
+    assert (calls_bad.scores < 0.15).all() and (calls_bad.genotypes == 0).all()
+
+
+def test_baf_lrr_interpolation() -> None:
+    n = 6
+    t = np.array([0.05, 0.5, 0.95, 0.275, 1.2, np.nan], np.float32)
+    r = np.array([1, 1, 2, 1, 1, np.nan], np.float32)
+    c = clusters(n)
+    io = np.array([0, 0, 0, 0, 0, 0], bool)
+    baf, lrr = baf_lrr(t, r, np.zeros(n, bool), c, io)
+    assert baf[:5].tolist() == pytest.approx([0.0, 0.5, 1.0, 0.25, 1.0])  # clamped at 1
+    assert lrr[:4].tolist() == pytest.approx([0.0, 0.0, 1.0, 0.0])
+    assert baf.view(np.uint32)[5] == 0xFFC00000 and lrr.view(np.uint32)[5] == 0xFFC00000
+    io[1] = True
+    c.mean_r[1] = 2.0
+    _, lrr2 = baf_lrr(t, r, np.zeros(n, bool), c, io)
+    assert lrr2[1] == pytest.approx(-1.0)  # intensity-only loci use the count-weighted mean R
+    assert np.isnan(NEG_NAN)
+
+
+def test_logr_dev_locus_set() -> None:
+    lrr = np.array([1, -1, 1, -1, 50, 50, np.nan], np.float32)
+    chroms = ["1", "2", "0", "22_KI270879v1_alt", "X", "MT", "3"]
+    auto = autosomal_mask(chroms)
+    assert auto.tolist() == [True, True, True, True, False, False, True]
+    assert logr_dev(lrr, auto, np.zeros(7, bool), np.ones(7, bool)) == pytest.approx(1.0)
+
+
+def test_sex_rules() -> None:
+    chroms = ["X"] * 30 + ["Y"] * 25 + ["1"] * 150
+    classes = chromosome_classes(chroms)
+    n = len(chroms)
+    gt = np.ones(n, np.uint8)
+    present = np.ones(n, bool)
+    gtrain = np.full(n, 0.8, np.float32)
+    r = np.ones(n, np.float32)
+    assert estimate_sex(gt, r, gtrain, present, classes) == "M"  # median Y R 1.0 > 0.3
+    r_f = r.copy()
+    r_f[30:55] = 0.05
+    assert estimate_sex(gt, r_f, gtrain, present, classes) == "F"  # X R 1.0 > 0.9
+    r_u = r_f.copy()
+    r_u[:30] = 0.5
+    assert estimate_sex(gt, r_u, gtrain, present, classes) == "U"
+    het = gt.copy()
+    het[:10] = 2  # X het rate 10/30 > 0.1
+    assert estimate_sex(het, r, gtrain, present, classes, SexParams(version=1)) == "F"
+    assert estimate_sex(gt, r, gtrain, present, classes, SexParams(version=1)) == "M"
+
+
+def test_base_calls_and_top_strand(tmp_path: Path) -> None:
+    loci = [
+        SynthLocus("a", snp="[A/G]", ilmn_strand="TOP"),
+        SynthLocus("b", address_a=2, snp="[T/C]", ilmn_strand="BOT"),
+        SynthLocus("c", address_a=3, snp="[N/A]", ilmn_strand="BOT"),
+    ]
+    (tmp_path / "m.bpm").write_bytes(encode_bpm(loci))
+    a, b = top_alleles(read_bpm(tmp_path / "m.bpm"))
+    assert bytes(a) == b"AA\x00" and bytes(b) == b"GGT"
+    bc = base_calls(np.array([1, 2, 0], np.uint8), a, b)
+    assert [bytes(x) for x in bc] == [b"AA", b"AG", b"--"]
+
+
+def test_invariant_datetime() -> None:
+    import datetime as dt
+
+    assert invariant_datetime(dt.datetime(2026, 1, 2, 0, 5, 9)) == "1/2/2026 12:05:09 AM"
+    assert invariant_datetime(dt.datetime(2026, 12, 31, 13, 0, 0)) == "12/31/2026 1:00:00 PM"
+
+
+def _synthetic_run(tmp_path: Path) -> dict[str, Path]:
+    """Two samples on a 600-locus synthetic array with a matching cluster file."""
+    rng = np.random.default_rng(4)
+    n = 600
+    chroms = ["1"] * 500 + ["X"] * 60 + ["Y"] * 40
+    loci = [SynthLocus(f"rs{i}", chrom=chroms[i], address_a=i + 1, pos=i + 1) for i in range(n)]
+    loci.append(SynthLocus("TGA_Ctrl_5716", address_a=n + 1))
+    (tmp_path / "m.bpm").write_bytes(encode_bpm(loci))
+    names = [lc.name for lc in loci]
+    tm = np.tile([0.03, 0.5, 0.97], (len(loci), 1))
+    rm = np.tile([1.0, 1.0, 1.0], (len(loci), 1))
+    (tmp_path / "c.egt").write_bytes(encode_egt(names, theta_mean=tm, r_mean=rm))
+    idats = tmp_path / "idats"
+    idats.mkdir()
+    addr = np.arange(1, len(loci) + 1, dtype=np.int32)
+    for pos in ("R01C01", "R02C01"):
+        g = rng.choice(3, len(loci))
+        inten = rng.uniform(0.6, 1.4, len(loci))
+        x = inten * np.select([g == 0, g == 1, g == 2], [1.0, 0.5, 0.02]) * 8000 + 300
+        y = inten * np.select([g == 0, g == 1, g == 2], [0.02, 0.5, 1.0]) * 9000 + 400
+        for ch, mean in (("Red", x), ("Grn", y)):
+            s = SynthIdat(
+                address_ids=addr,
+                mean=np.clip(mean, 0, 65535).astype(np.uint16),
+                n_beads=np.full(len(loci), 12, np.uint8),
+                barcode="200000000001",
+                position=pos,
+            )
+            (idats / f"200000000001_{pos}_{ch}.idat").write_bytes(encode_idat(s))
+    sheet = tmp_path / "s.csv"
+    sheet.write_text(
+        "[Data]\nSentrixBarcode_A,SentrixPosition_A,Sample_Name,Gender,Tissue\n"
+        "200000000001,R01C01,alpha,Male,blood\n200000000001,R02C01,beta,Female,saliva\n"
+        "200000000001,R09C01,ghost,Female,blood\n"
+    )
+    return {"bpm": tmp_path / "m.bpm", "egt": tmp_path / "c.egt", "idats": idats, "sheet": sheet}
+
+
+def test_genotype_cli_end_to_end(tmp_path: Path) -> None:
+    from simurg.cli.main import app
+
+    f = _synthetic_run(tmp_path)
+    runner = CliRunner()
+    args = [
+        "genotype",
+        "--quiet",
+        "--bpm",
+        str(f["bpm"]),
+        "--egt",
+        str(f["egt"]),
+        "--idats",
+        str(f["idats"]),
+        "--sheet",
+        str(f["sheet"]),
+        "--autocall-date",
+        "1/1/2026 12:00:00 PM",
+    ]
+    r1 = runner.invoke(app, [*args, "--out", str(tmp_path / "o1"), "--threads", "1"])
+    assert r1.exit_code == 3, r1.output  # 'ghost' has no IDATs
+    r2 = runner.invoke(app, [*args, "--out", str(tmp_path / "o2"), "--threads", "2", "--layout", "flat"])
+    assert r2.exit_code == 3
+    for stem in ("alpha", "beta"):
+        a = (tmp_path / "o1" / "samples" / stem / f"{stem}.gtc").read_bytes()
+        assert a == (tmp_path / "o2" / f"{stem}.gtc").read_bytes()  # thread count / layout never change bytes
+        g = read_gtc(tmp_path / "o1" / "samples" / stem / f"{stem}.gtc")
+        assert g.version == 5 and g.num_snps == 601 and g.autocall_version.startswith("SimurgArray ")
+        assert g.cluster_file == "c.egt" and g.snp_manifest == "m.bpm"
+        assert g.call_rate > 0.9
+    summary = json.loads((tmp_path / "o1" / "sample_summary.json").read_text())
+    assert summary["schema"] == "simurg.sample_summary/1"
+    cols = summary["columns"]
+    assert cols[:5] == ["sample_id", "sample_name", "sample_folder", "status", "fail_reason"]
+    assert "TGA_Ctrl_5716_norm_r" in cols and cols[-2:] == ["Gender", "Tissue"]
+    by = {s["sample_name"]: s for s in summary["samples"]}
+    assert by["ghost"]["status"] == "MISSING_INPUT" and by["alpha"]["status"] == "OK"
+    assert by["alpha"]["sex_reported"] == "M" and by["alpha"]["Tissue"] == "blood"
+    assert by["alpha"]["sex_check"] in ("match", "mismatch", "unknown")
+    tsv = (tmp_path / "o1" / "sample_summary.tsv").read_text().splitlines()
+    assert tsv[0].split("\t") == cols and len(tsv) == 4
+    r3 = runner.invoke(app, [*args, "--out", str(tmp_path / "o1")])
+    assert r3.exit_code == 2 and "E_OUTPUT_EXISTS" in r3.output
+
+
+def test_cluster_alignment_errors(tmp_path: Path) -> None:
+    f = _synthetic_run(tmp_path)
+    (tmp_path / "small.egt").write_bytes(encode_egt(["rs1"], theta_mean=[[0.1, 0.5, 0.9]], r_mean=[[1, 1, 1]]))
+    from simurg.core.errors import SimurgError
+
+    with pytest.raises(SimurgError) as e:
+        align_clusters(read_egt(tmp_path / "small.egt"), read_bpm(f["bpm"]))
+    assert e.value.code == "E_EGT_MANIFEST_MISMATCH"
+    c = align_clusters(read_egt(tmp_path / "small.egt"), read_bpm(f["bpm"]), allow_missing=True)
+    assert int(c.present.sum()) == 1