From 3c33ae5cb15bf158c86f0c27b1c714a0ecbdaa55 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Mon, 28 Sep 2026 17:00:37 +0100 Subject: [PATCH] P7: differential-testing harness (simurg compare) and nightly run on the Pi MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - compare/core: EXACT/NEAR/DISCREPANCY/UNKNOWN with the plan §16.3 tolerances (abs, float32 ULP), NaN = NaN, locus classes (assay, variant type, chromosome group, intensity-only); fields absent on both sides are not applicable - compare/records: GTC vs GTC (per-locus arrays, transforms per bin, controls, sample statistics; metadata listed, not counted) and VCF vs VCF (records keyed incl. ID; SNV/MNV/indel, merged, multi-allelic classes) - simurg compare A B (files or folders paired by name): summary.json, discrepancies.tsv.gz (never truncated), self-contained dashboard.html; --tolerance, --fail-on-discrepancy (exit 3) - scripts/nightly.py + install_nightly.sh (crontab 02:30): unit tests, GSA vs gtc2vcf (GTC + VCF), VCF regression vs golden, Omni2.5 subset vs vendor GTCs; first run PASS on every step - tests/unit/test_compare.py; docs/validation/P7.md, README --- README.md | 5 +- docs/validation/P7.md | 38 ++++ scripts/install_nightly.sh | 19 ++ scripts/nightly.py | 235 +++++++++++++++++++++ src/simurg/cli/main.py | 55 +++++ src/simurg/compare/__init__.py | 0 src/simurg/compare/core.py | 223 ++++++++++++++++++++ src/simurg/compare/records.py | 223 ++++++++++++++++++++ src/simurg/pipeline/compare.py | 359 +++++++++++++++++++++++++++++++++ tests/unit/test_compare.py | 112 ++++++++++ 10 files changed, 1268 insertions(+), 1 deletion(-) create mode 100644 docs/validation/P7.md create mode 100755 scripts/install_nightly.sh create mode 100644 scripts/nightly.py create mode 100644 src/simurg/compare/__init__.py create mode 100644 src/simurg/compare/core.py create mode 100644 src/simurg/compare/records.py create mode 100644 src/simurg/pipeline/compare.py create mode 100644 tests/unit/test_compare.py diff --git a/README.md b/README.md index 5328132..1a99fb6 100644 --- a/README.md +++ b/README.md @@ -21,9 +21,12 @@ sites, per-sample VCF.gz + tabix or a cohort BCF + CSI, optional INFO/intensity auxiliary loci) and `simurg export plink|bedgraph|finalreport|locus-summary` (PLINK 1 `.bed` and PLINK 2 `.pgen`). On the GSA demo, genotypes match gtc2vcf and GTCtoVCF with zero unexplained differences and PLINK files match plink2 (`docs/validation/P6.md`). +Phase **P7 (differential testing)** is complete: `simurg compare` (EXACT/NEAR/DISCREPANCY/UNKNOWN per field and +locus class, every discrepancy listed, HTML dashboard) and a nightly run on the Pi (GSA vs gtc2vcf, Omni2.5 subset vs +vendor GTCs, VCF regression) — `docs/validation/P7.md`. Available commands: `simurg version`, `simurg inspect [--json]`, `simurg ingest`, `simurg normalize`, `simurg genotype`, `simurg qc controls|samples|contamination|relatedness|loci|report|suggest-thresholds`, `simurg gtc-to-vcf`, -`simurg export plink|bedgraph|finalreport|locus-summary`. +`simurg export plink|bedgraph|finalreport|locus-summary`, `simurg compare`. ## Quick start diff --git a/docs/validation/P7.md b/docs/validation/P7.md new file mode 100644 index 0000000..beeb39b --- /dev/null +++ b/docs/validation/P7.md @@ -0,0 +1,38 @@ +# Phase P7 validation — differential-testing harness + +Date: 2026-09-28. All work on the Raspberry Pi (no cloud cost). + +## Delivered + +| Item | Module | +|---|---| +| Tolerance categories EXACT / NEAR / DISCREPANCY / UNKNOWN with the plan §16.3 tolerances; float32 ULP distance; NaN = NaN regardless of sign (K-014) | `compare/core.py` | +| Locus classes (assay, variant type, chromosome group incl. PAR/MT/unmapped, intensity-only) | `compare/core.py` | +| GTC-vs-GTC (every per-locus array, transforms per bin, controls, sample statistics; free-text metadata reported but not counted) and VCF-vs-VCF (per record and FORMAT field, records keyed by CHROM/POS/REF/ALT/ID, classes SNV/MNV/indel, merged sites, multi-allelic) | `compare/records.py` | +| `simurg compare A B` (files or folders paired by name; `--tolerance FIELD=…`, `--fail-on-discrepancy`) → `summary.json`, `discrepancies.tsv.gz` (every discrepancy and unmatched record, never truncated), self-contained `dashboard.html` | `pipeline/compare.py` | +| Nightly run: unit tests, GSA genotype vs gtc2vcf oracle, GSA VCF vs gtc2vcf (semantic) and vs golden (regression), Omni2.5 subset vs vendor GTCs; per-night page + history; 14 nights kept | `scripts/nightly.py`, `scripts/install_nightly.sh` (crontab, 02:30) | + +## Acceptance (plan §18 P7: nightly run on P0 + Omni subset produces report) + +First complete nightly run on the Pi (2026-09-28, 17 min): + +| Step | Compared | Result | +|---|---|---| +| unit | 206 unit tests | PASS | +| gsa_gtc | 36 GSA GTCs from IDATs vs gtc2vcf oracle GTCs, 20 fields, per locus class | **PASS**: 0 discrepancies, 0 unknown (free-text metadata listed separately) | +| gsa_vcf | 36 VCFs vs gtc2vcf (semantic, `validate_p6_vcf.py`) | **PASS**: 0 unexplained | +| gsa_vcf regression | 36 VCFs vs the golden set of the first run: 23 640 768 records | all records present, GT/GS/BAF/LRR EXACT (see note on BAF/LRR of merged sites) | +| omni_gtc | 6 Omni2.5 samples (`--preset autocall-1.6.2`) vs the vendor's AutoCall 1.6.2.2 GTCs, 18 fields | **PASS**: 0 discrepancies, 0 unknown | + +Installed as a crontab entry (every night 02:30, `flock` against overlapping runs; `bash scripts/install_nightly.sh +--remove` removes it). History page: `data/nightly/index.html`. + +## Notes + +* The harness found its own bug while being built: VCF records keyed without ID let duplicate designs overwrite + each other in `--unsquash-duplicates` output; keys now include the ID and repeated keys are counted, never merged. +* The first regression run reported 51 048 UNKNOWN values: BAF/LRR of multi-allelic merged records, which carry + `GT:GS` only on *both* sides. A field absent on both sides is now "not applicable" and not counted (GTC and VCF + alike); a field absent on one side stays UNKNOWN. +* Resource use on the Pi: one Omni2.5 sample genotypes in 93 s at 3.4 GB peak RSS (single process); the GSA VCF + step memory-maps the 3.1 GB genome cache, so its address-space cap is 12 GB (page cache, not RAM). diff --git a/scripts/install_nightly.sh b/scripts/install_nightly.sh new file mode 100755 index 0000000..26ed679 --- /dev/null +++ b/scripts/install_nightly.sh @@ -0,0 +1,19 @@ +#!/usr/bin/env bash +# Install (or remove) the nightly differential-testing run in the user's crontab (no root needed). +# bash scripts/install_nightly.sh # add: every night at 02:30 +# bash scripts/install_nightly.sh --remove # remove the entry +set -euo pipefail +cd "$(dirname "$0")/.." +REPO=$(pwd) +TAG="# simurg-nightly" +LINE="30 2 * * * cd $REPO && flock -n /tmp/simurg-nightly.lock .venv/bin/python scripts/nightly.py >> data/nightly/cron.log 2>&1 $TAG" +current=$(crontab -l 2>/dev/null | grep -v -- "$TAG" || true) +if [[ "${1:-}" == "--remove" ]]; then + printf '%s\n' "$current" | sed '/^$/d' | crontab - + echo "nightly entry removed" + exit 0 +fi +mkdir -p data/nightly +printf '%s\n%s\n' "$current" "$LINE" | sed '/^$/d' | crontab - +echo "installed: $LINE" +crontab -l | grep -- "$TAG" diff --git a/scripts/nightly.py b/scripts/nightly.py new file mode 100644 index 0000000..0bde8c4 --- /dev/null +++ b/scripts/nightly.py @@ -0,0 +1,235 @@ +#!/usr/bin/env python3 +"""Nightly differential-testing run (plan §18 P7): P0 GSA demo + Omni2.5 subset, every result compared with its +oracle by ``simurg compare``, one index page per night plus a history. + +Steps (each isolated; a failing step does not stop the others): + unit pytest tests/unit (memory-capped) + gsa_gtc simurg genotype on the 36 GSA IDATs → compare with the gtc2vcf oracle GTCs (expect PASS) + gsa_vcf simurg gtc-to-vcf --unsquash-duplicates → semantic check vs gtc2vcf (scripts/validate_p6_vcf.py, + expect 0 unexplained) and regression vs the golden VCFs of an earlier night (expect PASS) + omni_gtc simurg genotype --preset autocall-1.6.2 on the local Omni2.5 subset → compare with the vendor GTCs + +Output: data/nightly//{index.html, summary.json, /...}, data/nightly/index.html (history). +Usage: scripts/nightly.py [--root DIR] [--keep 14] [--steps unit,gsa_gtc,gsa_vcf,omni_gtc] [--update-golden] +""" + +from __future__ import annotations + +import argparse +import datetime as dt +import html +import json +import shutil +import subprocess +import sys +import time +from pathlib import Path +from typing import Any + +ROOT = Path(__file__).resolve().parents[1] +PY = ROOT / ".venv" / "bin" / "python" +SIMURG = ROOT / ".venv" / "bin" / "simurg" +GSA = ROOT / "data/p0/gsa_v4_pgx" +GSA_BPM = GSA / "GSA-PGx-48v4-0_20079540_E2_bpm/GSA-PGx-48v4-0_20079540_E2.bpm" +GSA_CSV = GSA / "GSA-PGx-48v4-0_20079540_E2_csv/GSA-PGx-48v4-0_20079540_E2.csv" +GSA_EGT = GSA / "GSA-PGx-48v4-0_20079540_E2_ClusterFile/GSA-PGx-48v4-0_20079540_E2_ClusterFile.egt" +GSA_IDATS = GSA / "idats_Demo_36/idats_Demo_36" +ORACLE_GTCS = ROOT / "data/p0/oracle/gsa_gt3_nofma" +ORACLE_VCF = ROOT / "data/p6/oracle/gsa36_gtc2vcf.vcf.gz" +FASTA = ROOT / "data/p6/GRCh38.primary.fa" +BCFTOOLS = ROOT / "refs/tools/bcftools-1.22/bcftools" +OMNI = ROOT / "data/p0/omni_head" +MEM_CAP = "--as=5500000000" # 5.5 GB address-space cap per step (the Pi has 8 GB) + + +def run(cmd: list[str], log: Path, timeout: int = 7200) -> tuple[int, float]: + t = time.time() + with log.open("a") as fh: + fh.write(f"$ {' '.join(cmd)}\n") + fh.flush() + try: + rc = subprocess.run(cmd, stdout=fh, stderr=subprocess.STDOUT, cwd=ROOT, timeout=timeout).returncode + except subprocess.TimeoutExpired: + fh.write("TIMEOUT\n") + rc = 124 + return rc, round(time.time() - t, 1) + + +def capped(cmd: list[str], cap: str = MEM_CAP) -> list[str]: + return ["prlimit", cap, *cmd] + + +# gtc-to-vcf memory-maps the whole genome cache (3.1 GB of file-backed address space that costs page cache, not RAM), +# so its address-space cap must be larger +VCF_CAP = "--as=12000000000" + + +def compare_verdict(folder: Path) -> dict[str, Any]: + p = folder / "summary.json" + if not p.exists(): + return {"verdict": "ERROR"} + s = json.loads(p.read_text()) + return {"verdict": s["verdict"], "discrepancies": s["discrepancies"], "unknown": s["unknown"], + "report": str((folder / "dashboard.html").relative_to(folder.parents[1]))} # fmt: skip + + +def step_unit(out: Path, log: Path) -> dict[str, Any]: + rc, sec = run(capped([str(PY), "-m", "pytest", "-q", "-m", "not slow", "tests/unit"]), log, 3600) + return {"verdict": "PASS" if rc == 0 else "FAIL", "seconds": sec} + + +def step_gsa_gtc(out: Path, log: Path) -> dict[str, Any]: + gt = out / "gt" + rc, sec = run(capped([str(SIMURG), "genotype", "--bpm", str(GSA_BPM), "--egt", str(GSA_EGT), "--idats", + str(GSA_IDATS), "--out", str(gt), "--threads", "2", "--quiet", "--force"]), log) # fmt: skip + if rc not in (0,): + return {"verdict": "ERROR", "seconds": sec, "exit": rc} + _, sec2 = run(capped([str(SIMURG), "compare", str(gt / "samples"), str(ORACLE_GTCS), "--out", str(out / "compare"), + "--bpm", str(GSA_BPM), "--threads", "2", "--title", "GSA GTCs vs gtc2vcf oracle", + "--quiet", "--force"]), log) # fmt: skip + return {**compare_verdict(out / "compare"), "seconds": sec + sec2} + + +def step_gsa_vcf(out: Path, log: Path, golden: Path, update_golden: bool) -> dict[str, Any]: + gt = out.parent / "gsa_gtc" / "gt" / "samples" + if not gt.exists(): + return {"verdict": "SKIPPED", "reason": "gsa_gtc did not run"} + vcf = out / "vcf" + rc, sec = run(capped([str(SIMURG), "gtc-to-vcf", "--bpm", str(GSA_BPM), "--csv", str(GSA_CSV), "--fasta", + str(FASTA), "--gtcs", str(gt), "--out", str(vcf), "--unsquash-duplicates", "--threads", "2", + "--quiet", "--force"], VCF_CAP), log) # fmt: skip + if rc != 0: + return {"verdict": "ERROR", "seconds": sec, "exit": rc} + res: dict[str, Any] = {"seconds": sec} + if ORACLE_VCF.exists() and BCFTOOLS.exists(): + _, sec2 = run(capped([str(PY), "scripts/validate_p6_vcf.py", "--ours", str(vcf), "--oracle", str(ORACLE_VCF), + "--fasta", str(FASTA), "--bcftools", str(BCFTOOLS), "--bpm", str(GSA_BPM), + "--out", str(out / "vs_gtc2vcf.json")]), log) # fmt: skip + d = json.loads((out / "vs_gtc2vcf.json").read_text()) if (out / "vs_gtc2vcf.json").exists() else {} + res["vs_gtc2vcf_unexplained"] = d.get("unexplained", "ERROR") + res["seconds"] += sec2 + if golden.exists() and not update_golden: + run(capped([str(SIMURG), "compare", str(golden), str(vcf), "--out", str(out / "regression"), "--threads", "2", + "--title", "GSA VCFs vs golden", "--quiet", "--force"]), log) # fmt: skip + res["regression"] = compare_verdict(out / "regression") + else: + if golden.exists(): + shutil.rmtree(golden) + golden.mkdir(parents=True) + for f in vcf.glob("*.snv.vcf.gz"): + shutil.copy2(f, golden / f.name) + res["regression"] = {"verdict": "GOLDEN_CREATED"} + ok = res.get("vs_gtc2vcf_unexplained", 0) == 0 and res["regression"]["verdict"] in ("PASS", "GOLDEN_CREATED") + res["verdict"] = "PASS" if ok else "FAIL" + return res + + +def step_omni_gtc(out: Path, log: Path) -> dict[str, Any]: + bpm, egt = OMNI / "HumanOmni2.5-4v1-Multi_B.bpm", OMNI / "HumanOmni2.5-4v1-Multi_B.egt" + if not bpm.exists(): + return {"verdict": "SKIPPED", "reason": "no local Omni subset"} + idats, vendor = out / "idats", out / "vendor" + for d in (idats, vendor): + d.mkdir(parents=True, exist_ok=True) + for f in sorted(OMNI.rglob("*.idat")): + (idats / f.name).unlink(missing_ok=True) + (idats / f.name).symlink_to(f) + for f in sorted(OMNI.rglob("*.gtc")): + (vendor / f.name).unlink(missing_ok=True) + (vendor / f.name).symlink_to(f) + gt = out / "gt" + rc, sec = run(capped([str(SIMURG), "genotype", "--bpm", str(bpm), "--egt", str(egt), "--idats", str(idats), + "--out", str(gt), "--preset", "autocall-1.6.2", "--threads", "1", "--no-cohort-qc", + "--quiet", "--force"]), log) # fmt: skip + if rc != 0: + return {"verdict": "ERROR", "seconds": sec, "exit": rc} + run(capped([str(SIMURG), "compare", str(gt / "samples"), str(vendor), "--out", str(out / "compare"), "--bpm", + str(bpm), "--threads", "1", "--title", "Omni2.5 subset vs vendor GTCs (AutoCall 1.6.2.2)", + "--quiet", "--force"]), log) # fmt: skip + return {**compare_verdict(out / "compare"), "seconds": sec} + + +CSS = ( + "body{font:14px system-ui;margin:16px} td,th{border:1px solid #ccc;padding:3px 8px}" + " table{border-collapse:collapse} .PASS td:nth-child(2){background:#d6e6f5}" + " .FAIL td:nth-child(2),.ERROR td:nth-child(2),.DISCREPANCIES td:nth-child(2){background:#f9d3ac}" +) + + +def render(night: dict[str, Any]) -> str: + esc = html.escape + rows = [] + for k, r in night["steps"].items(): + details = json.dumps({x: y for x, y in r.items() if x not in ("verdict", "seconds", "report")}) + link = f'dashboard' if r.get("report") else "" + verdict = esc(str(r.get("verdict"))) + rows.append( + f"{esc(k)}{verdict}{r.get('seconds', '')}" + f"{esc(details)}{link}" + ) + return ( + "" + f"Nightly {night['date']}" + f"

