P6 complete: optional VCF fields, SNV+indel merging, auxiliary MNV loci, PLINK 2 .pgen

- gtc-to-vcf --info-fields ASSAY_TYPE,NORM_BIN,GENTRAIN,CLUSTER_SEP,ILMN_STRAND,PROBE_IDS (--egt) and
  --extra-format X,Y,NORMX,NORMY,THETA,R, in per-sample VCFs and the --merge BCF (OUT-012, plan §8.3);
  equal to gtc2vcf's values (THETA/R missing where both normalized intensities are 0, K-023)
- --merge-snv-indel (OUT-021, EDGE-013): alleles extended to the longest REF; SNV probes read the base at POS
- --auxiliary-loci (OUT-035, EDGE-015): rules (a)-(e), E_AUX_INVALID; rule (d) uses the last plus-strand base
  for minus-strand designs (spec amended, as GTCtoVCF)
- export plink --format bed|pgen|both: PLINK 2 fixed-width .pgen/.pvar/.psam; plink2 reads it back to a .bed
  byte-identical to ours
- cohort store keeps raw/normalized intensities for the BCF path; tests (unit + gtc2vcf oracle); docs

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv

src/simurg/pipeline/vcf.py

@@ -57,6 +57,11 @@ class VcfOptions:
     ploidy_policy: str = "sex-aware"
     threads: int = 0
     merge: bool = False
+    merge_snv_indel: bool = False
+    auxiliary_loci: Path | None = None
+    info_fields: list[str] = field(default_factory=list)
+    extra_format: list[str] = field(default_factory=list)
+    egt: Path | None = None
     force: bool = False
     argv: list[str] = field(default_factory=list)
 
@@ -82,11 +87,33 @@ class SiteTable:
     site_id: list[str] = field(default_factory=list)
     site_alleles: list[list[str]] = field(default_factory=list)
     single_ab: npt.NDArray[np.int16] | None = None  # (n_single, 2) allele indices of the probe's A and B alleles
+    site_info: list[str] = field(default_factory=list)  # INFO text per site ('.' without --info-fields)
+    site_info_typed: list[list[tuple[str, str, list[Any]]]] = field(default_factory=list)  # (key, kind, values)
+    extra_format: list[str] = field(default_factory=list)
+    single_lookup: npt.NDArray[np.int32] | None = None  # normalization-transform index of single-probe loci
     inf1_nc_informative: bool = False
     ploidy_policy: str = "sex-aware"
     filter_gs: float | None = None
 
 
+INFO_DEFS = {  # OUT-012, one value per contributing probe (BPM order)
+    "ASSAY_TYPE": ("Integer", "Assay type per probe: 0 Infinium II, 1 Infinium I red, 2 Infinium I green"),
+    "NORM_BIN": ("Integer", "Normalization bin per probe"),
+    "GENTRAIN": ("Float", "GenTrain score per probe (cluster file)"),
+    "CLUSTER_SEP": ("Float", "Cluster separation score per probe (cluster file)"),
+    "ILMN_STRAND": ("String", "Design strand per probe (TOP, BOT, PLUS, MINUS)"),
+    "PROBE_IDS": ("String", "Names of the contributing probes"),
+}
+EXTRA_DEFS = {  # plan §8.3; single-probe records only
+    "X": ("Integer", "Raw intensity, X (red) channel"),
+    "Y": ("Integer", "Raw intensity, Y (green) channel"),
+    "NORMX": ("Float", "Normalized X intensity"),
+    "NORMY": ("Float", "Normalized Y intensity"),
+    "THETA": ("Float", "Normalized theta = (2/pi) atan2(NORMY, NORMX)"),
+    "R": ("Float", "Normalized R = NORMX + NORMY"),
+}
+
+
 def _site_class(contig: str, pos: int, par_flag: bool, build: str) -> int:
     c = contig[3:] if contig.lower().startswith("chr") else contig
     if c in ("MT", "M"):
