PGx star alleles: pgx star-allele call (PGX-011..021), probe masks vs 1000G, GeT-RM validation (K-029)

- simurg/pgx/star.py: solver-free diplotype enumeration incl. CN forms (deletion, xN), wildcards, array-equivalence
  classes, raw_score/score with CPIC frequency priors; calibrated error model (spec PGX-012 amended).
- pipeline/pgx_star.py: DB gene models, SNV + PGx CN VCF evidence, masks, diplotype_summary.tsv + pgx_calls.json.
- scripts/pgx: derive_masks.py (array vs 1000G 30x over HTTPS, GS calibration), getrm_truth.py,
  validate_star_getrm.py.
- GSA v4 PGx demo vs GeT-RM (16 CEPH): 100/109 gene calls; CYP2D6/2C19/2C9 41/45 (docs/validation/V2_pgx_star.md).

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

src/simurg/cli/main.py

@@ -838,6 +838,37 @@ def pgx_db_build_cmd(
     run_guarded(rep, lambda: run_db_build(o, rep))
 
 
+pgx_star_app = typer.Typer(help="Star-allele calling (PGX-011).")
+pgx_app.add_typer(pgx_star_app, name="star-allele")
+
+
+@pgx_star_app.command("call")
+def pgx_star_call_cmd(
+    vcfs: Annotated[list[Path], typer.Option("--vcfs", help="Folder(s) with SNV + PGx CN VCFs (repeatable).")],
+    db: Annotated[Path, typer.Option("--db", help="PGx DB from pgx db build (required).")],
+    out: Annotated[Path, typer.Option("--out", help="Output folder (required).")],
+    population: Annotated[str, typer.Option("--population", help="Frequency prior population.")] = "global",
+    max_candidates: Annotated[int, typer.Option("--max-candidates", min=1)] = 10,
+    feasible_ratio: Annotated[float, typer.Option("--feasible-ratio", min=0.0, max=1.0)] = 0.5,
+    low_gs: Annotated[float, typer.Option("--low-gs")] = 0.15,
+    no_cn: Annotated[bool, typer.Option("--no-cn", help="Run without PGx CN VCFs.")] = False,
+    allow_foreign_vcf: Annotated[bool, typer.Option("--allow-foreign-vcf")] = False,
+    gene: Annotated[list[str] | None, typer.Option("--gene", help="Restrict to gene(s).")] = None,
+    force: ForceOpt = False,
+    log_json: LogJson = False,
+    debug: Debug = False,
+    quiet: Quiet = False,
+) -> None:
+    """Call star-allele diplotypes (CLI-033, PGX-012 ... PGX-021)."""
+    from simurg.pipeline.pgx_star import StarOptions, run_star_call
+
+    rep = _reporter("pgx star-allele call", out, log_json, debug, quiet)
+    o = StarOptions(vcfs=vcfs, db=db, out=out, population=population, max_candidates=max_candidates,
+                    feasible_ratio=feasible_ratio, low_gs=low_gs, no_cn=no_cn, allow_foreign_vcf=allow_foreign_vcf,
+                    genes=gene or [], force=force)  # fmt: skip
+    run_guarded(rep, lambda: run_star_call(o, rep))
+
+
 @pgx_cn_app.command("train")
 def pgx_cn_train_cmd(
     gtcs: Annotated[Path, typer.Option("--gtcs", help="Folder searched recursively for training GTCs (required).")],

src/simurg/core/errors.py

@@ -171,6 +171,9 @@ CATALOGUE: dict[str, MessageDef] = dict(
         _m("W_PGX_NO_ALT_CONTIG", "warning", "FASTA lacks a contig carrying a PGx region; REF written as N", "IN-062"),
         _m("W_PGX_BATCH_COMPOSITION", "warning", "Batch spans several lots/batches or > 14 days of scans", "PCN-006"),
         _m("E_PGX_CN_FAILED", "error", "PGx CN calling failed for the sample", "PCN-002"),
+        _m("E_PGX_CN_MISSING", "error", "A sample has no PGx CN VCF (give --no-cn to run without)", "PGX-011"),
+        _m("E_PGX_VCF_UNINDEXED", "error", "SNV VCFs for PGx must be bgzipped and indexed", "PGX-011"),
+        _m("E_PGX_VCF_UNSUITABLE", "error", "SNV VCF unsuitable for PGx (squashing/locus filters/foreign)", "PGX-011"),
         _m("E_PGX_SOURCE", "error", "PGx source tables missing; run scripts/pgx/fetch_sources.sh", "PGX-001"),
         # --- methylation control QC (MET-001 .. MET-018) -------------------------------------------------
         _m("E_METH_NO_ASSAY", "error", "Probe section not found", "MET-002"),

src/simurg/pgx/star.py

@@ -0,0 +1,172 @@
+"""Star-allele / haplotype calling from array genotypes — spec PGX-012 … PGX-020. Pure functions, no I/O.
+
+Evidence model (PGX-012) [OURS]: each defining variant v of a gene is ``called`` (alt dosage g, GenCall score GS,
+BAF), ``no_call``, ``masked`` or ``not_covered``; only called variants contribute. With diploid CN 2 the genotype
+likelihood is P(g | d) = 1 − ε for the expected alt dosage d and ε/2 otherwise, with ε = 1e-3 for GS ≥ 0.15 and
+0.05 below (calibrated: after probe masking, 0 of 9740 GSA v4 PGx calls at DB variants disagreed with 1000G 30x in any
+GS bin — errors are probe-specific and handled by masks, not by GS; K-029). When the gene's copy number n ≠ 2
+(targeted CN), the BAF likelihood N(b; k/n, σ_b) replaces the
+genotype term (k = expected alt copies). Haplotype positions with an IUPAC value (CPIC core definitions) are wildcards:
+the best-fitting dosage is used.
+
+Candidates (PGX-013) [DIFF: solver-free]: every unordered pair of defined haplotypes, and for n ≠ 2 the structural
+forms of CYP2D6-like genes: n = 1 → one whole-gene deletion haplotype; n = 0 → two; n ≥ 3 → one haplotype duplicated
+(``*1x2``, n − 1 copies). Diplotypes that no called variant separates are array-equivalent (PGX-017) and form one
+solution, named after the member with the highest prior. Scores (PGX-014): raw_score = softmax(LL); score =
+softmax(LL + log π), π = Σ over members of f(h1)·f(h2) (×2 for h1 ≠ h2), f floored at 1e-5.
+"""
+
+from __future__ import annotations
+
+import math
+from dataclasses import dataclass, field
+
+import numpy as np
+import numpy.typing as npt
+
+F64 = npt.NDArray[np.float64]
+FREQ_FLOOR = 1e-5
+CALLED, NO_CALL, MASKED, NOT_COVERED = 0, 1, 2, 3
+STATE_NAMES = {CALLED: "called", NO_CALL: "no_call", MASKED: "masked", NOT_COVERED: "not_on_array"}
+
+
+@dataclass
+class GeneModel:
+    gene: str
+    variants: list[tuple[str, int, str, str]]  # (chrom, pos, ref, alt), one per defining variant
+    variant_ids: list[tuple[str, str | None]]  # (hgvs, rsid)
+    haplotypes: list[str]
+    carries: npt.NDArray[np.int8]  # (H, V): 1 carries alt, 0 reference, −1 wildcard
+    freq: F64  # (H,) population frequency of each haplotype (prior)
+    deletion_allele: str | None = None  # named whole-gene deletion (CYP2D6 *5), or "*0"
+    reference: str = ""
+
+
+@dataclass
+class Observations:
+    state: npt.NDArray[np.int8]  # (V,)
+    dosage: npt.NDArray[np.int8]  # alt dosage (0 … ploidy), −1 when not called
+    gs: F64
+    baf: F64
+
+
+@dataclass
+class Solution:
+    haplotypes: tuple[str, ...]  # representative, e.g. ("*1", "*4") or ("*2x2", "*4")
+    ll: float
+    prior: float
+    members: list[tuple[str, ...]] = field(default_factory=list)
+    raw_score: float = 0.0
+    score: float = 0.0
+    copy_number: int = 2
+
+    @property
+    def genotype(self) -> str:
+        return "/".join(self.haplotypes)
+
+
+ERR_CALLED = 1e-3
+ERR_LOW_GS = 0.05
+LOW_GS = 0.15
+
+
+def gs_error(gs: F64) -> F64:
+    """PGX-012: per-call error rate from the GenCall score (calibrated, see the module docstring)."""
+    return np.where(np.nan_to_num(gs, nan=0.0) >= LOW_GS, ERR_CALLED, ERR_LOW_GS)
+
+
+def _dosage_ll(obs: Observations, expected: npt.NDArray[np.int16], wild: npt.NDArray[np.bool_]) -> F64:
+    """(D, V) log-likelihood of each candidate's expected dosage (CN 2 genotype model)."""
+    eps = gs_error(obs.gs)[None, :]
+    g = obs.dosage[None, :].astype(np.int16)
+    hit = (expected == g) | wild
+    return np.where(hit, np.log(1 - eps), np.log(eps / 2))
+
+
+def _baf_ll(obs: Observations, frac_lo: F64, frac_hi: F64, sigma: float) -> F64:
+    """(D, V) BAF log-likelihood for expected alt fractions (a wildcard spans [lo, hi]: the closest value is used)."""
+    b = np.nan_to_num(obs.baf, nan=0.5)[None, :]
+    near = np.clip(b, frac_lo, frac_hi)
+    return np.log(0.98 * np.exp(-0.5 * ((b - near) / sigma) ** 2) / (sigma * math.sqrt(2 * math.pi)) + 0.02)
+
+
+def candidates(m: GeneModel, cn: int | None) -> list[tuple[tuple[int, ...], tuple[str, ...]]]:
+    """(haplotype index per gene copy, display names) for every candidate diplotype (PGX-013)."""
+    h = len(m.haplotypes)
+    out: list[tuple[tuple[int, ...], tuple[str, ...]]] = []
+    if cn is None or cn == 2:
+        for i in range(h):
+            for j in range(i, h):
+                out.append(((i, j), (m.haplotypes[i], m.haplotypes[j])))
+        return out
+    dele = m.deletion_allele or "*0"
+    if cn == 0:
+        return [((), (dele, dele))]
+    if cn == 1:
+        return [((i,), (m.haplotypes[i], dele)) for i in range(h)]
+    extra = cn - 2  # copies beyond two: one haplotype duplicated (e.g. CN 3 → *2x2/*4)
+    for i in range(h):
+        for j in range(h):
+            out.append(((i,) * (extra + 1) + (j,), (f"{m.haplotypes[i]}x{extra + 1}", m.haplotypes[j])))
+    return out
+
+
+def call_gene(
+    m: GeneModel, obs: Observations, cn: int | None = None, baf_sigma: float = 0.06, ploidy: int = 2,
+    max_candidates: int = 10, feasible_ratio: float = 0.5,
+) -> tuple[list[Solution], list[Solution]]:  # fmt: skip
+    """All scored solutions (ranked) and the feasible ones (PGX-014/015)."""
+    called = obs.state == CALLED
+    haploid = [((i,), (m.haplotypes[i],)) for i in range(len(m.haplotypes))]
+    cands = haploid if ploidy == 1 else candidates(m, cn)
+    if not cands:
+        return [], []
+    carries = m.carries.astype(np.int16)
+    lo = np.zeros((len(cands), carries.shape[1]))
+    hi = np.zeros_like(lo)
+    for k, (ix, _) in enumerate(cands):
+        for i in ix:
+            c = carries[i]
+            lo[k] += np.where(c < 0, 0, c)
+            hi[k] += np.where(c < 0, 1, c)
+    ncp = np.array([max(1, len(ix)) for ix, _ in cands], dtype=np.float64)[:, None]
+    if ploidy == 1 or cn is None or cn == 2:
+        wild = hi != lo
+        ll_v = _dosage_ll(obs, lo.astype(np.int16), wild)
+        if wild.any():  # a wildcard may also take the higher dosage
+            ll_v = np.maximum(ll_v, _dosage_ll(obs, hi.astype(np.int16), wild))
+    else:
+        ll_v = _baf_ll(obs, lo / ncp, hi / ncp, baf_sigma)
+    ll = (ll_v * called[None, :]).sum(axis=1)
+    f = np.maximum(m.freq, FREQ_FLOOR)
+    # array-equivalence classes: identical expected dosages at every called variant (PGX-017)
+    keys: dict[bytes, list[int]] = {}
+    for k in range(len(cands)):
+        sig = np.round(np.concatenate([lo[k, called], hi[k, called]]) / ncp[k, 0], 6).tobytes()
+        keys.setdefault(sig + str(len(cands[k][0])).encode(), []).append(k)
+    sols = []
+    for members in keys.values():
+        pri = []
+        for k in members:
+            ix = cands[k][0]
+            if not ix:
+                pri.append(1e-6)
+                continue
+            p = float(np.prod([f[i] for i in ix]))
+            if len(set(ix)) > 1:
+                p *= 2
+            pri.append(p)
+        best = members[int(np.argmax(pri))]
+        sols.append(Solution(cands[best][1], float(ll[best]), float(sum(pri)),
+                             [cands[k][1] for k in members], copy_number=cn if cn is not None else 2))  # fmt: skip
+    lls = np.array([s.ll for s in sols])
+    lp = lls + np.log([s.prior for s in sols])
+    raw = np.exp(lls - lls.max())
+    sc = np.exp(lp - lp.max())
+    raw, sc = raw / raw.sum(), sc / sc.sum()
+    for s, r, v in zip(sols, raw, sc, strict=True):
+        s.raw_score, s.score = float(r), float(v)
+    sols.sort(key=lambda s: (-s.score, -s.raw_score, s.genotype))
+    top = sols[0].score
+    feasible = [s for s in sols if s.score >= feasible_ratio * top]
+    return sols[: max(max_candidates, len(feasible))], feasible

src/simurg/pipeline/pgx_star.py

@@ -0,0 +1,256 @@
+"""``simurg pgx star-allele call`` — star alleles from SNV VCFs, PGx DB and PGx CN VCFs (spec PGX-011 … PGX-021,
+CLI-033). I/O and orchestration; the scoring lives in ``simurg.pgx.star``.
+
+Scope of this version: genes whose result basis is ``star_allele`` in the DB. Named-haplotype (G6PD), variant-based
+and tagging genes, sub-allele collapsing and hybrid CYP2D6 structures follow (K-029).
+"""
+
+from __future__ import annotations
+
+import gzip
+import json
+import math
+import sqlite3
+from collections import defaultdict
+from dataclasses import dataclass, field
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+from simurg import __version__
+from simurg.core.errors import ExitCode, SimurgError
+from simurg.core.log import Reporter
+from simurg.pgx.star import CALLED, MASKED, NO_CALL, NOT_COVERED, GeneModel, Observations, Solution, call_gene
+
+CN_REGION = {"CYP2D6": "CYP2D6.p5"}  # PGX-013: total gene copies from the 5′ region
+DELETION_ALLELE = {"CYP2D6": "*5"}
+
+
+@dataclass
+class StarOptions:
+    vcfs: list[Path]
+    db: Path
+    out: Path
+    population: str = "global"
+    max_candidates: int = 10
+    feasible_ratio: float = 0.5
+    low_gs: float = 0.15
+    no_cn: bool = False
+    allow_foreign_vcf: bool = False
+    genes: list[str] = field(default_factory=list)
+    force: bool = False
+
+
+def load_models(db: Path, population: str, genes: list[str] | None = None) -> dict[str, GeneModel]:
+    con = sqlite3.connect(db)
+    q = "SELECT gene FROM genes WHERE result_basis = 'star_allele' ORDER BY gene"
+    out = {}
+    for (gene,) in con.execute(q).fetchall():
+        if genes and gene not in genes:
+            continue
+        alle = con.execute("SELECT allele, kind FROM alleles WHERE gene = ? AND defined = 1 AND structural = 0 "
+                           "ORDER BY allele", (gene,)).fetchall()  # fmt: skip
+        names = [a for a, _ in alle]
+        ref = next((a for a, k in alle if k == "reference"), names[0] if names else "")
+        rows = con.execute("SELECT allele, chrom, pos, ref, alt, hgvs, rsid, location_id FROM allele_variants "
+                           "WHERE gene = ?", (gene,)).fetchall()  # fmt: skip
+        keys = sorted({(r[1], r[2], r[3], r[4]) for r in rows}, key=lambda k: (k[0], k[1], k[2], k[3]))
+        kidx = {k: i for i, k in enumerate(keys)}
+        ids: list[tuple[str, str | None]] = [("", None)] * len(keys)
+        loc_keys: dict[int, list[int]] = defaultdict(list)
+        carries = np.zeros((len(names), len(keys)), dtype=np.int8)
+        hidx = {n: i for i, n in enumerate(names)}
+        for allele, chrom, pos, r_, a_, hgvs, rsid, lid in rows:
+            k = kidx[(chrom, pos, r_, a_)]
+            ids[k] = (hgvs, rsid)
+            loc_keys[lid].append(k)
+            if allele in hidx:
+                carries[hidx[allele], k] = 1
+        for allele, lid in con.execute("SELECT allele, location_id FROM allele_wildcards WHERE gene = ?", (gene,)):
+            if allele in hidx:
+                for k in loc_keys.get(lid, []):
+                    carries[hidx[allele], k] = -1
+        freq = np.full(len(names), np.nan)
+        fr: dict[str, list[float]] = defaultdict(list)
+        for allele, pop, f in con.execute("SELECT allele, population, freq FROM frequencies WHERE gene = ?", (gene,)):
+            if population in ("global", pop):
+                fr[allele].append(float(f))
+        for i, n in enumerate(names):
+            if fr.get(n):
+                freq[i] = float(np.mean(fr[n]))
+        freq = np.nan_to_num(freq, nan=0.0)
+        out[gene] = GeneModel(gene, keys, ids, names, carries, freq, DELETION_ALLELE.get(gene), ref)
+    con.close()
+    return out
+
+
+def read_vcf_sites(
+    path: Path, wanted: dict[str, set[int]]
+) -> tuple[dict[tuple[str, int], list[dict[str, Any]]], dict[str, str]]:
+    """Records at the wanted positions, and the ``##`` header meta (``##key=value``)."""
+    recs: dict[tuple[str, int], list[dict[str, Any]]] = defaultdict(list)
+    meta: dict[str, str] = {}
+    with gzip.open(path, "rt") as fh:
+        for ln in fh:
+            if ln.startswith("##"):
+                k, _, v = ln[2:].rstrip("\n").partition("=")
+                meta.setdefault(k, v)
+                continue
+            if ln.startswith("#"):
+                continue
+            c, p, _ = ln.split("\t", 2)
+            ps = wanted.get(c)
+            if ps is None or int(p) not in ps:
+                continue
+            f = ln.rstrip("\n").split("\t")
+            fmt = f[8].split(":")
+            val = dict(zip(fmt, f[9].split(":"), strict=False))
+            recs[(c, int(p))].append({"id": f[2], "ref": f[3], "alts": f[4].split(","), "gt": val.get("GT", "./."),
+                                      "gs": val.get("GS"), "baf": val.get("BAF")})  # fmt: skip
+    return recs, meta
+
+
+def observations(
+    m: GeneModel, recs: dict[tuple[str, int], list[dict[str, Any]]], masked: set[str] | None = None
+) -> tuple[Observations, int]:
+    """PGX-012 evidence per defining variant; returns the observations and the ploidy seen in the GTs. A record with
+    any masked probe (PGX-010) makes the variant ``masked``."""
+    v = len(m.variants)
+    state = np.full(v, NOT_COVERED, dtype=np.int8)
+    dosage = np.full(v, -1, dtype=np.int8)
+    gs = np.full(v, np.nan)
+    baf = np.full(v, np.nan)
+    ploidy = 2
+    for k, (c, p, r_, a_) in enumerate(m.variants):
+        for rec in recs.get((c, p), []):
+            if rec["ref"] != r_ or a_ not in rec["alts"]:
+                continue
+            ai = rec["alts"].index(a_) + 1
+            if masked and any(pn in masked for pn in rec["id"].replace(";", ",").split(",")):
+                state[k] = MASKED
+                break
+            alleles = rec["gt"].replace("|", "/").split("/")
+            ploidy = len(alleles)
+            if any(a in (".", "") for a in alleles):
+                state[k] = NO_CALL
+                break
+            state[k] = CALLED
+            dosage[k] = sum(int(a) == ai for a in alleles)
+            gs[k] = float(rec["gs"]) if rec["gs"] not in (None, ".") else math.nan
+            baf[k] = float(rec["baf"]) if rec["baf"] not in (None, ".") else math.nan
+            if len(rec["alts"]) > 1 and not math.isnan(baf[k]):
+                baf[k] = math.nan  # BAF of a merged multi-allelic record is not per-ALT
+            break
+    return Observations(state, dosage, gs, baf), ploidy
+
+
+def read_cn_vcf(path: Path) -> dict[str, int | None]:
+    out: dict[str, int | None] = {}
+    with gzip.open(path, "rt") as fh:
+        for ln in fh:
+            if ln.startswith("#"):
+                continue
+            f = ln.rstrip("\n").split("\t")
+            region = f[7].split("REGION=")[1].split(";")[0]
+            cn = f[9].split(":")[0]
+            out[region] = None if cn == "." or "NoCall" in f[6] else int(cn)
+    return out
+
+
+def _discover(dirs: list[Path]) -> tuple[dict[str, Path], dict[str, Path]]:
+    snv, cn = {}, {}
+    for d in dirs:
+        for p in sorted(d.rglob("*.vcf.gz")):
+            if p.name.endswith(".pgx_cn.vcf.gz"):
+                cn[p.name[: -len(".pgx_cn.vcf.gz")]] = p
+            elif p.name.endswith(".snv.vcf.gz"):
+                snv[p.name[: -len(".snv.vcf.gz")]] = p
+    return snv, cn
+
+
+def run_star_call(o: StarOptions, rep: Reporter) -> int:
+    """``simurg pgx star-allele call`` (CLI-033)."""
+    from simurg.core.atomic import check_outputs_absent
+    from simurg.pipeline.qc import write_tsv
+
+    snv, cnv = _discover(o.vcfs)
+    if not snv:
+        raise SimurgError("E_NO_INPUT", f"no *.snv.vcf.gz under {', '.join(map(str, o.vcfs))}")
+    missing_cn = sorted(set(snv) - set(cnv))
+    if missing_cn and not o.no_cn:
+        raise SimurgError("E_PGX_CN_MISSING", f"{len(missing_cn)} sample(s) have no .pgx_cn.vcf.gz (e.g. "
+                          f"{missing_cn[0]})", hint="run pgx copy-number call, or give --no-cn")  # fmt: skip
+    for vp in snv.values():
+        if not Path(str(vp) + ".tbi").exists() and not Path(str(vp) + ".csi").exists():
+            raise SimurgError("E_PGX_VCF_UNINDEXED", f"{vp} is not indexed")
+    check_outputs_absent([o.out / "diplotype_summary.tsv"], force=o.force)
+    models = load_models(o.db, o.population, o.genes or None)
+    mcon = sqlite3.connect(o.db)
+    masks: dict[str, set[str]] = defaultdict(set)
+    for man, probe in mcon.execute("SELECT manifest, probe_name FROM masks"):
+        masks[man].add(probe)
+    mcon.close()
+    wanted: dict[str, set[int]] = defaultdict(set)
+    for m in models.values():
+        for c, p, _, _ in m.variants:
+            wanted[c].add(p)
+    rows = []
+    calls_json = []
+    for sid, path in sorted(snv.items()):
+        recs, meta = read_vcf_sites(path, wanted)
+        opts = meta.get("simurg_options", "")
+        if not opts and not o.allow_foreign_vcf:
+            raise SimurgError("E_PGX_VCF_UNSUITABLE", f"{path}: not a SimurgArray VCF (no ##simurg_options)",
+                              hint="give --allow-foreign-vcf to use it anyway")  # fmt: skip
+        kv = dict(x.partition("=")[::2] for x in opts.split(";") if x)
+        if (
+            kv.get("unsquash_duplicates") == "true"
+            or kv.get("exclude_loci", "none").lower() != "none"
+            or kv.get("include_loci", "none").lower() != "none"
+        ):  # PGX-011
+            raise SimurgError("E_PGX_VCF_UNSUITABLE", f"{path}: made with unsquashed duplicates or a locus filter")
+        cn = {} if o.no_cn or sid not in cnv else read_cn_vcf(cnv[sid])
+        masked = masks.get(meta.get("simurg_manifest", "").split(";")[0], set())
+        sample_name = sid
+        genes_out = []
+        for gene, m in sorted(models.items()):
+            obs, ploidy = observations(m, recs, masked)
+            region = CN_REGION.get(gene)
+            gcn = cn.get(region) if region else None
+            sols, feasible = call_gene(m, obs, cn=gcn if region else None, ploidy=ploidy,
+                                       max_candidates=o.max_candidates, feasible_ratio=o.feasible_ratio)  # fmt: skip
+            flags = []
+            low = [k for k in range(len(m.variants)) if obs.state[k] == CALLED and obs.gs[k] < o.low_gs]
+            if low:
+                flags.append("review_alternatives")
+            no_calls = [k for k in range(len(m.variants)) if obs.state[k] == NO_CALL]
+            genotype = ";".join(s.genotype for s in feasible)
+            for rank, s in enumerate(sols, start=1):
+                rows.append({"sample_id": sid, "sample_name": sample_name, "solution_id": rank, "gene": gene,
+                             "call_type": "star_allele", "solution": s.genotype,
+                             "array_equivalent_alleles_1": ";".join(sorted({mm[0] for mm in s.members})),
+                             "array_equivalent_alleles_2": ";".join(sorted({mm[-1] for mm in s.members})),
+                             "no_call_variants": ";".join(m.variant_ids[k][0] for k in no_calls),
+                             "alleles_impacted_by_no_call": "", "background_alleles": "", "variant_change": "",
+                             "score": f"{s.score:.4f}", "raw_score": f"{s.raw_score:.4f}",
+                             "copy_number_solution": f"{region}={gcn}" if region and gcn is not None else "",
+                             "flags": ";".join(flags + (["feasible"] if s in feasible else []))})  # fmt: skip
+            genes_out.append({"gene": gene, "genotype": genotype, "feasible_solutions": [s.genotype for s in feasible],
+                              "candidates": [_sol_json(s) for s in sols], "copy_number": gcn,
+                              "flags": flags})  # fmt: skip
+        calls_json.append({"sample_id": sid, "genes": genes_out})
+    cols = ["sample_id", "sample_name", "solution_id", "gene", "call_type", "solution", "array_equivalent_alleles_1",
+            "array_equivalent_alleles_2", "no_call_variants", "alleles_impacted_by_no_call", "background_alleles",
+            "variant_change", "score", "raw_score", "copy_number_solution", "flags"]  # fmt: skip
+    o.out.mkdir(parents=True, exist_ok=True)
+    write_tsv(o.out / "diplotype_summary.tsv", cols, rows)
+    (o.out / "pgx_calls.json").write_text(json.dumps({"software": f"SimurgArray {__version__}", "db": o.db.name,
+                                                      "samples": calls_json}, indent=1))  # fmt: skip
+    rep.info(f"star alleles: {len(calls_json)} sample(s), {len(models)} gene(s)")
+    return ExitCode.OK
+
+
+def _sol_json(s: Solution) -> dict[str, Any]:
+    return {"genotype": s.genotype, "score": round(s.score, 4), "raw_score": round(s.raw_score, 4),
+            "equivalents": ["/".join(m) for m in s.members][:20], "copy_number": s.copy_number}  # fmt: skip