Nightly {night['date']} — {esc(night['overall'])}

" + f"

{esc(night['version'])} · {night['seconds']} s

" + "" + + "".join(rows) + + "
stepverdictsecondsdetails
" + ) + + +def main() -> int: + ap = argparse.ArgumentParser() + ap.add_argument("--root", type=Path, default=ROOT / "data" / "nightly") + ap.add_argument("--keep", type=int, default=14, help="nights to keep") + ap.add_argument("--steps", default="unit,gsa_gtc,gsa_vcf,omni_gtc") + ap.add_argument("--update-golden", action="store_true") + a = ap.parse_args() + day = dt.date.today().isoformat() + base = a.root / day + base.mkdir(parents=True, exist_ok=True) + log = base / "run.log" + t0 = time.time() + version = subprocess.run([str(SIMURG), "version"], capture_output=True, text=True, cwd=ROOT).stdout.strip() + steps: dict[str, Any] = {} + for name in a.steps.split(","): + out = base / name + out.mkdir(exist_ok=True) + try: + if name == "unit": + steps[name] = step_unit(out, log) + elif name == "gsa_gtc": + steps[name] = step_gsa_gtc(out, log) + elif name == "gsa_vcf": + steps[name] = step_gsa_vcf(out, log, a.root / "golden" / "gsa_vcf", a.update_golden) + elif name == "omni_gtc": + steps[name] = step_omni_gtc(out, log) + except Exception as e: # a crashing step must not stop the others + steps[name] = {"verdict": "ERROR", "error": f"{type(e).__name__}: {e}"} + overall = "PASS" if all(s.get("verdict") in ("PASS", "SKIPPED") for s in steps.values()) else "FAIL" + night = {"date": day, "version": version, "overall": overall, "seconds": round(time.time() - t0), "steps": steps} + (base / "summary.json").write_text(json.dumps(night, indent=2) + "\n") + (base / "index.html").write_text(render(night)) + # large intermediates are not kept (reports and summaries are) + for heavy in ("gsa_gtc/gt/samples", "gsa_vcf/vcf", "omni_gtc/gt/samples"): + shutil.rmtree(base / heavy, ignore_errors=True) + nights = sorted(p for p in a.root.iterdir() if p.is_dir() and len(p.name) == 10 and p.name[4] == "-") + for old in nights[: max(0, len(nights) - a.keep)]: + shutil.rmtree(old, ignore_errors=True) + hist = [] + for p in sorted(a.root.glob("????-??-??/summary.json"), reverse=True)[:30]: + s = json.loads(p.read_text()) + cells = "".join(f"{html.escape(str(s['steps'].get(k, {}).get('verdict', '')))}" for k in steps) + hist.append(f"{s['date']}{s['overall']}{cells}") + head = "".join(f"{html.escape(k)}" for k in steps) + (a.root / "index.html").write_text( + f"SimurgArray nightly" + "

