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]