fix(genotype): novel-allele review round 4 — no unplaced leak, arg validation, docs reproducibility

Correctness:
- effective_dataconfig injects ONLY carried novel alleles (defined-but-unplaced
  novels no longer leak into the aligner reference / v_call). High.
- add_novel_allele tolerates missing/empty/mismatched gapped_seq (falls back to
  ungapped) instead of IndexError. Medium.
- delete_gene / duplicate_gene validate the haplotype argument (both/0/1) instead
  of silently no-op'ing or negative-indexing. Medium.

Docs:
- 'Reproduce it' now builds the EXACT 9-gene planted genotype behind the figure,
  writes reads.fasta (IgDiscover/partis) + germline + truth, and states the
  scoring + reported precision/recall. High.
- Added a method-signature table (segment/haplotype defaults) + a D/J example;
  a 'Supported loci and chains' note (VDJ vs VJ, BCR/TCR); 'at compile time'
  wording fix.
Tests: unplaced-novel non-leak, haplotype validation, missing-gapped-seq fallback.

Author: MuteJester

site_docs/guides/genotype.md

@@ -35,6 +35,15 @@ a D and a J **from the same chromosome**. That linkage is exactly the signal
 haplotype-inference methods exploit (e.g. the IGHJ6-anchor approach), and GenAIRR
 reproduces it.
 
+!!! note "Supported loci and chains"
+    Genotypes work on any GenAIRR reference cartridge — BCR **and** TCR, heavy
+    **and** light/α/β chains. On **VDJ** loci (IGH, TRB, TRD) the genotype spans
+    V, D and J and each rearrangement draws all three from one chromosome. On
+    **VJ** loci (IGK, IGL, TRA, TRG) there is no D segment: genotype V and J,
+    D rows are simply not required and are ignored. The examples below use the
+    human IGH cartridge, but the same API applies to every locus; just use that
+    cartridge's gene/allele names.
+
 ## Quick start
 
 ```python
@@ -100,7 +109,8 @@ the J drawn later (the phased choices are evaluated together, not independently)
 ### Strict vs permissive
 
 `Genotype.from_dataconfig(cfg)` is **strict**: any gene that could be used during
-recombination but was never specified is an error at attach time — you must define
+recombination but was never specified is an error when the experiment is compiled
+(`compile()` / `run_records()`) — you must define
 the whole genotype (use `complete_from_reference` to fill the genes you don't care
 about). This guarantees a genuine diploid complement, which is what you want for a
 ground-truth benchmark.
@@ -136,6 +146,27 @@ g.with_subject("DONOR01")                               # provenance label
 g.complete_from_reference("homozygous_first_reference") # fill every unspecified gene
 ```
 
+Every editing method takes a `segment` argument (`"V"` default, or `"D"` / `"J"`),
+so genotype the D and J loci too — important since J anchors and D/J usage drive
+haplotype-inference methods:
+
+| Method | Signature | Notes |
+|---|---|---|
+| `homozygous` | `(gene, allele, segment="V")` | one allele on both chromosomes |
+| `heterozygous` | `(gene, allele0, allele1, segment="V")` | one allele per chromosome |
+| `delete_gene` | `(gene, haplotype="both"\|0\|1, segment="V")` | whole-gene or one-chromosome (hemizygous) deletion |
+| `duplicate_gene` | `(gene, alleles=[...], haplotype=0\|1, segment="V")` | >1 copy on one chromosome |
+| `add_novel_allele` | `(name, *, base, mutations\|sequence, segment="V", allow_nonfunctional=False)` | define a private allele (see below) |
+| `chromosome_weights` | `(w0, w1)` | allelic-expression imbalance (default 0.5/0.5) |
+| `with_subject` | `(sid)` | provenance label stamped on every record |
+| `complete_from_reference` | `(policy="homozygous_first_reference"\|"heterozygous_first_two")` | fill unspecified genes |
+
+```python
+# Genotype the J locus too — e.g. heterozygous IGHJ6 + a homozygous IGHJ4:
+g.heterozygous("IGHJ6", "IGHJ6*02", "IGHJ6*03", segment="J")
+g.homozygous("IGHJ4", "IGHJ4*02", segment="J")
+```
+
 Notes and guard-rails:
 
 - **`delete_gene(..., haplotype=0|1)`** (one chromosome) requires the gene to be
