diff --git a/scripts/cloud/tasks/p5.sh b/scripts/cloud/tasks/p5.sh
new file mode 100755
index 0000000..59d1201
--- /dev/null
+++ b/scripts/cloud/tasks/p5.sh
@@ -0,0 +1,38 @@
+#!/bin/bash
+# Phase P5 bundle on the 1000 Genomes Omni2.5 set (2141 samples):
+# 1. P4 re-check with --preset autocall-1.6.2 (expect byte-identical GTCs except the K-009 samples);
+# 2. `simurg genotype` (legacy preset, GTC v5 so BAF exists) with control tables and cohort QC;
+# 3. P5 validation against the 1000G pedigree and VCF (scripts/validate_p5_large.py).
+# Args:
+set -uo pipefail
+D=$1; OUT=$2
+SRC=s3://1000genomes/release/20130502/supporting/hd_genotype_chip
+mkdir -p $OUT/p4 $OUT/p5 $D/egt
+cp $D/HumanOmni2.5-4v1-Multi_B.egt $D/egt/
+
+echo "=== P4 preset re-check $(date -u +%T) ==="
+/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/p4/summary.json --workers "$(nproc)" 2> $OUT/p4/time.txt \
+ || echo "P4 validation reported differences"
+tail -22 $OUT/p4/time.txt
+
+echo "=== P5 genotype + QC $(date -u +%T) ==="
+GT=/work/gt
+/usr/bin/time -v .venv/bin/simurg genotype --bpm $D/HumanOmni2.5-4v1-Multi_B.bpm --egt $D/HumanOmni2.5-4v1-Multi_B.egt \
+ --idats $D/data --out $GT --preset autocall-1.6.2 --gtc-version 5 --autocall-date "1/1/2026 12:00:00 AM" \
+ --qc-profile genotyping --threads "$(nproc)" --quiet 2> $OUT/p5/genotype_time.txt || echo "genotype exit $?"
+tail -22 $OUT/p5/genotype_time.txt
+cp $GT/sample_summary.tsv $GT/controls_qc.tsv $GT/relatedness.tsv $GT/run_manifest.json $OUT/p5/ 2>/dev/null
+cp $GT/mendel.tsv $OUT/p5/ 2>/dev/null
+cp -r $GT/logs $OUT/p5/ 2>/dev/null
+gzip -c $GT/controls_raw.tsv > $OUT/p5/controls_raw.tsv.gz 2>/dev/null
+df -h /work | tail -1
+
+echo "=== P5 validation $(date -u +%T) ==="
+aws s3 cp --no-sign-request "$SRC/broad_intensities/Omni25_genotypes_2141_samples.b37.v2.vcf.gz" /work/ --only-show-errors
+aws s3 cp --no-sign-request "s3://1000genomes/release/20130502/integrated_call_samples.20130502.ALL.ped" /work/ --only-show-errors
+/usr/bin/time -v .venv/bin/python scripts/validate_p5_large.py --gt-dir $GT --bpm $D/HumanOmni2.5-4v1-Multi_B.bpm \
+ --vcf /work/Omni25_genotypes_2141_samples.b37.v2.vcf.gz --ped /work/integrated_call_samples.20130502.ALL.ped \
+ --out $OUT/p5/summary.json --workers "$(nproc)" 2> $OUT/p5/validate_time.txt || echo "P5 validation reported differences"
+tail -22 $OUT/p5/validate_time.txt
+echo "=== task done $(date -u +%T) ==="
diff --git a/scripts/cloud/trial.sh b/scripts/cloud/trial.sh
index 41cfd31..5e31e5b 100755
--- a/scripts/cloud/trial.sh
+++ b/scripts/cloud/trial.sh
@@ -17,6 +17,8 @@ TYPE=${TYPE:-m7g.4xlarge}
MAX_PRICE=${MAX_PRICE:-0.30} # USD/hour spot cap (on-demand is 0.6528)
MAX_MINUTES=${MAX_MINUTES:-360} # hard cap: 6 hours
DISK_GB=${DISK_GB:-250}
+DISK_MBPS=${DISK_MBPS:-125} # gp3 throughput (125 = free baseline; each extra MB/s ~ $0.00005/h)
+DISK_IOPS=${DISK_IOPS:-3000} # gp3 IOPS (3000 = free baseline)
BUDGET_USD=${BUDGET_USD:-10}
TASK=${TASK:-p3} # scripts/cloud/tasks/.sh
TAG=simurg-trial
@@ -72,7 +74,7 @@ launch)
--iam-instance-profile Name=$ROLE \
--instance-initiated-shutdown-behavior terminate \
--instance-market-options "MarketType=spot,SpotOptions={MaxPrice=$MAX_PRICE,SpotInstanceType=one-time,InstanceInterruptionBehavior=terminate}" \
- --block-device-mappings "DeviceName=/dev/sda1,Ebs={VolumeSize=$DISK_GB,VolumeType=gp3,DeleteOnTermination=true}" \
+ --block-device-mappings "DeviceName=/dev/sda1,Ebs={VolumeSize=$DISK_GB,VolumeType=gp3,Throughput=$DISK_MBPS,Iops=$DISK_IOPS,DeleteOnTermination=true}" \
--metadata-options HttpTokens=required \
--user-data file:///tmp/simurg-userdata.sh \
--tag-specifications "ResourceType=instance,Tags=[{Key=Project,Value=$TAG},{Key=Run,Value=$RUN}]" \
diff --git a/scripts/validate_p5_large.py b/scripts/validate_p5_large.py
new file mode 100644
index 0000000..972af62
--- /dev/null
+++ b/scripts/validate_p5_large.py
@@ -0,0 +1,304 @@
+#!/usr/bin/env python3
+"""Phase P5 cohort-scale validation on the 1000 Genomes Omni2.5 set (2141 samples).
+
+Input is a ``simurg genotype`` output folder (GTC v5 + sample_summary.tsv + controls_qc.tsv + relatedness.tsv).
+The GTCs carry lab sample names, not 1000G IDs, so each sample is first identified by its heterozygosity pattern
+(orientation-free) against the public 1000G Omni2.5 VCF. Then:
+
+* sex estimate vs pedigree sex;
+* parent–offspring pairs of the pedigree are detected as first degree (KING-robust, relatedness.tsv);
+* Mendelian error rates of complete trios;
+* distributions of contamination, control metrics, statuses and flags.
+
+Usage: validate_p5_large.py --gt-dir DIR --bpm F --vcf F.vcf.gz --ped F.ped --out summary.json [--workers N]
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import gzip
+import json
+import multiprocessing as mp
+import os
+import sys
+import time
+from collections import Counter
+from concurrent.futures import ProcessPoolExecutor
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+_W: dict[str, Any] = {}
+
+
+def _init(idx: np.ndarray) -> None: # type: ignore[type-arg]
+ _W["idx"] = idx
+
+
+def _sel(path: str) -> np.ndarray: # type: ignore[type-arg]
+ from simurg.formats.gtc import read_gtc
+
+ return read_gtc(path).genotypes[_W["idx"]]
+
+
+def _trio(t: tuple[str, str, str, str, str, str]) -> dict[str, Any]:
+ from simurg.qc.relatedness import mendel_errors
+
+ k, f, m, pk, pf, pm = t
+ e, n = mendel_errors(_sel(pk), _sel(pf), _sel(pm))
+ return {"child": k, "father": f, "mother": m, "errors": e, "informative": n, "rate": e / n if n else None}
+
+
+def read_vcf_hets(
+ path: Path, names: dict[str, int], n_snps: int
+) -> tuple[list[str], list[int], np.ndarray, np.ndarray]: # type: ignore[type-arg]
+ """First ``n_snps`` common SNPs (0.2 <= AF <= 0.8) that the manifest has: het and called indicators."""
+ samples: list[str] = []
+ loci: list[int] = []
+ het_rows, call_rows = [], []
+ with gzip.open(path, "rt") as fh:
+ for line in fh:
+ if line.startswith("##"):
+ continue
+ f = line.rstrip("\n").split("\t")
+ if line.startswith("#"):
+ samples = f[9:]
+ continue
+ if f[2] not in names:
+ continue
+ info = dict(kv.split("=", 1) for kv in f[7].split(";") if "=" in kv)
+ try:
+ af = float(info.get("AF", "nan").split(",")[0])
+ except ValueError:
+ continue
+ if not 0.2 <= af <= 0.8:
+ continue
+ gts = f[9:]
+ het_rows.append(np.array([g[:3] in ("0/1", "1/0", "0|1", "1|0") for g in gts]))
+ call_rows.append(np.array([g[:1] != "." for g in gts]))
+ loci.append(names[f[2]])
+ if len(loci) >= n_snps:
+ break
+ return samples, loci, np.array(het_rows), np.array(call_rows)
+
+
+def main() -> int:
+ ap = argparse.ArgumentParser()
+ ap.add_argument("--gt-dir", required=True, type=Path)
+ ap.add_argument("--bpm", required=True, type=Path)
+ ap.add_argument("--vcf", required=True, type=Path)
+ ap.add_argument("--ped", required=True, 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("--n-snps", type=int, default=20000)
+ a = ap.parse_args()
+ t0 = time.time()
+
+ from simurg.formats.bpm import read_bpm
+ from simurg.formats.samplesheet import output_stem
+ from simurg.pipeline.qc import autosome_mask
+
+ manifest = read_bpm(a.bpm)
+ names = manifest.name_index()
+ summary = list(csv.DictReader((a.gt_dir / "sample_summary.tsv").open(), delimiter="\t"))
+ ok_rows = [r for r in summary if r["status"] == "OK"]
+ paths = [
+ str(a.gt_dir / "samples" / output_stem(r["sample_name"]) / f"{output_stem(r['sample_name'])}.gtc")
+ for r in ok_rows
+ ]
+ missing = [p for p in paths if not Path(p).exists()]
+ if missing:
+ print(f"GTC paths not found, e.g. {missing[:3]}", file=sys.stderr)
+ return 2
+
+ # 1. identify samples ------------------------------------------------------------------------------
+ vs, loci, vhet, vcall = read_vcf_hets(a.vcf, names, a.n_snps)
+ idx = np.array(loci, dtype=np.int64)
+ with ProcessPoolExecutor(a.workers, mp_context=mp.get_context("spawn"), initializer=_init, initargs=(idx,)) as ex:
+ ours = np.stack(list(ex.map(_sel, paths, chunksize=8))) # samples × snps
+ H = (ours == 2).astype(np.float32)
+ C = (ours > 0).astype(np.float32)
+ V = (vhet.T & vcall.T).astype(np.float32) # vcf samples × snps
+ W = vcall.T.astype(np.float32)
+ joint = C @ W.T
+ agree = H @ V.T + (C - H) @ (W - V).T
+ conc = agree / np.maximum(joint, 1)
+ best = conc.argmax(axis=1)
+ order = np.sort(conc, axis=1)
+ top, second = order[:, -1], order[:, -2]
+ sid2k = {}
+ ambiguous = []
+ for i, r in enumerate(ok_rows):
+ if top[i] >= 0.98 and top[i] - second[i] >= 0.05:
+ sid2k[r["sample_id"]] = vs[best[i]]
+ else:
+ ambiguous.append({"sample_id": r["sample_id"], "top": float(top[i]), "second": float(second[i])})
+ dup_targets = [k for k, c in Counter(sid2k.values()).items() if c > 1]
+
+ # 2. pedigree --------------------------------------------------------------------------------------
+ ped: dict[str, dict[str, str]] = {}
+ with a.ped.open() as fh:
+ for r in csv.DictReader(fh, delimiter="\t"):
+ ped[r["Individual ID"]] = r
+ k2sid = {k: s for s, k in sid2k.items()}
+
+ # 3. sex -------------------------------------------------------------------------------------------
+ sex_cmp = Counter()
+ sex_mismatch = []
+ for r in ok_rows:
+ k = sid2k.get(r["sample_id"])
+ if not k or k not in ped:
+ sex_cmp["no_pedigree"] += 1
+ continue
+ want = {"1": "M", "2": "F"}.get(ped[k]["Gender"], "U")
+ got = r["sex_est"]
+ sex_cmp["match" if got == want else "mismatch"] += 1
+ if got != want:
+ sex_mismatch.append({"sample_id": r["sample_id"], "id": k, "est": got, "ped": want})
+
+ # 4. relatedness -----------------------------------------------------------------------------------
+ rel = list(csv.DictReader((a.gt_dir / "relatedness.tsv").open(), delimiter="\t"))
+ found = {frozenset((p["sample_a"], p["sample_b"])): p for p in rel}
+ po_pairs = set()
+ for k, r in ped.items():
+ for par in (r["Paternal ID"], r["Maternal ID"]):
+ if k in k2sid and par in k2sid:
+ po_pairs.add(frozenset((k2sid[k], k2sid[par])))
+ po_detected = [found[p] for p in po_pairs if p in found and found[p]["degree"] == "first"]
+ po_kin = [float(found[p]["kinship"]) for p in po_pairs if p in found]
+ po_missed = [sorted(p) for p in po_pairs if p not in found or found[p]["degree"] != "first"]
+ sib = set()
+ for k, r in ped.items():
+ for s in (r.get("Siblings") or "").replace(",", " ").split():
+ if k in k2sid and s in k2sid:
+ sib.add(frozenset((k2sid[k], k2sid[s])))
+ close = [p for p in rel if p["degree"] in ("duplicate", "first", "second")]
+ unexplained = [
+ {**p, "id_a": sid2k.get(p["sample_a"], ""), "id_b": sid2k.get(p["sample_b"], "")}
+ for p in close
+ if frozenset((p["sample_a"], p["sample_b"])) not in po_pairs | sib
+ ]
+
+ # 5. Mendelian errors -------------------------------------------------------------------------------
+ auto = np.flatnonzero(autosome_mask(manifest.chroms) & ~manifest.intensity_only)
+ path_of = {r["sample_id"]: p for r, p in zip(ok_rows, paths, strict=True)}
+ trios = []
+ for k, r in ped.items():
+ f, m = r["Paternal ID"], r["Maternal ID"]
+ if k in k2sid and f in k2sid and m in k2sid:
+ trios.append((k, f, m))
+ tasks = [(k, f, m, path_of[k2sid[k]], path_of[k2sid[f]], path_of[k2sid[m]]) for k, f, m in trios]
+ with ProcessPoolExecutor(a.workers, mp_context=mp.get_context("spawn"), initializer=_init, initargs=(auto,)) as ex:
+ mendel = list(ex.map(_trio, tasks, chunksize=2))
+ rates = np.array([x["rate"] for x in mendel if x["rate"] is not None])
+
+ # 6. distributions ---------------------------------------------------------------------------------
+ def q(vals: list[float]) -> dict[str, float]:
+ v = np.array([x for x in vals if np.isfinite(x)])
+ if not len(v):
+ return {}
+ return {
+ k: float(np.percentile(v, p))
+ for k, p in (("min", 0), ("p01", 1), ("median", 50), ("p99", 99), ("max", 100))
+ } | {"n": len(v)}
+
+ def num(x: str) -> float:
+ try:
+ return float(x)
+ except ValueError:
+ return float("nan")
+
+ cqc = list(csv.DictReader((a.gt_dir / "controls_qc.tsv").open(), delimiter="\t"))
+ metric_cols = [
+ c
+ for c in (cqc[0] if cqc else {})
+ if c
+ not in (
+ "sample_id",
+ "array_type",
+ "imaging_date",
+ "sample_plate",
+ "sample_well",
+ "scanner_id",
+ "scanner_version",
+ "bg_offset",
+ "na_metrics",
+ )
+ ]
+ contam = [num(r["contamination"]) for r in ok_rows]
+ top_contam = sorted(
+ (
+ {"sample_id": r["sample_id"], "id": sid2k.get(r["sample_id"], ""), "c": num(r["contamination"])}
+ for r in ok_rows
+ ),
+ key=lambda x: -x["c"] if np.isfinite(x["c"]) else 0,
+ )[:15]
+ flags = Counter(f.split(":")[0] for r in summary for f in r["flags"].split(";") if f)
+ codes = Counter(c for r in summary for c in r["control_flag_codes"].split(";") if c)
+
+ out: dict[str, Any] = {
+ "n_samples": len(summary),
+ "n_ok": len(ok_rows),
+ "identification": {
+ "snps": len(loci),
+ "identified": len(sid2k),
+ "ambiguous": ambiguous[:20],
+ "n_ambiguous": len(ambiguous),
+ "ids_matched_twice": dup_targets,
+ "top_concordance": q(top.tolist()),
+ "margin": q((top - second).tolist()),
+ },
+ "sex_vs_pedigree": dict(sex_cmp),
+ "sex_mismatches": sex_mismatch[:30],
+ "parent_offspring": {
+ "pairs_in_pedigree": len(po_pairs),
+ "detected_first_degree": len(po_detected),
+ "kinship": q(po_kin),
+ "missed": po_missed[:20],
+ },
+ "close_pairs_total": len(close),
+ "close_pairs_unexplained_by_pedigree": len(unexplained),
+ "unexplained_examples": unexplained[:25],
+ "mendel": {
+ "trios": len(mendel),
+ "error_rate": q(rates.tolist()),
+ "worst": sorted(mendel, key=lambda x: -(x["rate"] or 0))[:10],
+ },
+ "contamination": q(contam) | {"above_0.02": int(sum(c > 0.02 for c in contam if np.isfinite(c)))},
+ "contamination_top": top_contam,
+ "functional_status": dict(Counter(r["functional_status"] for r in summary)),
+ "control_status": dict(Counter(r["control_status"] for r in summary)),
+ "control_flag_codes": dict(codes),
+ "flags": dict(flags),
+ "control_metrics": {c: q([num(r[c]) for r in cqc]) for c in metric_cols},
+ "het_rate": q([num(r["het_rate"]) for r in ok_rows]),
+ "seconds": round(time.time() - t0, 1),
+ }
+ po_rate = len(po_detected) / max(len(po_pairs), 1)
+ sex_rate = sex_cmp["match"] / max(sex_cmp["match"] + sex_cmp["mismatch"], 1)
+ out["checks"] = {
+ "identified_ge_99pct": len(sid2k) >= 0.99 * len(ok_rows),
+ "sex_agreement_ge_99pct": sex_rate >= 0.99,
+ "parent_offspring_detected_ge_99pct": po_rate >= 0.99,
+ "median_mendel_rate_lt_0.001": bool(len(rates) and float(np.median(rates)) < 0.001),
+ }
+ out["passed"] = all(out["checks"].values())
+ a.out.parent.mkdir(parents=True, exist_ok=True)
+ a.out.write_text(json.dumps(out, indent=2, default=str) + "\n")
+ print(json.dumps({k: out[k] for k in ("n_ok", "sex_vs_pedigree", "checks", "passed", "seconds")}, default=str))
+ print(
+ json.dumps(
+ {
+ "identification": {k: out["identification"][k] for k in ("identified", "n_ambiguous")},
+ "po": {k: out["parent_offspring"][k] for k in ("pairs_in_pedigree", "detected_first_degree")},
+ "mendel_trios": out["mendel"]["trios"],
+ }
+ )
+ )
+ return 0 if out["passed"] else 1
+
+
+if __name__ == "__main__":
+ sys.exit(main())
diff --git a/src/simurg/qc/contamination.py b/src/simurg/qc/contamination.py
index a8726b3..433932b 100644
--- a/src/simurg/qc/contamination.py
+++ b/src/simurg/qc/contamination.py
@@ -26,6 +26,7 @@ import numpy.typing as npt
F32 = np.float32
CHUNK = 4096
+MAP_BATCH = 32 # chunks in flight per map call
@dataclass
@@ -148,10 +149,19 @@ def baf_regress(
t = float(F32(threshold)) # the reference stores the threshold as float: 0.1f = 0.10000000149…
keep = np.flatnonzero(np.isfinite(af) & (af >= t) & (af <= 1.0 - t))
chunks = [keep[s : s + CHUNK] for s in range(0, len(keep), CHUNK)]
- tasks = ((np.asarray(genotypes[r]), np.asarray(baf[r], dtype=F32), adjust_min, median) for r in chunks)
+
+ def adjustments() -> Iterable[tuple[npt.NDArray[np.float32], npt.NDArray[np.float32]]]:
+ # bounded batches: an executor's map submits its whole input at once
+ for b in range(0, len(chunks), MAP_BATCH):
+ tasks = [
+ (np.asarray(genotypes[r]), np.asarray(baf[r], dtype=F32), adjust_min, median)
+ for r in chunks[b : b + MAP_BATCH]
+ ]
+ yield from map_fn(cluster_adjustments, tasks)
+
sx2, sxy, sx, sy = (np.zeros(n, F32) for _ in range(4))
cnt = np.zeros(n, np.int64)
- for rows, (adj_aa, adj_bb) in zip(chunks, map_fn(cluster_adjustments, tasks), strict=True):
+ for rows, (adj_aa, adj_bb) in zip(chunks, adjustments(), strict=True):
g = np.asarray(genotypes[rows])
b = np.asarray(baf[rows], dtype=F32)
p = af[rows]