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/genotype/baf_lrr.py

@@ -0,0 +1,77 @@
+"""B-allele frequency, log R ratio and LogR deviation — plan §7.6 (Peiffer et al. 2006).
+
+Exact mode follows gtc2vcf ``get_baf_lrr``/``calculate_baf_lrr`` (MIT): float32 interpolation between the cluster
+θ/R means, ``log2`` in double of a float32 ratio, NaN where θ is undefined, BAF clamped to [0, 1] with NaN kept.
+"""
+
+from __future__ import annotations
+
+import math
+
+import numpy as np
+import numpy.typing as npt
+
+from simurg.genotype.clusters import ClusterArrays
+from simurg.genotype.gencall import _libm
+
+F32 = np.float32
+F32Arr = npt.NDArray[np.float32]
+NEG_NAN = np.frombuffer(np.uint32(0xFFC00000).tobytes(), dtype=np.float32)[0]  # C's -NAN
+
+
+def baf_lrr(
+    theta: F32Arr, r: F32Arr, raw_zero: npt.NDArray[np.bool_], c: ClusterArrays, intensity_only: npt.NDArray[np.bool_]
+) -> tuple[F32Arr, F32Arr]:
+    t = np.asarray(theta, np.float32)
+    rr = np.asarray(r, np.float32)
+    at, bt, ct = c.theta_mean[:, 0], c.theta_mean[:, 1], c.theta_mean[:, 2]
+    ar, br, cr = c.r_mean[:, 0], c.r_mean[:, 1], c.r_mean[:, 2]
+    with np.errstate(invalid="ignore", divide="ignore", over="ignore"):
+        eq = t == bt
+        lo = t < bt
+        hi = t > bt
+        r_lo = (ar + ((t - at) * (ar - br)) / (at - bt)).astype(np.float32)
+        b_lo = (((t - at) / (bt - at)) * F32(0.5)).astype(np.float32)
+        r_hi = (br + ((t - bt) * (cr - br)) / (ct - bt)).astype(np.float32)
+        b_hi = (F32(0.5) + ((t - bt) / (ct - bt)) * F32(0.5)).astype(np.float32)
+        r_ref = np.where(eq, br, np.where(lo, r_lo, np.where(hi, r_hi, np.nan))).astype(np.float32)
+        baf = np.where(eq, F32(0.5), np.where(lo, b_lo, np.where(hi, b_hi, NEG_NAN))).astype(np.float32)
+        denom = np.where(intensity_only, c.mean_r, r_ref).astype(np.float32)
+        ratio = (rr / denom).astype(np.float32).astype(np.float64)
+        lrr = _libm(math.log2, ratio).astype(np.float32)
+        lrr = np.where(rr != F32(0.0), lrr, -np.finfo(np.float32).max).astype(np.float32)
+        undefined = ~(eq | lo | hi)
+        lrr[undefined] = NEG_NAN
+        baf = np.where(baf < 0.0, F32(0.0), np.where(baf > 1.0, F32(1.0), baf)).astype(np.float32)
+    skip = raw_zero | ~c.present
+    baf[skip] = NEG_NAN
+    lrr[skip] = NEG_NAN
+    # canonical NaN bit pattern (0xFFC00000), as written by the reference on every platform
+    baf[np.isnan(baf)] = NEG_NAN
+    lrr[np.isnan(lrr)] = NEG_NAN
+    return baf, lrr
+
+
+def chrom_core(chrom: str) -> str:
+    return chrom[3:] if chrom[:3].upper() == "CHR" else chrom
+
+
+def autosomal_mask(chroms: list[str]) -> npt.NDArray[np.bool_]:
+    """LogR-dev locus set: first character (after an optional 'chr') not X, Y or M (so '0' and alt contigs count)."""
+    return np.array([(chrom_core(ch)[:1].upper() not in ("X", "Y", "M")) for ch in chroms], dtype=np.bool_)
+
+
+def logr_dev(
+    lrr: F32Arr, autosomal: npt.NDArray[np.bool_], intensity_only: npt.NDArray[np.bool_], present: npt.NDArray[np.bool_]
+) -> np.float32:
+    """Population SD of LRR over autosomal, non-intensity-only, finite values (sequential double sums)."""
+    sel = autosomal & ~intensity_only & present & np.isfinite(lrr)
+    v = lrr[sel].astype(np.float64)
+    count = len(v)
+    if count == 0:
+        return F32(np.nan)
+    s = float(np.cumsum(v)[-1])
+    s2 = float(np.cumsum(v * v)[-1])
+    mean = s / count
+    var = s2 / count - mean * mean
+    return F32(math.sqrt(var)) if var >= 0 else F32(np.nan)

src/simurg/genotype/clusters.py