@@ -359,36 +390,57 @@ deletions correct.*
 
 ### Reproduce it
 
+This builds the **exact** genotype behind the figure — 3 heterozygous, 3
+homozygous, and 3 deleted study V genes, the rest filled from the reference —
+simulates 4,000 reads with light SHM at `seed=7`, and writes every input the two
+tools need plus the ground truth to score against:
+
 ```python
 import GenAIRR as ga
 import GenAIRR.data as gdata
 from GenAIRR.genotype import Genotype
 
 cfg = gdata.HUMAN_IGH_OGRDB
-g = (
-    Genotype.from_dataconfig(cfg)
-      .complete_from_reference("homozygous_first_reference")
-      .heterozygous("IGHVF1-G1", "IGHVF1-G1*01", "IGHVF1-G1*02")
-      .homozygous("IGHVF2-G4", "IGHVF2-G4*01")
-      .delete_gene("IGHVF3-G7", haplotype="both")
-      .with_subject("DONOR01")
-)
+HET = ["IGHVF1-G1", "IGHVF1-G2", "IGHVF1-G3"]   # 2 alleles each
+HOM = ["IGHVF2-G4", "IGHVF3-G5", "IGHVF3-G6"]   # 1 allele
+DEL = ["IGHVF3-G7", "IGHVF3-G8", "IGHVF3-G9"]   # deleted (both chromosomes)
+
+g = Genotype.from_dataconfig(cfg).complete_from_reference("homozygous_first_reference")
+for gene in HET:
+    a0, a1 = (a.name for a in cfg.v_alleles[gene][:2])
+    g.heterozygous(gene, a0, a1)
+for gene in HOM:
+    g.homozygous(gene, cfg.v_alleles[gene][0].name)
+for gene in DEL:
+    g.delete_gene(gene, haplotype="both")
+g.with_subject("DONOR01")
+
 res = (
     ga.Experiment.on(cfg).with_genotype(g).recombine()
       .mutate(rate=0.004)                       # light SHM, as in real data
-      .run_records(n=4000, seed=7)
+      .run_records(n=4000, seed=7, expose_provenance=True)
 )
+
 res.to_tsv("repertoire.tsv")                    # AIRR table → TIgGER
 g.to_tsv("truth_genotype.tsv")                  # ground truth to score against
 
-# export the cartridge V germline (names match v_call) for TIgGER's germline_db
-with open("germline_V.fasta", "w") as fh:
+with open("reads.fasta", "w") as fh:            # raw reads → IgDiscover / partis
+    for r in res:
+        fh.write(f">{r['sequence_id']}\n{r['sequence'].upper()}\n")
+
+with open("germline_V.fasta", "w") as fh:       # cartridge V germline (names match v_call)
     for gene, alleles in cfg.v_alleles.items():
         for a in alleles:
             fh.write(f">{a.name}\n{a.ungapped_seq.upper()}\n")
 ```
 
-Then run the R snippet above and compare `geno` against `truth_genotype.tsv`.
+**Score it.** Run TIgGER (R snippet above) on `repertoire.tsv`, or IgDiscover on
+`reads.fasta` with the cartridge as its starting database
+(`igdiscover init --database db/ --single-reads reads.fasta project/ && cd project
+&& igdiscover run`). Then compare each tool's per-gene allele set against
+`g.to_table()` (the planted truth): allele-presence precision/recall, zygosity,
+and deletion calls. With the genotype above this yields TIgGER precision/recall
+1.00 (52/52 genes) and IgDiscover precision 1.00 / recall 0.96 — the figure.
 
 ### Running other tools on the same data
 

