P4: unit tests for calling/BAF/LRR/sex/CLI; vendor validation script
- tests/unit/test_genotype.py: fuzzy memberships incl. NaN rules, calls at cluster centres, intensity threshold, BAF/LRR interpolation and NaN pattern, LogR-dev locus set, sex rules v1/v2, TOP-strand base calls, invariant dates, end-to-end `simurg genotype` (summary schema, pass-through columns, thread/layout determinism, overwrite refusal) - tests/oracle/test_p4_genotype.py: 36 GSA samples x 2 versions byte-identical to the oracle - scripts/validate_p4_large.py + cloud task p4: compare every vendor GTC field by field and as whole files; reports which sex rule matches
This commit is contained in:
kemal committed
1 parent
1f7b6d10c5
commit
e6cf273ab4
3 files changed
+496
No files matched your search
@@ -0,0 +1,235 @@
|
||||
"""Phase P4: calling, BAF/LRR, sex, base calls and `simurg genotype` on synthetic data (plan §7.5–§7.8)."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import json
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import pytest
|
||||
from typer.testing import CliRunner
|
||||
|
||||
from simurg.formats.bpm import read_bpm
|
||||
from simurg.formats.egt import read_egt
|
||||
from simurg.formats.gtc import read_gtc
|
||||
from simurg.genotype.baf_lrr import NEG_NAN, autosomal_mask, baf_lrr, logr_dev
|
||||
from simurg.genotype.clusters import ClusterArrays, align_clusters
|
||||
from simurg.genotype.gencall import call_genotypes, call_rate, smf, zmf
|
||||
from simurg.genotype.sample import base_calls, invariant_datetime, top_alleles
|
||||
from simurg.genotype.sex import SexParams, chromosome_classes, estimate_sex
|
||||
from simurg.synth.files import SynthIdat, SynthLocus, encode_bpm, encode_egt, encode_idat
|
||||
|
||||
F32 = np.float32
|
||||
|
||||
|
||||
def clusters(n: int, *, gentrain: float = 0.8, threshold: float = np.nan) -> ClusterArrays:
|
||||
tm = np.tile(np.array([0.05, 0.5, 0.95], np.float32), (n, 1))
|
||||
rm = np.ones((n, 3), np.float32)
|
||||
return ClusterArrays(
|
||||
present=np.ones(n, bool),
|
||||
n=np.full((n, 3), 100, np.int32),
|
||||
theta_mean=tm,
|
||||
theta_dev=np.full((n, 3), 0.03, np.float32),
|
||||
r_mean=rm,
|
||||
r_dev=np.full((n, 3), 0.1, np.float32),
|
||||
intensity_threshold=np.full(n, threshold, np.float32),
|
||||
gentrain=np.full(n, gentrain, np.float32),
|
||||
mean_r=np.ones(n, np.float32),
|
||||
)
|
||||
|
||||
|
||||
def test_fuzzy_memberships() -> None:
|
||||
x = np.array([-1, 0, 0.25, 0.5, 0.75, 1, 2, np.nan], np.float32)
|
||||
s = smf(x, F32(0), F32(1))
|
||||
assert s.tolist()[:7] == [0.0, 0.0, 0.125, 0.5, 0.875, 1.0, 1.0]
|
||||
assert s[7] == 1.0 # the reference returns 1 for NaN (all comparisons false)
|
||||
z = zmf(x, F32(0), F32(1))
|
||||
assert z.tolist()[:7] == [1.0, 1.0, 0.875, 0.5, 0.125, 0.0, 0.0] and z[7] == 0.0
|
||||
assert smf(np.array([0.4, 0.6], np.float32), F32(1), F32(0)).tolist() == [0.0, 1.0] # a >= b rule
|
||||
|
||||
|
||||
def test_calls_at_cluster_centres() -> None:
|
||||
n = 6
|
||||
t = np.array([0.05, 0.5, 0.95, 0.5, 0.5, 0.5], np.float32)
|
||||
r = np.ones(n, np.float32)
|
||||
io = np.array([0, 0, 0, 1, 0, 0], bool)
|
||||
raw_zero = np.array([0, 0, 0, 0, 1, 0], bool)
|
||||
c = clusters(n)
|
||||
c.intensity_threshold[5] = 2.0 # R below the cluster file's intensity threshold -> no-call
|
||||
calls = call_genotypes(t, r, raw_zero, c, io)
|
||||
assert calls.genotypes.tolist() == [1, 2, 3, 0, 0, 0]
|
||||
assert (calls.scores[:3] > 0.5).all() and calls.scores[4] == 0 and calls.scores[5] == 0
|
||||
assert calls.num_calls == 3 and calls.num_intensity_only == 1
|
||||
assert call_rate(calls, n) == F32(3 / 5)
|
||||
|
||||
|
||||
def test_low_gentrain_and_outliers_score_low() -> None:
|
||||
n = 2
|
||||
t = np.array([0.5, 0.3], np.float32) # second point between clusters
|
||||
calls_good = call_genotypes(t, np.ones(n, np.float32), np.zeros(n, bool), clusters(n), np.zeros(n, bool))
|
||||
calls_bad = call_genotypes(
|
||||
t, np.ones(n, np.float32), np.zeros(n, bool), clusters(n, gentrain=0.05), np.zeros(n, bool)
|
||||
)
|
||||
assert calls_good.scores[0] > calls_good.scores[1]
|
||||
assert (calls_bad.scores < 0.15).all() and (calls_bad.genotypes == 0).all()
|
||||
|
||||
|
||||
def test_baf_lrr_interpolation() -> None:
|
||||
n = 6
|
||||
t = np.array([0.05, 0.5, 0.95, 0.275, 1.2, np.nan], np.float32)
|
||||
r = np.array([1, 1, 2, 1, 1, np.nan], np.float32)
|
||||
c = clusters(n)
|
||||
io = np.array([0, 0, 0, 0, 0, 0], bool)
|
||||
baf, lrr = baf_lrr(t, r, np.zeros(n, bool), c, io)
|
||||
assert baf[:5].tolist() == pytest.approx([0.0, 0.5, 1.0, 0.25, 1.0]) # clamped at 1
|
||||
assert lrr[:4].tolist() == pytest.approx([0.0, 0.0, 1.0, 0.0])
|
||||
assert baf.view(np.uint32)[5] == 0xFFC00000 and lrr.view(np.uint32)[5] == 0xFFC00000
|
||||
io[1] = True
|
||||
c.mean_r[1] = 2.0
|
||||
_, lrr2 = baf_lrr(t, r, np.zeros(n, bool), c, io)
|
||||
assert lrr2[1] == pytest.approx(-1.0) # intensity-only loci use the count-weighted mean R
|
||||
assert np.isnan(NEG_NAN)
|
||||
|
||||
|
||||
def test_logr_dev_locus_set() -> None:
|
||||
lrr = np.array([1, -1, 1, -1, 50, 50, np.nan], np.float32)
|
||||
chroms = ["1", "2", "0", "22_KI270879v1_alt", "X", "MT", "3"]
|
||||
auto = autosomal_mask(chroms)
|
||||
assert auto.tolist() == [True, True, True, True, False, False, True]
|
||||
assert logr_dev(lrr, auto, np.zeros(7, bool), np.ones(7, bool)) == pytest.approx(1.0)
|
||||
|
||||
|
||||
def test_sex_rules() -> None:
|
||||
chroms = ["X"] * 30 + ["Y"] * 25 + ["1"] * 150
|
||||
classes = chromosome_classes(chroms)
|
||||
n = len(chroms)
|
||||
gt = np.ones(n, np.uint8)
|
||||
present = np.ones(n, bool)
|
||||
gtrain = np.full(n, 0.8, np.float32)
|
||||
r = np.ones(n, np.float32)
|
||||
assert estimate_sex(gt, r, gtrain, present, classes) == "M" # median Y R 1.0 > 0.3
|
||||
r_f = r.copy()
|
||||
r_f[30:55] = 0.05
|
||||
assert estimate_sex(gt, r_f, gtrain, present, classes) == "F" # X R 1.0 > 0.9
|
||||
r_u = r_f.copy()
|
||||
r_u[:30] = 0.5
|
||||
assert estimate_sex(gt, r_u, gtrain, present, classes) == "U"
|
||||
het = gt.copy()
|
||||
het[:10] = 2 # X het rate 10/30 > 0.1
|
||||
assert estimate_sex(het, r, gtrain, present, classes, SexParams(version=1)) == "F"
|
||||
assert estimate_sex(gt, r, gtrain, present, classes, SexParams(version=1)) == "M"
|
||||
|
||||
|
||||
def test_base_calls_and_top_strand(tmp_path: Path) -> None:
|
||||
loci = [
|
||||
SynthLocus("a", snp="[A/G]", ilmn_strand="TOP"),
|
||||
SynthLocus("b", address_a=2, snp="[T/C]", ilmn_strand="BOT"),
|
||||
SynthLocus("c", address_a=3, snp="[N/A]", ilmn_strand="BOT"),
|
||||
]
|
||||
(tmp_path / "m.bpm").write_bytes(encode_bpm(loci))
|
||||
a, b = top_alleles(read_bpm(tmp_path / "m.bpm"))
|
||||
assert bytes(a) == b"AA\x00" and bytes(b) == b"GGT"
|
||||
bc = base_calls(np.array([1, 2, 0], np.uint8), a, b)
|
||||
assert [bytes(x) for x in bc] == [b"AA", b"AG", b"--"]
|
||||
|
||||
|
||||
def test_invariant_datetime() -> None:
|
||||
import datetime as dt
|
||||
|
||||
assert invariant_datetime(dt.datetime(2026, 1, 2, 0, 5, 9)) == "1/2/2026 12:05:09 AM"
|
||||
assert invariant_datetime(dt.datetime(2026, 12, 31, 13, 0, 0)) == "12/31/2026 1:00:00 PM"
|
||||
|
||||
|
||||
def _synthetic_run(tmp_path: Path) -> dict[str, Path]:
|
||||
"""Two samples on a 600-locus synthetic array with a matching cluster file."""
|
||||
rng = np.random.default_rng(4)
|
||||
n = 600
|
||||
chroms = ["1"] * 500 + ["X"] * 60 + ["Y"] * 40
|
||||
loci = [SynthLocus(f"rs{i}", chrom=chroms[i], address_a=i + 1, pos=i + 1) for i in range(n)]
|
||||
loci.append(SynthLocus("TGA_Ctrl_5716", address_a=n + 1))
|
||||
(tmp_path / "m.bpm").write_bytes(encode_bpm(loci))
|
||||
names = [lc.name for lc in loci]
|
||||
tm = np.tile([0.03, 0.5, 0.97], (len(loci), 1))
|
||||
rm = np.tile([1.0, 1.0, 1.0], (len(loci), 1))
|
||||
(tmp_path / "c.egt").write_bytes(encode_egt(names, theta_mean=tm, r_mean=rm))
|
||||
idats = tmp_path / "idats"
|
||||
idats.mkdir()
|
||||
addr = np.arange(1, len(loci) + 1, dtype=np.int32)
|
||||
for pos in ("R01C01", "R02C01"):
|
||||
g = rng.choice(3, len(loci))
|
||||
inten = rng.uniform(0.6, 1.4, len(loci))
|
||||
x = inten * np.select([g == 0, g == 1, g == 2], [1.0, 0.5, 0.02]) * 8000 + 300
|
||||
y = inten * np.select([g == 0, g == 1, g == 2], [0.02, 0.5, 1.0]) * 9000 + 400
|
||||
for ch, mean in (("Red", x), ("Grn", y)):
|
||||
s = SynthIdat(
|
||||
address_ids=addr,
|
||||
mean=np.clip(mean, 0, 65535).astype(np.uint16),
|
||||
n_beads=np.full(len(loci), 12, np.uint8),
|
||||
barcode="200000000001",
|
||||
position=pos,
|
||||
)
|
||||
(idats / f"200000000001_{pos}_{ch}.idat").write_bytes(encode_idat(s))
|
||||
sheet = tmp_path / "s.csv"
|
||||
sheet.write_text(
|
||||
"[Data]\nSentrixBarcode_A,SentrixPosition_A,Sample_Name,Gender,Tissue\n"
|
||||
"200000000001,R01C01,alpha,Male,blood\n200000000001,R02C01,beta,Female,saliva\n"
|
||||
"200000000001,R09C01,ghost,Female,blood\n"
|
||||
)
|
||||
return {"bpm": tmp_path / "m.bpm", "egt": tmp_path / "c.egt", "idats": idats, "sheet": sheet}
|
||||
|
||||
|
||||
def test_genotype_cli_end_to_end(tmp_path: Path) -> None:
|
||||
from simurg.cli.main import app
|
||||
|
||||
f = _synthetic_run(tmp_path)
|
||||
runner = CliRunner()
|
||||
args = [
|
||||
"genotype",
|
||||
"--quiet",
|
||||
"--bpm",
|
||||
str(f["bpm"]),
|
||||
"--egt",
|
||||
str(f["egt"]),
|
||||
"--idats",
|
||||
str(f["idats"]),
|
||||
"--sheet",
|
||||
str(f["sheet"]),
|
||||
"--autocall-date",
|
||||
"1/1/2026 12:00:00 PM",
|
||||
]
|
||||
r1 = runner.invoke(app, [*args, "--out", str(tmp_path / "o1"), "--threads", "1"])
|
||||
assert r1.exit_code == 3, r1.output # 'ghost' has no IDATs
|
||||
r2 = runner.invoke(app, [*args, "--out", str(tmp_path / "o2"), "--threads", "2", "--layout", "flat"])
|
||||
assert r2.exit_code == 3
|
||||
for stem in ("alpha", "beta"):
|
||||
a = (tmp_path / "o1" / "samples" / stem / f"{stem}.gtc").read_bytes()
|
||||
assert a == (tmp_path / "o2" / f"{stem}.gtc").read_bytes() # thread count / layout never change bytes
|
||||
g = read_gtc(tmp_path / "o1" / "samples" / stem / f"{stem}.gtc")
|
||||
assert g.version == 5 and g.num_snps == 601 and g.autocall_version.startswith("SimurgArray ")
|
||||
assert g.cluster_file == "c.egt" and g.snp_manifest == "m.bpm"
|
||||
assert g.call_rate > 0.9
|
||||
summary = json.loads((tmp_path / "o1" / "sample_summary.json").read_text())
|
||||
assert summary["schema"] == "simurg.sample_summary/1"
|
||||
cols = summary["columns"]
|
||||
assert cols[:5] == ["sample_id", "sample_name", "sample_folder", "status", "fail_reason"]
|
||||
assert "TGA_Ctrl_5716_norm_r" in cols and cols[-2:] == ["Gender", "Tissue"]
|
||||
by = {s["sample_name"]: s for s in summary["samples"]}
|
||||
assert by["ghost"]["status"] == "MISSING_INPUT" and by["alpha"]["status"] == "OK"
|
||||
assert by["alpha"]["sex_reported"] == "M" and by["alpha"]["Tissue"] == "blood"
|
||||
assert by["alpha"]["sex_check"] in ("match", "mismatch", "unknown")
|
||||
tsv = (tmp_path / "o1" / "sample_summary.tsv").read_text().splitlines()
|
||||
assert tsv[0].split("\t") == cols and len(tsv) == 4
|
||||
r3 = runner.invoke(app, [*args, "--out", str(tmp_path / "o1")])
|
||||
assert r3.exit_code == 2 and "E_OUTPUT_EXISTS" in r3.output
|
||||
|
||||
|
||||
def test_cluster_alignment_errors(tmp_path: Path) -> None:
|
||||
f = _synthetic_run(tmp_path)
|
||||
(tmp_path / "small.egt").write_bytes(encode_egt(["rs1"], theta_mean=[[0.1, 0.5, 0.9]], r_mean=[[1, 1, 1]]))
|
||||
from simurg.core.errors import SimurgError
|
||||
|
||||
with pytest.raises(SimurgError) as e:
|
||||
align_clusters(read_egt(tmp_path / "small.egt"), read_bpm(f["bpm"]))
|
||||
assert e.value.code == "E_EGT_MANIFEST_MISMATCH"
|
||||
c = align_clusters(read_egt(tmp_path / "small.egt"), read_bpm(f["bpm"]), allow_missing=True)
|
||||
assert int(c.present.sum()) == 1
|
||||
Reference in new issue
Block a user