SimurgArray nightly history

" + f"{head}{''.join(hist)}
nightoverall
" + ) + print(json.dumps({"overall": overall, **{k: v.get("verdict") for k, v in steps.items()}})) + return 0 if overall == "PASS" else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index 9e9fe82..8d42c9f 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -444,6 +444,61 @@ def gtc_to_vcf_cmd( run_guarded(rep, lambda: run_gtc_to_vcf(opts, rep)) +# -------------------------------------------------------------------------------------------------------- +@app.command("compare") +def compare_cmd( + a: Annotated[Path, typer.Argument(help="GTC/VCF file or folder (A).")], + b: Annotated[Path, typer.Argument(help="GTC/VCF file or folder (B), paired with A by relative path.")], + out: Annotated[Path, typer.Option("--out", help="Output folder (summary.json, discrepancies, dashboard).")], + bpm: Annotated[Path | None, typer.Option("--bpm", help="Manifest, for per-locus-class counts of GTCs.")] = None, + tolerance: Annotated[ + list[str] | None, + typer.Option("--tolerance", help="FIELD=exact|abs=X|ulp=N (repeatable); VCF floats: float=abs=X"), + ] = None, + sample: Annotated[ + list[str] | None, typer.Option("--sample", help="VCF sample(s) to compare (default: all shared).") + ] = None, + title: Annotated[str, typer.Option("--title", help="Dashboard title.")] = "", + threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0, + fail_on_discrepancy: Annotated[ + bool, typer.Option("--fail-on-discrepancy", help="Exit 3 when any DISCREPANCY is found.") + ] = False, + force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, + log_json: Annotated[bool, typer.Option("--log-json", "--json-log")] = False, + debug: Annotated[bool, typer.Option("--debug")] = False, + quiet: Annotated[bool, typer.Option("--quiet", "--no-progress")] = False, +) -> None: + """Field-level differential test of two GTCs/VCFs (or folders): EXACT/NEAR/DISCREPANCY/UNKNOWN (plan §16.3).""" + from simurg.pipeline.compare import CompareOptions, run_compare + + rep = Reporter( + command="compare", level="debug" if debug else "info", json_stderr=log_json, debug=debug, quiet=quiet + ) + + def body() -> int: + tol: dict[str, str] = {} + for item in tolerance or []: + key, sep, val = item.partition("=") + if not sep: + raise SimurgError("E_CONFIG_KEY", f"--tolerance expects FIELD=VALUE, got {item!r}") + tol[key.strip()] = val.strip() + opts = CompareOptions( + a=a, + b=b, + out=out, + bpm=bpm, + tolerances=tol, + samples=list(sample) if sample else None, + threads=threads, + fail_on_discrepancy=fail_on_discrepancy, + title=title, + force=force, + ) + return run_compare(opts, rep) + + run_guarded(rep, body) + + # -------------------------------------------------------------------------------------------------------- export_app = typer.Typer(help="Exports from GTCs or VCFs: bedgraph, finalreport, plink, locus-summary.") app.add_typer(export_app, name="export") diff --git a/src/simurg/compare/__init__.py b/src/simurg/compare/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/simurg/compare/core.py b/src/simurg/compare/core.py new file mode 100644 index 0000000..8a2dc11 --- /dev/null +++ b/src/simurg/compare/core.py @@ -0,0 +1,223 @@ +"""Field-level comparison with tolerance categories (spec CLI-038; plan §16.3). + +Every compared value falls in exactly one category: + +* **EXACT** — bitwise equal (NaN equals NaN whatever its sign bit, K-014); +* **NEAR** — within the field's tolerance; +* **DISCREPANCY** — outside it (or NaN on one side only); +* **UNKNOWN** — one side lacks the field or the record. + +Counts are kept per locus class (assay, variant type, chromosome group, intensity-only), and every discrepancy is +kept with its identifiers — nothing is hidden (plan §16.3). +""" + +from __future__ import annotations + +from collections.abc import Mapping, Sequence +from dataclasses import dataclass +from dataclasses import field as dc_field +from typing import Any + +import numpy as np +import numpy.typing as npt + +EXACT, NEAR, DISCREPANCY, UNKNOWN = 0, 1, 2, 3 +CATEGORY_NAMES = ("EXACT", "NEAR", "DISCREPANCY", "UNKNOWN") + + +@dataclass(frozen=True) +class Tolerance: + abs: float | None = None # absolute difference allowed + ulp: int | None = None # float32 units in the last place allowed + + @property + def exact(self) -> bool: + return self.abs is None and self.ulp is None + + +EXACT_TOL = Tolerance() +# plan §16.3 (bit-exact mode targets) +GTC_TOLERANCES: dict[str, Tolerance] = { + "raw_x": EXACT_TOL, + "raw_y": EXACT_TOL, + "controls_x": EXACT_TOL, + "controls_y": EXACT_TOL, + "genotypes": EXACT_TOL, + "base_calls": EXACT_TOL, + "normalization_transforms": Tolerance(ulp=2), + "genotype_scores": Tolerance(abs=1e-5), + "b_allele_freqs": Tolerance(abs=1e-5), + "logr_ratios": Tolerance(abs=1e-5), + "call_rate": Tolerance(abs=1e-6), + "logr_dev": Tolerance(abs=1e-6), + "gc10": Tolerance(abs=1e-6), + "gc50": Tolerance(abs=1e-6), + "gender": EXACT_TOL, + "num_calls": EXACT_TOL, + "num_no_calls": EXACT_TOL, + "num_intensity_only": EXACT_TOL, + "percentiles_x": EXACT_TOL, + "percentiles_y": EXACT_TOL, +} +VCF_FLOAT_TOL = Tolerance(abs=1e-5) # GS, BAF, LRR, NORMX, NORMY, THETA, R +VCF_EXACT_FIELDS = {"GT", "X", "Y"} + + +def _float_ulps(a: npt.NDArray[np.float32], b: npt.NDArray[np.float32]) -> npt.NDArray[np.int64]: + """Distance in float32 ULPs (monotone integer mapping of the bit patterns).""" + + def key(x: npt.NDArray[np.float32]) -> npt.NDArray[np.int64]: + i = x.astype(np.float32).view(np.int32).astype(np.int64) + return np.where(i < 0, -(i & 0x7FFFFFFF), i) + + return np.abs(key(a) - key(b)) + + +def categorize(a: npt.ArrayLike, b: npt.ArrayLike, tol: Tolerance) -> npt.NDArray[np.int8]: + """Element-wise category of two equally shaped arrays (numbers or strings).""" + a, b = np.asarray(a), np.asarray(b) + if a.shape != b.shape: + raise ValueError(f"shape {a.shape} vs {b.shape}") + out = np.full(a.shape, DISCREPANCY, np.int8) + if a.dtype.kind in "fc" or b.dtype.kind in "fc": + fa, fb = a.astype(np.float64), b.astype(np.float64) + na, nb = np.isnan(fa), np.isnan(fb) + same = (fa == fb) | (na & nb) + if a.dtype == np.float32 and b.dtype == np.float32: + same &= (a.view(np.uint32) == b.view(np.uint32)) | (na & nb) | ((fa == 0) & (fb == 0)) + out[same] = EXACT + near = np.zeros(a.shape, np.bool_) + valid = ~(na | nb) & ~same + if tol.abs is not None: + near |= valid & (np.abs(fa - fb) <= tol.abs) + if tol.ulp is not None: + near |= valid & (_float_ulps(a.astype(np.float32), b.astype(np.float32)) <= tol.ulp) + out[near] = NEAR + return out + out[a == b] = EXACT + return out + + +@dataclass +class FieldResult: + """Category counts of one field: overall and per locus class, plus the discrepancies themselves.""" + + field: str + tolerance: str + counts: dict[str, list[int]] = dc_field(default_factory=dict) # class -> [EXACT, NEAR, DISCREPANCY, UNKNOWN] + discrepancies: list[dict[str, Any]] = dc_field(default_factory=list) + + def add(self, cls: str, cats: npt.NDArray[np.int8]) -> None: + c = self.counts.setdefault(cls, [0, 0, 0, 0]) + binc = np.bincount(cats.ravel().astype(np.int64), minlength=4) + for k in range(4): + c[k] += int(binc[k]) + + @property + def total(self) -> list[int]: + return self.counts.get("all", [0, 0, 0, 0]) + + +def tol_text(t: Tolerance) -> str: + if t.exact: + return "exact" + parts = [] + if t.abs is not None: + parts.append(f"abs<={t.abs:g}") + if t.ulp is not None: + parts.append(f"ulp<={t.ulp}") + return ",".join(parts) + + +def locus_classes( + chroms: Sequence[str], snps: Sequence[str], assay_types: npt.ArrayLike, intensity_only: npt.ArrayLike +) -> dict[str, npt.NDArray[np.bool_]]: + """Locus class masks (plan §16.3): assay, variant type, chromosome group, intensity-only.""" + at = np.asarray(assay_types) + io = np.asarray(intensity_only, dtype=np.bool_) + + def group(c: str) -> str: + c = c[3:] if c.lower().startswith("chr") else c + if c in ("0", ""): + return "unmapped" + if c == "XY": + return "PAR" + if c in ("X", "Y"): + return c + if c in ("MT", "M"): + return "MT" + return "autosome" if c.isdigit() else "other" + + groups = np.array([group(c) for c in chroms]) + indel = np.array([("D" in s or "I" in s) for s in snps]) + out = {"all": np.ones(len(groups), np.bool_), "assay=InfII": at == 0, "assay=InfI": at != 0, + "type=SNV": ~indel, "type=indel": indel, "intensity_only": io} # fmt: skip + for gname in ("autosome", "X", "Y", "PAR", "MT", "unmapped", "other"): + m = groups == gname + if m.any(): + out[f"chrom={gname}"] = m + return out + + +def compare_locus_field( + name: str, + a: npt.ArrayLike | None, + b: npt.ArrayLike | None, + tol: Tolerance, + classes: Mapping[str, npt.NDArray[np.bool_]], + ids: Sequence[str], + context: Mapping[str, str], +) -> FieldResult: + """Per-locus arrays (first axis = locus). Missing or empty on one side → UNKNOWN for every locus.""" + res = FieldResult(name, tol_text(tol)) + n = len(classes["all"]) + a_ok = a is not None and len(np.asarray(a)) == n + b_ok = b is not None and len(np.asarray(b)) == n + if not (a_ok and b_ok): + cats = np.full(n, UNKNOWN, np.int8) + else: + el = categorize(np.asarray(a), np.asarray(b), tol) + cats = el.max(axis=tuple(range(1, el.ndim))) if el.ndim > 1 else el # worst element of the locus + for i in np.flatnonzero(cats == DISCREPANCY).tolist(): + res.discrepancies.append( + {**context, "field": name, "locus": ids[i], "a": _fmt(np.asarray(a)[i]), "b": _fmt(np.asarray(b)[i])} + ) + for cls, mask in classes.items(): + res.add(cls, cats[mask]) + return res + + +def compare_scalar(name: str, a: Any, b: Any, tol: Tolerance, context: Mapping[str, str]) -> FieldResult: + res = FieldResult(name, tol_text(tol)) + if a is None or b is None: + cat = np.array([UNKNOWN], np.int8) + else: + aa, bb = np.asarray([a]), np.asarray([b]) + if aa.dtype.kind == "f" or bb.dtype.kind == "f": + aa, bb = aa.astype(np.float64), bb.astype(np.float64) + cat = categorize(aa, bb, tol) + res.add("sample", cat) + if cat[0] == DISCREPANCY: + res.discrepancies.append({**context, "field": name, "locus": "", "a": _fmt(a), "b": _fmt(b)}) + return res + + +def _fmt(v: Any) -> str: + arr = np.asarray(v) + if arr.ndim == 0: + x = arr.item() + return repr(x) if isinstance(x, float) else str(x) + return "[" + ",".join(_fmt(x) for x in arr.ravel().tolist()) + "]" + + +def merge_results(results: Sequence[FieldResult]) -> dict[str, FieldResult]: + """Sum results of the same field (e.g. over samples).""" + out: dict[str, FieldResult] = {} + for r in results: + m = out.setdefault(r.field, FieldResult(r.field, r.tolerance)) + for cls, c in r.counts.items(): + cur = m.counts.setdefault(cls, [0, 0, 0, 0]) + for k in range(4): + cur[k] += c[k] + m.discrepancies.extend(r.discrepancies) + return out diff --git a/src/simurg/compare/records.py b/src/simurg/compare/records.py new file mode 100644 index 0000000..60a02a5 --- /dev/null +++ b/src/simurg/compare/records.py @@ -0,0 +1,223 @@ +"""GTC-vs-GTC and VCF-vs-VCF comparison on parsed data (spec CLI-038; plan §16.3).""" + +from __future__ import annotations + +import math +from collections.abc import Mapping, Sequence +from dataclasses import dataclass, field +from typing import Any + +import numpy as np +import numpy.typing as npt +from numpy.lib import recfunctions as rfn + +from simurg.compare.core import ( + DISCREPANCY, + EXACT, + GTC_TOLERANCES, + UNKNOWN, + VCF_EXACT_FIELDS, + VCF_FLOAT_TOL, + FieldResult, + Tolerance, + categorize, + compare_locus_field, + compare_scalar, + tol_text, +) + +GTC_LOCUS_FIELDS = ("raw_x", "raw_y", "genotypes", "base_calls", "genotype_scores", "b_allele_freqs", "logr_ratios") +GTC_SCALARS = ("call_rate", "logr_dev", "gc10", "gc50", "gender", "num_calls", "num_no_calls", "num_intensity_only") +GTC_METADATA = ("sample_name", "sample_plate", "sample_well", "cluster_file", "snp_manifest", "imaging_date", + "autocall_date", "autocall_version", "slide_identifier") # fmt: skip + + +@dataclass +class Comparison: + """Results of one pair (one sample).""" + + label: str + kind: str # "gtc" | "vcf" + fields: list[FieldResult] = field(default_factory=list) + metadata: dict[str, list[str]] = field(default_factory=dict) # field -> [a, b] when they differ + notes: list[str] = field(default_factory=list) + + +def compare_gtc( + a: Any, + b: Any, + label: str, + classes: Mapping[str, npt.NDArray[np.bool_]], + names: Sequence[str], + tolerances: Mapping[str, Tolerance] = GTC_TOLERANCES, +) -> Comparison: + """Two parsed GTCs of the same manifest.""" + res = Comparison(label, "gtc") + ctx = {"pair": label} + if a.num_snps != b.num_snps: + res.notes.append(f"locus count differs: {a.num_snps} vs {b.num_snps}") + return res + for f in GTC_LOCUS_FIELDS: + va, vb = getattr(a, f, None), getattr(b, f, None) + if (va is None or len(va) == 0) and (vb is None or len(vb) == 0): + continue # absent on both sides (e.g. BAF/LRR in GTC v3): not applicable + res.fields.append(compare_locus_field(f, va, vb, tolerances[f], classes, names, ctx)) + # normalization transforms: one row per bin + ta, tb = a.normalization_transforms, b.normalization_transforms + fr = FieldResult("normalization_transforms", tol_text(tolerances["normalization_transforms"])) + if len(ta) != len(tb) or not len(ta): + fr.add("bins", np.full(max(len(ta), len(tb), 1), UNKNOWN, np.int8)) + else: + ua = rfn.structured_to_unstructured(ta).astype(np.float32) + ub = rfn.structured_to_unstructured(tb).astype(np.float32) + cats = categorize(ua, ub, tolerances["normalization_transforms"]).max(axis=1) + fr.add("bins", cats) + for i in np.flatnonzero(cats == DISCREPANCY).tolist(): + fr.discrepancies.append({**ctx, "field": "normalization_transforms", "locus": f"bin{i}", + "a": str(ua[i].tolist()), "b": str(ub[i].tolist())}) # fmt: skip + res.fields.append(fr) + # controls: one row per control bead + for f in ("controls_x", "controls_y"): + ca, cb = getattr(a, f), getattr(b, f) + fr = FieldResult(f, tol_text(tolerances[f])) + fr.add("controls", categorize(ca, cb, tolerances[f]) if len(ca) == len(cb) else np.array([UNKNOWN], np.int8)) + res.fields.append(fr) + for f in GTC_SCALARS: + res.fields.append(compare_scalar(f, _scalar(a, f), _scalar(b, f), tolerances[f], ctx)) + for f in ("percentiles_x", "percentiles_y"): + pa, pb = getattr(a, f), getattr(b, f) + fr = FieldResult(f, tol_text(tolerances[f])) + fr.add("sample", categorize(pa, pb, tolerances[f]) if len(pa) == len(pb) else np.array([UNKNOWN], np.int8)) + res.fields.append(fr) + for f in GTC_METADATA: + va, vb = str(getattr(a, f, "")), str(getattr(b, f, "")) + if va != vb: + res.metadata[f] = [va, vb] + if a.version != b.version: + res.notes.append(f"GTC versions differ: {a.version} vs {b.version} (fields absent in one are UNKNOWN)") + return res + + +def _scalar(g: Any, f: str) -> Any: + if f == "gender": + return g.gender_char + v = getattr(g, f, None) + if isinstance(v, float) and math.isnan(v): + return v + return v + + +# ------------------------------------------------------------------------------------------------ VCF +@dataclass +class VcfRecord: + chrom: str + pos: int + vid: str + ref: str + alt: str + fields: dict[str, str] # FORMAT key -> value of the compared sample + + +def vcf_classes(r: VcfRecord) -> list[str]: + c = r.chrom[3:] if r.chrom.lower().startswith("chr") else r.chrom + grp = "autosome" if c.isdigit() else {"X": "X", "Y": "Y", "MT": "MT", "M": "MT"}.get(c, "other") + alts = r.alt.split(",") + if all(len(x) == len(r.ref) == 1 for x in alts): + vtype = "SNV" + elif all(len(x) == len(r.ref) for x in alts): + vtype = "MNV" + else: + vtype = "indel" + out = ["all", f"chrom={grp}", f"type={vtype}"] + if "," in r.vid or ";" in r.vid: + out.append("merged_site") + if len(alts) > 1: + out.append("multiallelic") + return out + + +def _gt_norm(gt: str) -> str: + parts = gt.replace("|", "/").split("/") + if any(p == "." for p in parts): + return "." if len(parts) == 1 else "./." + return "/".join(sorted(parts, key=int)) + + +def _cat_scalar(va: str, vb: str, tol: Tolerance) -> int: + if va == vb: + return EXACT + x, y = _num(va), _num(vb) + if math.isnan(x) and math.isnan(y): + return EXACT + if math.isnan(x) or math.isnan(y): + return DISCREPANCY + if x == y: + return EXACT + if tol.abs is not None and abs(x - y) <= tol.abs: + return 1 # NEAR + return DISCREPANCY + + +def compare_vcf_records( + a: Mapping[tuple[str, int, str, str, str], VcfRecord], + b: Mapping[tuple[str, int, str, str, str], VcfRecord], + label: str, + float_tol: Tolerance = VCF_FLOAT_TOL, +) -> Comparison: + """Records keyed by (CHROM, POS, REF, ALT, ID); FORMAT fields of one sample on each side.""" + res = Comparison(label, "vcf") + ctx = {"pair": label} + keys = list(a) + [k for k in b if k not in a] + fields = sorted({f for r in a.values() for f in r.fields} | {f for r in b.values() for f in r.fields}) + counts: dict[str, dict[str, list[int]]] = {f: {} for f in ["record", *fields]} + disc: dict[str, list[dict[str, Any]]] = {f: [] for f in ["record", *fields]} + tols = {f: (Tolerance() if f in VCF_EXACT_FIELDS else float_tol) for f in fields} + + def bump(f: str, classes: list[str], cat: int) -> None: + per = counts[f] + for cls in classes: + c = per.get(cls) + if c is None: + c = per[cls] = [0, 0, 0, 0] + c[cat] += 1 + + for k in keys: + ra, rb = a.get(k), b.get(k) + r = ra or rb + assert r is not None + classes = vcf_classes(r) + where = f"{r.chrom}:{r.pos}:{r.ref}>{r.alt} {r.vid}" + if ra is None or rb is None: + bump("record", classes, UNKNOWN) + disc["record"].append({**ctx, "field": "record", "locus": where, + "a": "present" if ra else "absent", "b": "present" if rb else "absent"}) # fmt: skip + continue + bump("record", classes, EXACT) + for f in fields: + va, vb = ra.fields.get(f), rb.fields.get(f) + if va is None and vb is None: + continue # the field does not apply to this record on either side (e.g. BAF of merged sites) + if va is None or vb is None: + cat = UNKNOWN + elif f == "GT": + cat = EXACT if _gt_norm(va) == _gt_norm(vb) else DISCREPANCY + else: + cat = _cat_scalar(va, vb, tols[f]) + bump(f, classes, cat) + if cat == DISCREPANCY: + disc[f].append({**ctx, "field": f, "locus": where, "a": va, "b": vb}) + for f in ["record", *fields]: + fr = FieldResult( + f, "present in both" if f == "record" else ("exact" if f in VCF_EXACT_FIELDS else tol_text(float_tol)) + ) + fr.counts = counts[f] + fr.discrepancies = disc[f] + res.fields.append(fr) + return res + + +def _num(v: str) -> float: + try: + return float(v) + except ValueError: + return math.nan diff --git a/src/simurg/pipeline/compare.py b/src/simurg/pipeline/compare.py new file mode 100644 index 0000000..089364b --- /dev/null +++ b/src/simurg/pipeline/compare.py @@ -0,0 +1,359 @@ +"""``simurg compare A B``: field-level differential testing of GTCs or VCFs (spec CLI-038; plan §14, §16.3, §18 P7). + +A and B are two files or two folders; folders are paired by relative path. Outputs: ``summary.json``, +``discrepancies.tsv.gz`` (every discrepancy, never truncated) and ``dashboard.html``. +""" + +from __future__ import annotations + +import gzip +import html +import io +import json +import multiprocessing as mp +import os +from collections.abc import Iterator, Sequence +from concurrent.futures import ProcessPoolExecutor +from dataclasses import dataclass, field +from pathlib import Path +from typing import IO, Any, cast + +import numpy as np + +from simurg import PRODUCT_NAME, __version__ +from simurg.compare.core import CATEGORY_NAMES, GTC_TOLERANCES, FieldResult, Tolerance, locus_classes, merge_results +from simurg.compare.records import Comparison, VcfRecord, compare_gtc, compare_vcf_records +from simurg.core.atomic import atomic_open, atomic_write_text, check_outputs_absent +from simurg.core.errors import ExitCode, SimurgError +from simurg.core.log import Reporter + +GTC_EXT = (".gtc",) +VCF_EXT = (".vcf", ".vcf.gz") + + +@dataclass +class CompareOptions: + a: Path + b: Path + out: Path + bpm: Path | None = None + tolerances: dict[str, str] = field(default_factory=dict) # field -> "abs=1e-5" | "ulp=2" | "exact" + samples: list[str] | None = None + threads: int = 0 + fail_on_discrepancy: bool = False + title: str = "" + force: bool = False + + +def parse_tolerance(text: str) -> Tolerance: + text = text.strip().lower() + if text in ("exact", "0"): + return Tolerance() + kv = dict(part.split("=", 1) for part in text.split(",")) + try: + return Tolerance(abs=float(kv["abs"]) if "abs" in kv else None, ulp=int(kv["ulp"]) if "ulp" in kv else None) + except (KeyError, ValueError): + raise SimurgError("E_CONFIG_KEY", f"bad tolerance {text!r}; use exact, abs=X, ulp=N or both") from None + + +def _kind(p: Path) -> str: + n = p.name.lower() + if n.endswith(GTC_EXT): + return "gtc" + if n.endswith(VCF_EXT): + return "vcf" + raise SimurgError("E_FILE_UNREADABLE", f"{p}: compare accepts .gtc and .vcf(.gz) files (convert BCF with bcftools)") + + +def pair_inputs(a: Path, b: Path) -> tuple[list[tuple[str, Path, Path]], list[str]]: + """((label, file_a, file_b)…, unmatched notes). Folders are paired by file name (or relative path).""" + for p in (a, b): + if not p.exists(): + raise SimurgError("E_FILE_NOT_FOUND", f"not found: {p}", path=str(p)) + if a.is_file() and b.is_file(): + if _kind(a) != _kind(b): + raise SimurgError("E_FILE_UNREADABLE", f"cannot compare a {_kind(a)} with a {_kind(b)}") + return [(a.name, a, b)], [] + if not (a.is_dir() and b.is_dir()): + raise SimurgError("E_FILE_UNREADABLE", "give two files or two folders") + + def listing(root: Path) -> dict[str, Path]: + """Key files by name when names are unique in the tree (so nested and flat layouts pair up), + otherwise by relative path.""" + files = [p for p in sorted(root.rglob("*")) if p.is_file() and p.name.lower().endswith(GTC_EXT + VCF_EXT)] + names = [p.name for p in files] + if len(set(names)) == len(names): + return {p.name: p for p in files} + return {str(p.relative_to(root)): p for p in files} + + la, lb = listing(a), listing(b) + pairs = [(k, la[k], lb[k]) for k in sorted(la) if k in lb] + notes = [f"only in A: {k}" for k in sorted(set(la) - set(lb))] + [ + f"only in B: {k}" for k in sorted(set(lb) - set(la)) + ] + if not pairs: + raise SimurgError("E_NO_INPUT", f"no files with the same relative path under {a} and {b}") + return pairs, notes + + +# ------------------------------------------------------------------------------------------------ workers +_W: dict[str, Any] = {} + + +def _init(bpm: str | None, tolerances: dict[str, Tolerance], samples: list[str] | None) -> None: + _W.update(bpm=bpm, tol=tolerances, samples=samples) + + +def _manifest() -> tuple[Any, Any]: + if "classes" not in _W: + if _W["bpm"]: + from simurg.formats.bpm import read_bpm + + m = read_bpm(_W["bpm"]) + _W["classes"] = locus_classes(m.chroms, m.snps, m.assay_types, m.intensity_only) + _W["names"] = list(m.names) + else: + _W["classes"], _W["names"] = None, None + return _W["classes"], _W["names"] + + +VcfKey = tuple[str, int, str, str, str] + + +def read_vcf(path: Path, sample: str | None = None) -> tuple[list[str], dict[str, dict[VcfKey, VcfRecord]], int]: + """Records per sample (all samples, or one), keyed by (CHROM, POS, REF, ALT, ID). Repeated keys are kept under + a numbered ID so nothing is silently overwritten; their count is returned.""" + opener = gzip.open if str(path).endswith(".gz") else open + samples: list[str] = [] + out: dict[str, dict[VcfKey, VcfRecord]] = {} + seen: set[VcfKey] = set() + dups = 0 + with opener(path, "rt") as fh: + for line in fh: + if line.startswith("##"): + continue + f = line.rstrip("\n").split("\t") + if line.startswith("#"): + samples = f[9:] + keep = [i for i, s in enumerate(samples) if sample is None or s == sample] + out = {samples[i]: {} for i in keep} + continue + keys = f[8].split(":") + k = (f[0], int(f[1]), f[3], f[4], f[2]) + if k in seen: + dups += 1 + k = (f[0], int(f[1]), f[3], f[4], f"{f[2]}#{dups}") + seen.add(k) + for i in keep: + vals = f[9 + i].split(":") + out[samples[i]][k] = VcfRecord(f[0], int(f[1]), f[2], f[3], f[4], dict(zip(keys, vals, strict=False))) + return samples, out, dups + + +def _compare_pair(task: tuple[str, str, str]) -> list[Comparison]: + label, pa, pb = task + kind = _kind(Path(pa)) + try: + if kind == "gtc": + from simurg.formats.gtc import read_gtc + + ga, gb = read_gtc(pa, check_write_complete=False), read_gtc(pb, check_write_complete=False) + classes, names = _manifest() + if classes is None or len(names) != ga.num_snps: + classes = locus_classes(["1"] * ga.num_snps, ["[A/B]"] * ga.num_snps, np.zeros(ga.num_snps), + np.zeros(ga.num_snps, np.bool_)) # fmt: skip + classes = {"all": classes["all"]} + names = [str(i) for i in range(ga.num_snps)] + tol = {**GTC_TOLERANCES, **_W["tol"]} + return [compare_gtc(ga, gb, label, classes, names, tol)] + sa, ra, da = read_vcf(Path(pa)) + sb, rb, db = read_vcf(Path(pb)) + wanted = _W["samples"] or [s for s in sa if s in sb] + out = [] + for s in wanted: + if s not in ra or s not in rb: + c = Comparison(f"{label}:{s}", "vcf") + c.notes.append(f"sample {s} missing in {'A' if s not in ra else 'B'}") + out.append(c) + continue + ftol = _W["tol"].get("float") + if ftol is None: + out.append(compare_vcf_records(ra[s], rb[s], f"{label}:{s}")) + else: + out.append(compare_vcf_records(ra[s], rb[s], f"{label}:{s}", ftol)) + if not wanted: + c = Comparison(label, "vcf") + c.notes.append("no sample in common") + out.append(c) + if da or db: + out[0].notes.append(f"duplicate record keys: {da} in A, {db} in B") + return out + except SimurgError as e: + c = Comparison(label, kind) + c.notes.append(f"{e.code}: {e.message}") + return [c] + + +# ------------------------------------------------------------------------------------------------ outputs +def _totals(fields: dict[str, FieldResult]) -> dict[str, Any]: + out = {} + for f, r in fields.items(): + out[f] = { + "tolerance": r.tolerance, + "counts": {cls: dict(zip(CATEGORY_NAMES, c, strict=True)) for cls, c in sorted(r.counts.items())}, + "discrepancies": sum(c[2] for k, c in r.counts.items() if k in ("all", "sample", "bins", "controls")), + "listed": len(r.discrepancies), + } + return out + + +def write_discrepancies(path: Path, rows: Iterator[dict[str, Any]]) -> int: + n = 0 + with atomic_open(path, "wb") as raw: + gz = gzip.GzipFile(fileobj=cast(IO[bytes], raw), mode="wb", mtime=0) + fh = io.TextIOWrapper(gz, encoding="utf-8", newline="\n") + fh.write("#simurg_schema=simurg.compare_discrepancies/1\npair\tfield\tlocus\ta\tb\n") + for r in rows: + fh.write( + "\t".join(str(r.get(k, "")).replace("\t", " ") for k in ("pair", "field", "locus", "a", "b")) + "\n" + ) + n += 1 + fh.flush() + fh.detach() + gz.close() + return n + + +DASH_CSS = """ +:root{--bg:#fafafa;--fg:#222;--muted:#666;--line:#ddd;--ok:#d6e6f5;--near:#e8f2d6;--bad:#f9d3ac;--unk:#e6e6e6} +@media (prefers-color-scheme: dark){:root{--bg:#1d1f21;--fg:#e6e6e6;--muted:#aaa;--line:#444; + --ok:#23405c;--near:#35451f;--bad:#6b3d10;--unk:#3a3a3a}} +body{margin:0;padding:16px;background:var(--bg);color:var(--fg);font:14px/1.4 system-ui,sans-serif} +table{border-collapse:collapse;margin:6px 0} +td,th{border:1px solid var(--line);padding:3px 8px;text-align:right} +td:first-child,th:first-child{text-align:left} +.cEXACT{background:var(--ok)} .cNEAR{background:var(--near)} .cUNKNOWN{background:var(--unk)} +tr.bad .cDISCREPANCY{background:var(--bad)} +.verdict{font-size:20px;font-weight:600} .muted{color:var(--muted)} .wrap{overflow-x:auto} +""" + + +def _table(head: Sequence[str], rows: Sequence[Sequence[str]], cls: str = "") -> str: + th = "".join(f"{h}" for h in head) + body = "".join("" + "".join(f"{c}" for c in r) + "" for r in rows) + return f'
{th}{body}
' + + +def _counts_cells(c: dict[str, int] | list[int]) -> str: + vals = c if isinstance(c, list) else [c.get(k, 0) for k in CATEGORY_NAMES] + return "".join(f'{v:,}' for k, v in zip(CATEGORY_NAMES, vals, strict=True)) + + +def render_dashboard(summary: dict[str, Any], examples: dict[str, list[dict[str, Any]]]) -> str: + esc = html.escape + title = esc(summary["title"] or "Differential test") + verdict = summary["verdict"] + head = "".join(f"{k}" for k in CATEGORY_NAMES) + rows = [] + for f, d in summary["fields"].items(): + counts = d["counts"] + main: dict[str, int] = next((counts[k] for k in ("all", "sample", "bins", "controls") if k in counts), {}) + per = "".join(f"{esc(cls)}{_counts_cells(c)}" for cls, c in counts.items()) + inner = f"{head}{per}
class
" + bad = ' class="bad"' if d["discrepancies"] else "" + rows.append( + f"
{esc(f)}{inner}
" + f"{esc(d['tolerance'])}{_counts_cells(main)}" + ) + fields_html = ( + f'
{head}' + + "".join(rows) + + "
field (click for classes)tolerance
" + ) + parts = [ + f'

