"""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