SimurgArray/tests/oracle/test_p6_vcf.py

"""Phase P6 acceptance on P0 data (plan §18 P6): VCFs validate with bcftools and index correctly; genotypes agree
with gtc2vcf apart from documented categories; BAF is ALT-oriented (EDGE-017); bedGraph equals GTC values (EDGE-023).

Needs the GRCh38 FASTA (``scripts/fetch_data.sh p6``) and, for the comparison, the gtc2vcf oracle VCF
(``scripts/make_oracle_vcf.sh``).
"""

from __future__ import annotations

import gzip
import subprocess
from pathlib import Path

import numpy as np
import pytest
from tests.conftest import GSA_BPM, GSA_CSV, ORACLE_GTCS, REFS, ROOT, need

pytestmark = [pytest.mark.data, pytest.mark.slow]
FASTA = ROOT / "data" / "p6" / "GRCh38.primary.fa"
BCFTOOLS = REFS / "tools" / "bcftools-1.22" / "bcftools"
ORACLE_VCF = ROOT / "data" / "p6" / "oracle" / "gsa36_gtc2vcf.vcf.gz"


@pytest.fixture(scope="module")
def vcfs(tmp_path_factory: pytest.TempPathFactory) -> Path:
    need(GSA_BPM, GSA_CSV, FASTA, ORACLE_GTCS)
    from typer.testing import CliRunner

    from simurg.cli.main import app

    gtcs = tmp_path_factory.mktemp("gtc")
    for f in sorted(ORACLE_GTCS.glob("*.gtc"))[:3]:
        (gtcs / f.name).symlink_to(f)
    out = tmp_path_factory.mktemp("vcf")
    res = CliRunner().invoke(
        app,
        ["gtc-to-vcf", "--bpm", str(GSA_BPM), "--csv", str(GSA_CSV), "--fasta", str(FASTA), "--gtcs", str(gtcs),
         "--out", str(out), "--threads", "1", "--unsquash-duplicates", "--quiet", "--force"],
    )  # fmt: skip
    assert res.exit_code == 0, res.output
    return out


def test_bcftools_validates_and_index_works(vcfs: Path) -> None:  # OUT-038, OUT-039, EDGE-020, EDGE-021
    if not BCFTOOLS.exists():
        pytest.skip("bcftools not built (scripts/make_oracle_gtcs.sh)")
    f = next(vcfs.glob("*.snv.vcf.gz"))
    n = int(subprocess.run([str(BCFTOOLS), "view", "-H", str(f)], capture_output=True, check=True).stdout.count(b"\n"))
    nrec = subprocess.run([str(BCFTOOLS), "index", "--nrecords", str(f)], capture_output=True, text=True, check=True)
    assert int(nrec.stdout.strip()) == n > 600_000
    for region in ("1:1000000-2000000", "X:150000000-156000000", "MT", "Y"):
        idx = subprocess.run([str(BCFTOOLS), "view", "-H", "-r", region, str(f)], capture_output=True, check=True)
        chrom, _, span = region.partition(":")
        lo, _, hi = span.partition("-")
        scan = 0
        with gzip.open(f, "rt") as fh:
            for line in fh:
                if line.startswith("#"):
                    continue
                c, p = line.split("\t", 2)[:2]
                if c == chrom and (not lo or int(lo) <= int(p) <= int(hi)):
                    scan += 1
        assert idx.stdout.count(b"\n") == scan