{title}

{PRODUCT_NAME} {__version__} · {summary["pairs"]} pair(s)
', + f'
A = {esc(summary["a"])}
B = {esc(summary["b"])}
', + f'

{verdict}: {summary["discrepancies"]:,} discrepancies, {summary["unknown"]:,} unknown

', + "

Fields

" + fields_html, + ] + meta = [ + [esc(p), esc(k), esc(v[0]), esc(v[1])] + for p, md in summary["metadata_differences"].items() + for k, v in md.items() + ] + if meta: + parts.append("

Metadata differences (not counted)

" + _table(["pair", "field", "A", "B"], meta)) + if summary["notes"]: + parts.append("

Notes

    " + "".join(f"
  • {esc(n)}
    • " for n in summary["notes"]) + "
") + for f, lst in examples.items(): + if lst: + n = summary["fields"][f]["listed"] + body = [[esc(str(r.get(k, "")))[:200] for k in ("pair", "locus", "a", "b")] for r in lst] + parts.append( + f"

{esc(f)}: first {len(lst)} of {n:,} discrepancies

" + + _table(["pair", "locus", "A", "B"], body) + ) + parts.append('

Every discrepancy and every unmatched record is listed in discrepancies.tsv.gz.

') + return ( + '' + '' + f"{title} — {PRODUCT_NAME}" + + "\n".join(parts) + + "\n" + ) + + +def run_compare(opts: CompareOptions, rep: Reporter) -> int: + targets = [opts.out / "summary.json", opts.out / "discrepancies.tsv.gz", opts.out / "dashboard.html"] + check_outputs_absent(targets, force=opts.force) + tol = {k: parse_tolerance(v) for k, v in opts.tolerances.items()} + unknown = [k for k in tol if k not in GTC_TOLERANCES and k != "float"] + if unknown: + valid = ", ".join(sorted(GTC_TOLERANCES)) + raise SimurgError("E_CONFIG_KEY", f"unknown tolerance field(s) {unknown}; GTC fields: {valid}; VCF: float") + pairs, notes = pair_inputs(opts.a, opts.b) + tasks = [(label, str(pa), str(pb)) for label, pa, pb in pairs] + workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(tasks))) + rep.info(f"comparing {len(tasks)} pair(s) with {workers} worker(s)") + if workers == 1: + _init(str(opts.bpm) if opts.bpm else None, tol, opts.samples) + results = [c for t in tasks for c in _compare_pair(t)] + else: + with ProcessPoolExecutor(workers, mp_context=mp.get_context("spawn"), initializer=_init, + initargs=(str(opts.bpm) if opts.bpm else None, tol, opts.samples)) as ex: # fmt: skip + results = [c for batch in ex.map(_compare_pair, tasks) for c in batch] + merged = merge_results([f for c in results for f in c.fields]) + for c in results: + notes += [f"{c.label}: {n}" for n in c.notes] + + def main_counts(r: FieldResult) -> list[int]: + for k in ("all", "sample", "bins", "controls"): + if k in r.counts: + return r.counts[k] + return [0, 0, 0, 0] + + total_disc = sum(main_counts(r)[2] for r in merged.values()) + total_unknown = sum(main_counts(r)[3] for r in merged.values()) + summary = { + "schema": "simurg.compare/1", + "title": opts.title, + "a": str(opts.a), + "b": str(opts.b), + "pairs": len(results), + "discrepancies": total_disc, + "unknown": total_unknown, + "verdict": "DISCREPANCIES" if total_disc else ("UNKNOWNS" if total_unknown else "PASS"), + "fields": _totals(merged), + "metadata_differences": {c.label: c.metadata for c in results if c.metadata}, + "notes": notes, + } + opts.out.mkdir(parents=True, exist_ok=True) + write_discrepancies(opts.out / "discrepancies.tsv.gz", (d for r in merged.values() for d in r.discrepancies)) + atomic_write_text(opts.out / "summary.json", json.dumps(summary, indent=2) + "\n") + examples = {f: r.discrepancies[:50] for f, r in merged.items()} + atomic_write_text(opts.out / "dashboard.html", render_dashboard(summary, examples)) + rep.info(f"{summary['verdict']}: {total_disc} discrepancies, {total_unknown} unknown in {len(results)} pair(s); " + f"report {opts.out / 'dashboard.html'}") # fmt: skip + if total_disc and opts.fail_on_discrepancy: + return ExitCode.PARTIAL + return ExitCode.OK diff --git a/tests/unit/test_compare.py b/tests/unit/test_compare.py new file mode 100644 index 0000000..09b9780 --- /dev/null +++ b/tests/unit/test_compare.py @@ -0,0 +1,112 @@ +"""Phase P7: differential testing (plan §16.3, CLI-038) — categories, tolerances, classes, GTC/VCF pairs, CLI.""" + +from __future__ import annotations + +import gzip +import json +from pathlib import Path + +import numpy as np +from tests.unit.test_gtc import make_gtc +from typer.testing import CliRunner + +from simurg.cli.main import app +from simurg.compare.core import ( + DISCREPANCY, + EXACT, + NEAR, + Tolerance, + categorize, + compare_locus_field, + locus_classes, +) +from simurg.compare.records import VcfRecord, compare_gtc, compare_vcf_records +from simurg.formats.gtc import write_gtc + + +def test_categories_and_tolerances() -> None: + a = np.array([1.0, 1.0, 1.0, np.nan, np.nan, 0.0], np.float32) + b = np.array([1.0, 1.000001, 1.1, np.nan, 1.0, -0.0], np.float32) + assert categorize(a, b, Tolerance(abs=1e-5)).tolist() == [EXACT, NEAR, DISCREPANCY, EXACT, DISCREPANCY, EXACT] + assert categorize(a, b, Tolerance()).tolist()[:3] == [EXACT, DISCREPANCY, DISCREPANCY] + one_ulp = np.nextafter(np.float32(1.0), np.float32(2.0)) + assert categorize(np.array([1.0], np.float32), np.array([one_ulp], np.float32), Tolerance(ulp=2))[0] == NEAR + assert categorize(np.array([3, 4]), np.array([3, 5]), Tolerance()).tolist() == [EXACT, DISCREPANCY] + neg_nan = np.frombuffer(np.uint32(0xFFC00000).tobytes(), np.float32) + assert categorize(neg_nan, np.array([np.nan], np.float32), Tolerance())[0] == EXACT # K-014 + + +def test_locus_classes_and_unknown() -> None: + cls = locus_classes(["1", "X", "XY", "0", "MT"], ["[A/G]", "[I/D]", "[A/C]", "[A/G]", "[C/T]"], + [0, 1, 0, 2, 0], [False, False, False, True, False]) # fmt: skip + assert cls["type=indel"].tolist() == [False, True, False, False, False] + assert cls["chrom=PAR"].sum() == 1 and cls["chrom=unmapped"].sum() == 1 and cls["intensity_only"].sum() == 1 + r = compare_locus_field("baf", np.zeros(5), None, Tolerance(abs=1e-5), cls, list("abcde"), {}) + assert r.counts["all"] == [0, 0, 0, 5] # one side lacks the field → UNKNOWN + + +def test_gtc_pair_detects_every_difference() -> None: + a = make_gtc(40, 5, seed=1) + b = make_gtc(40, 5, seed=1) + b.genotypes = b.genotypes.copy() + b.genotypes[7] = (b.genotypes[7] + 1) % 4 + b.logr_ratios = b.logr_ratios.copy() + b.logr_ratios[3] += np.float32(2e-6) # within 1e-5 → NEAR + b.b_allele_freqs = b.b_allele_freqs.copy() + b.b_allele_freqs[5] += np.float32(0.01) # outside → DISCREPANCY + b.sample_name = "renamed" + cls = {"all": np.ones(40, np.bool_)} + names = [f"rs{i}" for i in range(40)] + res = compare_gtc(a, b, "pair", cls, names) + by = {f.field: f for f in res.fields} + assert by["genotypes"].counts["all"] == [39, 0, 1, 0] and by["genotypes"].discrepancies[0]["locus"] == "rs7" + assert by["logr_ratios"].counts["all"][NEAR] == 1 + assert by["b_allele_freqs"].counts["all"][DISCREPANCY] == 1 + assert by["raw_x"].counts["all"] == [40, 0, 0, 0] + assert res.metadata["sample_name"][1] == "renamed" # reported, not counted + + +def test_vcf_records() -> None: + def rec(gt: str, baf: str, vid: str = "rs1") -> VcfRecord: + return VcfRecord("1", 10, vid, "A", "G", {"GT": gt, "BAF": baf}) + + a = {("1", 10, "A", "G", "rs1"): rec("0/1", "0.5"), ("1", 20, "C", "T", "rs2"): rec("0/0", "0")} + b = {("1", 10, "A", "G", "rs1"): rec("1/0", "0.500001"), ("X", 5, "C", "T", "rs3"): rec("1", ".")} + res = compare_vcf_records(a, b, "p") + by = {f.field: f for f in res.fields} + assert by["GT"].counts["all"] == [1, 0, 0, 0] # 0/1 == 1/0 + assert by["BAF"].counts["all"] == [0, 1, 0, 0] + assert by["record"].counts["all"] == [1, 0, 0, 2] # one record only in A, one only in B + merged = {("1", 30, "A", "G,T", "a,b"): VcfRecord("1", 30, "a,b", "A", "G,T", {"GT": "0/1"})} + res = compare_vcf_records(merged, merged, "p") + assert {f.field for f in res.fields} == {"GT", "record"} # BAF absent on both sides: not applicable + + +def test_cli_end_to_end(tmp_path: Path) -> None: + a_dir, b_dir = tmp_path / "a" / "nested", tmp_path / "b" + a_dir.mkdir(parents=True) + b_dir.mkdir() + g = make_gtc(30, 5, seed=3) + write_gtc(g, a_dir / "s1.gtc") + g2 = make_gtc(30, 5, seed=3) + g2.genotype_scores = g2.genotype_scores.copy() + g2.genotype_scores[0] += np.float32(0.5) + write_gtc(g2, b_dir / "s1.gtc") + write_gtc(make_gtc(30, 5, seed=4), a_dir / "only_a.gtc") + out = tmp_path / "out" + res = CliRunner().invoke(app, ["compare", str(tmp_path / "a"), str(b_dir), "--out", str(out), "--threads", "1", + "--fail-on-discrepancy", "--quiet"]) # fmt: skip + assert res.exit_code == 3, res.output + s = json.loads((out / "summary.json").read_text()) + assert s["verdict"] == "DISCREPANCIES" and s["discrepancies"] == 1 + assert any("only in A: only_a.gtc" in n for n in s["notes"]) + rows = gzip.decompress((out / "discrepancies.tsv.gz").read_bytes()).decode().splitlines() + assert rows[2].split("\t")[1] == "genotype_scores" + html = (out / "dashboard.html").read_text() + assert "DISCREPANCIES" in html and "https://" not in html + res = CliRunner().invoke(app, ["compare", str(b_dir / "s1.gtc"), str(b_dir / "s1.gtc"), "--out", str(tmp_path / "o2"), + "--tolerance", "genotype_scores=abs=1", "--quiet"]) # fmt: skip + assert res.exit_code == 0 and json.loads((tmp_path / "o2" / "summary.json").read_text())["verdict"] == "PASS" + res = CliRunner().invoke(app, ["compare", str(b_dir), str(b_dir), "--out", str(tmp_path / "o3"), + "--tolerance", "nonsense=abs=1"]) # fmt: skip + assert res.exit_code == 2 and "E_CONFIG_KEY" in res.output