P7: differential-testing harness (simurg compare) and nightly run on the Pi
- 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
This commit is contained in:
kemal committed
1 parent
51747674ba
commit
3c33ae5cb1
10 files changed
+1268
-1
No files matched your search
@@ -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 <files…> [--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
|
||||
|
||||
|
||||
@@ -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).
|
||||
Executable
+19
@@ -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"
|
||||
@@ -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/<YYYY-MM-DD>/{index.html, summary.json, <step>/...}, 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'<a href="{esc(r["report"])}">dashboard</a>' if r.get("report") else ""
|
||||
verdict = esc(str(r.get("verdict")))
|
||||
rows.append(
|
||||
f"<tr class='{verdict}'><td>{esc(k)}</td><td>{verdict}</td><td>{r.get('seconds', '')}</td>"
|
||||
f"<td>{esc(details)}</td><td>{link}</td></tr>"
|
||||
)
|
||||
return (
|
||||
"<!doctype html><meta charset=utf-8><meta name=viewport content='width=device-width,initial-scale=1'>"
|
||||
f"<title>Nightly {night['date']}</title><style>{CSS}</style>"
|
||||
f"<h1>Nightly {night['date']} — {esc(night['overall'])}</h1>"
|
||||
f"<p>{esc(night['version'])} · {night['seconds']} s</p>"
|
||||
"<table><tr><th>step</th><th>verdict</th><th>seconds</th><th>details</th><th></th></tr>"
|
||||
+ "".join(rows)
|
||||
+ "</table>"
|
||||
)
|
||||
|
||||
|
||||
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"<td>{html.escape(str(s['steps'].get(k, {}).get('verdict', '')))}</td>" for k in steps)
|
||||
hist.append(f"<tr><td><a href='{s['date']}/index.html'>{s['date']}</a></td><td>{s['overall']}</td>{cells}</tr>")
|
||||
head = "".join(f"<th>{html.escape(k)}</th>" for k in steps)
|
||||
(a.root / "index.html").write_text(
|
||||
f"<!doctype html><meta charset=utf-8><title>SimurgArray nightly</title><style>{CSS}</style>"
|
||||
"<h1>SimurgArray nightly history</h1>"
|
||||
f"<table><tr><th>night</th><th>overall</th>{head}</tr>{''.join(hist)}</table>"
|
||||
)
|
||||
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())
|
||||
@@ -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")
|
||||
|
||||
Whitespace-only changes.
@@ -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
|
||||
@@ -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
|
||||
@@ -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"<th>{h}</th>" for h in head)
|
||||
body = "".join("<tr>" + "".join(f"<td>{c}</td>" for c in r) + "</tr>" for r in rows)
|
||||
return f'<div class="wrap"><table class="{cls}"><tr>{th}</tr>{body}</table></div>'
|
||||
|
||||
|
||||
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'<td class="c{k}">{v:,}</td>' 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"<th>{k}</th>" 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"<tr><td>{esc(cls)}</td>{_counts_cells(c)}</tr>" for cls, c in counts.items())
|
||||
inner = f"<table><tr><th>class</th>{head}</tr>{per}</table>"
|
||||
bad = ' class="bad"' if d["discrepancies"] else ""
|
||||
rows.append(
|
||||
f"<tr{bad}><td><details><summary>{esc(f)}</summary>{inner}</details></td>"
|
||||
f"<td>{esc(d['tolerance'])}</td>{_counts_cells(main)}</tr>"
|
||||
)
|
||||
fields_html = (
|
||||
f'<div class="wrap"><table><tr><th>field (click for classes)</th><th>tolerance</th>{head}</tr>'
|
||||
+ "".join(rows)
|
||||
+ "</table></div>"
|
||||
)
|
||||
parts = [
|
||||
f'<h1>{title}</h1><div class="muted">{PRODUCT_NAME} {__version__} · {summary["pairs"]} pair(s)</div>',
|
||||
f'<div class="muted">A = {esc(summary["a"])}<br>B = {esc(summary["b"])}</div>',
|
||||
f'<p class="verdict">{verdict}: {summary["discrepancies"]:,} discrepancies, {summary["unknown"]:,} unknown</p>',
|
||||
"<h2>Fields</h2>" + 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("<h2>Metadata differences (not counted)</h2>" + _table(["pair", "field", "A", "B"], meta))
|
||||
if summary["notes"]:
|
||||
parts.append("<h2>Notes</h2><ul>" + "".join(f"<li>{esc(n)}</li>" for n in summary["notes"]) + "</ul>")
|
||||
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"<h3>{esc(f)}: first {len(lst)} of {n:,} discrepancies</h3>"
|
||||
+ _table(["pair", "locus", "A", "B"], body)
|
||||
)
|
||||
parts.append('<p class="muted">Every discrepancy and every unmatched record is listed in discrepancies.tsv.gz.</p>')
|
||||
return (
|
||||
'<!doctype html><html lang="en"><head><meta charset="utf-8">'
|
||||
'<meta name="viewport" content="width=device-width, initial-scale=1">'
|
||||
f"<title>{title} — {PRODUCT_NAME}</title><style>{DASH_CSS}</style></head><body>"
|
||||
+ "\n".join(parts)
|
||||
+ "</body></html>\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
|
||||
@@ -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
|
||||
Reference in new issue
Block a user