@@ -0,0 +1,75 @@
+"""Cluster-file records aligned to manifest (GTC) locus order — input of genotype calling (plan §7.5–§7.6)."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+import numpy.typing as npt
+
+from simurg.core.errors import SimurgError
+from simurg.formats.bpm import Manifest
+from simurg.formats.egt import ClusterFile
+
+F32Arr = npt.NDArray[np.float32]
+
+
+@dataclass
+class ClusterArrays:
+    """Per-locus cluster parameters; (n, 3) arrays are ordered (AA, AB, BB)."""
+
+    present: npt.NDArray[np.bool_]  # locus has a record in the cluster file
+    n: npt.NDArray[np.int32]
+    theta_mean: F32Arr
+    theta_dev: F32Arr
+    r_mean: F32Arr
+    r_dev: F32Arr  # floored at 0.1 as the caller does (plan §7.5 step 2)
+    intensity_threshold: F32Arr  # NaN where the cluster file has none
+    gentrain: F32Arr
+    mean_r: F32Arr  # count-weighted mean R in float32 (expected R of intensity-only loci)
+
+    @property
+    def num_loci(self) -> int:
+        return len(self.present)
+
+
+def align_clusters(egt: ClusterFile, manifest: Manifest, *, allow_missing: bool = False) -> ClusterArrays:
+    """Join cluster records to the manifest by locus name (IN-043)."""
+    idx = egt.name_index()
+    rows = np.array([idx.get(nm, -1) for nm in manifest.names], dtype=np.int64)
+    present = rows >= 0
+    if not present.all() and not allow_missing:
+        missing = [manifest.names[i] for i in np.flatnonzero(~present)[:10]]
+        raise SimurgError(
+            "E_EGT_MANIFEST_MISMATCH",
+            f"{int((~present).sum())} manifest loci are missing from cluster file {egt.path}, "
+            f"e.g. {', '.join(missing)}",
+            hint="use the cluster file built for this manifest",
+        )
+    safe = np.where(present, rows, 0)
+
+    def take(a: npt.NDArray, fill: float = np.nan) -> npt.NDArray:  # type: ignore[type-arg]
+        out = a[safe].copy()
+        out[~present] = fill
+        return out
+
+    n = take(egt.n, 0).astype(np.int32)
+    r_mean = take(egt.r_mean).astype(np.float32)
+    r_dev = take(egt.r_dev).astype(np.float32)
+    r_dev = np.where(r_dev < np.float32(0.1), np.float32(0.1), r_dev).astype(np.float32)
+    # float32 weighted mean exactly as the reference: ((Naa*Raa + Nab*Rab) + Nbb*Rbb) / (Naa+Nab+Nbb)
+    nf = n.astype(np.float32)
+    with np.errstate(invalid="ignore", divide="ignore"):
+        num = (nf[:, 0] * r_mean[:, 0] + nf[:, 1] * r_mean[:, 1]) + nf[:, 2] * r_mean[:, 2]
+        mean_r = (num / (n.sum(axis=1)).astype(np.float32)).astype(np.float32)
+    return ClusterArrays(
+        present=present,
+        n=n,
+        theta_mean=take(egt.theta_mean).astype(np.float32),
+        theta_dev=take(egt.theta_dev).astype(np.float32),
+        r_mean=r_mean,
+        r_dev=r_dev,
+        intensity_threshold=take(egt.intensity_threshold).astype(np.float32),
+        gentrain=take(egt.gentrain_score).astype(np.float32),
+        mean_r=mean_r,
+    )

src/simurg/genotype/gencall.py

@@ -0,0 +1,270 @@
+"""GenCall-equivalent genotype calling from a cluster file — plan §7.5.
+
+Vectorised port of gtc2vcf ``gen_std_flair``/``modilik``/``compute_score_call_prelim``/``gencall_score_map``/
+``make_calls`` (MIT, © Giulio Genovese), which reproduce the vendor's GenTrain 6.x calling (concepts: Kermani
+patents, expired; Illumina GenCall technical note). Every expression keeps the reference's float/double types and
+evaluation order; scalar transcendental functions come from libm through ``math`` where NumPy's SIMD kernels could
+round differently.
+"""
+
+from __future__ import annotations
+
+import math
+from dataclasses import dataclass
+
+import numpy as np
+import numpy.typing as npt
+
+from simurg.genotype.clusters import ClusterArrays
+
+F32 = np.float32
+F32Arr = npt.NDArray[np.float32]
+F64Arr = npt.NDArray[np.float64]
+
+M_2_SQRTPI = 1.12837916709551257390
+M_SQRT1_2 = 0.707106781186547524401
+PEAK = 0.5 * M_2_SQRTPI * M_SQRT1_2  # constant sub-expression of the normal pdf (evaluated left to right)
+TINY_D = float(np.finfo(np.float64).tiny) * float(np.finfo(np.float64).eps)  # DBL_MIN * DBL_EPSILON
+DEFAULT_CUTOFF = 0.15
+
+
+def _c_call(fn, v: float) -> float:  # type: ignore[no-untyped-def]
+    """Call a ``math`` function with C semantics: domain/range errors give NaN/±inf instead of raising."""
+    if v != v:
+        return v
+    try:
+        return float(fn(v))
+    except ValueError:
+        if fn is math.log10 and v == 0.0:
+            return -math.inf
+        return math.nan
+    except OverflowError:
+        return math.inf
+
+
+def _libm(fn, x: F64Arr) -> F64Arr:  # type: ignore[no-untyped-def]
+    """Element-wise libm call (bit-identical to C)."""
+    return np.fromiter((_c_call(fn, v) for v in x.tolist()), dtype=np.float64, count=len(x))
+
+
+def _pow(e: float):  # type: ignore[no-untyped-def]
+    def f(v: float) -> float:
+        return math.pow(v, e)
+
+    return f
+
+
+def smf(x: F32Arr, a: F32Arr | np.float32, b: F32Arr | np.float32) -> F32Arr:
+    """MATLAB ``smf`` in float32, including the reference's ``a >= b`` rule."""
+    x = np.asarray(x, np.float32)
+    a = np.broadcast_to(np.asarray(a, np.float32), x.shape)
+    b = np.broadcast_to(np.asarray(b, np.float32), x.shape)
+    with np.errstate(invalid="ignore", divide="ignore", over="ignore"):
+        mid = ((a + b) / F32(2.0)).astype(np.float32)
+        u = ((x - a) / (b - a)).astype(np.float32)
+        v = ((x - b) / (b - a)).astype(np.float32)
+        rise = (F32(2.0) * (u * u)).astype(np.float32)
+        top = (F32(1.0) - F32(2.0) * (v * v)).astype(np.float32)
+        out = np.where(
+            x <= a, F32(0), np.where((a < x) & (x <= mid), rise, np.where((mid < x) & (x <= b), top, F32(1)))
+        ).astype(np.float32)
+        degenerate = a >= b
+        out = np.where(degenerate, np.where(x >= mid, F32(1), F32(0)), out).astype(np.float32)
+    return out
+
+
+def zmf(x: F32Arr, a: F32Arr | np.float32, b: F32Arr | np.float32) -> F32Arr:
+    """MATLAB ``zmf`` in float32 (callers guarantee a < b)."""
+    x = np.asarray(x, np.float32)
+    a = np.broadcast_to(np.asarray(a, np.float32), x.shape)
+    b = np.broadcast_to(np.asarray(b, np.float32), x.shape)
+    with np.errstate(invalid="ignore", divide="ignore", over="ignore"):
+        mid = ((a + b) / F32(2.0)).astype(np.float32)
+        u = ((x - a) / (b - a)).astype(np.float32)
+        v = ((x - b) / (b - a)).astype(np.float32)
+        fall = (F32(1.0) - F32(2.0) * (u * u)).astype(np.float32)
+        tail = (F32(2.0) * (v * v)).astype(np.float32)
+        return np.where(
+            x <= a, F32(1), np.where((a < x) & (x <= mid), fall, np.where((mid < x) & (x <= b), tail, F32(0)))
+        ).astype(np.float32)
+
+
+@dataclass
+class _Pdf:
+    """Normal pdf pieces for one cluster and one dimension: the full density and the flat peak."""
+
+    x: F32Arr
+    mu: F32Arr
+    full: F64Arr
+    peak: F64Arr
+    bad: npt.NDArray[np.bool_]
+
+    def get(self, side: str) -> F64Arr:
+        if side == "left":
+            out = np.where(self.x < self.mu, self.peak, self.full)
+        elif side == "right":
+            out = np.where(self.x > self.mu, self.peak, self.full)
+        else:
+            out = self.full
+        return np.where(self.bad, np.nan, out)
+
+
+def _pdf(x: F32Arr, mu: F32Arr, sigma: F32Arr) -> _Pdf:
+    """``exp(-0.5*((x-mu)/sigma)^2) * 0.5 * M_2_SQRTPI * M_SQRT1_2 / sigma`` with float32 z, double arithmetic."""
+    with np.errstate(invalid="ignore", divide="ignore", over="ignore"):
+        z = ((x - mu) / sigma).astype(np.float32).astype(np.float64)
+        e = _libm(math.exp, -0.5 * (z * z))
+        sd = sigma.astype(np.float64)
+        full = e * 0.5 * M_2_SQRTPI * M_SQRT1_2 / sd
+        peak = PEAK / sd
+    return _Pdf(x, mu, full, peak, sigma <= F32(0))
+
+
+@dataclass
+class Calls:
+    genotypes: npt.NDArray[np.uint8]  # 0 NC, 1 AA, 2 AB, 3 BB (after the cutoff)
+    scores: F32Arr  # GenCall scores (0 where not computable)
+    num_calls: int
+    num_intensity_only: int
+
+
+def call_genotypes(
+    theta: F32Arr,
+    r: F32Arr,
+    raw_zero: npt.NDArray[np.bool_],
+    c: ClusterArrays,
+    intensity_only: npt.NDArray[np.bool_],
+    cutoff: float = DEFAULT_CUTOFF,
+) -> Calls:
+    """Bayesian calls and GenCall scores for one sample (plan §7.5)."""
+    n = len(theta)
+    t = np.asarray(theta, np.float32)
+    rr = np.asarray(r, np.float32)
+    tm, td, rm, rd = c.theta_mean, c.theta_dev, c.r_mean, c.r_dev
+    with np.errstate(invalid="ignore", divide="ignore", over="ignore"):
+        # --- gen_std_flair ------------------------------------------------------------------------
+        z1 = (F32(0.5) * (tm[:, 1] - tm[:, 0]) / (td[:, 1] + td[:, 0])).astype(np.float32)
+        z2 = (F32(0.5) * (tm[:, 2] - tm[:, 1]) / (td[:, 2] + td[:, 1])).astype(np.float32)
+        mz = np.fmin(z1, z2)
+        third = F32(1.0) / F32(3.0)
+        alpha = np.fmax(third * mz, F32(1.0))
+        beta = np.fmin(third * mz, F32(1.0))
+        eta = (alpha * beta).astype(np.float32)
+        ftd = (td * eta[:, None]).astype(np.float32)
+        ftd = np.where(ftd < F32(0.02), F32(0.02), ftd).astype(np.float32)
+        med = np.fmax(np.fmin(td[:, 0], td[:, 1]), np.fmin(np.fmax(td[:, 0], td[:, 1]), td[:, 2]))
+        min_r = np.where(F32(0.2) < med, med, F32(0.2)).astype(np.float32)
+        min_rk = (rm / F32(7.0)).astype(np.float32)
+        min_rk = np.where(min_rk < min_r[:, None], min_r[:, None], min_rk).astype(np.float32)
+        frd = np.where(rd < min_rk, min_rk, rd).astype(np.float32)
+
+        # --- modilik: six ordered cases, later ones overwrite -----------------------------------------
+        d = [
+            100.0 * np.abs(t - tm[:, k]).astype(np.float64) + np.abs(rr - rm[:, k]).astype(np.float64) for k in range(3)
+        ]
+        daa, dab, dbb = d
+        ok_t = ~np.isnan(t)
+
+        pt = [_pdf(t, tm[:, k], ftd[:, k]) for k in range(3)]
+        pr = [_pdf(rr, rm[:, k], frd[:, k]) for k in range(3)]
+
+        def pdf2(k: int, st: str, sr: str) -> F64Arr:
+            return pt[k].get(st) * pr[k].get(sr)
+
+        L = np.zeros((n, 3), dtype=np.float64)
+        covered = np.zeros(n, dtype=np.bool_)
+        cases = [
+            ((daa <= dbb) & (dab <= dbb) & (rm[:, 0] <= rm[:, 1]), ((0, "left", "left"), (1, "right", "right"))),
+            ((daa <= dab) & (dbb <= dab) & (rm[:, 0] <= rm[:, 2]), ((0, "left", "left"), (2, "right", "right"))),
+            ((dab <= daa) & (dbb <= daa) & (rm[:, 1] <= rm[:, 2]), ((1, "left", "left"), (2, "right", "right"))),
+            ((daa <= dbb) & (dab <= dbb) & (rm[:, 0] > rm[:, 1]), ((0, "left", "right"), (1, "right", "left"))),
+            ((dab <= daa) & (dbb <= daa) & (rm[:, 1] > rm[:, 2]), ((1, "left", "right"), (2, "right", "left"))),
+            ((daa <= dab) & (dbb <= dab) & (rm[:, 0] > rm[:, 2]), ((0, "left", "right"), (2, "right", "left"))),
+        ]
+        for cond, assigns in cases:
+            m = cond & ok_t
+            if m.any():
+                for k, st, sr in assigns:
+                    L[m, k] = pdf2(k, st, sr)[m]
+                covered |= m
+        nc = ~covered
+        if nc.any():
+            for k in range(3):
+                L[nc, k] = pdf2(k, "full", "full")[nc]
+        L = np.where(np.isnan(L), 0.0, L)
+        L[~ok_t] = np.nan
+
+        # --- priors and posteriors -----------------------------------------------------------------
+        ntot = c.n.sum(axis=1)
+        P = (
+            (c.n.astype(np.float32) / ntot[:, None].astype(np.float32)).astype(np.float32).astype(np.float64) + 1.0
+        ) / 4.0
+        LP = L * P
+        ev = (LP[:, 0] + LP[:, 1]) + LP[:, 2]
+        ap = LP / ev[:, None]
+        aa, ab, bb = ap[:, 0], ap[:, 1], ap[:, 2]
+        iap = np.where(
+            (aa >= ab) & (aa >= bb), 1, np.where((ab >= aa) & (ab >= bb), 2, np.where((bb >= aa) & (bb >= ab), 3, 0))
+        ).astype(np.uint8)
+        mx = np.where(aa > ab, aa, ab)
+        scnd = np.where(aa > ab, ab, aa)
+        b_gt_mx = bb > mx
+        scnd = np.where(b_gt_mx, mx, np.where(bb > scnd, bb, scnd))
+        mx = np.where(b_gt_mx, bb, mx)
+        ap_lod = _libm(math.log10, mx / (TINY_D + scnd))
+        score_ap = smf(ap_lod.astype(np.float32), F32(0.0), F32(2.0)).astype(np.float32)
+
+        # --- sub-scores --------------------------------------------------------------------------
+        score_r1 = smf(rr, F32(0.0), F32(0.1))
+        nclus = (c.n > 0).sum(axis=1)
+        misclust = np.where(
+            (nclus == 1) & (c.n[:, 1] == 0), F32(0.7), np.where(nclus != 3, F32(0.95), F32(1.0))
+        ).astype(np.float32)
+        g = np.clip(iap.astype(np.int64) - 1, 0, 2)
+        gm_t, gd_t = tm[np.arange(n), g], ftd[np.arange(n), g]
+        gm_r, gd_r = rm[np.arange(n), g], frd[np.arange(n), g]
+        st_aa = zmf(t, gm_t + F32(2.0) * gd_t, gm_t + F32(6.0) * gd_t)
+        st_ab = zmf(np.abs((t - gm_t) / gd_t).astype(np.float32), F32(2.0), F32(6.0))
+        st_bb = smf(t, gm_t - F32(6.0) * gd_t, gm_t - F32(2.0) * gd_t)
+        score_t = np.where(iap == 1, st_aa, np.where(iap == 2, st_ab, np.where(iap == 3, st_bb, F32(1.0))))
+        r2 = smf(rr, F32(0.0), (gm_r / F32(10.0)).astype(np.float32))
+        r3 = smf(rr, gm_r - F32(6.0) * gd_r, gm_r - F32(2.0) * gd_r)
+        r4 = zmf(rr, gm_r + F32(6.0) * gd_r, gm_r + F32(12.0) * gd_r)
+        called = iap > 0
+        score_r = np.where(called, ((score_r1 * r2) * r3) * r4, F32(0.0)).astype(np.float32)
+        score_t = score_t.astype(np.float32)
+        score_n = np.ones(n, dtype=np.float32)
+        nan_gt = np.isnan(c.gentrain)
+        score_t[nan_gt] = np.nan
+        score_r[nan_gt] = np.nan
+        score_t[~ok_t] = np.nan
+        score_r[~ok_t] = np.nan
+        score_n[~ok_t] = np.nan
+        prelim = (((((score_ap * c.gentrain) * score_t) * score_r) * score_n) * misclust).astype(np.float32)
+
+        # intensity threshold: below it the call is NC with a NaN score
+        below = rr < c.intensity_threshold
+        prelim[below] = np.nan
+        iap[below] = 0
+
+        # --- gencall_score_map -------------------------------------------------------------------
+        x = prelim.astype(np.float64)
+        q = (x - 1.08) / 0.504
+        gbell = 1.0 / (1.0 + _libm(_pow(1.71), q * q))
+        score = (_libm(_pow(0.35), x) * gbell).astype(np.float32)
+    score = np.where(np.isnan(score), F32(0.0), score).astype(np.float32)
+
+    # loci without intensity or without a cluster record
+    no_call = raw_zero | ~c.present
+    score[no_call] = F32(0.0)
+    iap[no_call] = 0
+    geno = np.where(score < F32(cutoff), 0, iap).astype(np.uint8)
+    geno[intensity_only] = 0
+    return Calls(
+        genotypes=geno, scores=score, num_calls=int((geno > 0).sum()), num_intensity_only=int(intensity_only.sum())
+    )
+
+
+def call_rate(calls: Calls, num_snps: int) -> np.float32:
+    tiny = F32(np.finfo(np.float32).tiny) * F32(np.finfo(np.float32).eps)
+    return F32(F32(calls.num_calls) / ((F32(num_snps) - F32(calls.num_intensity_only)) + tiny))

