From fcc1f6c3e3e28153ab551db0e6986d5c1c7f5844 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Mon, 28 Sep 2026 00:07:37 +0100 Subject: [PATCH] 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) --- docs/messages.md | 1 + src/simurg/cli/main.py | 64 +++++- src/simurg/core/errors.py | 1 + src/simurg/genotype/baf_lrr.py | 77 +++++++ src/simurg/genotype/clusters.py | 75 +++++++ src/simurg/genotype/gencall.py | 270 ++++++++++++++++++++++++ src/simurg/genotype/sample.py | 196 ++++++++++++++++++ src/simurg/genotype/sex.py | 88 ++++++++ src/simurg/pipeline/genotype.py | 344 +++++++++++++++++++++++++++++++ src/simurg/pipeline/ingest.py | 61 +----- src/simurg/pipeline/inputs.py | 87 ++++++++ tests/oracle/test_p4_genotype.py | 59 ++++++ 12 files changed, 1270 insertions(+), 53 deletions(-) create mode 100644 src/simurg/genotype/baf_lrr.py create mode 100644 src/simurg/genotype/clusters.py create mode 100644 src/simurg/genotype/gencall.py create mode 100644 src/simurg/genotype/sample.py create mode 100644 src/simurg/genotype/sex.py create mode 100644 src/simurg/pipeline/genotype.py create mode 100644 src/simurg/pipeline/inputs.py create mode 100644 tests/oracle/test_p4_genotype.py diff --git a/docs/messages.md b/docs/messages.md index 87c3ff2..06f219a 100644 --- a/docs/messages.md +++ b/docs/messages.md @@ -12,6 +12,7 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `E_DUP_GTC` | error | Several GTC files resolve to the same sample | IN-030 | | `E_DUP_IDAT` | error | Several IDAT pairs resolve to the same sample | IN-022 | | `E_EGT_LOAD` | error | Cluster (EGT) file cannot be loaded | IN-044 | +| `E_EGT_MANIFEST_MISMATCH` | error | Cluster file lacks loci of the manifest | IN-043 | | `E_FILE_NOT_FOUND` | error | Input file or folder does not exist | CLI-016 | | `E_FILE_UNREADABLE` | error | Input file cannot be read | CLI-016 | | `E_GTC_CORRUPT` | warning | GTC file failed validation; sample skipped | IN-032 | diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index ef14a39..694e47b 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -1,6 +1,6 @@ """``simurg`` command-line interface (spec §2, plan §14). -Available: ``version``, ``inspect`` (P1), ``ingest`` (P2), ``normalize`` (P3). +Available: ``version``, ``inspect`` (P1), ``ingest`` (P2), ``normalize`` (P3), ``genotype`` (P4). Later phases add the other commands of CLI-001. """ @@ -274,5 +274,67 @@ def normalize_cmd( run_guarded(rep, lambda: run_normalize(opts, rep)) +# -------------------------------------------------------------------------------------------------------- +@app.command("genotype") +def genotype_cmd( + bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], + egt: Annotated[Path, typer.Option("--egt", help="EGT cluster file (required).")], + out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], + idats: Annotated[Path | None, typer.Option("--idats", help="Folder searched recursively for IDAT pairs.")] = None, + sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, + csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest, checked against the BPM.")] = None, + cutoff: Annotated[float, typer.Option("--cutoff", help="GenCall score below which a call is a no-call.")] = 0.15, + norm_version: Annotated[str, typer.Option("--norm-version", help="1.2.0 (current) or 1.1.2")] = "1.2.0", + sex_version: Annotated[int, typer.Option("--sex-version", min=1, max=2, help="Sex estimate rule (1 or 2).")] = 2, + gtc_version: Annotated[int, typer.Option("--gtc-version", help="GTC format version (3 or 5).")] = 5, + autocall_date: Annotated[ + str | None, typer.Option("--autocall-date", help="Fixed value for GTC field 201 (reproducible output).") + ] = None, + threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0, + force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, + layout: Annotated[str, typer.Option("--layout", help="nested|flat")] = "nested", + on_duplicate_idat: Annotated[str, typer.Option("--on-duplicate-idat", help="error|first|newest")] = "error", + log_json: Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")] = False, + debug: Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")] = False, + quiet: Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")] = False, +) -> None: + """IDAT pairs → GTC files and a sample summary (GenCall-equivalent calling; phase P4).""" + from simurg.pipeline.genotype import GenotypeOptions, run_genotype + + if layout not in ("nested", "flat"): + raise typer.BadParameter("--layout must be nested or flat") + if gtc_version not in (3, 5): + raise typer.BadParameter("--gtc-version must be 3 or 5") + if norm_version not in ("1.2.0", "1.1.2"): + raise typer.BadParameter("--norm-version must be 1.2.0 or 1.1.2") + rep = Reporter( + command="genotype", + level="debug" if debug else "info", + json_stderr=log_json, + debug=debug, + quiet=quiet, + log_file=out / "logs" / "genotype.log", + ) + opts = GenotypeOptions( + out=out, + bpm=bpm, + egt=egt, + idats=idats, + sheet=sheet, + csv=csv, + threads=threads, + force=force, + layout=layout, + on_duplicate_idat=on_duplicate_idat, + cutoff=cutoff, + norm_version=norm_version, + sex_version=sex_version, + gtc_version=gtc_version, + autocall_date=autocall_date, + argv=sys.argv[1:], + ) + run_guarded(rep, lambda: run_genotype(opts, rep)) + + if __name__ == "__main__": # pragma: no cover app() diff --git a/src/simurg/core/errors.py b/src/simurg/core/errors.py index e21bfd3..aa47d88 100644 --- a/src/simurg/core/errors.py +++ b/src/simurg/core/errors.py @@ -85,6 +85,7 @@ CATALOGUE: dict[str, MessageDef] = dict( _m("E_BPM_FORMAT", "error", "BPM manifest is malformed or unsupported", "IN-040"), _m("E_CSV_MANIFEST_FORMAT", "error", "CSV manifest is malformed", "IN-041"), _m("E_BPM_CSV_MISMATCH", "error", "BPM and CSV manifests describe different loci", "IN-042"), + _m("E_EGT_MANIFEST_MISMATCH", "error", "Cluster file lacks loci of the manifest", "IN-043"), _m("E_EGT_LOAD", "error", "Cluster (EGT) file cannot be loaded", "IN-044"), _m("W_EGT_NAME", "warning", "Cluster file name is not portable", "IN-045"), _m("W_MANIFEST_NO_REFSTRAND", "warning", "Manifest lacks RefStrand annotation", "IN-040"), diff --git a/src/simurg/genotype/baf_lrr.py b/src/simurg/genotype/baf_lrr.py new file mode 100644 index 0000000..f0646e1 --- /dev/null +++ b/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) diff --git a/src/simurg/genotype/clusters.py b/src/simurg/genotype/clusters.py new file mode 100644 index 0000000..8b54018 --- /dev/null +++ b/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, + ) diff --git a/src/simurg/genotype/gencall.py b/src/simurg/genotype/gencall.py new file mode 100644 index 0000000..ba96833 --- /dev/null +++ b/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)) diff --git a/src/simurg/genotype/sample.py b/src/simurg/genotype/sample.py new file mode 100644 index 0000000..1fdbe19 --- /dev/null +++ b/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 diff --git a/src/simurg/genotype/sex.py b/src/simurg/genotype/sex.py new file mode 100644 index 0000000..322ca42 --- /dev/null +++ b/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" diff --git a/src/simurg/pipeline/genotype.py b/src/simurg/pipeline/genotype.py new file mode 100644 index 0000000..1becc69 --- /dev/null +++ b/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 diff --git a/src/simurg/pipeline/ingest.py b/src/simurg/pipeline/ingest.py index 6cd3b10..2fd9e85 100644 --- a/src/simurg/pipeline/ingest.py +++ b/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) diff --git a/src/simurg/pipeline/inputs.py b/src/simurg/pipeline/inputs.py new file mode 100644 index 0000000..00bd82d --- /dev/null +++ b/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] diff --git a/tests/oracle/test_p4_genotype.py b/tests/oracle/test_p4_genotype.py new file mode 100644 index 0000000..d8c8257 --- /dev/null +++ b/tests/oracle/test_p4_genotype.py @@ -0,0 +1,59 @@ +"""Phase P4 acceptance on P0 data: GTCs produced from IDATs are byte-identical to oracle GTCs (gtc2vcf built without +FMA) once the four free-text metadata fields (cluster/manifest file names, AutoCall date and version) are aligned. +""" + +from __future__ import annotations + +import numpy as np +import pytest +from tests.conftest import GSA_BPM, GSA_EGT, GSA_IDATS, ORACLE_GTCS, ORACLE_GTCS_GT2, need + +from simurg.formats.bpm import read_bpm +from simurg.formats.egt import read_egt +from simurg.formats.gtc import encode_gtc, read_gtc +from simurg.formats.idat import read_idat +from simurg.genotype.clusters import align_clusters +from simurg.genotype.sample import GenotypeParams, genotype_sample, manifest_context + +pytestmark = pytest.mark.data +FREE_TEXT = ("cluster_file", "snp_manifest", "autocall_date", "autocall_version") +NEG_NAN = np.frombuffer(np.uint32(0xFFC00000).tobytes(), dtype=np.float32)[0] + + +def canonical_nans(g): # type: ignore[no-untyped-def] + """NaNs computed by libm carry the platform's default sign (+ on aarch64, − on x86-64); SimurgArray always + writes 0xFFC00000 like the vendor's x86 build. Compare everything else byte for byte (known_unknowns K-014).""" + for k in ("b_allele_freqs", "logr_ratios", "genotype_scores"): + a = getattr(g, k).copy() + a[np.isnan(a)] = NEG_NAN + setattr(g, k, a) + return g + + +@pytest.fixture(scope="module") +def ctx(): # type: ignore[no-untyped-def] + need(GSA_BPM, GSA_EGT) + m = read_bpm(GSA_BPM) + return manifest_context(m, align_clusters(read_egt(GSA_EGT), m)) + + +@pytest.mark.parametrize( + ("folder", "norm"), [(ORACLE_GTCS, "1.2.0"), (ORACLE_GTCS_GT2, "1.1.2")], ids=["1.2.0", "1.1.2"] +) +def test_gtc_byte_identical(ctx, folder, norm) -> None: # type: ignore[no-untyped-def] + need(folder, GSA_IDATS) + gtcs = sorted(folder.glob("*.gtc")) + assert len(gtcs) >= 18 + for f in gtcs: + sid = f.stem + chip = sid.split("_")[0] + ref = read_gtc(f) + g = genotype_sample( + read_idat(GSA_IDATS / chip / f"{sid}_Grn.idat"), + read_idat(GSA_IDATS / chip / f"{sid}_Red.idat"), + ctx, + GenotypeParams(norm_version=norm), + ) + for k in FREE_TEXT: + setattr(g, k, getattr(ref, k)) + assert encode_gtc(g) == encode_gtc(canonical_nans(ref)), sid