P6: gtc-to-vcf (VCF/BCF + indexes), harmonisation, exports (plink, bedgraph, finalreport)

- formats/fasta: plain/BGZF FASTA, .fai/.gzi built identical to htslib, index check, build detection, contig aliases, mmap genome cache (IN-060..064) - formats/bgzf, formats/bcf: pure-Python BGZF, tabix .tbi, BCF 2.2 and CSI (OUT-010, OUT-038, OUT-040) - harmonize/sites: plus-strand mapping with audit rows, indel resolution from SourceSeq (port of GTCtoVCF, Apache-2.0) and left normalisation (OUT-014..021, OUT-030, OUT-042) - harmonize/combine: OUT-031 genotype combination (ploidy-aware, Infinium I/II signal classes, OUT-037) - pipeline/vcf: `simurg gtc-to-vcf` (per-sample VCF.gz+tbi, --merge cohort BCF+csi, audit parquet) - pipeline/export: `simurg export bedgraph|finalreport|plink|locus-summary` (OUT-043, OUT-050..052, OUT-071) - validation: 36 GSA samples vs gtc2vcf: 23.66 M identical calls, 0 unexplained differences; PLINK identical to plink2 on 637,397 variants; bcftools validates VCF/BCF and our indexes answer region queries correctly - fetch_data.sh p6 (Ensembl GRCh38 with checksum), make_oracle_vcf.sh, validate_p6_vcf.py; K-022/K-023 Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
2026-09-28 12:17:53 +01:00
1 parent 60a05d5048
commit b56fdebaaf
@@ -0,0 +1,121 @@
+"""Genotypes of VCF sites from probe calls (spec OUT-018, OUT-031, OUT-032, OUT-037; plan §8.3).
+
+Single-probe sites map the A/B call straight to allele indices. Sites with several probes use the exact candidate
+algorithm of OUT-031: enumerate every genotype over {REF} ∪ ALT with the site's ploidy, predict each probe's call
+from the alleles' signal classes (Infinium II: dye group of the extended base; Infinium I and indels: exact allele,
+no signal otherwise) and keep the candidates consistent with every called probe.
+[OSS: Illumina GTCtoVCF ``GenotypeFormat.RecordCombiner`` (Apache-2.0) — the consistency idea; the class model and
+the ploidy-aware enumeration are ours.]
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from itertools import combinations_with_replacement
+
+from simurg.harmonize.sites import GREEN, RED, ProbeSite, Site
+
+NC, AA, AB, BB = 0, 1, 2, 3
+CALL_TEXT = {NC: "NC", AA: "AA", AB: "AB", BB: "BB"}
+
+
+def signal_class(p: ProbeSite, allele: str) -> str | None:
+    """OUT-031 step 1."""
+    a, b = p.alleles
+    if p.kind == "snv" and p.infinium == 2 and len(allele) == 1:
+        grp = RED if allele in RED else GREEN
+        if a in grp:
+            return "A"
+        if b in grp:
+            return "B"
+        return None
+    if allele == a:
+        return "A"
+    if allele == b:
+        return "B"
+    return None
+
+
+def predicted_call(p: ProbeSite, genotype: tuple[str, ...]) -> int:
+    classes = sorted(c for c in (signal_class(p, z) for z in genotype) if c is not None)
+    if not classes:
+        return NC
+    u = set(classes)
+    if u == {"A"}:
+        return AA
+    if u == {"B"}:
+        return BB
+    return AB
+
+
+@dataclass(frozen=True)
+class Outcome:
+    alleles: tuple[int, ...] | None  # allele indices (sorted) or None for missing
+    reason: str  # "" | "het_in_haploid" | W_MERGE_* | W_INF1_NC_CONFLICT
+
+
+def single(site: Site, call: int, ploidy: int) -> Outcome:
+    p = site.probes[0]
+    if call == NC:
+        return Outcome(None, "")
+    ia, ib = site.allele_index(p.alleles[0]), site.allele_index(p.alleles[1])
+    if ploidy == 1:
+        if call == AB:
+            return Outcome(None, "het_in_haploid")
+        return Outcome((ia if call == AA else ib,), "")
+    pair = {AA: (ia, ia), AB: (ia, ib), BB: (ib, ib)}[call]
+    return Outcome(tuple(sorted(pair)), "")
+
+
+def merged(site: Site, calls: tuple[int, ...], ploidy: int, inf1_nc_informative: bool = False) -> Outcome:
+    """OUT-031 steps 3–6."""
+    if all(c == NC for c in calls):
+        return Outcome(None, "W_MERGE_ALL_NC")
+    alleles = [site.ref, *site.alts]
+    cands = list(combinations_with_replacement(range(len(alleles)), ploidy))
+    keep = [
+        g
+        for g in cands
+        if all(
+            c == NC or predicted_call(p, tuple(alleles[i] for i in g)) == c
+            for p, c in zip(site.probes, calls, strict=True)
+        )
+    ]
+    reason = ""
+    if inf1_nc_informative:
+        for p, c in zip(site.probes, calls, strict=True):
+            if c == NC and p.infinium == 1:
+                sub = [g for g in keep if all(signal_class(p, alleles[i]) is None for i in g)]
+                if sub:
+                    keep = sub
+                else:
+                    reason = "W_INF1_NC_CONFLICT"
+    if len(keep) == 1:
+        return Outcome(keep[0], reason)
+    return Outcome(None, "W_MERGE_INCONSISTENT" if not keep else "W_MERGE_AMBIGUOUS")
+
+
+def genotype(site: Site, calls: tuple[int, ...], ploidy: int, inf1_nc_informative: bool = False) -> Outcome:
+    if len(site.probes) == 1:
+        return single(site, calls[0], ploidy)
+    key = (calls, ploidy, inf1_nc_informative)
+    hit = site.memo.get(key)
+    if hit is None:
+        hit = site.memo[key] = merged(site, calls, ploidy, inf1_nc_informative)
+    assert isinstance(hit, Outcome)
+    return hit
+
+
+def gt_text(o: Outcome, ploidy: int) -> str:
+    if o.alleles is None:
+        return "." if ploidy == 1 else "./."
+    return "/".join(str(i) for i in o.alleles)
+
+
+def failure_message(site: Site, calls: tuple[int, ...], reason: str) -> str:
+    what = {"W_MERGE_INCONSISTENT": "inconsistent", "W_MERGE_AMBIGUOUS": "ambiguous", "W_MERGE_ALL_NC": "all no-calls"}
+    probes = ", ".join(
+        f"{p.name} ({'InfiniumI' if p.infinium == 1 else 'InfiniumII'}): {CALL_TEXT[c]}"
+        for p, c in zip(site.probes, calls, strict=True)
+    )
+    return f"Failed to combine genotypes ({what.get(reason, reason)}) - {probes}"