src/simurg/genotype/sample.py

@@ -0,0 +1,196 @@
+"""One sample from IDAT pair to a complete GTC record — plan §7.1–§7.8, spec OUT-001/OUT-002.
+
+Order of operations follows the reference AutoCall pipeline (gtc2vcf ``gtc_init``, MIT): raw X/Y → normalization →
+calls → sex → BAF/LRR/LogR dev → GC10/GC50 → intensity percentiles → controls → sample metadata.
+"""
+
+from __future__ import annotations
+
+import datetime as _dt
+import re
+from dataclasses import dataclass, field
+from pathlib import Path
+
+import numpy as np
+import numpy.typing as npt
+
+from simurg import PRODUCT_NAME, __version__
+from simurg.formats.bpm import Manifest
+from simurg.formats.gtc import XFORM_DTYPE, Gtc, ScannerData
+from simurg.formats.idat import Idat
+from simurg.genotype.baf_lrr import autosomal_mask, baf_lrr, logr_dev
+from simurg.genotype.clusters import ClusterArrays
+from simurg.genotype.gencall import DEFAULT_CUTOFF, call_genotypes, call_rate
+from simurg.genotype.sex import SexParams, chromosome_classes, estimate_sex
+from simurg.signal.normalization.apply import normalize, polar
+from simurg.signal.normalization.fit import XForm, fit_transforms
+from simurg.signal.normalization.numerics import percentile
+from simurg.signal.raw_xy import LocusMapCache, assemble_raw_xy
+
+F32 = np.float32
+_COMPLEMENT = bytes.maketrans(b"ACGTDI", b"TGCADI")
+
+
+@dataclass
+class GenotypeParams:
+    cutoff: float = DEFAULT_CUTOFF
+    norm_version: str = "1.2.0"
+    sex: SexParams = field(default_factory=SexParams)
+    gtc_version: int = 5
+    autocall_version: str = f"{PRODUCT_NAME} {__version__}"
+    autocall_date: str | None = None  # None: run start time (invariant format, OUT-001)
+
+
+@dataclass
+class ManifestContext:
+    """Per-manifest data reused for every sample (computed once)."""
+
+    manifest: Manifest
+    clusters: ClusterArrays
+    norm_bins: npt.NDArray[np.uint8]
+    lookup: npt.NDArray[np.int32]
+    allele_a: npt.NDArray[np.uint8]
+    allele_b: npt.NDArray[np.uint8]
+    autosomal: npt.NDArray[np.bool_]
+    sex_classes: tuple[npt.NDArray[np.bool_], ...]
+    cache: LocusMapCache
+    cluster_file_name: str = ""
+
+
+def top_alleles(manifest: Manifest) -> tuple[npt.NDArray[np.uint8], npt.NDArray[np.uint8]]:
+    """TOP-strand A/B base per locus (gtc2vcf ``get_base_call``): SNP column on TOP designs, complemented otherwise.
+
+    Unknown bases complement to NUL, as in the reference's lookup table.
+    """
+    a = np.zeros(manifest.num_loci, dtype=np.uint8)
+    b = np.zeros(manifest.num_loci, dtype=np.uint8)
+    for i, (snp, strand) in enumerate(zip(manifest.snps, manifest.ilmn_strands, strict=True)):
+        raw = snp.encode("latin-1", errors="replace")
+        ca = raw[1:2] or b"\x00"
+        cb = raw[3:4] or b"\x00"
+        if strand[:1].upper() != "T":
+            ca = ca.translate(_COMPLEMENT) if ca in b"ACGTDI" else b"\x00"
+            cb = cb.translate(_COMPLEMENT) if cb in b"ACGTDI" else b"\x00"
+        a[i], b[i] = ca[0], cb[0]
+    return a, b
+
+
+def manifest_context(manifest: Manifest, clusters: ClusterArrays, cluster_file_name: str = "") -> ManifestContext:
+    a, b = top_alleles(manifest)
+    return ManifestContext(
+        manifest=manifest,
+        clusters=clusters,
+        norm_bins=manifest.norm_bin_ids(),
+        lookup=manifest.normalization_lookups(),
+        allele_a=a,
+        allele_b=b,
+        autosomal=autosomal_mask(manifest.chroms),
+        sex_classes=chromosome_classes(manifest.chroms),
+        cache=LocusMapCache(manifest),
+        cluster_file_name=cluster_file_name,
+    )
+
+
+def base_calls(
+    genotypes: npt.NDArray[np.uint8], a: npt.NDArray[np.uint8], b: npt.NDArray[np.uint8]
+) -> npt.NDArray[np.uint8]:
+    dash = np.uint8(ord("-"))
+    first = np.where(genotypes == 1, a, np.where(genotypes == 2, a, np.where(genotypes == 3, b, dash)))
+    second = np.where(genotypes == 1, a, np.where(genotypes == 2, b, np.where(genotypes == 3, b, dash)))
+    return np.stack([first, second], axis=1).astype(np.uint8)
+
+
+def invariant_datetime(t: _dt.datetime) -> str:
+    """``M/d/yyyy h:mm:ss tt`` with ASCII digits (OUT-001), independent of the process locale."""
+    hour = t.hour % 12 or 12
+    return f"{t.month}/{t.day}/{t.year} {hour}:{t.minute:02d}:{t.second:02d} {'AM' if t.hour < 12 else 'PM'}"
+
+
+def _param(pars: str, key: str) -> str:
+    m = re.search(rf"(?:^|\|){re.escape(key)}=([^|]*)", pars)
+    return m.group(1) if m else ""
+
+
+def xforms_to_records(xforms: list[XForm]) -> npt.NDArray[np.void]:
+    rec = np.zeros(len(xforms), dtype=XFORM_DTYPE)
+    for i, x in enumerate(xforms):
+        rec[i] = (
+            x.version,
+            x.offset_x,
+            x.offset_y,
+            x.scale_x,
+            x.scale_y,
+            x.shear,
+            x.theta,
+            (x.cvx, x.cvy, x.nn12, x.rr12, x.taa, x.tbb),
+        )
+    return rec
+
+
+def genotype_sample(
+    grn: Idat, red: Idat, ctx: ManifestContext, params: GenotypeParams, run_started: _dt.datetime | None = None
+) -> Gtc:
+    m, c = ctx.manifest, ctx.clusters
+    raw = assemble_raw_xy(grn, red, m, cache=ctx.cache)
+    xforms = fit_transforms(raw.raw_x, raw.raw_y, ctx.norm_bins, params.norm_version)
+    nx, ny = normalize(raw.raw_x, raw.raw_y, ctx.lookup, xforms)
+    theta, r = polar(nx, ny)
+    raw_zero = (raw.raw_x == 0) & (raw.raw_y == 0)
+    theta_c = np.where(raw_zero, np.float32(np.nan), theta).astype(np.float32)
+    r_c = np.where(raw_zero, np.float32(np.nan), r).astype(np.float32)
+
+    calls = call_genotypes(theta_c, r_c, raw_zero, c, m.intensity_only, params.cutoff)
+    sex = estimate_sex(calls.genotypes, r, c.gentrain, c.present, ctx.sex_classes, params.sex)
+    baf, lrr = baf_lrr(theta_c, r_c, raw_zero, c, m.intensity_only)
+    ldev = logr_dev(lrr, ctx.autosomal, m.intensity_only, c.present)
+
+    gs = np.sort(calls.scores[(calls.scores > F32(params.cutoff)) & ~m.intensity_only])
+    pct_x = [int(percentile(np.sort(raw.raw_x.astype(np.float32)), p)) for p in (5, 50, 95)]
+    pct_y = [int(percentile(np.sort(raw.raw_y.astype(np.float32)), p)) for p in (5, 50, 95)]
+
+    scan = red.scan_entry()
+    pars = scan.block_pars if scan else ""
+    started = run_started or _dt.datetime.now()
+    g = Gtc(
+        version=params.gtc_version,
+        num_snps=m.num_loci,
+        ploidy=2,
+        ploidy_type=1,
+        sample_name=red.sample_id or Path(grn.path).name.split("_Grn.idat")[0].split("_grn.idat")[0],
+        sample_plate=red.plate,
+        sample_well=red.well,
+        cluster_file=ctx.cluster_file_name,
+        snp_manifest=Path(m.path).name,
+        imaging_date=scan.run_time if scan else "",
+        autocall_date=params.autocall_date if params.autocall_date is not None else invariant_datetime(started),
+        autocall_version=params.autocall_version,
+        normalization_transforms=xforms_to_records(xforms),
+        controls_x=raw.controls_x,
+        controls_y=raw.controls_y,
+        raw_x=raw.raw_x,
+        raw_y=raw.raw_y,
+        genotypes=calls.genotypes,
+        base_calls=base_calls(calls.genotypes, ctx.allele_a, ctx.allele_b),
+        genotype_scores=calls.scores,
+        scanner_data=ScannerData(
+            name=_param(pars, "sherlockID"),
+            pmt_green=int(_param(pars, "PMTGainCY3") or 0).to_bytes(4, "little", signed=True),
+            pmt_red=int(_param(pars, "PMTGainCY5") or 0).to_bytes(4, "little", signed=True),
+            version=scan.code_version if scan else "",
+            user=_param(pars, "Username"),
+        ),
+        call_rate=float(call_rate(calls, m.num_loci)),
+        gender=sex.encode(),
+        logr_dev=float(ldev),
+        gc10=float(percentile(gs, 10)),
+        b_allele_freqs=baf,
+        logr_ratios=lrr,
+        percentiles_x=np.array(pct_x, dtype=np.uint16),
+        percentiles_y=np.array(pct_y, dtype=np.uint16),
+        slide_identifier=red.barcode,
+    )
+    g.sample_data["p50gc"] = percentile(gs, 50)
+    g.sample_data["num_calls"] = calls.num_calls
+    g.sample_data["num_no_calls"] = m.num_loci - calls.num_intensity_only - calls.num_calls
+    g.sample_data["num_intensity_only"] = calls.num_intensity_only
+    return g

