kemal
6bea6c4930
scripts/cloud/trial.sh launches one spot instance that validates the parsers on the 1000 Genomes Omni2.5 set (2141 samples, vendor GTCs) streamed from the public S3 bucket, uploads the summary and terminates itself. Cost guards: spot price cap, hard time limit (shutdown = terminate), EXIT-trap shutdown, CloudWatch idle alarm that terminates from outside, disk deleted with the instance, monthly budget alert, and an `audit` command listing anything still billable. scripts/validate_p1_large.py compares every GTC with BeadArrayFiles, checks byte-identical round-trips, IDAT pairs and BPM/EGT/CSV.
226 lines
9.3 KiB
Python
226 lines
9.3 KiB
Python
#!/usr/bin/env python3
|
|
"""Phase P1 validation on a large data set (e.g. 1000 Genomes Omni2.5, 2141 samples with vendor GTCs).
|
|
|
|
For every GTC: parse, compare every field with Illumina BeadArrayFiles (BSD-2) and test the byte-identical
|
|
round-trip. For every IDAT pair: parse, pair-check and compare the address set with the manifest. BPM/EGT/CSV are
|
|
compared with BeadArrayFiles once. Writes a JSON summary; exit code 0 only if everything matched.
|
|
|
|
Usage: validate_p1_large.py --root DIR --bpm F --egt F [--csv F] --out summary.json [--workers N] [--limit N]
|
|
"""
|
|
|
|
from __future__ import annotations
|
|
|
|
import argparse
|
|
import json
|
|
import os
|
|
import sys
|
|
import time
|
|
import traceback
|
|
from concurrent.futures import ProcessPoolExecutor
|
|
from pathlib import Path
|
|
from typing import Any
|
|
|
|
import numpy as np
|
|
|
|
ROOT = Path(__file__).resolve().parents[1]
|
|
sys.path.insert(0, str(ROOT / "refs" / "BeadArrayFiles"))
|
|
|
|
|
|
def check_gtc(path: str) -> dict[str, Any]:
|
|
from module.GenotypeCalls import GenotypeCalls # type: ignore[import-not-found]
|
|
|
|
from simurg.formats.gtc import encode_gtc, read_gtc
|
|
|
|
res: dict[str, Any] = {"file": path, "ok": True, "diffs": []}
|
|
try:
|
|
g = read_gtc(path, check_write_complete=False)
|
|
raw = Path(path).read_bytes()
|
|
res["version"] = g.version
|
|
res["autocall_version"] = g.autocall_version
|
|
res["roundtrip_identical"] = encode_gtc(g) == raw
|
|
if not res["roundtrip_identical"]:
|
|
res["roundtrip_canonical_identical"] = encode_gtc(g, canonical=True) == raw
|
|
res["size_ours_vs_file"] = [len(encode_gtc(g)), len(raw)]
|
|
b = GenotypeCalls(path, check_write_complete=False)
|
|
pairs = {
|
|
"raw_x": (g.raw_x, b.get_raw_x_intensities()),
|
|
"raw_y": (g.raw_y, b.get_raw_y_intensities()),
|
|
"genotypes": (g.genotypes, np.array(b.get_genotypes(), dtype=np.uint8)),
|
|
"scores": (g.genotype_scores, b.get_genotype_scores()),
|
|
}
|
|
if g.version >= 4:
|
|
pairs["baf"] = (g.b_allele_freqs, b.get_ballele_freqs())
|
|
pairs["lrr"] = (g.logr_ratios, b.get_logr_ratios())
|
|
for k, (a, c) in pairs.items():
|
|
if not np.array_equal(np.asarray(a), np.asarray(c), equal_nan=a.dtype.kind == "f"):
|
|
res["diffs"].append(k)
|
|
scalars = {
|
|
"call_rate": (g.call_rate, b.get_call_rate()),
|
|
"gender": (g.gender, b.get_gender()),
|
|
"logr_dev": (g.logr_dev, b.get_logr_dev()),
|
|
"gc10": (g.gc10, b.get_gc10()),
|
|
"gc50": (g.gc50, b.get_gc50()),
|
|
"num_calls": (g.num_calls, b.get_num_calls()),
|
|
"sample_name": (g.sample_name, b.get_sample_name()),
|
|
"base_calls": (g.base_call_strings(), [x.decode() for x in b.get_base_calls()]),
|
|
}
|
|
for k, (a, c) in scalars.items():
|
|
same = (a == c) or (isinstance(a, float) and np.isnan(a) and np.isnan(c))
|
|
if not same:
|
|
res["diffs"].append(k)
|
|
xf_b = [
|
|
(t.offset_x, t.offset_y, t.scale_x, t.scale_y, t.shear, t.theta) for t in b.get_normalization_transforms()
|
|
]
|
|
xf_g = [
|
|
tuple(float(x[k]) for k in ("offset_x", "offset_y", "scale_x", "scale_y", "shear", "theta"))
|
|
for x in g.normalization_transforms
|
|
]
|
|
if [tuple(float(v) for v in t) for t in xf_b] != xf_g:
|
|
res["diffs"].append("normalization_transforms")
|
|
res["ok"] = not res["diffs"]
|
|
except Exception as e: # noqa: BLE001 - report everything
|
|
res.update(ok=False, error=f"{type(e).__name__}: {e}", trace=traceback.format_exc(limit=3))
|
|
return res
|
|
|
|
|
|
def check_idat_pair(args: tuple[str, str]) -> dict[str, Any]:
|
|
from simurg.formats.idat import check_pair, read_idat
|
|
|
|
grn, red = args
|
|
res: dict[str, Any] = {"grn": grn, "ok": True}
|
|
try:
|
|
g, r = read_idat(grn), read_idat(red)
|
|
check_pair(g, r)
|
|
res.update(
|
|
n=g.n_addresses,
|
|
chip=g.chip_type,
|
|
midblock_itemsize=g.mid_block_itemsize,
|
|
scan=g.scan_time_iso,
|
|
name_ok=(g.barcode in Path(grn).name),
|
|
)
|
|
except Exception as e: # noqa: BLE001
|
|
res.update(ok=False, error=f"{type(e).__name__}: {e}")
|
|
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", required=True, type=Path)
|
|
ap.add_argument("--csv", type=Path)
|
|
ap.add_argument("--out", required=True, type=Path)
|
|
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.bpm import read_bpm
|
|
from simurg.formats.discovery import walk_files
|
|
from simurg.formats.egt import read_egt
|
|
|
|
t0 = time.time()
|
|
summary: dict[str, Any] = {"root": str(a.root), "started": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime())}
|
|
|
|
gtcs = [str(p) for p in walk_files(a.root, (".gtc",))]
|
|
idats = walk_files(a.root, (".idat", ".idat.gz"))
|
|
if a.limit:
|
|
gtcs = gtcs[: a.limit]
|
|
grn = {p.name.lower().replace("_grn.idat", ""): p for p in idats if "_grn.idat" in p.name.lower()}
|
|
red = {p.name.lower().replace("_red.idat", ""): p for p in idats if "_red.idat" in p.name.lower()}
|
|
pairs = [(str(grn[k]), str(red[k])) for k in sorted(set(grn) & set(red))]
|
|
if a.limit:
|
|
pairs = pairs[: a.limit]
|
|
summary["n_gtc"] = len(gtcs)
|
|
summary["n_idat_pairs"] = len(pairs)
|
|
summary["unpaired_idats"] = len(set(grn) ^ set(red))
|
|
|
|
with ProcessPoolExecutor(a.workers) as ex:
|
|
gres = list(ex.map(check_gtc, gtcs, chunksize=4))
|
|
ires = list(ex.map(check_idat_pair, pairs, chunksize=8))
|
|
summary["gtc_ok"] = sum(r["ok"] for r in gres)
|
|
summary["gtc_roundtrip_identical"] = sum(bool(r.get("roundtrip_identical")) for r in gres)
|
|
summary["gtc_versions"] = sorted({r.get("version") for r in gres if "version" in r})
|
|
summary["gtc_autocall_versions"] = sorted({r.get("autocall_version") for r in gres if "autocall_version" in r})
|
|
summary["gtc_failures"] = [r for r in gres if not r["ok"]][:50]
|
|
summary["gtc_roundtrip_failures"] = [r for r in gres if r.get("roundtrip_identical") is False][:20]
|
|
summary["idat_ok"] = sum(r["ok"] for r in ires)
|
|
summary["idat_failures"] = [r for r in ires if not r["ok"]][:50]
|
|
summary["idat_midblock_itemsizes"] = sorted({r.get("midblock_itemsize") for r in ires if r["ok"]})
|
|
summary["idat_chip_types"] = sorted({r.get("chip") for r in ires if r["ok"]})
|
|
|
|
# manifest / cluster file vs BeadArrayFiles
|
|
import module # type: ignore[import-not-found]
|
|
|
|
m = read_bpm(a.bpm)
|
|
b = module.BeadPoolManifest(str(a.bpm))
|
|
summary["bpm"] = {
|
|
"num_loci": m.num_loci,
|
|
"names": m.names == b.names,
|
|
"snps": m.snps == b.snps,
|
|
"chroms": m.chroms == b.chroms,
|
|
"map_infos": m.map_infos.tolist() == b.map_infos,
|
|
"address_a": m.address_a.tolist() == b.addresses,
|
|
"assay_types": m.assay_types.tolist() == b.assay_types,
|
|
"normalization_lookups": m.normalization_lookups().tolist() == b.normalization_lookups,
|
|
"n_bins": int(len(np.unique(m.norm_bin_ids()))),
|
|
}
|
|
e = read_egt(a.egt)
|
|
with open(a.egt, "rb") as fh:
|
|
be = module.ClusterFile.read_cluster_file(fh)
|
|
bad = 0
|
|
|
|
def same(a: float, b: float) -> bool:
|
|
return bool(np.float32(a) == np.float32(b) or (np.isnan(a) and np.isnan(b)))
|
|
|
|
for i, name in enumerate(e.names):
|
|
rec = be.get_record(name)
|
|
st = (rec.aa_cluster_stats, rec.ab_cluster_stats, rec.bb_cluster_stats)
|
|
if not all(
|
|
same(e.theta_mean[i, k], s.theta_mean)
|
|
and same(e.r_mean[i, k], s.r_mean)
|
|
and same(e.theta_dev[i, k], s.theta_dev)
|
|
and same(e.r_dev[i, k], s.r_dev)
|
|
and e.n[i, k] == s.N
|
|
for k, s in enumerate(st)
|
|
):
|
|
bad += 1
|
|
summary["egt"] = {
|
|
"num_records": e.num_records,
|
|
"dbv": e.data_block_version,
|
|
"extra_floats": e.extra_per_record,
|
|
"record_mismatches": bad,
|
|
}
|
|
if a.csv:
|
|
from simurg.formats.csv_manifest import compare_manifests, read_csv_manifest
|
|
|
|
summary["bpm_csv_problems"] = compare_manifests(m, read_csv_manifest(a.csv, load_sequences=False))
|
|
if gtcs:
|
|
g0 = __import__("simurg.formats.gtc", fromlist=["read_gtc"]).read_gtc(gtcs[0], check_write_complete=False)
|
|
summary["gtc_vs_bpm_num_snps"] = [g0.num_snps, m.num_loci]
|
|
summary["gtc_transforms_vs_bins"] = [len(g0.normalization_transforms), summary["bpm"]["n_bins"]]
|
|
|
|
summary["seconds"] = round(time.time() - t0, 1)
|
|
passed = (
|
|
summary["gtc_ok"] == len(gtcs)
|
|
and summary["idat_ok"] == len(pairs)
|
|
and all(v for k, v in summary["bpm"].items() if isinstance(v, bool))
|
|
and summary["egt"]["record_mismatches"] == 0
|
|
and not summary.get("bpm_csv_problems")
|
|
)
|
|
summary["passed"] = passed
|
|
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_gtc", "gtc_ok", "gtc_roundtrip_identical", "n_idat_pairs", "idat_ok", "passed", "seconds")
|
|
}
|
|
)
|
|
)
|
|
return 0 if passed else 1
|
|
|
|
|
|
if __name__ == "__main__":
|
|
sys.exit(main())
|