@@ -196,8 +223,21 @@ def build_sites(
             probes.append(p)
     for code, count in sorted(skips.items()):
         rep.warn(code, f"{count} locus/loci: {code}")
+    if opts.auxiliary_loci is not None:
+        from simurg.harmonize.sites import apply_auxiliary, read_auxiliary
+
+        aux = read_auxiliary(str(opts.auxiliary_loci))
+        strands = mv.ref_strands
+        probes, errors = apply_auxiliary(probes, aux, genome.fetch, cmap.get, lambda i: (strands[i] or "").strip())
+        if errors:
+            raise SimurgError("E_AUX_INVALID", f"{len(errors)} invalid auxiliary record(s):\n  " + "\n  ".join(errors))
+        rep.info(f"{len(aux)} auxiliary locus definition(s) applied")
+        names = {r.name for r in aux}
+        for a_row in audit:
+            if a_row.locus in names:
+                a_row.reason = "aux_locus"
     order = {e.name: k for k, e in enumerate(genome.entries)}
-    sites = group_sites(probes, order, unsquash=opts.unsquash_duplicates)
+    sites = group_sites(probes, order, unsquash=opts.unsquash_duplicates, merge_snv_indel=opts.merge_snv_indel)
     used = {s.contig for s in sites}
     contigs = [e.name for e in genome.entries if e.name in used]
     # static per-site text and single-probe genotype tables
@@ -222,6 +262,7 @@ def build_sites(
         else:
             merged_sites.append(s)
             merged_rows.append(k)
+    site_info, site_info_typed = _site_info(opts, bpm, sites)
     header = vcf_header(opts, genome, build, contigs, rm)
     table = SiteTable(
         prefix=prefix,
@@ -240,6 +281,10 @@ def build_sites(
         site_id=[s.id for s in sites],
         site_alleles=[[s.ref, *s.alts] for s in sites],
         single_ab=np.array(single_ab, np.int16).reshape(-1, 2),
+        site_info=site_info,
+        site_info_typed=site_info_typed,
+        extra_format=list(opts.extra_format),
+        single_lookup=bpm.normalization_lookups()[np.array(single_locus, np.int64)] if opts.extra_format else None,
         inf1_nc_informative=opts.inf1_nocall_informative,
         ploidy_policy=opts.ploidy_policy,
         filter_gs=opts.filter_gs,
@@ -251,6 +296,49 @@ def _sha(path: Path | None) -> str:
     return sha256_file(path) if path else ""
 
 
+def _site_info(
+    opts: VcfOptions, bpm: Any, sites: list[Site]
+) -> tuple[list[str], list[list[tuple[str, str, list[Any]]]]]:
+    """Static INFO per site for --info-fields (OUT-012)."""
+    if not opts.info_fields:
+        return ["."] * len(sites), [[] for _ in sites]
+    gentrain = sep = None
+    if {"GENTRAIN", "CLUSTER_SEP"} & set(opts.info_fields):
+        from simurg.formats.egt import read_egt
+
+        assert opts.egt is not None
+        egt = read_egt(opts.egt)
+        eidx = egt.name_index()
+        rows = np.array([eidx.get(n, -1) for n in bpm.names], np.int64)
+        gentrain = np.where(rows >= 0, egt.gentrain_score[np.maximum(rows, 0)], np.nan)
+        sep = np.where(rows >= 0, egt.cluster_sep[np.maximum(rows, 0)], np.nan)
+    bins = bpm.norm_bin_ids()
+    texts, typed = [], []
+    for s in sites:
+        idx = [p.index for p in s.probes]
+        vals: dict[str, list[Any]] = {
+            "ASSAY_TYPE": [int(bpm.assay_types[i]) for i in idx],
+            "NORM_BIN": [int(bins[i]) for i in idx],
+            "GENTRAIN": [float(gentrain[i]) for i in idx] if gentrain is not None else [],
+            "CLUSTER_SEP": [float(sep[i]) for i in idx] if sep is not None else [],
+            "ILMN_STRAND": [bpm.ilmn_strands[i] or "." for i in idx],
+            "PROBE_IDS": [p.name for p in s.probes],
+        }
+        parts, tv = [], []
+        for key in opts.info_fields:
+            v = vals[key]
+            kind = INFO_DEFS[key][0]
+            text = ",".join(
+                "." if (isinstance(x, float) and math.isnan(x)) else (f"{x:.6g}" if isinstance(x, float) else str(x))
+                for x in v
+            )
+            parts.append(f"{key}={text}")
+            tv.append((key, kind, v))
+        texts.append(";".join(parts))
+        typed.append(tv)
+    return texts, typed
+
+
 def vcf_header(opts: VcfOptions, genome: Any, build: str, contigs: list[str], rm: RunManifest) -> str:
     """OUT-011 (without the ##SAMPLE and #CHROM lines, which are per sample)."""
 
@@ -271,14 +359,18 @@ def vcf_header(opts: VcfOptions, genome: Any, build: str, contigs: list[str], rm
         "##simurg_options="
         f"unsquash_duplicates={str(opts.unsquash_duplicates).lower()};exclude_loci={lst(opts.exclude_loci)};"
         f"include_loci={lst(opts.include_loci)};inf1_nocall_informative={str(opts.inf1_nocall_informative).lower()};"
-        f"auxiliary_loci=none;ploidy_policy={opts.ploidy_policy};contig_alias={opts.contig_alias}"
+        f"auxiliary_loci={lst(opts.auxiliary_loci)};ploidy_policy={opts.ploidy_policy};contig_alias={opts.contig_alias}"
         + (";drop_ambiguous=true" if opts.drop_ambiguous else "")
+        + (";merge_snv_indel=true" if opts.merge_snv_indel else "")
         + (";unknown_sex_policy=diploid_x" if opts.ploidy_policy == "sex-aware" else "")
     )
     lines.append(f"##simurg_run={rm.run_id}")
     for c in contigs:
         lines.append(f"##contig=<ID={c},length={genome.contig_length(c)}>")
     lines.append('##FILTER=<ID=PASS,Description="All filters passed">')
+    for key in opts.info_fields:
+        kind, desc = INFO_DEFS[key]
+        lines.append(f'##INFO=<ID={key},Number=.,Type={kind},Description="{desc}">')
     if opts.filter_gs is not None:
         lines.append(f'##FILTER=<ID=LowGS,Description="GenCall score below {opts.filter_gs}">')
     lines += [
@@ -289,6 +381,9 @@ def vcf_header(opts: VcfOptions, genome: Any, build: str, contigs: list[str], rm
         'ALT (REF=A, ALT=B relative to the reference genome)">',
         '##FORMAT=<ID=LRR,Number=1,Type=Float,Description="Log R ratio">',
     ]
+    for key in opts.extra_format:
+        kind, desc = EXTRA_DEFS[key]
+        lines.append(f'##FORMAT=<ID={key},Number=1,Type={kind},Description="{desc}">')
     return "\n".join(lines) + "\n"
 
 
@@ -338,6 +433,22 @@ def write_sample(task: tuple[str, str, str, str, bool]) -> dict[str, Any]:
             b = np.where(t.single_flip, (np.float32(1.0) - b).astype(np.float32), b)
             baf[rows] = b.astype(np.float64)
             lrr[rows] = g.logr_ratios[loc].astype(np.float64)
+        extra_cols: list[list[str]] = []
+        if t.extra_format:
+            from simurg.signal.normalization.apply import normalize, polar
+
+            rx, ry = g.raw_x[loc], g.raw_y[loc]
+            nx, ny = normalize(rx, ry, t.single_lookup, g.normalization_transforms)  # type: ignore[arg-type]
+            zero = (rx == 0) & (ry == 0)  # no intensity: normalized values are missing (as in the GTC semantics)
+            nx = np.where(zero, np.float32(np.nan), nx).astype(np.float32)
+            ny = np.where(zero, np.float32(np.nan), ny).astype(np.float32)
+            th, rr = polar(nx, ny)
+            src: dict[str, npt.NDArray[Any]] = {"X": rx, "Y": ry, "NORMX": nx, "NORMY": ny, "THETA": th, "R": rr}
+            for key in t.extra_format:
+                col = np.full(n, "", dtype=object)
+                vals = src[key]
+                col[rows] = [str(int(v)) for v in vals] if key in ("X", "Y") else [_g8(float(v)) for v in vals]
+                extra_cols.append(col.tolist())
         het_haploid = int(((calls == 2) & (pl == 1)).sum())
         if het_haploid:
             res["reasons"]["het_in_haploid"] += het_haploid
@@ -373,10 +484,16 @@ def write_sample(task: tuple[str, str, str, str, bool]) -> dict[str, Any]:
                 raw.write(header.encode())  # type: ignore[arg-type]
             for k in range(n):
                 filt = "LowGS" if low is not None and not math.isnan(gsl[k]) and gsl[k] < low else "PASS"
+                info = t.site_info[k] if t.site_info else "."
                 if merged_mask[k]:
-                    body = f"{filt}\t.\tGT:GS\t{gt[k]}:{_f4(gsl[k])}\n"
+                    body = f"{filt}\t{info}\tGT:GS\t{gt[k]}:{_f4(gsl[k])}\n"
                 else:
-                    body = f"{filt}\t.\tGT:GS:BAF:LRR\t{gt[k]}:{_f4(gsl[k])}:{_g8(bafl[k])}:{_g8(lrrl[k])}\n"
+                    extra_keys = "".join(":" + key for key in t.extra_format)
+                    extra_vals = "".join(":" + c[k] for c in extra_cols)
+                    body = (
+                        f"{filt}\t{info}\tGT:GS:BAF:LRR{extra_keys}\t"
+                        f"{gt[k]}:{_f4(gsl[k])}:{_g8(bafl[k])}:{_g8(lrrl[k])}{extra_vals}\n"
+                    )
                 line = t.prefix[k] + body.encode()
                 if bgzip:
                     vs, ve = w.write_record(line)
@@ -413,10 +530,20 @@ def run_gtc_to_vcf(opts: VcfOptions, rep: Reporter) -> int:
             (opts.sheet, "sample sheet"),
             (opts.exclude_loci, "locus list"),
             (opts.include_loci, "locus list"),
+            (opts.auxiliary_loci, "auxiliary loci VCF"),
         ]
     )
     if opts.ploidy_policy not in ("sex-aware", "diploid"):
         raise SimurgError("E_CONFIG_KEY", "--ploidy-policy must be sex-aware or diploid")
+    bad = [x for x in opts.info_fields if x not in INFO_DEFS] + [x for x in opts.extra_format if x not in EXTRA_DEFS]
+    if bad:
+        raise SimurgError(
+            "E_CONFIG_KEY",
+            f"unknown field(s) {', '.join(bad)}; --info-fields takes {','.join(INFO_DEFS)}; "
+            f"--extra-format takes {','.join(EXTRA_DEFS)}",
+        )
+    if {"GENTRAIN", "CLUSTER_SEP"} & set(opts.info_fields) and opts.egt is None:
+        raise SimurgError("E_CONFIG_KEY", "--info-fields GENTRAIN/CLUSTER_SEP need the cluster file (--egt)")
     sheet = parse_samplesheet(opts.sheet) if opts.sheet else None
     root = opts.gtcs or (sheet.root_folder if sheet else None)
     if root is None:
@@ -534,25 +661,31 @@ def write_cohort_bcf(
     from simurg.formats.bcf import (
         FLOAT,
         INT8,
+        INT32,
         BcfWriter,
         CsiBuilder,
         _type_byte,
         encode_record,
         header_dictionaries,
+        typed_floats,
+        typed_ints,
+        typed_string,
     )
     from simurg.formats.bpm import read_bpm
     from simurg.pipeline.qc import default_tmp, load_cohort
+    from simurg.signal.normalization.apply import polar
 
-    n_loci = read_bpm(opts.bpm).num_loci
+    bpm = read_bpm(opts.bpm)
     workers = max(1, opts.threads or (os.cpu_count() or 1))
     store, infos = load_cohort(
         [(sid, path) for sid, _, path in samples],
-        n_loci,
+        bpm.num_loci,
         default_tmp(out.parent),
         baf=True,
         scores=True,
         lrr=True,
         workers=workers,
+        lookup=bpm.normalization_lookups() if t.extra_format else None,
     )
     try:
         for s in infos:
@@ -591,6 +724,9 @@ def write_cohort_bcf(
                 sc = np.asarray(store.score[:, sel]) if len(sel) else np.zeros((n, 0), np.float32)  # type: ignore[index]
                 bf = np.asarray(store.baf[:, sel]) if len(sel) else np.zeros((n, 0), np.float32)  # type: ignore[index]
                 lr = np.asarray(store.lrr[:, sel]) if len(sel) else np.zeros((n, 0), np.float32)  # type: ignore[index]
+                if t.extra_format and len(sel):
+                    ex_rx, ex_ry = np.asarray(store.rx[:, sel]), np.asarray(store.ry[:, sel])  # type: ignore[index]
+                    ex_nx, ex_ny = np.asarray(store.nx[:, sel]), np.asarray(store.ny[:, sel])  # type: ignore[index]
                 for k in rows:
                     pl = ploidy[k]
                     if single_pos[k] >= 0:
@@ -616,6 +752,25 @@ def write_cohort_bcf(
                         b = np.where(t.single_flip[j], (np.float32(1.0) - b).astype(np.float32), b)
                         fmt = [(k_gt, _type_byte(width, INT8), cells.tobytes()), (k_gs, fdesc, _f32_payload(gs)),
                                (k_baf, fdesc, _f32_payload(b)), (k_lrr, fdesc, _f32_payload(lr[:, c]))]  # fmt: skip
+                        if t.extra_format:
+                            zero = (ex_rx[:, c] == 0) & (ex_ry[:, c] == 0)
+                            enx = np.where(zero, np.float32(np.nan), ex_nx[:, c]).astype(np.float32)
+                            eny = np.where(zero, np.float32(np.nan), ex_ny[:, c]).astype(np.float32)
+                            th, rr = polar(enx, eny)
+                            xsrc = {
+                                "X": ex_rx[:, c],
+                                "Y": ex_ry[:, c],
+                                "NORMX": enx,
+                                "NORMY": eny,
+                                "THETA": th,
+                                "R": rr,
+                            }
+                            for key in t.extra_format:
+                                if key in ("X", "Y"):
+                                    payload = xsrc[key].astype("<i4").tobytes()
+                                    fmt.append((strings[key], _type_byte(1, INT32), payload))
+                                else:
+                                    fmt.append((strings[key], fdesc, _f32_payload(xsrc[key].astype(np.float32))))
                         gs_min = gs
                     else:
                         site = t.merged_sites[merged_pos[k]]
@@ -637,9 +792,18 @@ def write_cohort_bcf(
                                     cells[si, q] = (a + 1) << 1
                         fmt = [(k_gt, _type_byte(width, INT8), cells.tobytes()), (k_gs, fdesc, _f32_payload(gs_min))]
                     low = opts.filter_gs is not None and bool(np.nanmin(gs_min) < opts.filter_gs) if n else False
+                    info = []
+                    for key, kind, vals in t.site_info_typed[k] if t.site_info_typed else []:
+                        if kind == "Integer":
+                            enc = typed_ints([int(v) for v in vals])
+                        elif kind == "Float":
+                            enc = typed_floats([float(v) for v in vals])
+                        else:
+                            enc = typed_string(",".join(str(v) for v in vals))
+                        info.append((strings[key], enc))
                     rec = encode_record(contigs[t.site_contig[k]], int(t.site_beg[k]),
                                         int(t.site_end[k] - t.site_beg[k]), t.site_id[k], t.site_alleles[k],
-                                        f_low if low and f_low >= 0 else f_pass, n, fmt)  # fmt: skip
+                                        f_low if low and f_low >= 0 else f_pass, n, fmt, info)  # fmt: skip
                     vs, ve = w.write(rec)
                     csi.add(contigs[t.site_contig[k]], int(t.site_beg[k]), int(t.site_end[k]), vs, ve)
                     written += 1