src/simurg/genotype/sex.py

@@ -0,0 +1,88 @@
+"""Sex estimate from X heterozygosity and X/Y intensities — plan §7.8.
+
+Port of gtc2vcf ``estimate_gender`` (MIT), which reproduces the vendor's two generations:
+
+* version 1 (AutoConvert 1.x): male unless the X heterozygosity rate among called X loci exceeds 0.1;
+* version 2 (AutoConvert 2.0, IAAP, current CLI): median normalized R of Y loci > 0.3 → M; otherwise median R of
+  X loci > 0.9 → F (or F when there are too few X loci), else U. Falls back to version 1 when there are too few Y or
+  autosomal loci.
+
+Only loci whose cluster record has a non-zero GenTrain score count.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+import numpy.typing as npt
+
+from simurg.genotype.baf_lrr import chrom_core
+from simurg.signal.normalization.numerics import median_sorted
+
+F32 = np.float32
+
+
+@dataclass(frozen=True)
+class SexParams:
+    version: int = 2
+    min_loci: int = 100
+    min_x_loci: int = 20
+    min_y_loci: int = 20
+    call_rate_threshold: float = 0.0
+    y_threshold: float = 0.3
+    x_threshold: float = 0.9
+    x_het_rate_threshold: float = 0.1
+
+
+def chromosome_classes(chroms: list[str]) -> tuple[npt.NDArray[np.bool_], ...]:
+    core = [chrom_core(c) for c in chroms]
+    is_x = np.array([c == "X" for c in core])
+    is_y = np.array([c == "Y" for c in core])
+    is_auto = np.array([c not in ("X", "Y", "XY", "MT") for c in core])
+    return is_x, is_y, is_auto
+
+
+def estimate_sex(
+    genotypes: npt.NDArray[np.uint8],
+    r_all: npt.NDArray[np.float32],
+    gentrain: npt.NDArray[np.float32],
+    present: npt.NDArray[np.bool_],
+    classes: tuple[npt.NDArray[np.bool_], ...],
+    p: SexParams | None = None,
+) -> str:
+    """``r_all``: normalized R of every locus, computed even where raw intensities are zero (as the reference)."""
+    p = p or SexParams()
+    is_x, is_y, is_auto = classes
+
+    def thr(v: float) -> float:  # thresholds are C floats promoted to double in the comparisons
+        return float(F32(v))
+
+    use = present & (gentrain != F32(0.0))  # NaN != 0 is True, as in C
+    x = use & is_x
+    y = use & is_y
+    auto = use & is_auto & ~is_x & ~is_y
+    x_count, y_count, auto_count = int(x.sum()), int(y.sum()), int(auto.sum())
+    called = genotypes > 0
+    x_non_missing = int((x & called).sum())
+    x_hets = int((x & (genotypes == 2)).sum())
+    auto_non_missing = int((auto & called).sum())
+
+    if p.version == 1 or y_count < p.min_y_loci or auto_count < p.min_loci:
+        if x_non_missing > p.min_x_loci:
+            rate = float(F32(F32(x_hets) / F32(x_non_missing)))
+            return "M" if not rate > thr(p.x_het_rate_threshold) else "F"
+        return "U"
+    if auto_count > 0 and auto_non_missing / auto_count > thr(p.call_rate_threshold):
+        ry = r_all[y].astype(np.float32)
+        ry = np.where(np.isfinite(ry), ry, F32(0.0)).astype(np.float32)
+        y_med = median_sorted(np.sort(ry))
+        if float(y_med) > thr(p.y_threshold):
+            return "M"
+        if x_count < p.min_x_loci:
+            return "F"
+        rx = r_all[x].astype(np.float32)
+        rx = np.where(np.isfinite(rx), rx, F32(0.0)).astype(np.float32)
+        x_med = median_sorted(np.sort(rx))
+        return "F" if float(x_med) > thr(p.x_threshold) else "U"
+    return "U"