fix(genotype): novel-allele review round — functional validation, gene identity, gapped projection

Addresses critic findings on the novel-allele slice:
- #1 functional validation: synthesized V/J coding sequence is checked for an
  intact conserved anchor codon (Cys/Trp|Phe) and stop-free coding frame; a
  broken variant is rejected unless allow_nonfunctional=True (then kept + marked
  non-functional). Closes the 'nonfunctional emitted as productive' gap.
- #2 gene identity: the novel allele's gene is taken from its NAME and must equal
  the base allele's gene; dropped the gene= override that left allele.gene stale.
- #3 anchor: inherited from base (correct for same-length/substitution-only
  variants — the conserved residue does not move) and validated to remain intact;
  no reliance on the unavailable _native anchor resolver.
- #4 name uniqueness enforced across all segments + novel set (prevents truth-table
  mislabeling).
- #5 to_tsv now emits the 'novel' column.
- #6 substitutions projected onto the gapped sequence (no stale gapped_seq).
- #7 mutation positions/bases type-checked with clean ValueErrors.
Tests expanded: stop/anchor rejection + allow_nonfunctional, gene-mismatch,
cross-segment collision, tsv export, type-check.

Author: MuteJester

site_docs/guides/genotype.md

@@ -250,7 +250,7 @@ so it flows through alignment and AIRR output as a genuine allele:
 g = (
     Genotype.from_dataconfig(cfg)
       .add_novel_allele("IGHVF1-G1*i01", base="IGHVF1-G1*01",
-                        mutations=[(120, "T"), (250, "G")])   # two point variants
+                        mutations=[(38, "C"), (41, "A")])   # two point variants
       .complete_from_reference()
       .heterozygous("IGHVF1-G1", "IGHVF1-G1*01", "IGHVF1-G1*i01")  # one reference + one private
       .with_subject("DONOR_N")
