SimurgArray/scripts/validate_p6_vcf.py

#!/usr/bin/env python3
"""Phase P6 acceptance: SimurgArray VCFs vs a reference converter's multi-sample VCF (gtc2vcf or GTCtoVCF).

Both sides are left-aligned with ``bcftools norm -f`` and compared per probe (ID) and sample as unordered allele
sets (a haploid call ``a`` equals ``a/a``). Every difference is put into a category; the acceptance criterion
(plan §18 P6) is zero *unexplained* genotype differences.

Usage: validate_p6_vcf.py --ours DIR --oracle F.vcf.gz --fasta F --bcftools PATH --out summary.json
"""

from __future__ import annotations

import argparse
import json
import subprocess
import sys
from collections import Counter, defaultdict
from pathlib import Path
from typing import Any


def query(bcftools: str, fasta: str, path: Path) -> list[list[str]]:
    norm = subprocess.Popen([bcftools, "norm", "-f", fasta, "-c", "w", "-Ou", str(path)], stdout=subprocess.PIPE,
                            stderr=subprocess.DEVNULL)  # fmt: skip
    out = subprocess.run(
        [bcftools, "query", "-f", "%ID\t%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n", "-"],
        stdin=norm.stdout, capture_output=True, text=True, check=True,
    ).stdout  # fmt: skip
    norm.wait()
    return [ln.split("\t") for ln in out.splitlines()]


def alleles_of(gt: str, ref: str, alt: str) -> tuple[str, ...] | None:
    if gt in (".", "./.", ".|."):
        return None
    al = [ref, *alt.split(",")]
    idx = gt.replace("|", "/").split("/")
    if any(i == "." for i in idx):
        return None
    got = tuple(sorted(al[int(i)] for i in idx))
    return got * 2 if len(got) == 1 else got


def main() -> int:
    ap = argparse.ArgumentParser()
    ap.add_argument("--ours", required=True, type=Path)
    ap.add_argument("--oracle", required=True, type=Path)
    ap.add_argument("--fasta", required=True)
    ap.add_argument("--bcftools", required=True)
    ap.add_argument("--out", required=True, type=Path)
    ap.add_argument("--name", default="oracle")
    ap.add_argument("--bpm", type=Path, help="manifest, to recognise opposite-strand indel designs (EDGE-012)")
    a = ap.parse_args()

    samples = subprocess.run([a.bcftools, "query", "-l", str(a.oracle)], capture_output=True, text=True,
                             check=True).stdout.split()  # fmt: skip
    oracle: dict[str, tuple[str, str, str, str, list[str]]] = {}
    oracle_group: dict[str, frozenset[str]] = {}
    for f in query(a.bcftools, a.fasta, a.oracle):
        ids = f[0].replace(";", ",").split(",")
        for pid in ids:
            oracle[pid] = (f[1], f[2], f[3], f[4], f[5:])
            oracle_group[pid] = frozenset(ids)
    snp: dict[str, str] = {}
    if a.bpm:
        sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
        from simurg.formats.bpm import read_bpm

        m = read_bpm(a.bpm)
        snp = dict(zip(m.names, m.snps, strict=True))
    first = a.ours / f"{samples[0]}.snv.vcf.gz"
    ours_all = {pid for f in query(a.bcftools, a.fasta, first) for pid in f[0].replace(";", ",").split(",")}
    cats: Counter[str] = Counter()
    examples: dict[str, list[Any]] = defaultdict(list)
    ours_ids: set[str] = set()
    probes_in_merged = 0
    for si, sample in enumerate(samples):
        vcf = a.ours / f"{sample}.snv.vcf.gz"
        if not vcf.exists():
            cats["sample_missing"] += 1
            continue
        for f in query(a.bcftools, a.fasta, vcf):
            ids = f[0].replace(";", ",").split(",")
            merged = len(ids) > 1
            ours_al = alleles_of(f[5], f[3], f[4])
            haploid = "/" not in f[5] and f[5] != "."
            group = frozenset(ids)
            mixed_indel = {snp.get(x) for x in ids} >= {"[I/D]", "[D/I]"}
            for pid in ids:
                if si == 0:
                    ours_ids.add(pid)
                    probes_in_merged += merged
                o = oracle.get(pid)
                if o is None:
                    cats["only_ours"] += 1
                    continue
                if (o[0], o[1], o[2]) != (f[1], f[2], f[3]):
                    cats["representation_diff"] += 1
                    if len(examples["representation_diff"]) < 20:
                        examples["representation_diff"].append([pid, f[1:5], list(o[:4])])
                    continue
                theirs = alleles_of(o[4][si], o[2], o[3])
                if ours_al == theirs:
                    cats["match_called" if ours_al else "match_missing"] += 1
                    continue
                og = oracle_group[pid]
                if og - ours_all:
                    key = "oracle_site_includes_probes_we_skip"  # e.g. intensity-only CNV probes (OUT-019)
                elif og != group and (merged or len(og) > 1):
                    key = "grouping_differs"  # same normalised variant, different manifest positions (OUT-030)
                elif ours_al is not None and theirs is None and mixed_indel:
                    key = "oracle_opposite_strand_indel_bug"  # EDGE-012: [I/D] + [D/I] both het → 0/1
                elif ours_al is None:
                    if theirs is not None and theirs[0] != theirs[1] and f[5] == ".":
                        key = "ours_missing_het_in_haploid"  # OUT-018
                    elif f[1] in ("Y", "chrY") and f[5] == ".":
                        key = "ours_missing_haploid_Y"  # OUT-018: Y of female samples
                    elif merged:
                        key = "ours_missing_merged_site"
                    else:
                        key = "ours_missing_UNEXPLAINED"
                elif theirs is None:
                    if merged and haploid:
                        key = "oracle_diploid_candidates_in_haploid_context"  # OUT-031 enumerates with the ploidy
                    else:
                        key = "oracle_missing_ours_called_merged" if merged else "oracle_missing_ours_called"
                else:
                    key = "diff_merged_site" if merged else "diff_UNEXPLAINED"
                cats[key] += 1
                if len(examples[key]) < 20:
                    examples[key].append([sample, pid, f[1:6], o[4][si], list(o[:4])])
    only_oracle = sorted(set(oracle) - ours_ids)
    summary = {
        "oracle": a.name,
        "samples": len(samples),
        "probes_oracle": len(oracle),
        "probes_ours": len(ours_ids),
        "probes_in_merged_records": probes_in_merged,
        "probes_only_in_oracle": len(only_oracle),
        "only_in_oracle_examples": only_oracle[:20],
        "categories": dict(sorted(cats.items())),
        "examples": examples,
    }
    unexplained = sum(v for k, v in cats.items() if "UNEXPLAINED" in k)
    summary["unexplained"] = unexplained
    a.out.parent.mkdir(parents=True, exist_ok=True)
    a.out.write_text(json.dumps(summary, indent=2) + "\n")
    print(json.dumps({k: summary[k] for k in ("probes_oracle", "probes_ours", "probes_only_in_oracle", "categories",
                                                "unexplained")}))  # fmt: skip
    return 0 if unexplained == 0 else 1


if __name__ == "__main__":
    sys.exit(main())