diff --git a/site_docs/guides/genotype.md b/site_docs/guides/genotype.md index 19c8957..63b1b2a 100644 --- a/site_docs/guides/genotype.md +++ b/site_docs/guides/genotype.md @@ -316,6 +316,99 @@ excludes an allele), and unspecified genes fall back to uniform. explicit opt-in that reuses the cartridge's recombination `allele_usage` as a frequency proxy — convenient but biologically approximate. +## Population genotype models on a cartridge + +The `allele_frequencies` / `haplotype_deletion_prob` you pass to +`Genotype.sample` can instead be **authored once on the cartridge** as a +*population genotype model* — a donor-population germline prior. It is a distinct +plane from `reference_models.allele_usage`: `allele_usage` weights how often each +allele is *expressed* during recombination, whereas a genotype prior describes +which alleles a *donor population carries* (frequencies, gene-deletion rates, and +population novel alleles). It lives on `DataConfig.genotype_priors`. + +### Authoring and attaching a model + +```python +from GenAIRR.genotype_priors import PopulationGenotypeModel, PopulationNovelAllele + +model = PopulationGenotypeModel( + model_id="IGH-toy-1", source="hand-authored", # identity is required + allele_frequencies={"V": {"IGHVF1-G1": {"IGHVF1-G1*01": 3.0, "IGHVF1-G1*02": 1.0}}}, + haplotype_deletion_prob={"V": {"IGHVF1-G1": 0.1}}, +) +# Attach via the cartridge builder (validated against the chain type + catalogue): +# cfg = builder.set_genotype_priors(model).build() +``` + +`set_genotype_priors` validates the model against the cartridge: unknown +genes/alleles raise, and any population novel allele goes through the same +functional validation as `add_novel_allele` (conserved anchor, stop-free frame, +base-allele match). A non-`None` plane becomes part of cartridge identity (it +folds into `compute_checksum()` and the manifest). + +### Estimating a model from observed genotypes + +Given a list of observed `Genotype` objects (e.g. one per donor), estimate a +prior directly: + +```python +model = PopulationGenotypeModel.from_genotypes( + [g_donor1, g_donor2, g_donor3], + cfg=cfg, model_id="cohort-est", source="my-cohort", +) +# or, attaching to a builder in one chained step: +# builder.estimate_genotype_priors([g_donor1, g_donor2, g_donor3], source="my-cohort") +``` + +Estimator conventions: allele frequencies are counted **per carried chromosome** +(homozygous contributes 2, hemizygous 1, deleted 0); gene deletion probability is +`deleted_haplotypes / (2 × n_subjects)`. `pseudocount` smooths the per-gene +catalogue-allele counts only (deletion gets none). Genotypes carrying a +duplicated gene are rejected — the plane is deletion-only. + +### Drawing from the cartridge plane + +When a cartridge carries a plane, `Genotype.sample(cfg)` **auto-uses it** and +records where every input came from: + +```python +g = Genotype.sample(cfg, seed=7) # plane supplies freqs / deletion / weights +print(g.prior_provenance) +# {'allele_frequencies': 'cartridge', 'haplotype_deletion_prob': 'cartridge', +# 'chromosome_weights': 'cartridge', 'novel_alleles': 'none', # 'cartridge' if the model has novels +# 'model_id': 'IGH-toy-1', 'model_checksum': '…'} +print(g.to_metadata()) # subject_id + provenance + source/effective refdata hashes +``` + +Pass `use_cartridge_priors=False` for a clean uniform, catalogue-only draw (all +plane consumption — including novels — disabled). Each input is sourced +**independently**, so you can mix explicit and cartridge values: + +```python +g = Genotype.sample( + cfg, seed=7, + allele_frequencies={"V": {"IGHVF1-G1": {"IGHVF1-G1*01": 1.0}}}, # explicit + # haplotype_deletion_prob and chromosome_weights left to the plane +) +print(g.prior_provenance["allele_frequencies"]) # 'explicit' +print(g.prior_provenance["haplotype_deletion_prob"]) # 'cartridge' +print(g.prior_provenance["chromosome_weights"]) # 'cartridge' +``` + +**Population novel alleles** on the plane are *candidate* alleles for the draw: a +novel that gets sampled is carried (and flows into `v_call` / reads / truth like +any allele); a novel that isn't drawn never pollutes the output reference. By +default (`include_cartridge_novel_alleles="auto"`) novels are injected only when +the allele frequencies are cartridge- or uniform-sourced — supplying an +**explicit** frequency table keeps it explicit. Pass +`include_cartridge_novel_alleles=True` to inject them anyway, or `False` to never. + +### Auditing the plane + +`cfg.cartridge_manifest()["models"]["genotype_priors"]` reports availability, +`model_id` / `source` / `version`, the plane's `model_checksum`, per-segment gene +counts, the novel-allele count, and `source_field` (`"DataConfig.genotype_priors"`). + ## Novel / private alleles Individuals carry germline alleles that aren't in any reference — *private* or @@ -527,8 +620,6 @@ The genotype foundation is deliberately scoped. Deferred to later work: (`with_genotype` is single-subject; `result.genotypes` is a one-element list). - **External loaders** — importing genotypes from VDJbase / TIgGER / IgDiscover / partis output. -- **Cartridge genotype plane** — persisting a population genotype model in a - cartridge. - **Same-haplotype receptor revision** — `receptor_revision` with a genotype is rejected for now. diff --git a/src/GenAIRR/cartridge_builder.py b/src/GenAIRR/cartridge_builder.py index 2feff30..460da76 100644 --- a/src/GenAIRR/cartridge_builder.py +++ b/src/GenAIRR/cartridge_builder.py @@ -90,6 +90,7 @@ class _BuilderDAllele(DAllele): from .dataconfig.config_info import ConfigInfo from .dataconfig.data_config import DataConfig from .dataconfig.enums import ChainType, Species +from .genotype_priors import PopulationGenotypeModel from .reference_models import ( AlleleUsageSpec, EmpiricalDistributionSpec, @@ -346,6 +347,7 @@ def __init__(self, chain_type: ChainType) -> None: self._metadata: Optional[ConfigInfo] = None self._reference_rules: Optional[ReferenceRulesSpec] = None self._reference_models: Optional[ReferenceEmpiricalModels] = None + self._genotype_priors: Optional[PopulationGenotypeModel] = None self._report = CartridgeBuildReport() # ────────────────────────────────────────────────────────── @@ -714,6 +716,81 @@ def with_models( # Estimators # ────────────────────────────────────────────────────────── + def _build_working_cfg(self): + """A lightweight DataConfig carrying only the parsed catalogues — used by + genotype-prior validation/estimation to resolve gene/allele names and + construct throwaway Genotypes for novel functional validation.""" + return DataConfig( + name=self._name or "WORKING", + metadata=self._metadata, + v_alleles=self._v_alleles, + d_alleles=self._d_alleles or None, + j_alleles=self._j_alleles, + c_alleles=None, + ) + + def set_genotype_priors( + self, model: PopulationGenotypeModel + ) -> "ReferenceCartridgeBuilder": + """Attach a hand-authored population genotype prior, validated against + this cartridge's chain type and catalogue. Chainable.""" + from GenAIRR.genotype import Genotype + + if not isinstance(model, PopulationGenotypeModel): + raise TypeError( + f"set_genotype_priors expects a PopulationGenotypeModel, " + f"got {type(model).__name__}") + chain_label = "vdj" if self._chain_type.has_d else "vj" + model.validate(chain_type=chain_label) # catalogue-free + + # Catalogue-aware checks against the builder's pools. + cfg = self._build_working_cfg() + by_seg = {"V": cfg.v_alleles, "D": cfg.d_alleles or {}, "J": cfg.j_alleles} + for table_name, table in (("allele_frequencies", model.allele_frequencies), + ("haplotype_deletion_prob", model.haplotype_deletion_prob)): + for seg, genes in (table or {}).items(): + catalogue = by_seg[seg] + for gene, payload in genes.items(): + if gene not in catalogue: + raise ValueError( + f"genotype_priors.{table_name}: {seg} gene {gene!r} is not " + f"in the cartridge") + if table_name == "allele_frequencies": + names = {a.name for a in catalogue[gene]} + for allele in payload: + if allele not in names: + raise ValueError( + f"genotype_priors.allele_frequencies: {allele!r} is " + f"not a known allele of {gene!r}") + # Novels: reuse Genotype.add_novel_allele for full functional validation. + helper = Genotype.from_dataconfig(cfg) + for nv in model.novel_alleles: + helper.add_novel_allele( + nv.name, base=nv.base_allele, sequence=nv.sequence.upper(), + segment=nv.segment, allow_nonfunctional=nv.allow_nonfunctional) + + self._genotype_priors = model + self._report.stages.append({ + "stage": "set_genotype_priors", + "inputs": {"model_id": model.model_id, "source": model.source, + "chain_type_label": chain_label}, + "inferred": {"model_checksum": model.content_checksum(), + "novel_allele_count": len(model.novel_alleles)}, + "warnings": [], + }) + return self + + def estimate_genotype_priors( + self, genotypes, **kwargs + ) -> "ReferenceCartridgeBuilder": + """Estimate a population genotype prior from observed ``Genotype`` objects + and attach it (chainable). Thin wrapper over + :meth:`PopulationGenotypeModel.from_genotypes` followed by + :meth:`set_genotype_priors`.""" + model = PopulationGenotypeModel.from_genotypes( + genotypes, cfg=self._build_working_cfg(), **kwargs) + return self.set_genotype_priors(model) + def estimate_allele_usage( self, rearrangements: Any, @@ -2276,6 +2353,7 @@ def build(self) -> DataConfig: c_alleles=None, # v1 boundary reference_rules=self._reference_rules, reference_models=self._reference_models, + genotype_priors=self._genotype_priors, ) # Pull a manifest snapshot + checksum. The manifest call # runs before verify_integrity so the report carries the diff --git a/src/GenAIRR/dataconfig/data_config.py b/src/GenAIRR/dataconfig/data_config.py index 086b362..0a1c3c3 100644 --- a/src/GenAIRR/dataconfig/data_config.py +++ b/src/GenAIRR/dataconfig/data_config.py @@ -15,6 +15,7 @@ from GenAIRR.alleles.allele import Allele from GenAIRR.reference_models import ReferenceEmpiricalModels from GenAIRR.reference_rules import ReferenceRulesSpec +from GenAIRR.genotype_priors import PopulationGenotypeModel DEFAULT_P_NUCLEOTIDE_LENGTH_PROBS = {0: 0.50, 1: 0.25, 2: 0.15, 3: 0.07, 4: 0.03} @@ -136,6 +137,84 @@ def _allele_usage_manifest_block(cfg): } +def _genotype_priors_manifest_block(cfg): + """Build the ``models.genotype_priors`` manifest block (Slice — Cartridge + genotype plane). Audit-sized: counts and identity, never the full tables. + Reads the top-level ``DataConfig.genotype_priors`` plane (independent of + ``reference_models``).""" + def _empty(available, valid, error=None): + b = { + "available": available, + "valid": valid, + "model_id": None, + "source": None, + "version": None, + "model_checksum": None, + "segments_with_frequencies": [], + "freq_gene_counts": {"V": 0, "D": 0, "J": 0}, + "deletion_gene_counts": {"V": 0, "D": 0, "J": 0}, + "novel_allele_count": 0, + "chromosome_weights": None, + "source_field": "DataConfig.genotype_priors", + } + if error is not None: + b["validation_error"] = error + return b + + model = getattr(cfg, "genotype_priors", None) + if model is None: + return _empty(False, None) + if not isinstance(model, PopulationGenotypeModel): + # Field is typed Optional[PopulationGenotypeModel] but Python won't enforce + # it; a garbage value must not crash the manifest (called from build()). + return _empty(True, False, + f"genotype_priors is not a PopulationGenotypeModel " + f"(got {type(model).__name__})") + + # A plane may have been attached directly (bypassing builder validation). + # Report validity rather than leaking non-JSON-clean numerics (e.g. NaN + # chromosome weights) into the manifest. + try: + model.validate( + chain_type=getattr(getattr(cfg, "metadata", None), "chain_type", None)) + valid, validation_error = True, None + except ValueError as exc: + valid, validation_error = False, str(exc) + + # content_checksum / float() can themselves raise on a malformed model; never + # let that crash the manifest — report None and the validity flag instead. + try: + checksum = model.content_checksum() + except Exception: + checksum = None + cw = None + if valid: + try: + cw = [float(model.chromosome_weights[0]), float(model.chromosome_weights[1])] + except Exception: + cw = None + freq = model.allele_frequencies if isinstance(model.allele_frequencies, dict) else {} + dele = model.haplotype_deletion_prob if isinstance(model.haplotype_deletion_prob, dict) else {} + novels = model.novel_alleles if isinstance(model.novel_alleles, (list, tuple)) else [] + block = { + "available": True, + "valid": valid, + "model_id": (model.model_id or None) if isinstance(model.model_id, str) else None, + "source": (model.source or None) if isinstance(model.source, str) else None, + "version": (model.version or None) if isinstance(model.version, str) else None, + "model_checksum": checksum, + "segments_with_frequencies": [s for s in ("V", "D", "J") if freq.get(s)], + "freq_gene_counts": {s: len(freq.get(s, {})) for s in ("V", "D", "J")}, + "deletion_gene_counts": {s: len(dele.get(s, {})) for s in ("V", "D", "J")}, + "novel_allele_count": len(novels), + "chromosome_weights": cw, + "source_field": "DataConfig.genotype_priors", + } + if validation_error is not None: + block["validation_error"] = validation_error + return block + + def _np_length_models_manifest_block(cfg): """Build the ``models.np_length_models`` manifest block per the NP Length Distribution Estimation v1 audit @@ -301,6 +380,15 @@ class DataConfig: # ``reference_rules``. reference_models: Optional[ReferenceEmpiricalModels] = None + # Donor-population germline prior plane (Slice — Cartridge genotype plane). + # ``None`` means no population prior; ``Genotype.sample(cfg)`` then falls + # back to a uniform synthetic prior. A non-``None`` plane is cartridge + # identity (folds into compute_checksum). See + # ``site_docs/guides/genotype.md`` ("Population genotype models on a + # cartridge"). Same soft-transition checksum policy as ``reference_rules`` / + # ``reference_models``. + genotype_priors: Optional[PopulationGenotypeModel] = None + def __getattr__(self, name): # Backward-compat shim for pickled DataConfigs missing post-v1 # fields. Note: schema_version / schema_sha256 fall through to @@ -319,6 +407,8 @@ def __getattr__(self, name): return None if name == 'reference_models': return None + if name == 'genotype_priors': + return None raise AttributeError(f"'{type(self).__name__}' object has no attribute '{name}'") def compute_checksum(self) -> str: @@ -355,6 +445,8 @@ def compute_checksum(self) -> str: pop_rules = 'reference_rules' in self.__dict__ and rr_value is None rm_value = self.__dict__.get('reference_models') pop_models = 'reference_models' in self.__dict__ and rm_value is None + gp_value = self.__dict__.get('genotype_priors') + pop_priors = 'genotype_priors' in self.__dict__ and gp_value is None self.schema_sha256 = "" if had_report: @@ -363,6 +455,8 @@ def compute_checksum(self) -> str: del self.__dict__['reference_rules'] if pop_models: del self.__dict__['reference_models'] + if pop_priors: + del self.__dict__['genotype_priors'] try: blob = pickle.dumps(self, protocol=4) return hashlib.sha256(blob).hexdigest() @@ -374,6 +468,8 @@ def compute_checksum(self) -> str: self.__dict__['reference_rules'] = rr_value if pop_models: self.__dict__['reference_models'] = rm_value + if pop_priors: + self.__dict__['genotype_priors'] = gp_value def verify_integrity(self) -> None: """Validate schema_version and schema_sha256. @@ -917,6 +1013,10 @@ def cartridge_manifest( # ``np_length_keys`` / ``legacy_np_lengths_present`` # entries above. "np_length_models": _np_length_models_manifest_block(self), + # Donor-population germline prior plane (Slice — Cartridge genotype + # plane). Read from the top-level ``DataConfig.genotype_priors`` + # field (NOT ``reference_models``); ``source_field`` records that. + "genotype_priors": _genotype_priors_manifest_block(self), } # Bridge once (or accept the provided refdata) to read the diff --git a/src/GenAIRR/genotype.py b/src/GenAIRR/genotype.py index 47d3cfc..a896126 100644 --- a/src/GenAIRR/genotype.py +++ b/src/GenAIRR/genotype.py @@ -36,6 +36,10 @@ def __init__(self, cfg, *, permissive: bool = False): self._cfg = cfg self._permissive = bool(permissive) self.subject_id: Optional[str] = None # plain attribute (read anywhere) + # Provenance of any population prior used to build this genotype. For + # builder-constructed genotypes every source is "manual" / None; for + # Genotype.sample(...) it records explicit/cartridge/uniform per input. + self.prior_provenance: Dict = self._manual_provenance() self._chromosome_weights: Tuple[float, float] = (0.5, 0.5) # segment -> gene -> [hap0 list[(allele, copies, weight)], hap1 ...] self._slots: Dict[str, Dict[str, List[List[Tuple[str, int, float]]]]] = { @@ -66,6 +70,17 @@ def with_subject(self, sid: str) -> "Genotype": self.subject_id = str(sid) return self + @staticmethod + def _manual_provenance() -> Dict: + return { + "allele_frequencies": "manual", + "haplotype_deletion_prob": "manual", + "chromosome_weights": "manual", + "novel_alleles": "manual", + "model_id": None, + "model_checksum": None, + } + @staticmethod def _check_chromosome_weights(w0, w1) -> Tuple[float, float]: for w in (w0, w1): @@ -391,6 +406,7 @@ def _snapshot(self) -> "Genotype": g._slots = _copy.deepcopy(self._slots) g._novel = _copy.deepcopy(self._novel) g._source_hash = self._source_hash + g.prior_provenance = _copy.deepcopy(self.prior_provenance) return g # ── population sampling ─────────────────────────────────────── @@ -544,7 +560,10 @@ def _resolve_allele_frequencies(cls, cfg, spec, segs): continue pairs: List[Tuple[str, float]] = [] total = 0.0 - for nm, w in supplied.items(): + # Sort by allele name so the draw is independent of the spec's + # dict insertion order — two content-equal models (same + # content_checksum) then produce identical draws at a given seed. + for nm, w in sorted(supplied.items()): if nm not in names: raise ValueError(f"{seg} gene {gene}: {nm!r} is not a known allele") if isinstance(w, bool) or not isinstance(w, (int, float)) or not math.isfinite(w) or w < 0: @@ -616,26 +635,30 @@ def sample( *, seed: int = 0, allele_frequencies=None, - haplotype_deletion_prob=0.0, + haplotype_deletion_prob=None, segments_to_sample=None, - chromosome_weights: Tuple[float, float] = (0.5, 0.5), + chromosome_weights: Optional[Tuple[float, float]] = None, subject_id: Optional[str] = None, ensure_viable: bool = True, max_resamples: int = 1000, + use_cartridge_priors: bool = True, + include_cartridge_novel_alleles="auto", ) -> "Genotype": """Sample a fully-specified diploid genotype from population priors. Independent per-gene, per-chromosome Hardy-Weinberg model: each gene on - each chromosome is independently deleted (prob - ``haplotype_deletion_prob``) or assigned one allele drawn from the gene's - allele frequencies. Homozygous/heterozygous/hemizygous/deleted states - emerge at Hardy-Weinberg rates. + each chromosome is independently deleted or assigned one allele drawn + from the gene's allele frequencies. Homozygous/heterozygous/hemizygous/ + deleted states emerge at Hardy-Weinberg rates. + + When ``cfg`` carries a ``genotype_priors`` plane and an argument is left + at its default, the plane supplies it (``use_cartridge_priors=False`` + disables ALL plane consumption — a clean uniform catalogue-only draw). + Each input is sourced independently and recorded in + ``g.prior_provenance`` (explicit / cartridge / uniform / default). This is NOT a population haplotype model — no linkage disequilibrium, gene - co-deletion blocks, ancestry, or donor-specific haplotype structure. It - samples catalogue alleles only (no novel alleles) and deletion only (no - copy-number duplication). The default prior is uniform within each gene; - supply ``allele_frequencies`` for realistic per-gene frequencies. + co-deletion blocks, ancestry, or donor-specific haplotype structure. With ``ensure_viable=True`` (default), the draw is repeated (with a deterministic sub-seed) up to ``max_resamples`` times until at least one @@ -651,10 +674,71 @@ def sample( if isinstance(max_resamples, bool) or not isinstance(max_resamples, int) or max_resamples < 1: raise ValueError(f"max_resamples must be an int >= 1, got {max_resamples!r}") - cw = cls._check_chromosome_weights(*chromosome_weights) + if not isinstance(use_cartridge_priors, bool): + raise ValueError( + f"use_cartridge_priors must be a bool, got {use_cartridge_priors!r}") + # Identity checks (not `in (...)`): Python's `1 == True` / `0 == False` + # would otherwise let integers slip through and silently act like False. + _icna = include_cartridge_novel_alleles + if not (_icna is True or _icna is False or _icna == "auto"): + raise ValueError( + "include_cartridge_novel_alleles must be 'auto', True, or False, " + f"got {include_cartridge_novel_alleles!r}") + + plane = cls._resolve_plane(cfg, use_cartridge_priors) + if plane is not None: + # A plane attached via the builder is already validated, but a plane + # set directly on the DataConfig bypasses that — validate before use so + # a malformed prior (empty model_id, NaN weights, D-on-VJ) fails loudly + # rather than being silently sampled and stamped into provenance. + plane.validate(chain_type=getattr(getattr(cfg, "metadata", None), "chain_type", None)) + + # Per-input source resolution. + if allele_frequencies is not None: + freq_spec, freq_src = allele_frequencies, "explicit" + elif plane is not None and plane.allele_frequencies: + freq_spec, freq_src = plane.allele_frequencies, "cartridge" + else: + freq_spec, freq_src = None, "uniform" + + if haplotype_deletion_prob is not None: + del_spec, del_src = haplotype_deletion_prob, "explicit" + elif plane is not None and plane.haplotype_deletion_prob: + del_spec, del_src = plane.haplotype_deletion_prob, "cartridge" + else: + del_spec, del_src = 0.0, "default" + + if chromosome_weights is not None: + cw_in, cw_src = chromosome_weights, "explicit" + elif plane is not None: + cw_in, cw_src = plane.chromosome_weights, "cartridge" + else: + cw_in, cw_src = (0.5, 0.5), "default" + + cw = cls._check_chromosome_weights(*cw_in) segs = cls._resolve_sample_segments(cfg, segments_to_sample) - freqs = cls._resolve_allele_frequencies(cfg, allele_frequencies, segs) - delp = cls._resolve_haplotype_deletion(cfg, haplotype_deletion_prob, segs) + + # Candidate-novel injection. Plane novels are draw CANDIDATES (injected + # into the sampling cfg) only when their source is active; the returned + # genotype still exports only CARRIED novels via effective_dataconfig(). + inject_novels = False + if plane is not None and plane.novel_alleles: + if include_cartridge_novel_alleles is True: + inject_novels = True + elif include_cartridge_novel_alleles == "auto": + inject_novels = freq_src in ("cartridge", "uniform") + # False -> never + novel_src = "none" + sampling_cfg = cfg + novel_helper = None + if inject_novels: + novel_helper = cls._register_plane_novels(cfg, plane) + sampling_cfg = cls._dataconfig_injecting_all_novels(novel_helper) + freq_spec = cls._augment_freqs_with_novels(sampling_cfg, freq_spec, plane, segs) + novel_src = "cartridge" + + freqs = cls._resolve_allele_frequencies(sampling_cfg, freq_spec, segs) + delp = cls._resolve_haplotype_deletion(sampling_cfg, del_spec, segs) # Compute the cartridge content hash ONCE (it rebuilds refdata + hashes); # reuse it across all draws instead of recomputing per attempt. source_hash = cfg.cartridge_manifest()["hashes"]["refdata_content_hash"] @@ -662,11 +746,22 @@ def sample( # `seed` cannot collide with a direct draw at `seed + 1`. base_rng = random.Random(seed) + provenance = { + "allele_frequencies": freq_src, + "haplotype_deletion_prob": del_src, + "chromosome_weights": cw_src, + "novel_alleles": novel_src, + "model_id": plane.model_id if plane is not None else None, + "model_checksum": plane.content_checksum() if plane is not None else None, + } + attempts = max_resamples if ensure_viable else 1 for _attempt in range(attempts): sub_seed = base_rng.getrandbits(63) - g = cls._draw_one(cfg, sub_seed, segs, freqs, delp, cw, subject_id, source_hash) + g = cls._draw_one(sampling_cfg, sub_seed, segs, freqs, delp, cw, subject_id, source_hash) if not ensure_viable or g._is_viable(cfg, cw): + cls._rebind_to_base(g, cfg, novel_helper) + g.prior_provenance = dict(provenance) return g raise ValueError( f"could not sample a viable genotype after {max_resamples} attempts; " @@ -675,6 +770,82 @@ def sample( f"(chromosome_weights={cw})" ) + @staticmethod + def _resolve_plane(cfg, use_cartridge_priors): + if not use_cartridge_priors: + return None + return getattr(cfg, "genotype_priors", None) + + @classmethod + def _register_plane_novels(cls, cfg, plane): + """Build a throwaway Genotype carrying every plane novel as a registered + novel allele. This is where catalogue-aware + functional validation of + plane novels happens (via add_novel_allele). Returns the helper.""" + helper = cls.from_dataconfig(cfg) + for nv in plane.novel_alleles: + helper.add_novel_allele( + nv.name, base=nv.base_allele, sequence=nv.sequence.upper(), + segment=nv.segment, allow_nonfunctional=nv.allow_nonfunctional) + return helper + + @staticmethod + def _dataconfig_injecting_all_novels(helper): + """A cfg copy with ALL of the helper's registered novels injected as + catalogue alleles — the *sampling* reference (candidates), distinct from + effective_dataconfig() which injects only carried novels.""" + import copy as _copy + cfg = _copy.deepcopy(helper._cfg) + by_seg = {"V": cfg.v_alleles, "D": cfg.d_alleles, "J": cfg.j_alleles} + for name, info in helper._novel.items(): + d = by_seg[info["segment"]] + existing = list(d.get(info["gene"], [])) + existing.append(_copy.deepcopy(info["allele"])) + d[info["gene"]] = existing + return cfg + + @classmethod + def _augment_freqs_with_novels(cls, sampling_cfg, freq_spec, plane, segs): + """Return a nested freq spec that includes each plane novel in its gene's + table per the synthesis rule: authored table -> preserve + add novel; + no table -> catalogue alleles 1.0 + novel frequency. Collisions raise.""" + if freq_spec == "usage_as_prior": + nested = cls._usage_frequencies(sampling_cfg, segs) + nested = {seg: {g: dict(al) for g, al in genes.items()} + for seg, genes in nested.items()} + else: + nested = {seg: {g: dict(al) for g, al in genes.items()} + for seg, genes in cls._normalize_freq_spec(sampling_cfg, freq_spec, segs).items()} + for nv in plane.novel_alleles: + seg = nv.segment + if seg not in segs: + continue + gene = nv.name.split("*")[0] + gene_tbl = nested.setdefault(seg, {}).get(gene) + if gene_tbl is None: + # no authored table for this gene: catalogue 1.0 + novel + catalogue = _alleles_by_gene(sampling_cfg, seg).get(gene, []) + gene_tbl = {a.name: 1.0 for a in catalogue if a.name != nv.name} + nested[seg][gene] = gene_tbl + if nv.name in gene_tbl: + raise ValueError( + f"plane novel {nv.name!r} collides with an existing allele weight") + gene_tbl[nv.name] = float(nv.frequency) + return nested + + @staticmethod + def _rebind_to_base(g, base_cfg, novel_helper): + """Point a drawn genotype back at the base cfg and register only the + novels it actually carries, so effective_dataconfig() injects carried + novels (and nothing else). Task 7 supplies ``novel_helper``.""" + import copy as _copy + g._cfg = base_cfg + if novel_helper is None: + return + carried = g._carried_allele_names() + for name, info in novel_helper._novel.items(): + if name in carried: + g._novel[name] = _copy.deepcopy(info) + @classmethod def _draw_one(cls, cfg, seed, segs, freqs, delp, cw, subject_id, source_hash): import random @@ -690,6 +861,7 @@ def _draw_one(cls, cfg, seed, segs, freqs, delp, cw, subject_id, source_hash): g._slots = {s: {} for s in _SEGMENTS} g._novel = {} g._source_hash = source_hash + g.prior_provenance = cls._manual_provenance() # sample() overwrites with real sources for seg in segs: for gene in _alleles_by_gene(cfg, seg): pdel = delp[seg][gene] @@ -744,6 +916,27 @@ def _zygosity(h0: List, h1: List) -> str: return "homozygous" return "heterozygous" + def to_metadata(self) -> Dict: + """Flat genotype-level metadata (subject + prior provenance + refdata + hashes) for sidecar export / benchmark tooling. Distinct from + ``to_table`` (one row per gene), which carries no provenance. + + ``source_refdata_hash`` is the base cartridge's content hash; + ``effective_refdata_hash`` is the hash the engine actually runs against — + equal to the source hash unless carried novel alleles are injected, in + which case it reflects the effective (novel-augmented) reference.""" + import copy as _copy + effective = self._source_hash + if self.has_novel(): + effective = self.effective_dataconfig().cartridge_manifest()[ + "hashes"]["refdata_content_hash"] + return { + "subject_id": self.subject_id, + "source_refdata_hash": self._source_hash, + "effective_refdata_hash": effective, + "prior_provenance": _copy.deepcopy(self.prior_provenance), + } + def to_table(self) -> List[Dict]: """One row per (segment, gene) with full diploid truth: zygosity (incl. ``hemizygous`` / ``deleted``), the carried alleles per diff --git a/src/GenAIRR/genotype_priors.py b/src/GenAIRR/genotype_priors.py new file mode 100644 index 0000000..8dfdb73 --- /dev/null +++ b/src/GenAIRR/genotype_priors.py @@ -0,0 +1,434 @@ +"""Population genotype model — a donor-population germline prior carried as a +top-level cartridge plane (``DataConfig.genotype_priors``). + +This is NOT an empirical recombination model (those live on +``ReferenceEmpiricalModels``); it is a per-gene carriage/deletion prior plus a +catalogue of population novel/private alleles. ``Genotype.sample`` consumes it to +draw a per-individual diploid genotype. See +``.private/specs/2026-06-17-genotype-cartridge-plane-design.md``. + +The plane stays decoupled from a specific catalogue: ``validate()`` does shape + +numeric/DNA sanity only. Catalogue-aware checks (gene/allele existence, novel +functional validation, viability) run at attach time +(``ReferenceCartridgeBuilder.set_genotype_priors``) and at sample time +(``Genotype.sample``). +""" +from __future__ import annotations + +import hashlib +import json +import math +from dataclasses import dataclass, field +from typing import Dict, List, Optional, Tuple + +_SEGMENTS = ("V", "D", "J") +_DNA = set("ACGT") +SCHEMA_TAG = "population_genotype_model/1" + + +@dataclass +class PopulationNovelAllele: + """A population novel/private allele carried on a cartridge prior. + + ``name``'s gene token (text before ``"*"``) must equal ``base_allele``'s + gene (validated catalogue-aware); ``sequence`` is substitution-only (same + length as the base, validated catalogue-aware). ``frequency`` is the + population weight the allele competes with inside its base gene's draw. + """ + + name: str + segment: str + base_allele: str + sequence: str + frequency: float + allow_nonfunctional: bool = False + + +def _check_weight(w, where, *, strict_positive): + if isinstance(w, bool) or not isinstance(w, (int, float)): + raise ValueError(f"{where}: weight must be a finite number, got {type(w).__name__}") + wf = float(w) + if not math.isfinite(wf): + raise ValueError(f"{where}: weight must be finite, got {w!r}") + if strict_positive and wf <= 0.0: + raise ValueError(f"{where}: weight must be > 0, got {wf}") + if not strict_positive and wf < 0.0: + raise ValueError(f"{where}: weight must be >= 0, got {wf}") + return wf + + +@dataclass +class PopulationGenotypeModel: + """A donor-population germline prior (see module docstring).""" + + allele_frequencies: Dict[str, Dict[str, Dict[str, float]]] = field(default_factory=dict) + haplotype_deletion_prob: Dict[str, Dict[str, float]] = field(default_factory=dict) + chromosome_weights: Tuple[float, float] = (0.5, 0.5) + novel_alleles: List[PopulationNovelAllele] = field(default_factory=list) + model_id: str = "" + source: str = "" + description: str = "" + version: str = "" + + def validate(self, chain_type=None, *, name: str = "genotype_priors") -> None: + """Catalogue-free shape + numeric/DNA sanity. Raises ``ValueError`` on + the first violation. ``chain_type`` rejects any D entry on a VJ chain; + it accepts ``"vj"`` / ``"vdj"`` strings or a ``ChainType``-like object + exposing ``has_d``.""" + ct = self._normalize_chain_type(chain_type) + if not isinstance(self.model_id, str) or not self.model_id: + raise ValueError(f"{name}.model_id must be a non-empty string") + if not isinstance(self.source, str) or not self.source: + raise ValueError(f"{name}.source must be a non-empty string") + for label in ("description", "version"): + if not isinstance(getattr(self, label), str): + raise ValueError(f"{name}.{label} must be a string") + + # allele_frequencies: seg -> gene -> {allele: weight >= 0}, >=1 positive + if not isinstance(self.allele_frequencies, dict): + raise ValueError(f"{name}.allele_frequencies must be a dict") + for seg, genes in self.allele_frequencies.items(): + self._check_segment(seg, ct, f"{name}.allele_frequencies") + if not isinstance(genes, dict): + raise ValueError(f"{name}.allele_frequencies[{seg!r}] must be a dict") + for gene, alleles in genes.items(): + if not isinstance(gene, str) or not gene: + raise ValueError( + f"{name}.allele_frequencies[{seg!r}]: gene names must be " + f"non-empty strings, got {gene!r}") + if not isinstance(alleles, dict) or not alleles: + raise ValueError( + f"{name}.allele_frequencies[{seg!r}][{gene!r}] must be a " + f"non-empty mapping of allele -> weight") + positive = 0 + for allele, w in alleles.items(): + if not isinstance(allele, str) or not allele: + raise ValueError( + f"{name}.allele_frequencies[{seg!r}][{gene!r}]: allele names " + f"must be non-empty strings, got {allele!r}") + wf = _check_weight(w, f"{name}.allele_frequencies[{seg}][{gene}][{allele}]", + strict_positive=False) + if wf > 0: + positive += 1 + if positive == 0: + raise ValueError( + f"{name}.allele_frequencies[{seg!r}][{gene!r}]: at least one " + f"allele weight must be > 0") + + # haplotype_deletion_prob: seg -> gene -> prob in [0, 1] + if not isinstance(self.haplotype_deletion_prob, dict): + raise ValueError(f"{name}.haplotype_deletion_prob must be a dict") + for seg, genes in self.haplotype_deletion_prob.items(): + self._check_segment(seg, ct, f"{name}.haplotype_deletion_prob") + if not isinstance(genes, dict): + raise ValueError(f"{name}.haplotype_deletion_prob[{seg!r}] must be a dict") + for gene, p in genes.items(): + if not isinstance(gene, str) or not gene: + raise ValueError( + f"{name}.haplotype_deletion_prob[{seg!r}]: gene names must be " + f"non-empty strings, got {gene!r}") + if isinstance(p, bool) or not isinstance(p, (int, float)) or not math.isfinite(p): + raise ValueError( + f"{name}.haplotype_deletion_prob[{seg}][{gene}] must be finite, got {p!r}") + if not (0.0 <= float(p) <= 1.0): + raise ValueError( + f"{name}.haplotype_deletion_prob[{seg}][{gene}] must be in [0, 1], got {p}") + + # chromosome_weights: 2-tuple, finite, non-negative, >=1 positive + cw = self.chromosome_weights + if not (isinstance(cw, (tuple, list)) and len(cw) == 2): + raise ValueError(f"{name}.chromosome_weights must be a 2-tuple, got {cw!r}") + for w in cw: + if isinstance(w, bool) or not isinstance(w, (int, float)) or not math.isfinite(w): + raise ValueError(f"{name}.chromosome_weights must be finite numbers, got {cw!r}") + if cw[0] < 0 or cw[1] < 0 or (cw[0] + cw[1]) <= 0: + raise ValueError( + f"{name}.chromosome_weights must be non-negative and sum>0, got {tuple(cw)}") + + # novel_alleles + seen = set() + for nv in self.novel_alleles: + if not isinstance(nv, PopulationNovelAllele): + raise ValueError(f"{name}.novel_alleles entries must be PopulationNovelAllele") + for attr in ("name", "base_allele", "sequence"): + v = getattr(nv, attr) + if not isinstance(v, str) or not v: + raise ValueError(f"{name}.novel_alleles: {attr} must be a non-empty string") + self._check_segment(nv.segment, ct, f"{name}.novel_alleles[{nv.name!r}]") + if "*" not in nv.name or not nv.name.split("*")[0]: + raise ValueError( + f"{name}.novel_alleles: name {nv.name!r} must contain a gene token before '*'") + if any(b not in _DNA for b in nv.sequence.upper()): + raise ValueError( + f"{name}.novel_alleles[{nv.name!r}]: sequence must be DNA (A/C/G/T only)") + _check_weight(nv.frequency, f"{name}.novel_alleles[{nv.name!r}].frequency", + strict_positive=True) + if not isinstance(nv.allow_nonfunctional, bool): + raise ValueError( + f"{name}.novel_alleles[{nv.name!r}].allow_nonfunctional must be a " + f"bool, got {type(nv.allow_nonfunctional).__name__}") + if nv.name in seen: + raise ValueError(f"{name}.novel_alleles: duplicate novel name {nv.name!r}") + seen.add(nv.name) + + @classmethod + def from_genotypes(cls, genotypes, *, cfg=None, pseudocount=0.0, + include_novel=True, min_subjects=None, + subject_id_policy="require_unique", segments=None, + model_id="estimated", source="from_genotypes", + description="", version="") -> "PopulationGenotypeModel": + """Estimate a population genotype model from observed ``Genotype`` objects. + + Allele frequencies are counted **per carried chromosome** (homozygous=2, + hemizygous=1, deleted=0); ``pseudocount`` is added to every catalogue + allele per gene (and to a novel only when ``include_novel`` and the novel + belongs to the gene). Deletion probability is per haplotype: + ``deleted_haplotypes / (2 * n_subjects)`` (NO pseudocount). Genotypes + carrying a duplicated gene (copy-number > 1) are rejected — the plane is + deletion-only. ``subject_id_policy`` ('require_unique' | 'allow_duplicates') + guards against double-counting; ``min_subjects`` guards tiny estimates. + """ + gts = list(genotypes) + if subject_id_policy not in ("require_unique", "allow_duplicates"): + raise ValueError( + f"subject_id_policy must be 'require_unique' or 'allow_duplicates', " + f"got {subject_id_policy!r}") + if subject_id_policy == "require_unique": + ids = [g.subject_id for g in gts] + none_count = sum(1 for i in ids if i is None) + if 0 < none_count < len(ids): + raise ValueError( + "subject_id_policy='require_unique': some genotypes have a " + "subject_id and others don't") + if none_count == len(ids) and len(ids) > 1: + # all-None: uniqueness cannot be verified, so double-counting + # cannot be ruled out — fail loudly rather than silently allow it. + raise ValueError( + "subject_id_policy='require_unique': none of the genotypes have a " + "subject_id, so uniqueness cannot be verified; set with_subject(...) " + "or pass subject_id_policy='allow_duplicates'") + present = [i for i in ids if i is not None] + if len(present) != len(set(present)): + raise ValueError( + "subject_id_policy='require_unique': duplicate subject_id found") + n = len(gts) + if n == 0: + raise ValueError("from_genotypes: need at least one genotype") + if min_subjects is not None: + if isinstance(min_subjects, bool) or not isinstance(min_subjects, int) or min_subjects < 1: + raise ValueError( + f"from_genotypes: min_subjects must be an int >= 1 (or None), " + f"got {min_subjects!r}") + if n < min_subjects: + raise ValueError( + f"from_genotypes: min_subjects={min_subjects} but only {n} given") + if (isinstance(pseudocount, bool) or not isinstance(pseudocount, (int, float)) + or not math.isfinite(pseudocount) or pseudocount < 0): + raise ValueError( + f"from_genotypes: pseudocount must be a finite number >= 0, " + f"got {pseudocount!r}") + if not isinstance(include_novel, bool): + raise ValueError( + f"from_genotypes: include_novel must be a bool, got {include_novel!r}") + + if cfg is None: + cfg = gts[0]._cfg + seg_list = cls._resolve_estimate_segments(cfg, segments) + + by_seg = {"V": cfg.v_alleles, "D": cfg.d_alleles, "J": cfg.j_alleles} + + # reject duplicated genes (copy-number > 1 on any slot) — both multi-entry + # slots AND a single entry whose copy_count > 1. The plane is deletion-only. + for g in gts: + for seg in _SEGMENTS: + for gene, haps in g._slots[seg].items(): + for hap in haps: + if len(hap) > 1 or any(c > 1 for (_a, c, _w) in hap): + raise ValueError( + f"from_genotypes: genotype {g.subject_id!r} carries a " + f"duplicated gene ({seg} {gene}, copy-number > 1); the " + f"population plane is deletion-only — remove duplications " + f"or collapse them before estimating") + + # Require completeness: a gene absent from _slots is UNSPECIFIED (unknown), + # NOT deleted. Counting absence as deletion silently inflates p_del, so we + # require every estimated-segment catalogue gene to be specified (call + # complete_from_reference() to fill, or delete_gene() to mark absence). + for g in gts: + for seg in seg_list: + for gene in (by_seg[seg] or {}): + if gene not in g._slots[seg]: + raise ValueError( + f"from_genotypes: genotype {g.subject_id!r} does not specify " + f"{seg} gene {gene!r}; an unspecified gene is unknown, not " + f"deleted. Call complete_from_reference() (or delete_gene to " + f"mark it absent) before estimating") + + freqs: Dict[str, Dict[str, Dict[str, float]]] = {} + dele: Dict[str, Dict[str, float]] = {} + novel_freq: Dict[str, float] = {} + novel_spec: Dict[str, "PopulationNovelAllele"] = {} + + for seg in seg_list: + catalogue = by_seg[seg] or {} + freqs[seg] = {} + dele[seg] = {} + for gene, alleles in catalogue.items(): + counts = {a.name: 0.0 for a in alleles} + deleted_haps = 0 + for g in gts: + haps = g._slots[seg].get(gene, [[], []]) + for hap in haps: + names = {a for (a, _c, _w) in hap} + if not names: + deleted_haps += 1 + continue + for nm in names: + if nm in counts: + counts[nm] += 1.0 + continue + # Not a catalogue allele of this gene: it must be a + # registered novel of THIS genotype matching gene + + # segment, otherwise it is unknown/corrupt and we must + # not silently drop it. + info = getattr(g, "_novel", {}).get(nm) + if (info is None or info.get("gene") != gene + or info.get("segment") != seg): + raise ValueError( + f"from_genotypes: genotype {g.subject_id!r} carries " + f"allele {nm!r} in {seg} {gene} that is neither a " + f"catalogue allele nor a registered novel of that " + f"gene/segment") + if include_novel: + novel_freq[nm] = novel_freq.get(nm, 0.0) + 1.0 + # include_novel=False -> registered novel intentionally excluded + if pseudocount: + for nm in counts: + counts[nm] += float(pseudocount) + kept = {nm: c for nm, c in counts.items() if c > 0} + if kept: + freqs[seg][gene] = kept + dele[seg][gene] = deleted_haps / (2.0 * n) + + novels: List[PopulationNovelAllele] = [] + if include_novel: + for g in gts: + for name, info in getattr(g, "_novel", {}).items(): + if name in novel_freq and name not in novel_spec: + novel_spec[name] = PopulationNovelAllele( + name=name, segment=info["segment"], + base_allele=info["base"], + sequence=info["allele"].ungapped_seq.upper(), + frequency=novel_freq[name], + allow_nonfunctional=not info.get("functional", True)) + # Novel frequency is carried on the PopulationNovelAllele itself, NOT + # duplicated into allele_frequencies — putting a non-catalogue name + # there would (a) fail set_genotype_priors' catalogue-aware check and + # (b) collide with sample-time novel synthesis. sample() combines the + # catalogue table with novel frequencies via _augment_freqs_with_novels. + novels.extend(novel_spec.values()) + + model = cls(allele_frequencies=freqs, haplotype_deletion_prob=dele, + chromosome_weights=(0.5, 0.5), novel_alleles=novels, + model_id=model_id, source=source, description=description, + version=version) + # Required-field + shape validation still applies to estimator output + # (empty model_id/source, non-string description/version, etc.). + model.validate(chain_type=getattr(getattr(cfg, "metadata", None), "chain_type", None)) + return model + + @staticmethod + def _segments_for(cfg): + segs = ["V"] + if getattr(cfg, "d_alleles", None): + segs.append("D") + segs.append("J") + return segs + + @classmethod + def _resolve_estimate_segments(cls, cfg, segments): + """Validate / canonicalize the ``segments`` argument: a non-empty + iterable of unique V/D/J labels each present in the cartridge. A bare + string is rejected (it would iterate as characters).""" + if segments is None: + return cls._segments_for(cfg) + if isinstance(segments, str): + raise ValueError( + f"from_genotypes: segments must be a list/tuple of segment labels, " + f"not a string ({segments!r})") + seq = list(segments) + if not seq: + raise ValueError("from_genotypes: segments must be a non-empty iterable") + by_seg = {"V": cfg.v_alleles, "D": cfg.d_alleles, "J": cfg.j_alleles} + seen = [] + for s in seq: + if s not in _SEGMENTS: + raise ValueError( + f"from_genotypes: unknown segment {s!r}; expected one of {_SEGMENTS}") + if s in seen: + raise ValueError(f"from_genotypes: duplicate segment {s!r}") + if not by_seg[s]: + raise ValueError(f"from_genotypes: cartridge has no {s} segment") + seen.append(s) + # canonical V/D/J order + return [s for s in _SEGMENTS if s in seen] + + @staticmethod + def _normalize_chain_type(chain_type): + """Return ``"vj"`` / ``"vdj"`` / ``None`` from a string or a + ``ChainType``-like object (one exposing ``has_d``).""" + if chain_type is None: + return None + if isinstance(chain_type, str): + return chain_type.lower() + has_d = getattr(chain_type, "has_d", None) + if has_d is None: + return None + return "vdj" if has_d else "vj" + + @staticmethod + def _check_segment(seg, ct, where): + if seg not in _SEGMENTS: + raise ValueError(f"{where}: segment {seg!r} must be one of {_SEGMENTS}") + if ct == "vj" and seg == "D": + raise ValueError( + f"{where}: D-segment prior on a VJ chain is meaningless (no D pool). " + f"Drop the D entries or use a VDJ cartridge.") + + def content_checksum(self) -> str: + """Canonical sha256 of the plane's semantic content — stable across dict + insertion order, int-vs-float spelling, and DNA case. This (not the + pickle-based DataConfig checksum) is what manifest/provenance report.""" + + def _num(x): + return repr(float(x)) + + def _freqs(d): + return {seg: {gene: {al: _num(w) for al, w in sorted(alleles.items())} + for gene, alleles in sorted(genes.items())} + for seg, genes in sorted(d.items())} + + def _del(d): + return {seg: {gene: _num(p) for gene, p in sorted(genes.items())} + for seg, genes in sorted(d.items())} + + novels = sorted( + ({"name": nv.name, "segment": nv.segment, "base_allele": nv.base_allele, + "sequence": nv.sequence.upper(), "frequency": _num(nv.frequency), + "allow_nonfunctional": bool(nv.allow_nonfunctional)} + for nv in self.novel_alleles), + key=lambda r: r["name"], + ) + canon = { + "schema": SCHEMA_TAG, + "model_id": self.model_id, + "source": self.source, + "description": self.description, + "version": self.version, + "allele_frequencies": _freqs(self.allele_frequencies), + "haplotype_deletion_prob": _del(self.haplotype_deletion_prob), + "chromosome_weights": [_num(self.chromosome_weights[0]), _num(self.chromosome_weights[1])], + "novel_alleles": novels, + } + blob = json.dumps(canon, sort_keys=True, separators=(",", ":")).encode("utf-8") + return hashlib.sha256(blob).hexdigest() diff --git a/tests/test_genotype_cartridge_plane.py b/tests/test_genotype_cartridge_plane.py new file mode 100644 index 0000000..410078d --- /dev/null +++ b/tests/test_genotype_cartridge_plane.py @@ -0,0 +1,686 @@ +"""Cartridge genotype plane: PopulationGenotypeModel + Genotype.sample consumption.""" +import pickle + +import pytest + +import GenAIRR as ga +import GenAIRR.data as gdata +from GenAIRR.genotype import Genotype +from GenAIRR.genotype_priors import PopulationGenotypeModel, PopulationNovelAllele + + +def _cfg(): + return gdata.HUMAN_IGH_OGRDB + + +def _vg_two_alleles(cfg): + """A V gene with >= 2 alleles, plus its first two allele names.""" + vg = next(g for g, al in cfg.v_alleles.items() if len(al) >= 2) + a0, a1 = (a.name for a in cfg.v_alleles[vg][:2]) + return vg, a0, a1 + + +def test_model_requires_identity(): + m = PopulationGenotypeModel(model_id="", source="x") + with pytest.raises(ValueError, match="model_id"): + m.validate() + m = PopulationGenotypeModel(model_id="x", source="") + with pytest.raises(ValueError, match="source"): + m.validate() + + +def test_model_validates_shapes(): + cfg = _cfg() + vg, a0, a1 = _vg_two_alleles(cfg) + ok = PopulationGenotypeModel( + model_id="toy", source="unit-test", + allele_frequencies={"V": {vg: {a0: 2.0, a1: 1.0}}}, + haplotype_deletion_prob={"V": {vg: 0.1}}, + chromosome_weights=(0.5, 0.5), + ) + ok.validate(chain_type="vdj") # no raise + + with pytest.raises(ValueError, match="finite"): + PopulationGenotypeModel(model_id="t", source="s", + allele_frequencies={"V": {vg: {a0: float("nan")}}}).validate() + with pytest.raises(ValueError, match=r"\[0, 1\]"): + PopulationGenotypeModel(model_id="t", source="s", + haplotype_deletion_prob={"V": {vg: 1.5}}).validate() + with pytest.raises(ValueError, match="non-negative"): + PopulationGenotypeModel(model_id="t", source="s", + chromosome_weights=(-1.0, 1.0)).validate() + with pytest.raises(ValueError, match="at least one"): + PopulationGenotypeModel(model_id="t", source="s", + allele_frequencies={"V": {vg: {a0: 0.0, a1: 0.0}}}).validate() + + +def test_novel_shape_validation(): + base = PopulationNovelAllele(name="IGHV1-2*99", segment="V", + base_allele="IGHV1-2*02", sequence="ACGT", frequency=1.0) + PopulationGenotypeModel(model_id="t", source="s", novel_alleles=[base]).validate() + + with pytest.raises(ValueError, match="must be > 0"): + PopulationGenotypeModel(model_id="t", source="s", novel_alleles=[ + PopulationNovelAllele(name="IGHV1-2*99", segment="V", + base_allele="IGHV1-2*02", sequence="ACGT", frequency=0.0)]).validate() + with pytest.raises(ValueError, match="A/C/G/T|DNA"): + PopulationGenotypeModel(model_id="t", source="s", novel_alleles=[ + PopulationNovelAllele(name="IGHV1-2*99", segment="V", + base_allele="IGHV1-2*02", sequence="ACGX", frequency=1.0)]).validate() + with pytest.raises(ValueError, match="duplicate|unique"): + PopulationGenotypeModel(model_id="t", source="s", novel_alleles=[ + PopulationNovelAllele(name="IGHV1-2*99", segment="V", + base_allele="IGHV1-2*02", sequence="ACGT", frequency=1.0), + PopulationNovelAllele(name="IGHV1-2*99", segment="V", + base_allele="IGHV1-2*02", sequence="ACGA", frequency=1.0)]).validate() + + +def test_d_on_vj_rejected(): + m = PopulationGenotypeModel(model_id="t", source="s", + allele_frequencies={"D": {"IGHD1-1": {"IGHD1-1*01": 1.0}}}) + with pytest.raises(ValueError, match="VJ|D-segment|D segment"): + m.validate(chain_type="vj") + + +def test_d_on_vj_rejected_via_chaintype_object(): + class _FakeChainType: # ChainType-like: exposes has_d + has_d = False + m = PopulationGenotypeModel(model_id="t", source="s", + haplotype_deletion_prob={"D": {"IGHD1-1": 0.5}}) + with pytest.raises(ValueError, match="VJ|D-segment|D segment"): + m.validate(chain_type=_FakeChainType()) + + +def test_gene_keys_must_be_strings(): + with pytest.raises(ValueError, match="gene"): + PopulationGenotypeModel(model_id="t", source="s", + allele_frequencies={"V": {123: {"IGHV1-2*02": 1.0}}}).validate() + with pytest.raises(ValueError, match="gene"): + PopulationGenotypeModel(model_id="t", source="s", + haplotype_deletion_prob={"V": {123: 0.1}}).validate() + + +def test_allow_nonfunctional_must_be_bool(): + with pytest.raises(ValueError, match="allow_nonfunctional"): + PopulationGenotypeModel(model_id="t", source="s", novel_alleles=[ + PopulationNovelAllele(name="IGHV1-2*99", segment="V", + base_allele="IGHV1-2*02", sequence="ACGT", frequency=1.0, + allow_nonfunctional="yes")]).validate() + + +def test_content_checksum_is_canonical(): + vg = "IGHV1-2" + m1 = PopulationGenotypeModel(model_id="m", source="s", + allele_frequencies={"V": {vg: {"IGHV1-2*02": 2, "IGHV1-2*04": 1}}}, + novel_alleles=[PopulationNovelAllele(name="IGHV1-2*99", segment="V", + base_allele="IGHV1-2*02", sequence="acgt", frequency=1.0)]) + # same content, different dict insertion order + int-vs-float + DNA case + m2 = PopulationGenotypeModel(model_id="m", source="s", + allele_frequencies={"V": {vg: {"IGHV1-2*04": 1.0, "IGHV1-2*02": 2.0}}}, + novel_alleles=[PopulationNovelAllele(name="IGHV1-2*99", segment="V", + base_allele="IGHV1-2*02", sequence="ACGT", frequency=1.0)]) + assert m1.content_checksum() == m2.content_checksum() + + m3 = PopulationGenotypeModel(model_id="m", source="s", + allele_frequencies={"V": {vg: {"IGHV1-2*02": 3.0, "IGHV1-2*04": 1.0}}}) + assert m3.content_checksum() != m1.content_checksum() + + +def test_genotype_priors_field_default_and_checksum_invariant(): + cfg = _cfg() + # default-new and bundled both have no plane + assert getattr(cfg, "genotype_priors", "MISSING") is None + base_checksum = cfg.compute_checksum() + + import copy + cfg2 = copy.deepcopy(cfg) + cfg2.genotype_priors = None # explicitly None must not change the checksum + assert cfg2.compute_checksum() == base_checksum + + cfg3 = copy.deepcopy(cfg) + cfg3.genotype_priors = PopulationGenotypeModel(model_id="m", source="s") + assert cfg3.compute_checksum() != base_checksum # a real plane is cartridge identity + + +def test_genotype_priors_pickle_round_trip(): + import copy + cfg = copy.deepcopy(_cfg()) + m = PopulationGenotypeModel(model_id="m", source="s", + haplotype_deletion_prob={"V": {next(iter(cfg.v_alleles)): 0.2}}) + cfg.genotype_priors = m + back = pickle.loads(pickle.dumps(cfg, protocol=4)) + assert back.genotype_priors.model_id == "m" + assert back.genotype_priors.content_checksum() == m.content_checksum() + + +def test_manifest_genotype_priors_block(): + import copy + cfg = copy.deepcopy(_cfg()) + block = cfg.cartridge_manifest()["models"]["genotype_priors"] + assert block["available"] is False + + vg = next(iter(cfg.v_alleles)) + cfg.genotype_priors = PopulationGenotypeModel( + model_id="m1", source="VDJbase-toy", version="1", + allele_frequencies={"V": {vg: {cfg.v_alleles[vg][0].name: 1.0}}}, + haplotype_deletion_prob={"V": {vg: 0.1}}, + novel_alleles=[PopulationNovelAllele(name="IGHV1-2*99", segment="V", + base_allele="IGHV1-2*02", sequence="ACGT", frequency=1.0)], + ) + block = cfg.cartridge_manifest()["models"]["genotype_priors"] + assert block["available"] is True + assert block["model_id"] == "m1" + assert block["source"] == "VDJbase-toy" + assert block["model_checksum"] == cfg.genotype_priors.content_checksum() + assert block["freq_gene_counts"]["V"] == 1 + assert block["deletion_gene_counts"]["V"] == 1 + assert block["novel_allele_count"] == 1 + assert block["chromosome_weights"] == [0.5, 0.5] + assert block["source_field"] == "DataConfig.genotype_priors" + + +# ── Task 5: provenance scaffolding ─────────────────────────────── + + +def test_prior_provenance_default_and_metadata(): + cfg = _cfg() + vg, a0, _a1 = _vg_two_alleles(cfg) + g = Genotype.from_dataconfig(cfg).homozygous(vg, a0).with_subject("S1") + # builder-path genotype: every source non-cartridge, no model id + prov = g.prior_provenance + assert prov["allele_frequencies"] == "manual" + assert prov["model_id"] is None + md = g.to_metadata() + assert md["subject_id"] == "S1" + assert md["prior_provenance"] == prov + assert "source_refdata_hash" in md + + +def test_prior_provenance_survives_snapshot(): + cfg = _cfg() + g = Genotype.sample(cfg, seed=1) + snap = g._snapshot() + assert snap.prior_provenance == g.prior_provenance + + +# ── Task 6: sample plane consumption (catalogue alleles) ───────── + + +def _planed_cfg(): + import copy + cfg = copy.deepcopy(_cfg()) + vg, a0, a1 = _vg_two_alleles(cfg) + cfg.genotype_priors = PopulationGenotypeModel( + model_id="m1", source="toy", + allele_frequencies={"V": {vg: {a0: 100.0, a1: 1.0}}}, + haplotype_deletion_prob={"V": {vg: 0.0}}, + chromosome_weights=(0.5, 0.5), + ) + return cfg, vg, a0, a1 + + +def test_sample_auto_uses_plane_with_provenance(): + cfg, vg, a0, a1 = _planed_cfg() + homo_a0 = sum(1 for s in range(60) + if Genotype.sample(cfg, seed=s).carried_alleles("V", vg) == {a0}) + assert homo_a0 > 40 # dominant plane allele -> usually homozygous-common + g = Genotype.sample(cfg, seed=1) + assert g.prior_provenance["allele_frequencies"] == "cartridge" + assert g.prior_provenance["haplotype_deletion_prob"] == "cartridge" + assert g.prior_provenance["chromosome_weights"] == "cartridge" + assert g.prior_provenance["model_id"] == "m1" + assert g.prior_provenance["model_checksum"] == cfg.genotype_priors.content_checksum() + + +def test_sample_opt_out_is_uniform(): + cfg, vg, a0, a1 = _planed_cfg() + g = Genotype.sample(cfg, seed=1, use_cartridge_priors=False) + assert g.prior_provenance["allele_frequencies"] == "uniform" + assert g.prior_provenance["haplotype_deletion_prob"] == "default" + assert g.prior_provenance["chromosome_weights"] == "default" + assert g.prior_provenance["model_id"] is None + # uniform: a1 should appear materially more than under the biased plane + a1_seen = sum(1 for s in range(60) + if a1 in Genotype.sample(cfg, seed=s, use_cartridge_priors=False) + .carried_alleles("V", vg)) + assert a1_seen > 5 + + +def test_sample_explicit_overrides_plane(): + cfg, vg, a0, a1 = _planed_cfg() + g = Genotype.sample(cfg, seed=1, allele_frequencies={"V": {vg: {a1: 1.0}}}) + assert g.prior_provenance["allele_frequencies"] == "explicit" + assert g.carried_alleles("V", vg) <= {a1} + + +def test_sample_mixed_sourcing(): + cfg, vg, a0, a1 = _planed_cfg() + g = Genotype.sample(cfg, seed=2, allele_frequencies={"V": {vg: {a0: 1.0}}}) + # explicit freq, cartridge deletion, cartridge chromosome weights + assert g.prior_provenance["allele_frequencies"] == "explicit" + assert g.prior_provenance["haplotype_deletion_prob"] == "cartridge" + assert g.prior_provenance["chromosome_weights"] == "cartridge" + + +def test_include_cartridge_novel_alleles_validated(): + cfg, vg, a0, a1 = _planed_cfg() + with pytest.raises(ValueError, match="include_cartridge_novel_alleles"): + Genotype.sample(cfg, seed=1, include_cartridge_novel_alleles="yes") + # integers must NOT slip through (1 == True / 0 == False in Python) + for bad in (0, 1, 2): + with pytest.raises(ValueError, match="include_cartridge_novel_alleles"): + Genotype.sample(cfg, seed=1, include_cartridge_novel_alleles=bad) + + +def test_sample_no_plane_unchanged(): + cfg = _cfg() # no plane + g = Genotype.sample(cfg, seed=3) + assert g.prior_provenance["allele_frequencies"] == "uniform" + assert g.prior_provenance["haplotype_deletion_prob"] == "default" + + +# ── Task 7: candidate-vs-carried novel injection ───────────────── + + +def _functional_novel_seq(cfg, base_name): + """Find a single-base substitution off ``base_name`` that Genotype accepts + as functional (avoids stop codons / anchor breakage), returning the seq.""" + base_allele = next(a for g in cfg.v_alleles.values() for a in g if a.name == base_name) + base_seq = base_allele.ungapped_seq.upper() + for pos in range(len(base_seq)): + for nt in "ACGT": + if nt == base_seq[pos]: + continue + cand = base_seq[:pos] + nt + base_seq[pos + 1:] + try: + (Genotype.from_dataconfig(cfg) + .add_novel_allele(f"{base_name.split('*')[0]}*97", base=base_name, + sequence=cand, segment="V")) + return cand + except ValueError: + continue + raise AssertionError("no functional single-base novel found") + + +def _planed_cfg_with_novel(freq_novel=1000.0): + import copy + cfg = copy.deepcopy(_cfg()) + vg = next(g for g, al in cfg.v_alleles.items() if len(al) >= 1) + base = cfg.v_alleles[vg][0].name + novel_seq = _functional_novel_seq(cfg, base) + novel = f"{vg}*97" + cfg.genotype_priors = PopulationGenotypeModel( + model_id="mN", source="toy", + allele_frequencies={"V": {vg: {base: 1.0}}}, + haplotype_deletion_prob={"V": {vg: 0.0}}, + novel_alleles=[PopulationNovelAllele(name=novel, segment="V", + base_allele=base, sequence=novel_seq, frequency=freq_novel)], + ) + return cfg, vg, base, novel + + +def test_plane_novel_is_drawn_and_carried(): + cfg, vg, base, novel = _planed_cfg_with_novel() + g = Genotype.sample(cfg, seed=1) # auto -> cartridge freqs -> novels injected + assert g.prior_provenance["novel_alleles"] == "cartridge" + assert novel in g.carried_alleles("V", vg) # dominant novel frequency + eff = g.effective_dataconfig() + assert any(a.name == novel for a in eff.v_alleles[vg]) + + +def test_uncarried_plane_novel_absent_from_export(): + # rare novel: base dominates -> novel essentially never carried + cfg, vg, base, novel = _planed_cfg_with_novel(freq_novel=1e-9) + g = Genotype.sample(cfg, seed=7) + assert novel not in g.carried_alleles("V", vg) + eff = g.effective_dataconfig() + assert all(a.name != novel for a in eff.v_alleles[vg]) + + +def test_novel_skipped_with_explicit_freqs_auto(): + cfg, vg, base, novel = _planed_cfg_with_novel() + g = Genotype.sample(cfg, seed=1, allele_frequencies={"V": {vg: {base: 1.0}}}) + assert g.prior_provenance["novel_alleles"] == "none" + assert novel not in g.carried_alleles("V", vg) + + +def test_novel_forced_with_explicit_freqs_true(): + cfg, vg, base, novel = _planed_cfg_with_novel() + g = Genotype.sample(cfg, seed=1, allele_frequencies={"V": {vg: {base: 1.0}}}, + include_cartridge_novel_alleles=True) + assert g.prior_provenance["novel_alleles"] == "cartridge" + # synthesized table: explicit base at 1.0 + novel at huge frequency -> novel wins + assert novel in g.carried_alleles("V", vg) + + +def test_novel_disabled(): + cfg, vg, base, novel = _planed_cfg_with_novel() + g = Genotype.sample(cfg, seed=1, include_cartridge_novel_alleles=False) + assert g.prior_provenance["novel_alleles"] == "none" + assert novel not in g.carried_alleles("V", vg) + + +# ── Task 8: from_genotypes pure estimator ──────────────────────── + + +def test_from_genotypes_counts_chromosomes_and_deletions(): + cfg = _cfg() + vg, a0, a1 = _vg_two_alleles(cfg) + # subject A: homozygous a0 (2 chromosomes of a0) + gA = (Genotype.from_dataconfig(cfg).homozygous(vg, a0) + .complete_from_reference().with_subject("A")) + # subject B: heterozygous a0/a1 (1 each) + gB = (Genotype.from_dataconfig(cfg).heterozygous(vg, a0, a1) + .complete_from_reference().with_subject("B")) + # subject C: a0 on hap0, deleted on hap1 (hemizygous) -> 1 deleted haplotype + gC = (Genotype.from_dataconfig(cfg).homozygous(vg, a0) + .delete_gene(vg, haplotype=1).complete_from_reference().with_subject("C")) + + m = PopulationGenotypeModel.from_genotypes([gA, gB, gC], cfg=cfg, + model_id="est", source="unit") + # a0 chromosomes: A=2, B=1, C=1 => 4 ; a1: B=1 => 1 + assert m.allele_frequencies["V"][vg][a0] == 4.0 + assert m.allele_frequencies["V"][vg][a1] == 1.0 + # deletions for vg: only C hap1 => 1 / (2*3) subjects + assert m.haplotype_deletion_prob["V"][vg] == pytest.approx(1.0 / 6.0) + m.validate(chain_type="vdj") # estimator output is shape-valid + + +def test_from_genotypes_duplicate_subject_id_rejected(): + cfg = _cfg() + g1 = Genotype.from_dataconfig(cfg).complete_from_reference().with_subject("X") + g2 = Genotype.from_dataconfig(cfg).complete_from_reference().with_subject("X") + with pytest.raises(ValueError, match="duplicate subject"): + PopulationGenotypeModel.from_genotypes([g1, g2], cfg=cfg, + model_id="e", source="u") + + +def test_from_genotypes_min_subjects(): + cfg = _cfg() + g1 = Genotype.from_dataconfig(cfg).complete_from_reference().with_subject("X") + with pytest.raises(ValueError, match="min_subjects"): + PopulationGenotypeModel.from_genotypes([g1], cfg=cfg, min_subjects=5, + model_id="e", source="u") + + +def test_from_genotypes_rejects_duplicated_gene(): + cfg = _cfg() + vg, a0, a1 = _vg_two_alleles(cfg) + g = (Genotype.from_dataconfig(cfg).homozygous(vg, a0) + .duplicate_gene(vg, [a0, a1], haplotype=0) # copy-number > 1 on a slot + .complete_from_reference().with_subject("D")) + with pytest.raises(ValueError, match="copy-number|duplicated gene"): + PopulationGenotypeModel.from_genotypes([g], cfg=cfg, model_id="e", source="u") + + +def test_from_genotypes_pseudocount_scope(): + cfg = _cfg() + vg, a0, a1 = _vg_two_alleles(cfg) + gA = (Genotype.from_dataconfig(cfg).homozygous(vg, a0) + .complete_from_reference().with_subject("A")) + m = PopulationGenotypeModel.from_genotypes([gA], cfg=cfg, pseudocount=0.5, + subject_id_policy="allow_duplicates", + model_id="e", source="u") + # a0 observed twice + 0.5 ; a1 unobserved but catalogue -> 0.5 (pseudocount only) + assert m.allele_frequencies["V"][vg][a0] == pytest.approx(2.5) + assert m.allele_frequencies["V"][vg][a1] == pytest.approx(0.5) + + +# ── Task 9: builder set/estimate_genotype_priors ───────────────── + + +def _builder_from_bundled(): + """A ReferenceCartridgeBuilder seeded from the bundled IGH catalogue.""" + from GenAIRR.cartridge_builder import ReferenceCartridgeBuilder + from GenAIRR.dataconfig.enums import ChainType + cfg = _cfg() + b = ReferenceCartridgeBuilder(ChainType.BCR_HEAVY) + b._v_alleles = {g: list(a) for g, a in cfg.v_alleles.items()} + b._d_alleles = {g: list(a) for g, a in (cfg.d_alleles or {}).items()} + b._j_alleles = {g: list(a) for g, a in cfg.j_alleles.items()} + b._metadata = cfg.metadata + return b, cfg + + +def test_set_genotype_priors_catalogue_aware(): + b, cfg = _builder_from_bundled() + vg = next(iter(cfg.v_alleles)) + good = PopulationGenotypeModel(model_id="m", source="s", + allele_frequencies={"V": {vg: {cfg.v_alleles[vg][0].name: 1.0}}}) + assert b.set_genotype_priors(good) is b # chainable + assert b._genotype_priors is good + + with pytest.raises(ValueError, match="not in the cartridge|unknown"): + b.set_genotype_priors(PopulationGenotypeModel(model_id="m", source="s", + allele_frequencies={"V": {"NOPE-GENE": {"NOPE*01": 1.0}}})) + with pytest.raises(ValueError, match="not a known allele|allele"): + b.set_genotype_priors(PopulationGenotypeModel(model_id="m", source="s", + allele_frequencies={"V": {vg: {"IGHVNOPE*99": 1.0}}})) + + +def test_set_genotype_priors_validates_novel(): + b, cfg = _builder_from_bundled() + with pytest.raises(ValueError, match="base allele|not found"): + b.set_genotype_priors(PopulationGenotypeModel(model_id="m", source="s", + novel_alleles=[PopulationNovelAllele(name="IGHV1-2*99", segment="V", + base_allele="IGHVNOPE*01", sequence="ACGT", frequency=1.0)])) + + +def test_estimate_genotype_priors_chainable_and_attaches(): + b, cfg = _builder_from_bundled() + vg, a0, a1 = _vg_two_alleles(cfg) + gA = Genotype.from_dataconfig(cfg).homozygous(vg, a0).complete_from_reference().with_subject("A") + gB = Genotype.from_dataconfig(cfg).heterozygous(vg, a0, a1).complete_from_reference().with_subject("B") + out = b.estimate_genotype_priors([gA, gB], model_id="est", source="cohort") + assert out is b + assert b._genotype_priors.allele_frequencies["V"][vg][a0] == 3.0 + + +def test_from_genotypes_requires_complete_genotypes(): + # A genotype missing a catalogue gene (no complete_from_reference) must NOT + # be silently scored as a double deletion — it should raise. + cfg = _cfg() + vg, a0, _a1 = _vg_two_alleles(cfg) + g = Genotype.from_dataconfig(cfg).homozygous(vg, a0).with_subject("P") # partial + with pytest.raises(ValueError, match="does not specify|complete_from_reference"): + PopulationGenotypeModel.from_genotypes([g], cfg=cfg, model_id="x", source="y") + + +def test_from_genotypes_all_none_subject_ids_rejected_under_require_unique(): + cfg = _cfg() + g1 = Genotype.from_dataconfig(cfg).complete_from_reference() # no subject id + g2 = Genotype.from_dataconfig(cfg).complete_from_reference() + with pytest.raises(ValueError, match="subject_id"): + PopulationGenotypeModel.from_genotypes([g1, g2], cfg=cfg, + model_id="x", source="y") + + +def test_estimate_with_novel_round_trips_and_samples(): + cfg, vg, base, novel = _planed_cfg_with_novel() # gives us a functional novel seq + novel_seq = next(a for a in cfg.genotype_priors.novel_alleles).sequence + base_cfg = _cfg() + gn = (Genotype.from_dataconfig(base_cfg) + .add_novel_allele(novel, base=base, sequence=novel_seq, segment="V") + .homozygous(vg, novel).complete_from_reference().with_subject("N")) + from GenAIRR.cartridge_builder import ReferenceCartridgeBuilder + from GenAIRR.dataconfig.enums import ChainType + b = ReferenceCartridgeBuilder(ChainType.BCR_HEAVY) + b._v_alleles = {g: list(a) for g, a in base_cfg.v_alleles.items()} + b._d_alleles = {g: list(a) for g, a in (base_cfg.d_alleles or {}).items()} + b._j_alleles = {g: list(a) for g, a in base_cfg.j_alleles.items()} + b._metadata = base_cfg.metadata + # round-trip must not raise (novel lives in novel_alleles, not allele_frequencies) + b.estimate_genotype_priors([gn], model_id="e", source="c", + subject_id_policy="allow_duplicates") + model = b._genotype_priors + assert any(nv.name == novel for nv in model.novel_alleles) + assert novel not in model.allele_frequencies.get("V", {}).get(vg, {}) + + +def test_from_genotypes_pseudocount_validation(): + cfg = _cfg() + g = Genotype.from_dataconfig(cfg).complete_from_reference().with_subject("A") + for bad in (-1.0, float("nan"), float("inf"), True): + with pytest.raises(ValueError, match="pseudocount"): + PopulationGenotypeModel.from_genotypes([g], cfg=cfg, pseudocount=bad, + model_id="x", source="y") + + +def test_from_genotypes_rejects_unknown_carried_allele(): + cfg = _cfg() + vg, a0, _a1 = _vg_two_alleles(cfg) + g = Genotype.from_dataconfig(cfg).homozygous(vg, a0).complete_from_reference().with_subject("A") + # forge an unregistered, non-catalogue allele directly into a slot + g._slots["V"][vg] = [[("BOGUS*01", 1, 1.0)], [(a0, 1, 1.0)]] + with pytest.raises(ValueError, match="neither a catalogue allele nor a registered novel"): + PopulationGenotypeModel.from_genotypes([g], cfg=cfg, model_id="x", source="y") + + +def test_from_genotypes_rejects_copy_count_gt_one(): + cfg = _cfg() + vg, a0, _a1 = _vg_two_alleles(cfg) + g = Genotype.from_dataconfig(cfg).homozygous(vg, a0).complete_from_reference().with_subject("A") + g._slots["V"][vg] = [[(a0, 2, 1.0)], [(a0, 1, 1.0)]] # single entry, copies=2 + with pytest.raises(ValueError, match="copy-number|duplicated"): + PopulationGenotypeModel.from_genotypes([g], cfg=cfg, model_id="x", source="y") + + +def test_from_genotypes_validates_identity_before_return(): + cfg = _cfg() + g = Genotype.from_dataconfig(cfg).complete_from_reference().with_subject("A") + with pytest.raises(ValueError, match="model_id"): + PopulationGenotypeModel.from_genotypes([g], cfg=cfg, model_id="", source="y") + + +def test_sample_validates_directly_attached_plane(): + import copy + cfg = copy.deepcopy(_cfg()) + vg, a0, _a1 = _vg_two_alleles(cfg) + # invalid plane attached directly (bypassing the builder) + cfg.genotype_priors = PopulationGenotypeModel(model_id="", source="", + allele_frequencies={"V": {vg: {a0: 1.0}}}) + with pytest.raises(ValueError, match="model_id|source"): + Genotype.sample(cfg, seed=1) + + +def test_manifest_invalid_plane_marked_not_valid(): + import copy, math + cfg = copy.deepcopy(_cfg()) + cfg.genotype_priors = PopulationGenotypeModel(model_id="m", source="s", + chromosome_weights=(float("nan"), 1.0)) + block = cfg.cartridge_manifest()["models"]["genotype_priors"] + assert block["available"] is True + assert block["valid"] is False + # no non-JSON-clean (NaN) numerics leak through + cw = block["chromosome_weights"] + assert cw is None or all(math.isfinite(x) for x in cw) + + +def test_to_metadata_effective_refdata_hash_for_novel(): + cfg, vg, base, novel = _planed_cfg_with_novel() + g = Genotype.sample(cfg, seed=1) + assert novel in g.carried_alleles("V", vg) + md = g.to_metadata() + assert "effective_refdata_hash" in md + assert md["effective_refdata_hash"] != md["source_refdata_hash"] + # a no-novel genotype: effective == source + g2 = Genotype.from_dataconfig(_cfg()).complete_from_reference() + md2 = g2.to_metadata() + assert md2["effective_refdata_hash"] == md2["source_refdata_hash"] + + +def test_manifest_never_crashes_on_garbage_plane(): + import copy, math + cfg = copy.deepcopy(_cfg()) + # non-finite frequency weight + bad chromosome weight type -> content_checksum + # and float() would raise; the manifest must stay JSON-clean and not crash. + cfg.genotype_priors = PopulationGenotypeModel(model_id="m", source="s", + chromosome_weights=("x", 1.0)) + block = cfg.cartridge_manifest()["models"]["genotype_priors"] + assert block["available"] is True and block["valid"] is False + assert block["chromosome_weights"] is None + assert block["model_checksum"] is None or isinstance(block["model_checksum"], str) + + # genotype_priors set to a non-model object must not crash the manifest either + cfg.genotype_priors = object() + block = cfg.cartridge_manifest()["models"]["genotype_priors"] + assert block["available"] is True and block["valid"] is False + + +def test_use_cartridge_priors_must_be_bool(): + cfg, vg, a0, a1 = _planed_cfg() + for bad in ("False", 1, 0, None): + with pytest.raises(ValueError, match="use_cartridge_priors"): + Genotype.sample(cfg, seed=1, use_cartridge_priors=bad) + + +def test_from_genotypes_include_novel_must_be_bool(): + cfg = _cfg() + g = Genotype.from_dataconfig(cfg).complete_from_reference().with_subject("A") + for bad in ("yes", 1, None): + with pytest.raises(ValueError, match="include_novel"): + PopulationGenotypeModel.from_genotypes([g], cfg=cfg, include_novel=bad, + model_id="x", source="y") + + +def test_from_genotypes_min_subjects_and_segments_validation(): + cfg = _cfg() + g = Genotype.from_dataconfig(cfg).complete_from_reference().with_subject("A") + for bad in (-1, 0, True, "2"): + with pytest.raises(ValueError, match="min_subjects"): + PopulationGenotypeModel.from_genotypes([g], cfg=cfg, min_subjects=bad, + model_id="x", source="y") + with pytest.raises(ValueError, match="segment"): + PopulationGenotypeModel.from_genotypes([g], cfg=cfg, segments=("X",), + model_id="x", source="y") + with pytest.raises(ValueError, match="segment"): + PopulationGenotypeModel.from_genotypes([g], cfg=cfg, segments="V", + model_id="x", source="y") + with pytest.raises(ValueError, match="segment"): + PopulationGenotypeModel.from_genotypes([g], cfg=cfg, segments=("V", "V"), + model_id="x", source="y") + + +def test_sample_draw_independent_of_freq_dict_order(): + import copy + cfg = _cfg() + vg, a0, a1 = _vg_two_alleles(cfg) + c1 = copy.deepcopy(cfg) + c1.genotype_priors = PopulationGenotypeModel(model_id="m", source="s", + allele_frequencies={"V": {vg: {a0: 2.0, a1: 3.0}}}) + c2 = copy.deepcopy(cfg) + c2.genotype_priors = PopulationGenotypeModel(model_id="m", source="s", + allele_frequencies={"V": {vg: {a1: 3.0, a0: 2.0}}}) # reversed insertion order + assert c1.genotype_priors.content_checksum() == c2.genotype_priors.content_checksum() + for s in range(15): + assert (Genotype.sample(c1, seed=s).to_table() + == Genotype.sample(c2, seed=s).to_table()) + + +# ── Task 10: end-to-end ────────────────────────────────────────── + + +def test_end_to_end_planed_cartridge_truth_calls_carried(): + cfg, vg, a0, a1 = _planed_cfg() + g = Genotype.sample(cfg, seed=11) + assert g.prior_provenance["allele_frequencies"] == "cartridge" + res = ga.Experiment.on(cfg).with_genotype(g).recombine().run_records( + n=200, seed=2, expose_provenance=True) + assert len(res) == 200 + for r in res: + for seg, col in (("V", "truth_v_call"), ("D", "truth_d_call"), ("J", "truth_j_call")): + call = r[col] + if not call: + continue + gene = call.split("*")[0] + assert call in g.carried_alleles(seg, gene), (seg, call) + + +def test_end_to_end_planed_novel_flows_to_records(): + cfg, vg, base, novel = _planed_cfg_with_novel() + g = Genotype.sample(cfg, seed=1) + assert novel in g.carried_alleles("V", vg) + res = ga.Experiment.on(cfg).with_genotype(g).recombine().run_records( + n=50, seed=3, expose_provenance=True) + # the novel can legitimately appear as a V truth call + assert any(r["truth_v_call"] == novel for r in res) diff --git a/tests/test_reference_cartridge_authoring_contract.py b/tests/test_reference_cartridge_authoring_contract.py index 62c78d9..447ef88 100644 --- a/tests/test_reference_cartridge_authoring_contract.py +++ b/tests/test_reference_cartridge_authoring_contract.py @@ -434,6 +434,9 @@ def test_pin_estimate_step_method_boundary_post_p_nucleotide_length_slice() -> N "estimate_np_length_distributions", "estimate_np_base_model", "estimate_p_nucleotide_lengths", + # Cartridge genotype plane slice — estimates a population genotype + # prior from observed Genotype objects and attaches it. + "estimate_genotype_priors", } for owner_name in ("DataConfig", "RefDataConfig", "ReferenceCartridgeBuilder"): owner = getattr(ga, owner_name, None)