65 lines
2.9 KiB
Python
65 lines
2.9 KiB
Python
#!/usr/bin/env python3
|
|
"""One-line-per-call-type summary of a tuning round: calls, fraction in NYGC truth, fraction seen in the other
|
|
replicate; plus truth sensitivity by size bin (QUAL >= 20, PASS or SampleQuality)."""
|
|
|
|
from __future__ import annotations
|
|
|
|
import csv
|
|
import json
|
|
import sys
|
|
from collections import defaultdict
|
|
from pathlib import Path
|
|
|
|
|
|
def ro(a: tuple[int, int], b: tuple[int, int]) -> float:
|
|
ov = min(a[1], b[1]) - max(a[0], b[0]) + 1
|
|
return 0.0 if ov <= 0 else min(ov / (a[1] - a[0] + 1), ov / (b[1] - b[0] + 1))
|
|
|
|
|
|
def main() -> int:
|
|
run, mapf = Path(sys.argv[1]), Path(sys.argv[2])
|
|
with mapf.open() as fh:
|
|
smap = {r["sample_id"]: r["individual"] for r in csv.DictReader(fh, delimiter="\t")}
|
|
truth: dict[str, list[tuple[str, int, int, str]]] = defaultdict(list)
|
|
with open("data/v2cn/cloud1/truth.tsv") as fh:
|
|
for line in fh:
|
|
c, s, e, t, i = line.rstrip("\n").split("\t")
|
|
truth[i].append((c.removeprefix("chr"), int(s), int(e), t))
|
|
side = json.loads((run / "cn_call_manifest.json").read_text())
|
|
calls = {s["sample_id"]: [e for e in s["events"] if e["quality"] >= 20 and set(e["filters"]) <= {"PASS",
|
|
"SampleQuality"}] for s in side["samples"] if s.get("status") == "OK"} # fmt: skip
|
|
cls = {"DEL": "DEL", "MDEL": "DEL", "DUP": "DUP", "MDUP": "DUP"}
|
|
by_ind: dict[str, list[str]] = defaultdict(list)
|
|
for sid in calls:
|
|
by_ind[smap[sid]].append(sid)
|
|
st: dict[str, list[int]] = defaultdict(lambda: [0, 0, 0, 0])
|
|
for sid, evs in calls.items():
|
|
ind = smap[sid]
|
|
others = [calls[x] for x in by_ind[ind] if x != sid]
|
|
for e in evs:
|
|
k = e["call"]
|
|
v = st[k]
|
|
v[0] += 1
|
|
if k in cls:
|
|
v[1] += any(t == cls[k] and c == e["chrom"] and ro((s, en), (e["start"], e["end"])) >= 0.5
|
|
for c, s, en, t in truth[ind]) # fmt: skip
|
|
for oc in others:
|
|
v[2] += any(o["chrom"] == e["chrom"] and o["call"] == k and ro((o["start"], o["end"]),
|
|
(e["start"], e["end"])) >= 0.5 for o in oc) # fmt: skip
|
|
v[3] += 1
|
|
print(f"{run.name}: {len(calls)} samples")
|
|
for k, v in sorted(st.items()):
|
|
tr = f"{v[1] / v[0]:.2f}" if k in cls and v[0] else " - "
|
|
rp = f"{v[2] / v[3]:.2f}" if v[3] else " - "
|
|
print(f" {k:8s} calls {v[0]:5d} in-truth {tr} replicated {rp}")
|
|
val = json.loads((run / "val" / "summary.json").read_text())
|
|
for t, bins in val["sensitivity"].items():
|
|
print(f" sensitivity {t}: " + ", ".join(f"{b} {x['found']}/{x['truth']}" for b, x in sorted(bins.items())))
|
|
for t, bins in val.get("sensitivity_by_carrier_frequency", {}).items():
|
|
print(f" by frequency {t}: " + ", ".join(f"{b} {x['found']}/{x['truth']}" for b, x in sorted(bins.items())))
|
|
return 0
|
|
|
|
|
|
if __name__ == "__main__":
|
|
raise SystemExit(main())
|