From e6cf273ab4bdaa1e094476cee7fbba7e8d2967d3 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Mon, 28 Sep 2026 00:10:47 +0100 Subject: [PATCH] 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 --- scripts/cloud/tasks/p4.sh | 13 ++ scripts/validate_p4_large.py | 248 +++++++++++++++++++++++++++++++++++ tests/unit/test_genotype.py | 235 +++++++++++++++++++++++++++++++++ 3 files changed, 496 insertions(+) create mode 100755 scripts/cloud/tasks/p4.sh create mode 100644 scripts/validate_p4_large.py create mode 100644 tests/unit/test_genotype.py diff --git a/scripts/cloud/tasks/p4.sh b/scripts/cloud/tasks/p4.sh new file mode 100755 index 0000000..1e67025 --- /dev/null +++ b/scripts/cloud/tasks/p4.sh @@ -0,0 +1,13 @@ +#!/bin/bash +# Phase P4: genotype every Omni2.5 IDAT pair and compare the GTC with the vendor's (fields + whole file). +# Args: +set -uo pipefail +D=$1; OUT=$2 +SRC=s3://1000genomes/release/20130502/supporting/hd_genotype_chip/broad_intensities +mkdir -p $D/egt +aws s3 cp --no-sign-request "$SRC/HumanOmni2.5-4v1-Multi_H.egt" $D/egt/ --only-show-errors +cp $D/HumanOmni2.5-4v1-Multi_B.egt $D/egt/ +/usr/bin/time -v .venv/bin/python scripts/validate_p4_large.py --root $D/data --bpm $D/HumanOmni2.5-4v1-Multi_B.bpm \ + --egt-dir $D/egt --out $OUT/summary.json --norm-version 1.1.2 --sex-version 1 --workers "$(nproc)" \ + 2> $OUT/time.txt || echo "P4 validation reported differences" +tail -25 $OUT/time.txt diff --git a/scripts/validate_p4_large.py b/scripts/validate_p4_large.py new file mode 100644 index 0000000..f2ccf63 --- /dev/null +++ b/scripts/validate_p4_large.py @@ -0,0 +1,248 @@ +#!/usr/bin/env python3 +"""Phase P4 validation against vendor GTCs: genotype every IDAT pair and compare the produced GTC with the vendor's, +field by field and as whole files (after aligning free-text metadata that cannot match). + +Usage: validate_p4_large.py --root DIR --bpm F --egt-dir DIR --out summary.json + [--norm-version 1.1.2] [--sex-version 1] [--workers N] [--limit N] +""" + +from __future__ import annotations + +import argparse +import json +import multiprocessing as mp +import os +import re +import sys +import time +import traceback +from concurrent.futures import ProcessPoolExecutor +from pathlib import Path +from typing import Any + +import numpy as np + +NEG_NAN = np.frombuffer(np.uint32(0xFFC00000).tobytes(), dtype=np.float32)[0] +ALIGN = ( + "cluster_file", + "snp_manifest", + "autocall_date", + "autocall_version", + "imaging_date", + "sample_name", + "sample_plate", + "sample_well", + "scanner_data", +) +_W: dict[str, Any] = {} + + +def _init(bpm: str, egt_dir: str, norm: str, sex_version: int) -> None: + _W.update(bpm=bpm, egt_dir=egt_dir, norm=norm, sex_version=sex_version, ctx={}) + + +def _ctx(cluster_field: str) -> Any: + from simurg.formats.bpm import read_bpm + from simurg.formats.egt import read_egt + from simurg.genotype.clusters import align_clusters + from simurg.genotype.sample import manifest_context + + name = re.sub(r"\(MD5Checksum=.*\)$", "", cluster_field).strip() + if name not in _W["ctx"]: + if "manifest" not in _W: + _W["manifest"] = read_bpm(_W["bpm"]) + path = Path(_W["egt_dir"]) / name + if not path.exists(): + raise FileNotFoundError(f"cluster file {name!r} named in the GTC is not in {_W['egt_dir']}") + m = _W["manifest"] + _W["ctx"][name] = manifest_context(m, align_clusters(read_egt(path), m), cluster_file_name=name) + return _W["ctx"][name] + + +def _canon(a: np.ndarray) -> np.ndarray: # type: ignore[type-arg] + a = np.asarray(a).copy() + if a.dtype.kind == "f": + a[np.isnan(a)] = NEG_NAN + return a + + +def check(args: tuple[str, str, str]) -> dict[str, Any]: + from simurg.formats.gtc import encode_gtc, read_gtc + from simurg.formats.idat import read_idat + from simurg.genotype.sample import GenotypeParams, genotype_sample + from simurg.genotype.sex import SexParams + + gtc, grn, red = args + res: dict[str, Any] = {"gtc": gtc} + try: + ref = read_gtc(gtc, check_write_complete=False) + ctx = _ctx(ref.cluster_file) + params = GenotypeParams( + norm_version=_W["norm"], sex=SexParams(version=_W["sex_version"]), gtc_version=ref.version + ) + g = genotype_sample(read_idat(grn), read_idat(red), ctx, params) + res["cluster_file"] = ctx.cluster_file_name + diffs: dict[str, Any] = {} + for k in ("raw_x", "raw_y", "genotypes", "genotype_scores", "base_calls", "controls_x", "controls_y"): + a, b = _canon(getattr(g, k)), _canon(getattr(ref, k)) + if a.shape != b.shape: + diffs[k] = f"shape {a.shape} vs {b.shape}" + elif a.tobytes() != b.tobytes(): + n = int((a.reshape(len(a), -1) != b.reshape(len(b), -1)).any(axis=1).sum()) + diffs[k] = n + if g.normalization_transforms.tobytes() != ref.normalization_transforms.tobytes(): + diffs["normalization_transforms"] = True + if ref.version >= 4: + for k in ("b_allele_freqs", "logr_ratios"): + if _canon(getattr(g, k)).tobytes() != _canon(getattr(ref, k)).tobytes(): + diffs[k] = int( + (_canon(getattr(g, k)).view(np.uint32) != _canon(getattr(ref, k)).view(np.uint32)).sum() + ) + for k in ( + "call_rate", + "gender", + "logr_dev", + "gc10", + "gc50", + "num_calls", + "num_no_calls", + "num_intensity_only", + "ploidy", + "ploidy_type", + "dx", + "slide_identifier", + ): + a, b = getattr(g, k), getattr(ref, k) + same = a == b or (isinstance(a, float) and np.isnan(a) and np.isnan(b)) + if not same: + diffs[k] = [repr(a), repr(b)] + if "genotype_scores" in diffs: + d = np.abs(g.genotype_scores.astype(np.float64) - ref.genotype_scores.astype(np.float64)) + res["max_score_abs_diff"] = float(np.nanmax(d)) + res["field_diffs"] = diffs + # also evaluate the other sex-estimate generation on the same calls (cheap), to identify the vendor's rule + from simurg.genotype.sex import estimate_sex + from simurg.signal.normalization.apply import normalize, polar + + nx, ny = normalize(g.raw_x, g.raw_y, ctx.lookup, g.normalization_transforms) + _, r_all = polar(nx, ny) + res["sex_by_version"] = { + v: estimate_sex( + g.genotypes, r_all, ctx.clusters.gentrain, ctx.clusters.present, ctx.sex_classes, SexParams(version=v) + ) + for v in (1, 2) + } + res["sex_ref"] = ref.gender_char + for k in ALIGN: + setattr(g, k, getattr(ref, k)) + for k in ("b_allele_freqs", "logr_ratios", "genotype_scores"): + setattr(g, k, _canon(getattr(g, k))) + setattr(ref, k, _canon(getattr(ref, k))) + res["byte_identical"] = encode_gtc(g) == encode_gtc(ref) + res["ref_layout_canonical"] = encode_gtc(ref) == Path(gtc).read_bytes() or None + res["metadata"] = {k: str(getattr(ref, k))[:80] for k in ("autocall_version", "cluster_file", "imaging_date")} + except Exception as e: + res["error"] = f"{type(e).__name__}: {e}" + res["trace"] = traceback.format_exc(limit=4) + return res + + +def main() -> int: + ap = argparse.ArgumentParser() + ap.add_argument("--root", required=True, type=Path) + ap.add_argument("--bpm", required=True, type=Path) + ap.add_argument("--egt-dir", required=True, type=Path) + ap.add_argument("--out", required=True, type=Path) + ap.add_argument("--norm-version", default="1.1.2") + ap.add_argument("--sex-version", type=int, default=1) + ap.add_argument("--workers", type=int, default=os.cpu_count() or 1) + ap.add_argument("--limit", type=int, default=0) + a = ap.parse_args() + + from simurg.formats.discovery import walk_files + + t0 = time.time() + files = walk_files(a.root, (".gtc", ".idat")) + key_re = re.compile(r"(\d+_r\d{2}c\d{2})", re.IGNORECASE) + grn, red, gtcs = {}, {}, [] + for p in files: + low = p.name.lower() + m = key_re.search(p.name) + if not m: + continue + k = m.group(1).lower() + if low.endswith("_grn.idat"): + grn[k] = str(p) + elif low.endswith("_red.idat"): + red[k] = str(p) + elif low.endswith(".gtc"): + gtcs.append((k, str(p))) + triples = [(g, grn[k], red[k]) for k, g in gtcs if k in grn and k in red] + if a.limit: + triples = triples[: a.limit] + with ProcessPoolExecutor( + a.workers, + mp_context=mp.get_context("spawn"), + initializer=_init, + initargs=(str(a.bpm), str(a.egt_dir), a.norm_version, a.sex_version), + ) as ex: + results = list(ex.map(check, triples, chunksize=2)) + + field_counts: dict[str, int] = {} + for r in results: + for k in r.get("field_diffs", {}): + field_counts[k] = field_counts.get(k, 0) + 1 + summary: dict[str, Any] = { + "n_gtc": len(gtcs), + "n_compared": len(results), + "errors": [r for r in results if "error" in r][:20], + "byte_identical": sum(bool(r.get("byte_identical")) for r in results), + "samples_with_no_field_diff": sum(not r.get("field_diffs") and "error" not in r for r in results), + "samples_with_field_diff": field_counts, + "genotype_diff_loci_total": sum( + int(r.get("field_diffs", {}).get("genotypes", 0) or 0) + if isinstance(r.get("field_diffs", {}).get("genotypes", 0), int) + else 0 + for r in results + ), + "score_diff_loci_total": sum( + int(r.get("field_diffs", {}).get("genotype_scores", 0) or 0) + if isinstance(r.get("field_diffs", {}).get("genotype_scores", 0), int) + else 0 + for r in results + ), + "max_score_abs_diff": max((r.get("max_score_abs_diff", 0.0) for r in results), default=0.0), + "cluster_files": sorted({r.get("cluster_file", "") for r in results}), + "sex_agreement_by_version": { + str(v): sum(r.get("sex_by_version", {}).get(v) == r.get("sex_ref") for r in results) for v in (1, 2) + }, + "ref_sex_counts": {s: sum(r.get("sex_ref") == s for r in results) for s in ("M", "F", "U")}, + "examples": [r for r in results if r.get("field_diffs")][:15], + "params": {"norm_version": a.norm_version, "sex_version": a.sex_version}, + "seconds": round(time.time() - t0, 1), + } + summary["passed"] = summary["byte_identical"] == len(results) > 0 + a.out.parent.mkdir(parents=True, exist_ok=True) + a.out.write_text(json.dumps(summary, indent=2, default=str) + "\n") + print( + json.dumps( + { + k: summary[k] + for k in ( + "n_compared", + "byte_identical", + "samples_with_no_field_diff", + "samples_with_field_diff", + "genotype_diff_loci_total", + "score_diff_loci_total", + "passed", + "seconds", + ) + } + ) + ) + return 0 if summary["passed"] else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/unit/test_genotype.py b/tests/unit/test_genotype.py new file mode 100644 index 0000000..3e3ac16 --- /dev/null +++ b/tests/unit/test_genotype.py @@ -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