def test_matches_gtc2vcf(vcfs: Path, tmp_path: Path) -> None:  # plan §18 P6 gate
    if not ORACLE_VCF.exists() or not BCFTOOLS.exists():
        pytest.skip("gtc2vcf oracle VCF missing (scripts/make_oracle_vcf.sh)")
    names = ",".join(p.name.removesuffix(".snv.vcf.gz") for p in sorted(vcfs.glob("*.snv.vcf.gz")))
    sub = tmp_path / "oracle.vcf.gz"
    subprocess.run([str(BCFTOOLS), "view", "-s", names, "-Oz", "-o", str(sub), str(ORACLE_VCF)], check=True)
    out = tmp_path / "cmp.json"
    r = subprocess.run(
        [str(ROOT / ".venv" / "bin" / "python"), str(ROOT / "scripts" / "validate_p6_vcf.py"), "--ours", str(vcfs),
         "--oracle", str(sub), "--fasta", str(FASTA), "--bcftools", str(BCFTOOLS), "--out", str(out)],
        capture_output=True, text=True,
    )  # fmt: skip
    assert r.returncode == 0, r.stdout + r.stderr
    import json

    d = json.loads(out.read_text())
    assert d["unexplained"] == 0 and d["categories"]["match_called"] > 1_500_000


def test_baf_is_alt_oriented(vcfs: Path) -> None:  # OUT-016, EDGE-017
    from simurg.formats.bpm import read_bpm
    from simurg.formats.gtc import read_gtc

    f = sorted(vcfs.glob("*.snv.vcf.gz"))[0]
    g = read_gtc(ORACLE_GTCS / f.name.replace(".snv.vcf.gz", ".gtc"))
    m = read_bpm(GSA_BPM)
    idx = m.name_index()
    flipped = same = 0
    with gzip.open(f, "rt") as fh:
        for line in fh:
            if line.startswith("#"):
                continue
            c = line.rstrip("\n").split("\t")
            if "," in c[2] or ";" in c[2] or c[8] != "GT:GS:BAF:LRR" or "," in c[4]:
                continue
            baf = c[9].split(":")[2]
            if baf == ".":
                continue
            i = idx[c[2]]
            gtc_baf = float(g.b_allele_freqs[i])
            gt = c[9].split(":")[0]
            call = int(g.genotypes[i])
            if gt in ("0/0", "1/1") and call in (1, 3):
                a_is_ref = (gt == "0/0") == (call == 1)
                want = gtc_baf if a_is_ref else float(np.float32(1.0) - np.float32(gtc_baf))
                assert abs(float(baf) - want) < 1e-6, (c[2], baf, gtc_baf)
                flipped += not a_is_ref
                same += a_is_ref
            if flipped > 500 and same > 500:
                break
    assert flipped > 500 and same > 500  # both orientations exercised


def test_bedgraph_equals_gtc(tmp_path: Path) -> None:  # OUT-050, OUT-052, EDGE-023
    need(GSA_BPM, ORACLE_GTCS)
    from simurg.core.log import Reporter
    from simurg.formats.bpm import read_bpm
    from simurg.formats.gtc import read_gtc
    from simurg.pipeline.export import BedgraphOptions, run_export_bedgraph

    gtcs = tmp_path / "g"
    gtcs.mkdir()
    src = sorted(ORACLE_GTCS.glob("*.gtc"))[0]
    (gtcs / src.name).symlink_to(src)
    out = tmp_path / "bg"
    assert run_export_bedgraph(BedgraphOptions(out=out, bpm=GSA_BPM, gtcs=gtcs), Reporter(command="t", quiet=True)) == 0
    m = read_bpm(GSA_BPM)
    g = read_gtc(src)
    single = {}
    seen: dict[tuple[str, int], int] = {}
    for i, (c, p) in enumerate(zip(m.chroms, m.map_infos.tolist(), strict=True)):
        key = ("X" if c == "XY" else c, p)
        seen[key] = seen.get(key, 0) + 1
        single[key] = i
    lines = gzip.decompress((out / f"{src.stem}.lrr.bedgraph.gz").read_bytes()).decode().splitlines()
    assert lines[0].startswith('track type=bedGraph name="') and "viewLimits=-2:2" in lines[0]
    checked = 0
    for ln in lines[1:]:
        c, s, e, v = ln.split("\t")
        key = (c, int(e))
        assert int(e) - int(s) == 1
        if seen[key] == 1:
            assert v == f"{float(g.logr_ratios[single[key]]):.6g}"
            checked += 1
    assert checked > 500_000