PGx: CYP2D6 structures from sub-region CN — *68+*4, *36+*10, *13 (PGX-013/018)

Sub-region CNs (p5, intron 2, exon 9) resolve full copies and *36/*68/*13 units; candidates carry the
tandem names, background alleles, and CPIC activity sums. GeT-RM: CYP2D6 13/15 top (NA12878 *3/*68+*4 and
NA11832 *1/*68+*4 exact), 15/15 within feasible; all genes 132/144.

src/simurg/pgx/star.py

@@ -115,14 +115,81 @@ def candidates(m: GeneModel, cn: int | None) -> list[tuple[tuple[int, ...], tupl
     return out
 
 
+@dataclass(frozen=True)
+class Cyp2d6Structure:
+    """CYP2D6 copies from the sub-region CNs (PGX-013) [OURS; DOC: PharmVar CYP2D6 structural variation]:
+    ``full`` complete CYP2D6 copies, ``u36`` *36-type units (CYP2D7 exon 9), ``u68`` *68-type units (CYP2D6 exon 1
+    only), ``u13`` *13-type units (CYP2D7 5′ part). CN(p5) = full + u36 + u68, CN(intron 2) = full + u36 + u13,
+    CN(exon 9) = full + u13; with u68 · u13 = 0 this has one solution."""
+
+    full: int
+    u36: int
+    u68: int
+    u13: int
+
+
+def cyp2d6_structure(p5: int | None, intron2: int | None, exon9: int | None) -> Cyp2d6Structure | None:
+    if p5 is None:
+        return None
+    ni = p5 if intron2 is None else intron2
+    n9 = ni if exon9 is None else exon9
+    u36 = max(0, ni - n9)
+    x = p5 - u36 - n9
+    if x >= 0:
+        full, u68, u13 = n9, x, 0
+    else:
+        full, u68, u13 = p5 - u36, 0, -x
+    if full < 0:
+        return Cyp2d6Structure(max(0, p5), 0, 0, 0)
+    return Cyp2d6Structure(full, u36, u68, u13)
+
+
+def structure_candidates(m: GeneModel, st: Cyp2d6Structure) -> list[tuple[tuple[int, ...], tuple[str, ...]]]:
+    """Candidates for a CYP2D6 structure: SNV copies = full + u36 (a *36 unit carries the *10 sequence up to exon 8);
+    *36 units attach to a *10 copy (``*36+*10``), *68 units to *4 when present (``*68+*4``), *13 units are separate
+    haplotypes."""
+    snv_copies = st.full + st.u36
+    base = candidates(m, snv_copies)
+    i10 = m.haplotypes.index("*10") if "*10" in m.haplotypes else None
+    out = []
+    for ix, names in base:
+        names_l = list(names)
+        if st.u36:
+            if i10 is None or ix.count(i10) < st.u36 + 1:  # the *36 unit(s) plus the *10 they precede
+                continue
+            k = next((n for n, x in enumerate(names_l) if x.startswith("*10")), None)
+            if k is None:
+                continue
+            unit = "*36" if st.u36 == 1 else f"*36x{st.u36}"
+            names_l[k] = f"{unit}+*10"
+        if st.u68:
+            unit = "*68" if st.u68 == 1 else f"*68x{st.u68}"
+            k = next((n for n, x in enumerate(names_l) if x.split("+")[-1].split("x")[0] == "*4"), 0)
+            names_l[k] = f"{unit}+{names_l[k]}"
+        if st.u13:
+            if len(names_l) == 2 and names_l[1] == (m.deletion_allele or "*0"):
+                names_l[1] = "*13" if st.u13 == 1 else f"*13x{st.u13}"
+            else:
+                names_l[0] = f"*13+{names_l[0]}"
+        out.append((ix, tuple(names_l)))
+    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,
+    max_candidates: int = 10, feasible_ratio: float = 0.5, structure: Cyp2d6Structure | None = None,
 ) -> 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)
+    haploid: list[tuple[tuple[int, ...], tuple[str, ...]]] = [((i,), (m.haplotypes[i],)) for i in
+                                                               range(len(m.haplotypes))]  # fmt: skip
+    if ploidy == 1:
+        cands = haploid
+    elif structure is not None and (structure.u36 or structure.u68 or structure.u13):
+        cands = structure_candidates(m, structure)
+        cn = structure.full + structure.u36
+    else:
+        cands = candidates(m, cn)
     if not cands:
         return [], []
     carries = m.carries.astype(np.int16)