From bc24148e31afcf70f70e31ac044f18eb8d2b8408 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 01:39:49 +0100 Subject: [PATCH] v2 CN: truth validation script (NYGC 30x SV freeze V3) and cloud task cnv --- scripts/cloud/tasks/cnv.sh | 60 ++++++++++++ scripts/v2/cn_validate_truth.py | 166 ++++++++++++++++++++++++++++++++ 2 files changed, 226 insertions(+) create mode 100644 scripts/cloud/tasks/cnv.sh create mode 100644 scripts/v2/cn_validate_truth.py diff --git a/scripts/cloud/tasks/cnv.sh b/scripts/cloud/tasks/cnv.sh new file mode 100644 index 0000000..7eaf0fe --- /dev/null +++ b/scripts/cloud/tasks/cnv.sh @@ -0,0 +1,60 @@ +#!/bin/bash +# v2 CN/LOH validation on GSE165845 (360 GSA arrays, 120 1000 Genomes individuals x 3): simurg cn call vs the NYGC +# 30x structural-variant callset (1KGP_3202 GATK-SV/svtools freeze V3) and between replicates. +# Args: . Run with OMNI=0. +set -uo pipefail +OUT=$2 +W=/work/cnv; R=$OUT/cnv; mkdir -p $W $R +N=$(nproc) +export DEBIAN_FRONTEND=noninteractive +apt-get install -y -q bcftools tabix gcc > /dev/null +curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal > /dev/null 2>&1 +export PATH="$HOME/.cargo/bin:$PATH" +VIRTUAL_ENV=$PWD/.venv uv pip install -q ./rust/simurg-core || echo "!! kernels not built (pure Python)" + +SV=https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20210124.SV_Illumina_Integration/1KGP_3202.gatksv_svtools_novelins.freeze_V3.wAF.vcf.gz +( curl -sSfL --retry 5 -o $W/sv.vcf.gz $SV && curl -sSfL --retry 5 -o $W/sv.vcf.gz.tbi $SV.tbi \ + && echo "SV callset $(du -h $W/sv.vcf.gz | cut -f1) at $(date -u +%T)" ) > $W/sv.log 2>&1 & +SV_PID=$! +( bash scripts/fetch_data.sh p6 ) > $W/fasta.log 2>&1 & +FASTA_PID=$! + +echo "=== data $(date -u +%T) ===" +GSA=https://webdata.illumina.com/downloads/productfiles/global-screening-array/v1-0 +mkdir -p $W/gsa && cd $W/gsa +for z in infinium-global-screening-array-24-v1-0-c2-manifest-file-bpm-build38.zip \ + infinium-global-screening-array-v1-0-c1-cluster-file.zip; do + curl -sSfL --retry 5 -O "$GSA/$z" && unzip -oq $z +done +cd - > /dev/null +BPM=$W/gsa/GSA-24v1-0_C2.bpm; EGT=$W/gsa/GSA-24v1-0_C1_ClusterFile.egt +.venv/bin/python scripts/p8/prepare.py --tsv scripts/p8/gse165845_samples.tsv --dl $W/geo --out $W/idats --workers 24 \ + | tail -2 || exit 1 + +echo "=== genotype $(date -u +%T) ===" +.venv/bin/simurg genotype --bpm $BPM --egt $EGT --idats $W/idats --out $W/gt --threads $N --quiet +cp $W/idats/samples_map.tsv $W/gt/sample_summary.tsv $R/ + +echo "=== model $(date -u +%T) ===" +wait $FASTA_PID; tail -1 $W/fasta.log +/usr/bin/time -v .venv/bin/simurg cn build-model --bpm $BPM --fasta data/p6/GRCh38.primary.fa --egt $EGT \ + --reference-gtcs $W/gt --max-reference 120 --out $W/gsa1.simurg-cyto 2> $R/model_time.txt | grep -v "GC windows" +grep -E "Elapsed|Maximum resident" $R/model_time.txt + +echo "=== cn call $(date -u +%T) ===" +/usr/bin/time -v .venv/bin/simurg cn call --bpm $BPM --model $W/gsa1.simurg-cyto --gtcs $W/gt \ + --fasta data/p6/GRCh38.primary.fa --out $W/cn --threads $N --quiet 2> $R/call_time.txt +grep -E "Elapsed|Maximum resident" $R/call_time.txt +.venv/bin/simurg cn annotate --vcfs $W/cn --quiet || echo "annotate failed" +cp $W/cn/cn_sample_summary.tsv $W/cn/cn_call_manifest.json $R/ + +echo "=== truth $(date -u +%T) ===" +wait $SV_PID; cat $W/sv.log +cut -f3 $W/idats/samples_map.tsv | tail -n +2 | sort -u > $W/inds.txt +bcftools view -S $W/inds.txt --force-samples -i 'INFO/SVTYPE="DEL" || INFO/SVTYPE="DUP"' -Ou $W/sv.vcf.gz \ + | bcftools query -i 'GT="alt"' -f '[%CHROM\t%POS\t%INFO/END\t%INFO/SVTYPE\t%SAMPLE\n]' > $W/truth.tsv +wc -l < $W/truth.tsv +.venv/bin/python scripts/v2/cn_validate_truth.py --calls $W/cn --map $W/idats/samples_map.tsv --truth $W/truth.tsv \ + --bpm $BPM --out $R/validation +gzip -c $W/truth.tsv > $R/truth.tsv.gz +echo "=== task done $(date -u +%T) ===" diff --git a/scripts/v2/cn_validate_truth.py b/scripts/v2/cn_validate_truth.py new file mode 100644 index 0000000..0ddaf2b --- /dev/null +++ b/scripts/v2/cn_validate_truth.py @@ -0,0 +1,166 @@ +#!/usr/bin/env python3 +"""Validate `simurg cn call` on GSE165845 (120 1000 Genomes individuals × 3 GSA arrays) against the NYGC 30x +structural-variant callset (1KGP_3202 GATK-SV/svtools freeze V3) and between replicates. + + python scripts/v2/cn_validate_truth.py --calls CN_OUT --map samples_map.tsv --truth truth.tsv --bpm M.bpm --out DIR + +``truth.tsv``: ``chrom start end svtype sample`` for non-reference DEL/DUP genotypes of the 120 individuals +(``bcftools query``; chr prefix removed). Matching: reciprocal overlap ≥ 50 % and the same class (DEL ↔ losses +incl. mosaic, DUP ↔ gains). Only truth events with ≥ 10 array probes are counted as detectable. +""" + +from __future__ import annotations + +import argparse +import csv +import json +import statistics +from collections import defaultdict +from pathlib import Path + +import numpy as np + +from simurg.formats.bpm import read_bpm + +BINS = [(25_000, 50_000), (50_000, 100_000), (100_000, 500_000), (500_000, 1_000_000), (1_000_000, 10**9)] +CLASS = {"DEL": "DEL", "MDEL": "DEL", "DUP": "DUP", "MDUP": "DUP"} + + +def ro(a: tuple[int, int], b: tuple[int, int]) -> float: + ov = min(a[1], b[1]) - max(a[0], b[0]) + 1 + if ov <= 0: + return 0.0 + return min(ov / (a[1] - a[0] + 1), ov / (b[1] - b[0] + 1)) + + +def bin_of(size: int) -> str: + for lo, hi in BINS: + if lo <= size < hi: + return f"{lo // 1000}-{hi // 1000 if hi < 10**9 else 'inf'}kb" + return "<25kb" + + +def main() -> int: + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--calls", required=True, type=Path) + ap.add_argument("--map", required=True, type=Path) + ap.add_argument("--truth", required=True, type=Path) + ap.add_argument("--bpm", required=True) + ap.add_argument("--min-qual", type=float, default=20) + ap.add_argument("--out", required=True, type=Path) + a = ap.parse_args() + a.out.mkdir(parents=True, exist_ok=True) + man = read_bpm(a.bpm) + chroms = np.array(["X" if c == "XY" else c for c in man.chroms]) + pos = man.map_infos.astype(np.int64) + probe_pos = {c: np.sort(pos[(chroms == c) & (pos > 0)]) for c in set(chroms.tolist())} + + def n_probes(c: str, s: int, e: int) -> int: + p = probe_pos.get(c) + return 0 if p is None else int(np.searchsorted(p, e, "right") - np.searchsorted(p, s, "left")) + + with a.map.open() as fh: + smap = {r["sample_id"]: r for r in csv.DictReader(fh, delimiter="\t")} + side = json.loads((a.calls / "cn_call_manifest.json").read_text()) + calls: dict[str, list[dict]] = {} + for s in side["samples"]: + if s.get("status") != "OK": + continue + calls[s["sample_id"]] = [e for e in s["events"] if e["quality"] >= a.min_qual and e["filters"] == ["PASS"]] + truth: dict[str, list[tuple[str, int, int, str]]] = defaultdict(list) + with a.truth.open() as fh: + for line in fh: + c, s_, e_, t, ind = line.rstrip("\n").split("\t") + c = c[3:] if c.startswith("chr") else c + if t in ("DEL", "DUP") and int(e_) - int(s_) + 1 >= 25_000: + truth[ind].append((c, int(s_), int(e_), t)) + + # sensitivity per class and size bin + sens: dict[str, dict[str, list[int]]] = defaultdict(lambda: defaultdict(lambda: [0, 0])) + prec = defaultdict(lambda: [0, 0]) + per_sample = [] + for sid, evs in calls.items(): + ind = smap.get(sid, {}).get("individual") + if ind is None: + continue + tr = [t for t in truth.get(ind, []) if n_probes(t[0], t[1], t[2]) >= 10] + hit_t = 0 + for c, s_, e_, t in tr: + ok = any(CLASS.get(e["call"]) == t and e["chrom"] == c and ro((s_, e_), (e["start"], e["end"])) >= 0.5 + for e in evs) # fmt: skip + cell = sens[t][bin_of(e_ - s_ + 1)] + cell[0] += ok + cell[1] += 1 + hit_t += ok + hit_c = 0 + cnv_calls = [e for e in evs if e["call"] in CLASS] + for e in cnv_calls: + ok = any(t == CLASS[e["call"]] and c == e["chrom"] and ro((s_, e_), (e["start"], e["end"])) >= 0.5 + for c, s_, e_, t in truth.get(ind, [])) # fmt: skip + p = prec[CLASS[e["call"]]] + p[0] += ok + p[1] += 1 + hit_c += ok + per_sample.append({"sample": sid, "individual": ind, "truth_detectable": len(tr), "truth_found": hit_t, + "cnv_calls": len(cnv_calls), "calls_in_truth": hit_c, + "loh_calls": sum(e["call_type"] == "LOH" for e in evs)}) # fmt: skip + + # replicate reproducibility of CNV and LOH calls + by_ind: dict[str, list[str]] = defaultdict(list) + for sid in calls: + if sid in smap: + by_ind[smap[sid]["individual"]].append(sid) + rep = {"CNV": [0, 0], "LOH": [0, 0]} + + def kind_of(e: dict) -> str | None: + if e["call_type"] == "LOH": + return "LOH" + return "CNV" if e["call"] in CLASS else None + + def matched(e: dict, others: list[dict], kind: str) -> bool: + return any( + o["chrom"] == e["chrom"] + and kind_of(o) == kind + and ro((o["start"], o["end"]), (e["start"], e["end"])) >= 0.5 + for o in others + ) + + for sids in by_ind.values(): + for i, s1 in enumerate(sids): + for e in calls[s1]: + kind = kind_of(e) + if kind is None: + continue + for j, s2 in enumerate(sids): + if j != i: + rep[kind][0] += matched(e, calls[s2], kind) + rep[kind][1] += 1 + + def rate(v: list[int]) -> float | None: + return round(v[0] / v[1], 4) if v[1] else None + + summary = { + "samples": len(per_sample), + "min_qual": a.min_qual, + "sensitivity": { + t: {b: {"found": v[0], "truth": v[1], "rate": rate(v)} for b, v in sorted(d.items())} + for t, d in sens.items() + }, + "precision": {t: {"in_truth": v[0], "calls": v[1], "rate": rate(v)} for t, v in prec.items()}, + "replicate_reproducibility": { + k: {"matched": v[0], "comparisons": v[1], "rate": rate(v)} for k, v in rep.items() + }, + "calls_per_sample_median": statistics.median([p["cnv_calls"] for p in per_sample]) if per_sample else None, + "loh_calls_per_sample_median": statistics.median([p["loh_calls"] for p in per_sample]) if per_sample else None, + } + (a.out / "summary.json").write_text(json.dumps(summary, indent=1)) + with (a.out / "per_sample.tsv").open("w") as fh: + w = csv.DictWriter(fh, fieldnames=list(per_sample[0]) if per_sample else ["sample"], delimiter="\t") + w.writeheader() + w.writerows(per_sample) + print(json.dumps(summary, indent=1)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main())