src/GenAIRR/genotype.py

@@ -109,6 +109,10 @@ def heterozygous(
         return self
 
     def delete_gene(self, gene: str, haplotype="both", segment: str = "V") -> "Genotype":
+        if haplotype not in ("both", 0, 1):
+            raise ValueError(
+                f"haplotype must be 'both', 0, or 1, got {haplotype!r}"
+            )
         # One-haplotype (hemizygous) deletion requires the gene to be
         # specified first, otherwise the *other* haplotype would also be
         # empty — silently producing a full deletion that
@@ -133,11 +137,13 @@ def delete_gene(self, gene: str, haplotype="both", segment: str = "V") -> "Genot
     def duplicate_gene(
         self, gene: str, alleles: List[str], haplotype: int, segment: str = "V"
     ) -> "Genotype":
+        if haplotype not in (0, 1):
+            raise ValueError(f"haplotype must be 0 or 1, got {haplotype!r}")
         for a in alleles:
             self._check_allele(segment, gene, a)
         cur = self._slots[segment].get(gene, [[], []])
         cur = [list(cur[0]), list(cur[1])]
-        cur[int(haplotype)] = [(a, 1, 1.0) for a in alleles]
+        cur[haplotype] = [(a, 1, 1.0) for a in alleles]
         self._slots[segment][gene] = cur
         return self
 
@@ -210,9 +216,13 @@ def add_novel_allele(
             raise ValueError("provide exactly one of `mutations` or `sequence`")
 
         base_ungapped = base_allele.ungapped_seq.upper()
-        gapped = list(base_allele.gapped_seq)
-        # ungapped index -> gapped index (positions of non-gap characters)
-        ung_to_gap = [i for i, ch in enumerate(base_allele.gapped_seq) if ch != "."]
+        gapped = list(base_allele.gapped_seq or "")
+        # ungapped index -> gapped index (positions of non-gap characters).
+        # Some custom cartridges carry no (or inconsistent) gapped sequence;
+        # in that case we can't project onto gaps, so fall back to an
+        # ungapped novel sequence (no gap-derived metadata).
+        ung_to_gap = [i for i, ch in enumerate(base_allele.gapped_seq or "") if ch != "."]
+        project_gaps = len(ung_to_gap) == len(base_ungapped)
         seq = list(base_ungapped)
         if sequence is not None:
             sequence = sequence.upper()
@@ -241,15 +251,20 @@ def add_novel_allele(
         if new_ungapped == base_ungapped:
             raise ValueError("novel allele is identical to its base allele")
         # Project the substitutions onto the gapped sequence too, so
-        # gap-dependent metadata stays consistent.
-        for k, b in enumerate(seq):
-            gapped[ung_to_gap[k]] = b
+        # gap-dependent metadata stays consistent. If the base has no
+        # usable gapped sequence, fall back to the ungapped form.
+        if project_gaps:
+            for k, b in enumerate(seq):
+                gapped[ung_to_gap[k]] = b
+            new_gapped = "".join(gapped)
+        else:
+            new_gapped = new_ungapped
 
         novel = _copy.deepcopy(base_allele)
         novel.name = name
         novel.gene = gene
         novel.ungapped_seq = new_ungapped
-        novel.gapped_seq = "".join(gapped)
+        novel.gapped_seq = new_gapped
         if hasattr(novel, "ungapped_len"):
             novel.ungapped_len = len(new_ungapped)
 
@@ -296,15 +311,31 @@ def has_novel(self) -> bool:
     def novel_allele_names(self) -> Set[str]:
         return set(self._novel)
 
+    def _carried_allele_names(self) -> Set[str]:
+        """Allele names actually placed on a haplotype (across segments)."""
+        names: Set[str] = set()
+        for seg in _SEGMENTS:
+            for haps in self._slots[seg].values():
+                for hap in haps:
+                    names.update(a for (a, _c, _w) in hap)
+        return names
+
     def effective_dataconfig(self):
         """Return a copy of the source ``DataConfig`` with this genotype's
-        novel alleles appended to their genes' allele lists — the reference
-        the engine actually runs against when novel alleles are present."""
+        **carried** novel alleles appended to their genes' allele lists —
+        the reference the engine actually runs against when novel alleles
+        are present. A novel allele that was defined but never placed on a
+        haplotype is NOT injected (it would otherwise pollute the aligner
+        reference and could surface in ``v_call`` despite being absent from
+        the ground truth)."""
         import copy as _copy
 
         cfg = _copy.deepcopy(self._cfg)
         by_seg = {"V": cfg.v_alleles, "D": cfg.d_alleles, "J": cfg.j_alleles}
-        for info in self._novel.values():
+        carried = self._carried_allele_names()
+        for name, info in self._novel.items():
+            if name not in carried:
+                continue
             d = by_seg[info["segment"]]
             existing = list(d.get(info["gene"], []))
             existing.append(_copy.deepcopy(info["allele"]))

tests/test_genotype_novel.py

@@ -136,6 +136,57 @@ def test_to_table_and_tsv_expose_novel(tmp_path):
     assert "novel" in header
 
 
+def test_unplaced_novel_allele_not_injected_or_emitted():
+    cfg = _cfg()
+    gene = "IGHVF1-G1"
+    # Define a novel allele but NEVER place it on a haplotype.
+    g = (
+        Genotype.from_dataconfig(cfg)
+        .add_novel_allele(f"{gene}*unplaced", base=f"{gene}*01", mutations=_SAFE)
+        .complete_from_reference()
+        .with_subject("S1")
+    )
+    # Effective reference must not contain the unplaced novel allele.
+    eff = g.effective_dataconfig()
+    eff_names = {a.name for alleles in eff.v_alleles.values() for a in alleles}
+    assert f"{gene}*unplaced" not in eff_names
+    res = ga.Experiment.on(cfg).with_genotype(g).recombine().run_records(
+        n=200, seed=4, expose_provenance=True
+    )
+    assert all(f"{gene}*unplaced" not in r["v_call"] for r in res)
+    assert all(f"{gene}*unplaced" not in r["truth_v_call"] for r in res)
+
+
+def test_haplotype_argument_validation():
+    cfg = _cfg()
+    gene = "IGHVF1-G1"
+    a0 = cfg.v_alleles[gene][0].name
+    g = Genotype.from_dataconfig(cfg).homozygous(gene, a0)
+    with pytest.raises(ValueError, match="haplotype must be"):
+        g.delete_gene(gene, haplotype=2)
+    with pytest.raises(ValueError, match="haplotype must be"):
+        g.delete_gene(gene, haplotype="x")
+    with pytest.raises(ValueError, match="haplotype must be 0 or 1"):
+        Genotype.from_dataconfig(cfg).duplicate_gene(gene, [a0], haplotype=-1)
+    with pytest.raises(ValueError, match="haplotype must be 0 or 1"):
+        Genotype.from_dataconfig(cfg).duplicate_gene(gene, [a0], haplotype=2)
+
+
+def test_novel_synthesis_tolerates_missing_gapped_seq():
+    import copy
+
+    cfg = copy.deepcopy(_cfg())
+    # Simulate a custom cartridge whose base allele has no gapped sequence.
+    cfg.v_alleles["IGHVF1-G1"][0].gapped_seq = ""
+    g = Genotype.from_dataconfig(cfg).add_novel_allele(
+        "IGHVF1-G1*i01", base="IGHVF1-G1*01", mutations=_SAFE
+    )
+    nv = g._novel["IGHVF1-G1*i01"]["allele"]
+    # falls back to ungapped form (no crash)
+    assert nv.gapped_seq == nv.ungapped_seq
+    assert nv.ungapped_seq[38] == "C"
+
+
 def test_genotype_without_novel_is_unaffected():
     cfg = _cfg()
     g = Genotype.from_dataconfig(cfg).complete_from_reference().with_subject("S1")