@@ -264,8 +264,14 @@ result = (
 # its name appears in v_call / truth_v_call and the reads carry its variants.
 ```
 
-Novel alleles are flagged in the ground truth: each `to_table()` row carries a
-`novel` list of the private alleles carried at that gene.
+The novel allele's **gene is taken from its name** and must match the base
+allele's gene; it must be a same-length (substitution-only) variant. The
+synthesized coding sequence is **validated** — for V/J the conserved anchor codon
+must still encode the conserved residue (Cys for V, Trp/Phe for J) and the coding
+frame must be stop-free. A variant that breaks either is rejected unless you pass
+`allow_nonfunctional=True` (then it is kept and marked non-functional). Novel
+alleles are flagged in the ground truth: each `to_table()`/`to_tsv()` row carries
+a `novel` list of the private alleles carried at that gene.
 
 **Benchmarking novel-allele discovery.** Plant a novel allele, simulate, then run
 the discovery tool against the **base** germline (the cartridge *without* your

src/GenAIRR/genotype.py

@@ -157,71 +157,136 @@ def add_novel_allele(
         mutations: Optional[List[Tuple[int, str]]] = None,
         sequence: Optional[str] = None,
         segment: str = "V",
-        gene: Optional[str] = None,
+        allow_nonfunctional: bool = False,
     ) -> "Genotype":
         """Define a private/novel allele not present in the reference.
 
-        Derive it from a reference ``base`` allele by either applying
-        point ``mutations`` (a list of ``(0-based position, base)``) or by
-        supplying an explicit ``sequence`` (same length as the base, so the
-        inherited anchor/sub-regions stay valid). Gene, anchor, functional
-        status and V sub-regions are inherited from the base allele.
+        Derive it from a reference ``base`` allele by either applying point
+        ``mutations`` (a list of ``(0-based position, base)``) or supplying
+        an explicit same-length ``sequence`` (substitutions only — no
+        indels — so the gene's reading frame and the conserved-anchor
+        position are preserved). The novel allele's **gene is taken from
+        its name** and must equal the base allele's gene; sub-regions and
+        the anchor position are inherited (valid for same-length variants).
+
+        The synthesized coding sequence is **validated**: for V/J the
+        conserved anchor codon must still encode the conserved residue
+        (Cys for V, Trp/Phe for J) and the coding frame must contain no
+        internal stop codon. A variant that breaks either is rejected
+        unless ``allow_nonfunctional=True`` (in which case it is kept and
+        marked non-functional).
 
         Registers the novel allele under ``name`` so it can then be placed
-        with :meth:`homozygous` / :meth:`heterozygous` / :meth:`duplicate_gene`
-        like any allele. At compile time it is injected as a real entry in
-        an *effective* reference, so it flows through alignment and AIRR
-        output exactly like a catalogue allele.
+        with :meth:`homozygous` / :meth:`heterozygous` / :meth:`duplicate_gene`.
+        At compile time it is injected as a real entry in an *effective*
+        reference, so it flows through alignment and AIRR output like a
+        catalogue allele.
         """
         import copy as _copy
 
+        from .utilities.misc import translate
+
         base_allele = self._find_ref_allele(segment, base)
         if base_allele is None:
             raise ValueError(f"base allele {base!r} not found in {segment} reference")
-        gene = gene or base_allele.gene
-        if gene not in _alleles_by_gene(self._cfg, segment):
-            raise ValueError(f"{gene!r} is not a known {segment} gene in this cartridge")
-        # name must not collide with a catalogue allele or another novel.
-        if self._find_ref_allele(segment, name) is not None or name in self._novel:
-            raise ValueError(f"novel allele name {name!r} collides with an existing allele")
+        # Gene identity comes from the name; it must match the base's gene
+        # (no cross-gene synthesis — that would corrupt gene identity).
+        name_gene = name.split("*")[0]
+        if name_gene != base_allele.gene:
+            raise ValueError(
+                f"novel name {name!r} implies gene {name_gene!r} but base {base!r} "
+                f"belongs to gene {base_allele.gene!r}; a novel allele must belong "
+                f"to its base allele's gene"
+            )
+        gene = base_allele.gene
+        # Name must be unique across the WHOLE catalogue (all segments) and
+        # all previously-defined novel alleles.
+        if name in self._novel:
+            raise ValueError(f"novel allele name {name!r} already defined")
+        for seg in _SEGMENTS:
+            if self._find_ref_allele(seg, name) is not None:
+                raise ValueError(f"novel allele name {name!r} collides with a catalogue allele")
         if (mutations is None) == (sequence is None):
             raise ValueError("provide exactly one of `mutations` or `sequence`")
 
-        seq = list(base_allele.ungapped_seq.upper())
+        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 != "."]
+        seq = list(base_ungapped)
         if sequence is not None:
             sequence = sequence.upper()
             if len(sequence) != len(seq):
                 raise ValueError(
                     f"explicit sequence length {len(sequence)} != base length {len(seq)}; "
-                    "use `mutations` for indels-free variants or match the base length"
+                    "novel alleles are substitution-only (same length as the base)"
                 )
             if any(b not in "ACGT" for b in sequence):
                 raise ValueError("sequence must contain only A/C/G/T")
-            new_seq = sequence
+            seq = list(sequence)
         else:
             if not mutations:
                 raise ValueError("`mutations` must be a non-empty list of (position, base)")
             for pos, b in mutations:
+                if not isinstance(pos, int) or isinstance(pos, bool):
+                    raise ValueError(f"mutation position must be an int, got {pos!r}")
+                if not (isinstance(b, str) and len(b) == 1):
+                    raise ValueError(f"mutation base must be a single character, got {b!r}")
                 if not (0 <= pos < len(seq)):
                     raise ValueError(f"mutation position {pos} out of range [0,{len(seq)})")
                 if b.upper() not in "ACGT":
                     raise ValueError(f"mutation base {b!r} must be A/C/G/T")
                 seq[pos] = b.upper()
-            new_seq = "".join(seq)
-        if new_seq == base_allele.ungapped_seq.upper():
+        new_ungapped = "".join(seq)
+        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
 
         novel = _copy.deepcopy(base_allele)
         novel.name = name
-        novel.ungapped_seq = new_seq
+        novel.gene = gene
+        novel.ungapped_seq = new_ungapped
+        novel.gapped_seq = "".join(gapped)
         if hasattr(novel, "ungapped_len"):
-            novel.ungapped_len = len(new_seq)
+            novel.ungapped_len = len(new_ungapped)
+
+        # Functional validation (V/J have a conserved coding frame).
+        functional, reason = True, None
+        anchor = getattr(novel, "anchor", None)
+        if segment in ("V", "J") and anchor is not None:
+            conserved = {"V": {"C"}, "J": {"W", "F"}}[segment]
+            anchor_aa = translate(new_ungapped[anchor : anchor + 3])
+            if anchor_aa not in conserved:
+                functional = False
+                reason = (
+                    f"conserved anchor codon now encodes {anchor_aa!r}, "
+                    f"expected one of {sorted(conserved)}"
+                )
+            coding = new_ungapped[:anchor] if segment == "V" else new_ungapped[anchor:]
+            if "*" in translate(coding):
+                reason = (reason + "; " if reason else "") + "internal stop codon in coding frame"
+                functional = False
+        if not functional and not allow_nonfunctional:
+            raise ValueError(
+                f"novel allele {name!r} is non-functional ({reason}); pass "
+                f"allow_nonfunctional=True to keep it anyway"
+            )
+        if not functional:
+            try:
+                novel.functional_status = "pseudogene"
+            except Exception:
+                pass
+
         self._novel[name] = {
             "allele": novel,
             "gene": gene,
             "segment": segment,
             "base": base,
             "mutations": list(mutations) if mutations else None,
+            "functional": functional,
         }
         return self
 
@@ -366,6 +431,7 @@ def _fmt(detail):
                     "zygosity",
                     "haplotype_0",
                     "haplotype_1",
+                    "novel",
                     "permissive",
                 ]
             )
@@ -378,6 +444,7 @@ def _fmt(detail):
                         r["zygosity"],
                         _fmt(r["haplotype_0_detail"]),
                         _fmt(r["haplotype_1_detail"]),
+                        ";".join(r["novel"]),
                         r["permissive"],
                     ]
                 )

tests/test_genotype_novel.py

@@ -1,4 +1,4 @@
-"""Novel / private allele support on genotypes (PR additions)."""
+"""Novel / private allele support on genotypes."""
 import pytest
 
 import GenAIRR as ga
@@ -10,46 +10,97 @@ def _cfg():
     return gdata.HUMAN_IGH_OGRDB
 
 
+# Guaranteed-safe SNPs for IGHVF1-G1*01: wobble of non-T-starting codons in
+# the framework — cannot create a stop, anchor untouched.
+_SAFE = [(38, "C"), (41, "A")]
+
+
 def test_add_novel_allele_synthesizes_from_base_and_mutations():
     cfg = _cfg()
     base = cfg.v_alleles["IGHVF1-G1"][0]
     g = Genotype.from_dataconfig(cfg).add_novel_allele(
-        "IGHVF1-G1*i01", base="IGHVF1-G1*01", mutations=[(120, "T"), (130, "A")]
+        "IGHVF1-G1*i01", base="IGHVF1-G1*01", mutations=_SAFE
     )
     assert g.has_novel()
     assert "IGHVF1-G1*i01" in g.novel_allele_names()
     nv = g._novel["IGHVF1-G1*i01"]["allele"]
-    assert nv.ungapped_seq[120] == "T" and nv.ungapped_seq[130] == "A"
-    assert len(nv.ungapped_seq) == len(base.ungapped_seq)  # point mutations: no length change
-    assert nv.gene == "IGHVF1-G1" and nv.anchor == base.anchor  # metadata inherited
+    assert nv.ungapped_seq[38] == "C" and nv.ungapped_seq[41] == "A"
+    assert len(nv.ungapped_seq) == len(base.ungapped_seq)  # substitution-only
+    assert nv.gene == "IGHVF1-G1" and nv.anchor == base.anchor  # gene/anchor inherited
+    # gapped sequence projected consistently (ungapped derived from it)
+    assert nv.gapped_seq.replace(".", "") == nv.ungapped_seq
+    assert g._novel["IGHVF1-G1*i01"]["functional"] is True
 
 
-def test_novel_allele_validation():
+def test_novel_name_gene_must_match_base_gene():
+    cfg = _cfg()
+    with pytest.raises(ValueError, match="implies gene"):
+        Genotype.from_dataconfig(cfg).add_novel_allele(
+            "IGHVF1-G2*i01", base="IGHVF1-G1*01", mutations=_SAFE  # name gene != base gene
+        )
+
+
+def test_novel_allele_basic_validation():
     cfg = _cfg()
     G = lambda: Genotype.from_dataconfig(cfg)
     with pytest.raises(ValueError, match="base allele"):
-        G().add_novel_allele("X*i01", base="NOPE*01", mutations=[(1, "T")])
-    with pytest.raises(ValueError, match="collides"):
-        G().add_novel_allele("IGHVF1-G1*01", base="IGHVF1-G1*01", mutations=[(1, "T")])
+        G().add_novel_allele("NOPE*i01", base="NOPE*01", mutations=[(1, "T")])
+    with pytest.raises(ValueError, match="collides with a catalogue"):
+        G().add_novel_allele("IGHVF1-G1*01", base="IGHVF1-G1*01", mutations=_SAFE)
     with pytest.raises(ValueError, match="out of range"):
-        G().add_novel_allele("X*i01", base="IGHVF1-G1*01", mutations=[(99999, "T")])
+        G().add_novel_allele("IGHVF1-G1*i01", base="IGHVF1-G1*01", mutations=[(99999, "T")])
     with pytest.raises(ValueError, match="exactly one"):
-        G().add_novel_allele("X*i01", base="IGHVF1-G1*01")  # neither mutations nor sequence
-    with pytest.raises(ValueError, match="identical"):
-        # mutate to the same base it already is
-        b = cfg.v_alleles["IGHVF1-G1"][0].ungapped_seq.upper()
-        G().add_novel_allele("X*i01", base="IGHVF1-G1*01", mutations=[(0, b[0])])
+        G().add_novel_allele("IGHVF1-G1*i01", base="IGHVF1-G1*01")
+    with pytest.raises(ValueError, match="must be an int"):
+        G().add_novel_allele("IGHVF1-G1*i01", base="IGHVF1-G1*01", mutations=[(1.5, "A")])
+
+
+def test_cross_segment_name_collision_rejected():
+    # A novel can't be named for a different gene/segment than its base, so
+    # naming a V-derived novel after a real J allele is rejected outright —
+    # the truth table can never mislabel the real J row as novel.
+    cfg = _cfg()
+    j_name = next(iter(cfg.j_alleles[next(iter(cfg.j_alleles))])).name
+    with pytest.raises(ValueError):
+        Genotype.from_dataconfig(cfg).add_novel_allele(
+            j_name, base="IGHVF1-G1*01", mutations=_SAFE
+        )
+
+
+def test_nonfunctional_novel_rejected_by_default_and_allowed_explicitly():
+    cfg = _cfg()
+    base = cfg.v_alleles["IGHVF1-G1"][0]
+    # Force a stop codon at framework codon 13 (positions 39,40,41 -> TAA).
+    stop = [(39, "T"), (40, "A"), (41, "A")]
+    with pytest.raises(ValueError, match="non-functional.*stop codon"):
+        Genotype.from_dataconfig(cfg).add_novel_allele(
+            "IGHVF1-G1*i01", base="IGHVF1-G1*01", mutations=stop
+        )
+    # Explicit override keeps it, marked non-functional.
+    g = Genotype.from_dataconfig(cfg).add_novel_allele(
+        "IGHVF1-G1*i01", base="IGHVF1-G1*01", mutations=stop, allow_nonfunctional=True
+    )
+    assert g._novel["IGHVF1-G1*i01"]["functional"] is False
+
+
+def test_broken_anchor_codon_rejected():
+    cfg = _cfg()
+    base = cfg.v_alleles["IGHVF1-G1"][0]
+    a = base.anchor
+    # rewrite the conserved Cys anchor codon to GGG (Gly)
+    with pytest.raises(ValueError, match="conserved anchor codon"):
+        Genotype.from_dataconfig(cfg).add_novel_allele(
+            "IGHVF1-G1*i01", base="IGHVF1-G1*01",
+            mutations=[(a, "G"), (a + 1, "G"), (a + 2, "G")],
+        )
 
 
 def test_novel_allele_placed_and_sampled_appears_in_airr():
     cfg = _cfg()
     gene = "IGHVF1-G1"
-    base_seq = cfg.v_alleles[gene][0].ungapped_seq.upper()
-    pos = 120
-    new = "T" if base_seq[pos] != "T" else "A"
     g = (
         Genotype.from_dataconfig(cfg)
-        .add_novel_allele(f"{gene}*i01", base=f"{gene}*01", mutations=[(pos, new)])
+        .add_novel_allele(f"{gene}*i01", base=f"{gene}*01", mutations=_SAFE)
         .complete_from_reference()
         .homozygous(gene, f"{gene}*i01")   # carry ONLY the novel allele for this gene
         .with_subject("DONOR_N")
@@ -60,21 +111,34 @@ def test_novel_allele_placed_and_sampled_appears_in_airr():
         .recombine()
         .run_records(n=300, seed=3, expose_provenance=True)
     )
-    # every read assigned to this gene must be the novel allele (truth)
     truth_for_gene = {
         r["truth_v_call"] for r in res if r["truth_v_call"].startswith(gene + "*")
     }
     assert truth_for_gene == {f"{gene}*i01"}, truth_for_gene
-    # and the engine actually produced the mutated base in those reads
     novel_reads = [r for r in res if r["truth_v_call"] == f"{gene}*i01"]
-    assert novel_reads, "expected some reads from the novel allele"
-    assert any(r["sequence"][pos].upper() == new for r in novel_reads)
+    assert novel_reads
+    assert any(r["sequence"][38].upper() == "C" for r in novel_reads)
+
+
+def test_to_table_and_tsv_expose_novel(tmp_path):
+    cfg = _cfg()
+    gene = "IGHVF1-G1"
+    g = (
+        Genotype.from_dataconfig(cfg)
+        .add_novel_allele(f"{gene}*i01", base=f"{gene}*01", mutations=_SAFE)
+        .heterozygous(gene, f"{gene}*01", f"{gene}*i01")
+    )
+    row = next(r for r in g.to_table() if r["gene"] == gene)
+    assert row["novel"] == [f"{gene}*i01"]
+    p = tmp_path / "truth.tsv"
+    g.to_tsv(str(p))
+    header = p.read_text().splitlines()[0].split("\t")
+    assert "novel" in header
 
 
 def test_genotype_without_novel_is_unaffected():
     cfg = _cfg()
     g = Genotype.from_dataconfig(cfg).complete_from_reference().with_subject("S1")
     assert g.has_novel() is False
-    # still runs (uses base refdata)
     res = ga.Experiment.on(cfg).with_genotype(g).recombine().run_records(n=10, seed=1)
     assert len(res) == 10