diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 08ed8f6..5461123 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -17,7 +17,7 @@ jobs: python-version: "3.12" - name: Install the released downstream consumer run: | - bash scripts/pip_install.sh pytest 'vaxrank==3.35.0' -e . + bash scripts/pip_install.sh pytest 'vaxrank==3.36.0' -e . python -m pip check - name: Verify table features through DSL scoring and vaccine construction run: python -m pytest scripts/check_vaxrank_candidates.py -q diff --git a/CHANGELOG.md b/CHANGELOG.md index 3919777..1cd51a7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,19 @@ between published tags; older pre-5.0 notes are retained below. For current interfaces, see the [consumer guide](docs/consumer-guide.md). +## 5.92.0 + +- Add all-candidate self matching within an explicit same-length Hamming radius + and exact matching across complete vaccine windows. Retain every origin, + including shared CTA/non-CTA sequences, with caller-resolved exclusions (#455). +- Preserve supplied presentation observations, allele attribution, conflicting + predictions and same-allele coverage separately from search completeness. + Typed CSV/TSV retains evidence and search identity; missing evidence remains + unknown. Neither presentation nor similarity establishes TCR recognition. +- Test against Vaxrank 3.36.0's published `builtin:openvax-v1` bundle. A synthetic + consumer workflow changes the selected window using self evidence while + retaining both intended epitopes and alleles; Topiary adds no selection rule. + ## 5.91.0 - Score explicitly identified peptide occurrences with their original flanks, diff --git a/RELEASING.md b/RELEASING.md index 666ebf9..b2c1a8c 100644 --- a/RELEASING.md +++ b/RELEASING.md @@ -57,7 +57,7 @@ CI installs a pinned published Vaxrank release and runs the real candidate scoring and vaccine-construction workflow. To repeat it locally: ```bash -"$PYTHON" -m pip install 'vaxrank==3.35.0' +"$PYTHON" -m pip install 'vaxrank==3.36.0' "$PYTHON" -m pip check "$PYTHON" -m pytest scripts/check_vaxrank_candidates.py -q ``` diff --git a/docs/api.md b/docs/api.md index d0d6a01..530883d 100644 --- a/docs/api.md +++ b/docs/api.md @@ -70,6 +70,18 @@ symmetric conservative transformation, and O/U/X/* or unrecognized characters receive the symmetric worst-case canonical distance (15). Both functions return fresh read-only int8 arrays, so callers never share mutable matrix state. +## Self-match evidence + +`match_self_peptides(reference, peptides, *, max_mismatches=0, alleles=None, +excluded_gene_ids=(), observations=None, predictions=None)` returns every +same-length Hamming match within the radius, with all source origins and +supplied presentation/prediction records. `reference.match_candidates(...)` +delegates to it. `self_matches_in_windows(windows, reference, *, peptide_lengths, +...)` searches exact peptides at every requested window position. Results use +ordinary typed CSV/TSV, with explicit coverage and evidence identity. +See [self evidence](self_proteome.md#all-match-evidence) for the input schemas, +policy composition and interpretation of missing results. + ## TopiaryPredictor | Parameter | Type | Description | diff --git a/docs/combined-sources.md b/docs/combined-sources.md index c15680e..7fdb211 100644 --- a/docs/combined-sources.md +++ b/docs/combined-sources.md @@ -54,9 +54,11 @@ version-qualified expression for this simple case. `SelectionPolicy` saves the exact candidate filter, score, model/version selections, ranking direction, duplicate policy and strata under a stable name. -It has no implicit scientific recipe. Topiary does not currently ship an -official `openvax-v1` preset; that name can identify a reviewed definition you -freeze, and later definitions can use new names. +It has no implicit scientific recipe. **Vaxrank owns the published +`builtin:openvax-v1` bundle**, shipped with experimental overlays in +[Vaxrank 3.36.0 / PR #565](https://github.com/openvax/vaxrank/pull/565). +Topiary evaluates its policy subtree and tests against that released bundle; +it does not maintain a second definition. Give changed policies new names. ```python from topiary import ( @@ -115,7 +117,7 @@ site-level evidence remains [#288](https://github.com/openvax/topiary/issues/288 ### Composed consumer configuration Vaxrank owns YAML composition and window/construct configuration. Its repeated -`--config openvax-v1.yaml --config overrides.yaml` workflow merges mappings +`--config builtin:openvax-v1 --config overrides.yaml` workflow merges mappings left to right and replaces lists/scalars. The Topiary integration boundary is the **already-composed policy subtree**: @@ -142,16 +144,13 @@ the old expression. Model selections are mappings; version selections are lists of `{kind, method, version}` records, so an override replaces that entire list. Topiary does not merge files or parse unrelated consumer settings. -Vaxrank adoption of this subtree remains -[Vaxrank #497](https://github.com/openvax/vaxrank/issues/497); the example above -describes the consumer boundary, not a newly supported Vaxrank configuration -key. Its current occurrence-level scorer can consume `policy.score_by`, -`policy.filter_by`, and the selection mappings through the existing public -`EvalContext` / `apply_filter` APIs, retaining its explicit grouping and -per-occurrence `alleles` lookup. Vaxrank's post-score minimum gate, missing-score -fill, window rules, RNA weighting and construct settings must also be frozen in -its bundle. A cutoff inside a score expression must not become a destructive -pre-filter when capturing that baseline. +Vaxrank 3.36.0 consumes this subtree in its shipped configuration bundle. +Its post-score minimum gate, missing-score fill, window rules, RNA weighting +and construct settings belong to that bundle too. The consumer fixture in +`scripts/check_vaxrank_candidates.py` loads `builtin:openvax-v1` through +Vaxrank's public config loader and compares its numeric scores with +`evaluate_selection_policy`. A cutoff inside a score expression must not become +a destructive pre-filter when replaying that baseline. The representative ranking above is one view. `evaluate_selection_policy` retains every input row and evaluates explicit occurrence identities before @@ -175,8 +174,8 @@ evaluated.evidence.to_tsv("complete-evaluation.tsv") replayed = replay_selection_policy(read_tsv("complete-evaluation.tsv")) ``` -This example freezes a scoring convention; it does not establish an official -`openvax-v1` recipe or a calibrated probability of immunogenicity. The cutoff +This example illustrates a scoring convention; use Vaxrank's published bundle +for `openvax-v1`. Neither expression is a calibrated probability of immunogenicity. The cutoff remains inside the score, with a separate inclusive minimum-score gate. `raw_score` preserves missingness even when `score_fill=0` supplies an effective zero. Pre-filtered groups are not scored and cannot be restored by filling. diff --git a/docs/self_proteome.md b/docs/self_proteome.md index a8ffe42..81479ba 100644 --- a/docs/self_proteome.md +++ b/docs/self_proteome.md @@ -1,21 +1,107 @@ -# Cross-reactivity: `self_nearest_*` via `SelfProteome` +# Self-sequence and presentation evidence -Given a mutant peptide, find its closest match in a reference proteome -of healthy self peptides. The distance (and the source gene of the -match) is a cross-reactivity risk signal: mutant neoantigens that look -a lot like a real self-peptide may trigger T-cell cross-reactivity -against healthy tissue. +`SelfProteome` searches an explicit reference corpus. Sequence similarity, +predicted binding and observed presentation are separate facts; none establishes +TCR recognition. Peptide-loaded or minigene recognition can also differ from +recognition of native protein processing and presentation ([primary study](https://www.nature.com/articles/s41541-023-00713-y)). -This page covers the `SelfProteome` class and the `self_nearest_*` -columns it adds to `TopiaryPredictor` output. +Use [all-match evidence](#all-match-evidence) for complete vaccine-window searches +or multiple similar candidates. The existing `nearest()` and `self_nearest_*` +predictor columns remain a single-nearest-sequence view. > **Status.** `SelfProteome` currently exposes the sequence-nearest > axis: substitutions plus 1aa indel neighbors against a scoped self > proteome (`include="all"`, `"non_cta"`, `"protected_tissues"`, or a > callable). Binding-aware axes (`self_mimic_*`, `self_strongest_nearby_*`) -> and the structured full-candidate column (`self_nearest_candidates`) -> require MHC prediction on candidate peptides and are tracked under -> [#124](https://github.com/openvax/topiary/issues/124). +> remain [#412](https://github.com/openvax/topiary/issues/412). +> All-match evidence accepts supplied observations and predictions; it does not +> run models, rank by binding or infer TCR recognition. + +## All-match evidence + +```python +from topiary import SelfProteome, match_self_peptides, self_matches_in_windows + +# Synthetic reference. In production, retain an unfiltered corpus so a shared +# sequence keeps both CTA and non-CTA origins. oncoref owns CTA membership. +reference = SelfProteome.from_peptides( + {"CTA": "GILGFVFTL", "healthy": "GILGFVFTL", "near": "SIINFEKM"}, + peptide_lengths=[8, 9], +) +windows = {"full": "GILGFVFTLGGGSIINFEKLELAGIGILT", "trimmed": "SIINFEKLELAGIGILT"} +exact = self_matches_in_windows( + windows, reference, peptide_lengths=[8, 9], + alleles={name: ["A0201", "B0702"] for name in windows}, + excluded_gene_ids={"CTA"}, # illustrative caller-resolved exclusion +) +similar = match_self_peptides(reference, ["SIINFEKL"], max_mismatches=1, alleles=["A0201"]) +exact.to_tsv("window-self-evidence.tsv") +``` + +`reference.match_candidates(...)` delegates to `match_self_peptides`. Queries +retain order and duplicates through `query_index`; every matched sequence and +gene/transcript/reference-offset origin has a row. `self_match_id` identifies the +reference occurrence across queries. Window output also carries `window_id`, +`window_sequence` and every zero-based `peptide_offset`. The exact adapter checks +all requested lengths and offsets, not only mutation-overlapping peptides. + +Exclusions flag `self_in_scope=False`; they never delete origins. A reference +already filtered with `include="non_cta"` cannot recover removed CTA origins. +For an Ensembl corpus use `include="all"` and resolve the exclusion set with +`cta_gene_ids(source="oncoref")` for the appropriate species/tier. A sequence +match alone is not evidence that the peptide was presented. + +Both functions accept `observations=` and `predictions=` as tables or record +lists. They preserve multiple records in typed list/dict columns: + +| Input | Required fields | Interpretation | +|---|---|---| +| Observations | `evidence_id`, `peptide`, `source`, `evidence_kind`, `allele_assignment` | `evidence_kind` is `observed` or `predicted`; allele assignment is `confirmed`, `predicted` or `unknown`. Confirmed/predicted assignments require an explicit `allele`. | +| Predictions | `peptide`, `allele`, `kind`, `value`, `prediction_method_name`, `predictor_version` | Preserve supplied model facts, including unknown versions, missing values and conflicts; no averaging or model choice. | + +Observation records may carry tissue, assay, source version and other +JSON-compatible provenance. An optional `gene_id` restricts attribution to that +origin; its absence does not identify the producing gene. An unresolved +`allele_set` remains a genotype, not confirmed restriction. Alleles are parsed +with mhcgnomes, including non-human alleles. + +`self_observations` and `self_predictions` retain the records. +`self_same_allele_predictions` contains only finite numeric `value` measurements +at the query allele with per-allele MHC scope. Haplotype presenter labels do not +become per-allele evidence. This reports any supplied measurement, not complete +coverage of a desired model panel. Other-allele records remain in +`self_predictions`. Observation flags distinguish confirmed and inferred +same-allele assignments; unreported observations remain unknown. + +Coverage is explicit: + +| Field | Meaning | +|---|---| +| `self_search_status` | `matched`, `no_match_in_scope`, or `unassessed`. | +| `self_search_complete`, `self_search_reason` | Whether the requested canonical, same-length search is complete; unavailable lengths and unsupported query/reference sequences have explicit reasons. A match can coexist with incomplete coverage. | +| `self_observation_status` | Observation input absent, a record reported, or no record reported for this match/origin. | +| `self_prediction_status` | Prediction input absent, query allele unknown, same-allele measurement reported, or not reported. | + +`no_match_in_scope` is negative evidence only within the declared reference, +lengths and Hamming radius. It does not mean no presentable peptide or no risk. +Indels, substitutions outside the radius and unsearched reference data are not +assessed. Comparisons use chunks of 65,536 reference rows; exact queries reuse +the reference index. Metadata under `extra['self_search']` records the reference +identity, parameters, coverage counts and hashes of normalized supplied evidence. +`extra['self_window_coverage']` also records requests too short for a peptide. + +CSV/TSV preserves these fields and nested records. This is an evidence table, +not a fabricated prediction table: join it to candidate measurements at the +explicit query/allele identity before using a prediction policy. Preserve all +matches and their IDs when there are several per candidate. The composed test +in `tests/test_consumer_workflows.py` checks named-criterion unknown handling and +exact replay after joining and saving the evidence. + +Vaxrank owns window selection. `scripts/check_vaxrank_candidates.py` consumes its +released `builtin:openvax-v1` bundle and shows that self evidence changes a window +choice while retaining both intended epitopes and their alleles. That synthetic +fixture uses 100% target-score retention and checks every required target; it +does not establish that a general aggregate-score threshold protects each target. ## Basic usage @@ -24,7 +110,7 @@ from topiary import SelfProteome, TopiaryPredictor from mhctools import NetMHCpan ref = SelfProteome.from_ensembl(species="human") -# Default include="non_cta" strips CTAs via pirlygenes. +# Default include="non_cta" strips CTAs using oncoref membership via pirlygenes. predictor = TopiaryPredictor( models=NetMHCpan, @@ -64,8 +150,11 @@ Three construction modes: |---|---|---| | `"all"` | Whole proteome, no filter | — | | `"non_cta"` (default for human Ensembl) | Remove CTA genes | `cta_source="pirlygenes"` default; set / callable accepted | -| `"protected_tissues"` | Keep only genes expressed in named tissues | `tissues=[…]`, `tissue_source="hpa"`/`"gtex"`, `min_expression=…` | -| callable | Arbitrary `gene → bool` filter | — | +| `"protected_tissues"` | Keep only genes expressed in named tissues | `tissues=[…]`, `min_tissue_ntpm=…` for human HPA data; or explicit `tissue_gene_ids={…}` for any species | +| callable | Arbitrary `gene_id → bool` filter | — | + +Tissue expression selects the reference genes; it is not evidence of peptide +presentation. Supply presentation observations separately when available. **Human users** get zero-config `include="non_cta"` via pirlygenes: @@ -74,7 +163,7 @@ ref = SelfProteome.from_ensembl(species="human") ``` **Non-human users** must either use `include="all"` or supply their own -CTA source, because pirlygenes is human-only today: +CTA source, because oncoref's CTA membership is human-only today: ```python ref = SelfProteome.from_ensembl(species="mouse", release=102, include="all") @@ -89,8 +178,7 @@ ref = SelfProteome.from_ensembl( ``` A non-human `include="non_cta"` call without `cta_source=` raises at -construction — silent unfiltered results would be a misleading -cross-reactivity signal. +construction — silent unfiltered results would misstate the reference scope. ## Non-Ensembl sources @@ -115,12 +203,12 @@ ref = SelfProteome.from_peptides( ## Reference version Every row of the output carries a `self_nearest_reference_version` -string. Two runs produce interchangeable `self_nearest_*` values iff the -strings match. Its form is +string. Matching strings identify the same indexed reference; comparing results +also requires the same queries, metric and search settings. Its form is `{source}-{species}[-{release}]+include-{scope}+sha256:{digest}`: ``` -ensembl-human-115+include-non_cta+cta-pirlygenes-6.0.4+sha256:3f2a9c1e7b40 +ensembl-human-115+include-non_cta+cta-oncoref-VERSION+sha256:3f2a9c1e7b40 ensembl-mouse-102+include-all+sha256:91d0c4e2a8b7 fasta-fasta+include-callable-keep_named+sha256:0be51d2c9f63 peptides-synthetic+include-all+sha256:bd8a4e854d1c @@ -153,12 +241,10 @@ when enabled. - Construction: one pass over the reference proteome extracts every L-mer for each configured length, dedupes per length, and encodes - into a `(M, L) int8` array. For human non-CTA Ensembl × length 9, - expect ~200k rows. -- Lookup: per query, the full reference array is compared in one SIMD - operation. Chunked to bound peak memory. Typical throughput for - ~200k reference × ~10k queries is seconds. - -Seed-and-extend indexing and 1aa indel candidates are queued in -[#124](https://github.com/openvax/topiary/issues/124) — benchmark -decides the default algorithm. + into a `(M, L) int8` array. Corpus size depends on the release, scope and lengths. +- Lookup: vectorized comparisons are chunked to bound temporary memory. + Runtime depends on the corpus and query count; exact matching reuses a hash index. + +The existing `nearest()` path checks one-residue indel neighbors. The new +all-match API is limited to same-length Hamming candidates; broader candidate +search and binding-ranked axes remain [#412](https://github.com/openvax/topiary/issues/412). diff --git a/scripts/check_vaxrank_candidates.py b/scripts/check_vaxrank_candidates.py index 12c1922..87c604d 100644 --- a/scripts/check_vaxrank_candidates.py +++ b/scripts/check_vaxrank_candidates.py @@ -22,6 +22,91 @@ from vaxrank.vaccine_peptide import VaccinePeptide +def test_self_evidence_changes_vaxrank_window_without_losing_required_targets(tmp_path): + from topiary import ( + SelfProteome, self_matches_in_windows, read_tsv, from_predictions, + TopiaryResult, evaluate_selection_policy, + ) + from vaxrank.config.loader import load_vaxrank_config, extract_epitope_config_kwargs + from vaxrank.epitope_dsl import score_predictions, prediction_group_columns + from vaxrank.vaccine_antigen import SelfReferenceMatch, SelfReferenceSource + from vaxrank.window_selection import WindowSelection, select_windows + + # Consume Vaxrank's actual frozen bundle; Topiary keeps no second recipe. + bundle = load_vaxrank_config(config_path="builtin:openvax-v1") + cfg = EpitopeConfig(**extract_epitope_config_kwargs(bundle)) + baseline = SelectionPolicy.from_dict(cfg.selection_policy) + assert baseline.name == "openvax-v1" + full = "GILGFVFTLGGGSIINFEKLELAGIGILT" + windows = {"full": full, "trimmed": full[12:]} + alleles = ["HLA-A*02:01", "HLA-B*07:02"] + targets = {"SIINFEKL", "ELAGIGILT"} + ref = SelfProteome.from_peptides({"CTA": "GILGFVFTL", "healthy": "GILGFVFTL", "near": "SIINFEKM"}, + peptide_lengths=[8, 9]) + evidence = self_matches_in_windows(windows, ref, peptide_lengths=[8, 9], + alleles={name: alleles for name in windows}, excluded_gene_ids={"CTA"}) + evidence.to_tsv(tmp_path / "self-evidence.tsv") + evidence = read_tsv(tmp_path / "self-evidence.tsv") + hits = evidence.df[evidence.df.self_peptide.eq("GILGFVFTL")] + assert set(hits.self_gene_id) == {"CTA", "healthy"} + assert hits.groupby("allele").self_in_scope.sum().eq(1).all() + + def vaccine(name, use_evidence): + sequence = windows[name] + rows, measured = [], [] + for peptide in ["GILGFVFTL", "SIINFEKL", "ELAGIGILT"]: + if peptide not in sequence: + continue + offset = sequence.index(peptide) + origins = evidence.df[(evidence.df.window_id == name) & (evidence.df.peptide == peptide) + & evidence.df.self_in_scope.eq(True)] + sources = tuple(SelfReferenceSource(gene_id=gene, transcript_id=transcript) + for gene, transcript in origins[["self_gene_id", "self_transcript_id"]] + .drop_duplicates().itertuples(index=False, name=None)) if use_evidence else () + match = SelfReferenceMatch(peptide, bool(sources), "mutation", excluded_gene_ids=("CTA",), + sources=sources, source_provenance_complete=use_evidence, + genome_release=ref.reference_version) + for allele in alleles: + prediction = Prediction(peptide=peptide, allele=allele, kind="pMHC_affinity", value=50., + score=.9, predictor_name="synthetic", predictor_version="1") + rows.append(dict(peptide=peptide, source=sequence, offset=offset, prediction_id=name, + mutant=prediction, source_class="mutation", overlaps_targetable=peptide in targets, + occurs_in_reference=bool(sources), occurs_in_non_CTA_reference=bool(sources), + self_reference_match=match, patient_alleles=alleles)) + measured.append(dict(peptide=peptide, allele=allele, kind="pMHC_affinity", value=50., score=.9, + prediction_method_name="synthetic", predictor_version="1", prediction_id=name, + peptide_offset=offset, n_flank="", c_flank="")) + epitopes = candidate_epitopes_from_rows(rows) + frame = from_predictions(pd.DataFrame(measured)) + scored = attach_per_allele_scores(epitopes, cfg, topiary_df=frame) + expected = score_predictions(epitopes, cfg, topiary_df=frame) + evaluated = evaluate_selection_policy(TopiaryResult(frame), baseline, group_keys=prediction_group_columns(frame)) + actual = evaluated.selected.set_index(prediction_group_columns(frame)).score.sort_index() + pd.testing.assert_series_equal(actual, expected.sort_index(), check_names=False, check_exact=True) + antigen = VaccineAntigen(kind="mutation", amino_acids=sequence, + targetable_mask=TargetableMask((AminoAcidInterval(sequence.index("SIINFEKL"), len(sequence)),)), + tumor_specificity=TumorSpecificityAttestation(status="admitted", evidence_kind="synthetic_fixture", + evidence_source="check_vaxrank_candidates", patient_specific=True, rationale_code="test_only"), + source_identifier=name, self_reference_excluded_gene_ids=("CTA",)) + return VaccinePeptide(antigen=antigen, epitopes=scored, combined_score_expr="target_epitope_score", + ranking_rules=("target_epitope_score",)) + + policy = WindowSelection(self_weight=1., min_target_fraction=1., serum_weight=0.) + unchanged = select_windows([vaccine(name, False) for name in windows], policy, preferred_length=len(full), limit=1) + candidates = [vaccine(name, True) for name in windows] + changed = select_windows(candidates, policy, preferred_length=len(full), limit=1) + assert unchanged[0].amino_acids == full + assert changed[0].amino_acids == windows["trimmed"] + required = {(peptide, allele) for peptide in targets for allele in alleles} + for candidate in (candidates[0], changed[0]): + retained = {(epitope.sequence, allele) for epitope in candidate.target_epitopes + for allele, score in epitope.per_allele_scores.items() if score > 0} + assert retained == required + assert changed[0].target_epitope_score == candidates[0].target_epitope_score + assert changed[0].window_selection_audit["non_cta_self_score"] == 0. + assert candidates[0].window_selection_audit["non_cta_self_score"] > 0. + + def test_explicit_windows_keep_repeated_source_ids_through_vaxrank_scoring(): from topiary import TopiaryPredictor, TopiaryResult, evaluate_selection_policy from vaxrank.epitope_dsl import prediction_group_columns diff --git a/tests/test_consumer_workflows.py b/tests/test_consumer_workflows.py index 583e784..f5210b4 100644 --- a/tests/test_consumer_workflows.py +++ b/tests/test_consumer_workflows.py @@ -277,6 +277,47 @@ def test_peptide_level_measurement_conflicts_are_rejected_before_projection(): ) +def test_self_evidence_unknown_decisions_survive_policy_replay(tmp_path, pandas_string_inference): + from dataclasses import replace + from topiary import (SelfProteome, SelectionCriterion, SelectionPolicy, TopiaryResult, match_self_peptides, + evaluate_selection_policy, replay_selection_policy, from_predictions) + from .test_twin_conformance import DELIMITED_IO_TWINS + + reference = SelfProteome.from_peptides({"one": "SIINFEKL", "two": "SIINFEKM"}, peptide_lengths=[8]) + evidence = match_self_peptides(reference, ["SIINFEKL", "SIINFEKM", "AAAAAAAA"], + alleles=["A0201"] * 3, observations=[dict(evidence_id="observed", peptide="SIINFEKL", + source="synthetic", tissue="heart", evidence_kind="observed", + allele_assignment="confirmed", allele="A0201")]) + # Attach the evidence to explicit synthetic candidate measurements. The raw + # self-search table carries observations, not invented MHC measurements. + predictions = from_predictions(pd.DataFrame(dict( + peptide=["SIINFEKL", "SIINFEKM", "AAAAAAAA"], allele="HLA-A*02:01", + query_index=[0, 1, 2], kind="pMHC_affinity", value=[10., 20., 30.], + prediction_method_name="synthetic", predictor_version="1"))) + joined = predictions.merge(evidence.df, on=["query_index", "peptide", "allele", "peptide_length"], + validate="one_to_one") + candidates = TopiaryResult(joined, extra=evidence.extra) + # This selects reported observations for review; it is not a vaccine-safety rule. + policy = SelectionPolicy("reported-self", "affinity.value", filter_by='criterion("observed")', strata=(), + criteria=(SelectionCriterion("observed", "self_observed_confirmed_same_allele == 1", "eligibility"),)) + keys = ["query_index", "peptide", "allele"] + evaluated = evaluate_selection_policy(candidates, policy, group_keys=keys) + assert evaluated.selected.query_index.tolist() == [0] + assert evaluated.audit.status.tolist() == ["pass", "unknown", "unknown"] + included = evaluate_selection_policy(candidates, replace(policy, unknown="include"), group_keys=keys) + assert included.selected.query_index.tolist() == [0, 1, 2] + assert included.audit.status.tolist() == ["pass", "unknown", "unknown"] + assert included.evidence.df.self_search_status.tolist() == ["matched", "matched", "no_match_in_scope"] + for suffix, writer, _, reader in DELIMITED_IO_TWINS: + path = tmp_path / ("self-policy." + suffix) + writer(included.evidence, path) + replay = replay_selection_policy(reader(path)) + pd.testing.assert_frame_equal(replay.occurrences, included.occurrences, check_exact=True) + pd.testing.assert_frame_equal(replay.audit, included.audit, check_exact=True) + assert replay.evidence.df.self_observations.tolist() == evidence.df.self_observations.tolist() + assert replay.evidence.extra["self_search"] == evidence.extra["self_search"] + + @pytest.mark.parametrize("extension", ["csv", "tsv"]) def test_explicit_occurrences_from_mixed_sources_rescore_and_replay(tmp_path, monkeypatch, extension): from dataclasses import replace diff --git a/tests/test_self_evidence.py b/tests/test_self_evidence.py new file mode 100644 index 0000000..3096108 --- /dev/null +++ b/tests/test_self_evidence.py @@ -0,0 +1,158 @@ +"""All-match reference searches preserve origins, scope and evidence facts.""" + +import pandas as pd +import pytest +import numpy as np + +from topiary import SelfProteome, match_self_peptides, self_matches_in_windows, read_tsv +from .test_twin_conformance import SELF_MATCH_TWINS, SELF_WINDOW_TWINS + + +def reference(): + return SelfProteome.from_peptides({"CTA": "SIINFEKL", "healthy": "SIINFEKL", + "paralog": "SIINFEKM"}, peptide_lengths=[8]) + + +def observations(): + return [dict(evidence_id="ms-confirmed", peptide="SIINFEKL", source="synthetic-ms", tissue="heart", + evidence_kind="observed", allele_assignment="confirmed", allele="A0201"), + dict(evidence_id="ms-predicted", peptide="SIINFEKL", source="synthetic-ms", tissue="liver", + evidence_kind="observed", allele_assignment="predicted", allele="B0702"), + dict(evidence_id="unresolved", peptide="SIINFEKM", source="synthetic-ms", tissue="lung", + evidence_kind="observed", allele_assignment="unknown", allele_set="A0201,B0702")] + + +def predictions(): + return [dict(peptide="SIINFEKL", allele="A0201", kind="pMHC_affinity", value=value, + prediction_method_name="synthetic", predictor_version=None) for value in [0., 50.]] + + +@pytest.mark.parametrize("search", SELF_MATCH_TWINS) +def test_all_origins_and_candidate_sequences_are_retained(search): + result = search(reference(), ["SIINFEKL", "SIINFEKL"], max_mismatches=1, excluded_gene_ids={"CTA"}) + assert result.df.self_peptide.tolist() == ["SIINFEKL", "SIINFEKL", "SIINFEKM"] * 2 + assert result.df.self_gene_id.tolist() == ["CTA", "healthy", "paralog"] * 2 + assert result.df.self_in_scope.tolist() == [False, True, True] * 2 + assert result.df.self_mismatches.tolist() == [0, 0, 1] * 2 + assert result.df.query_index.tolist() == [0, 0, 0, 1, 1, 1] + assert result.extra["self_search"]["include_indels"] is False + assert result.df.self_match_id.nunique() == 3 + + +@pytest.mark.parametrize("search", SELF_WINDOW_TWINS) +def test_window_and_direct_exact_queries_agree(search, pandas_string_inference): + result = search(reference(), ["SIINFEKL", "AAAAAAAA", "SIINFEKL"], excluded_gene_ids={"CTA"}) + assert result.df.self_gene_id.isna().tolist() == [False, False, True, False, False] + assert result.df.self_gene_id.dropna().tolist() == ["CTA", "healthy", "CTA", "healthy"] + assert result.df.self_search_status.tolist() == ["matched", "matched", "no_match_in_scope", "matched", "matched"] + assert result.df.query_index.tolist() == [0, 0, 1, 2, 2] + + +def test_every_window_position_and_explicit_alleles_survive(): + result = self_matches_in_windows({"repeated": "SIINFEKLSIINFEKL", "short": "AC"}, reference(), + peptide_lengths=[8, 9], alleles={"repeated": ["A0201", "B0702"]}) + matched = result.df[result.df.self_peptide.notna()] + assert matched.peptide_offset.tolist() == [0] * 4 + [8] * 4 + assert set(matched.allele) == {"HLA-A*02:01", "HLA-B*07:02"} + unavailable = result.df[result.df.peptide.str.len().eq(9)] + assert unavailable.self_search_reason.eq("reference_length_unavailable").all() + assert result.extra["self_window_coverage"][-1] == dict( + window_id="short", peptide_length=9, n_occurrences=0, reason="window_shorter_than_peptide_length") + + +def test_observed_presented_and_predicted_facts_do_not_become_recognition(): + obs, pred = observations(), predictions() + before = repr((obs, pred)) + result = match_self_peptides(reference(), ["SIINFEKL", "SIINFEKL", "SIINFEKM", "SIINFEKL"], + alleles=["A0201", "B0702", "A0201", None], observations=obs, predictions=pred) + by_query = result.df.groupby("query_index", sort=False).first() + assert by_query.self_observed.tolist() == [True] * 4 + assert bool(by_query.self_observed_confirmed_same_allele.iloc[0]) + assert not bool(by_query.self_observed_confirmed_same_allele.iloc[1]) + assert bool(by_query.self_observed_predicted_same_allele.iloc[1]) + # A candidate genotype does not establish restriction to this allele. + assert not bool(by_query.self_observed_confirmed_same_allele.iloc[2]) + assert pd.isna(by_query.self_observed_confirmed_same_allele.iloc[3]) + assert by_query.self_prediction_status.tolist() == [ + "reported_for_allele", "not_reported_for_allele", "not_reported_for_allele", "query_allele_unknown"] + assert [row["value"] for row in by_query.self_same_allele_predictions.iloc[0]] == [0., 50.] + assert by_query.self_same_allele_predictions.iloc[1] == [] + assert all(row["predictor_version"] is None for row in by_query.self_predictions.iloc[0]) + assert "gene_id" not in by_query.self_observations.iloc[0][0] # no inferred gene attribution + assert "recognition" not in " ".join(result.df.columns) + assert repr((obs, pred)) == before + + +def test_gene_attributed_observations_do_not_leak_to_another_origin(): + obs = [dict(observations()[0], gene_id="healthy")] + result = match_self_peptides(reference(), ["SIINFEKL"], alleles=["A0201"], observations=obs) + assert result.df.self_observation_status.tolist() == ["not_reported", "reported"] + assert pd.isna(result.df.self_observed.iloc[0]) + + +def test_haplotype_presenter_is_not_same_allele_binding_evidence(): + pred = [dict(predictions()[0], kind="pMHC_presentation", allele_set="A0201,B0702", + prediction_mhc_dependence="haplotype")] + result = match_self_peptides(reference(), ["SIINFEKL"], alleles=["A0201"], predictions=pred) + assert result.df.self_prediction_status.eq("not_reported_for_allele").all() + assert len(result.df.self_predictions.iloc[0]) == 1 + assert result.df.self_predictions.iloc[0][0]["allele_set"] == "HLA-A*02:01,HLA-B*07:02" + + +def test_mouse_alleles_and_missing_values_are_explicit(): + obs = [dict(observations()[0], allele="H-2-Kb")] + result = match_self_peptides(reference(), ["SIINFEKL", "AAAAAAAA", "SIINFEKLL", "SIINFEKX"], + alleles=["H-2-K*b", None, None, None], observations=obs, predictions=[]) + assert result.df.self_observed_confirmed_same_allele.iloc[0] + assert result.df.self_search_status.tolist() == ["matched", "matched", "no_match_in_scope", "unassessed", "unassessed"] + assert result.df.self_search_reason.iloc[-1] == "query_sequence_unsupported" + assert result.df.self_observation_status.iloc[-1] == "not_reported" + + +def test_incomplete_reference_cannot_produce_a_complete_negative(): + ref = SelfProteome.from_peptides({"unknown": "XIINFEKL"}, peptide_lengths=[8]) + for radius in (0, 1): + result = match_self_peptides(ref, ["SIINFEKL"], max_mismatches=radius) + assert result.df.self_search_status.tolist() == ["unassessed"] + assert result.df.self_search_reason.tolist() == ["unsupported_reference_sequences"] + assert result.extra["self_search"]["reference_coverage"]["8"]["n_unsupported"] == 1 + + +def test_search_finds_a_candidate_beyond_the_first_reference_chunk(): + ref = SelfProteome.from_peptides({"background": "AAAAAAAA", "hit": "SIINFEKL"}, peptide_lengths=[8]) + ref._reference_arrays[8] = np.concatenate([np.repeat(ref._reference_arrays[8][:1], 65536, axis=0), + ref._reference_arrays[8][1:]]) + ref._reference_peptides[8] = ["AAAAAAAA"] * 65536 + ["SIINFEKL"] + result = match_self_peptides(ref, ["SIINFEKM"], max_mismatches=1) + assert result.df.self_peptide.tolist() == ["SIINFEKL"] + assert result.df.self_mismatches.tolist() == [1] + + +def test_search_evidence_roundtrip_and_identity(tmp_path, pandas_string_inference): + result = match_self_peptides(reference(), ["SIINFEKL", "SIINFEKM", "AAAAAAAA"], max_mismatches=1, + alleles=["A0201"] * 3, observations=observations(), predictions=predictions()) + path = tmp_path / "self.tsv" + result.to_tsv(path) + restored = read_tsv(path) + # Readers add the input filename as source; every original fact survives. + assert restored.df.source.eq("self.tsv").all() + pd.testing.assert_frame_equal(restored.df[result.df.columns].convert_dtypes(), + result.df.convert_dtypes(), check_dtype=False) + assert restored.extra == result.extra + changed = match_self_peptides(reference(), ["SIINFEKL"], observations=[dict(observations()[0], tissue="kidney")]) + assert changed.extra["self_search"]["observation_sha256"] != result.extra["self_search"]["observation_sha256"] + assert match_self_peptides(reference(), [], observations=[]).empty + assert self_matches_in_windows({}, reference(), peptide_lengths=[8]).empty + + +@pytest.mark.parametrize("options,match", [ + ({"max_mismatches": -1}, "nonnegative"), + ({"max_mismatches": True}, "nonnegative"), + ({"alleles": []}, "one allele"), + ({"observations": observations() * 2}, "unique"), + ({"observations": [dict(observations()[0], allele=None)]}, "assigned"), + ({"observations": [dict(observations()[0], allele_assignment="certain") ]}, "Unknown"), +]) +def test_invalid_evidence_and_search_settings_fail(options, match): + with pytest.raises(ValueError, match=match): + match_self_peptides(reference(), ["SIINFEKL"], **options) diff --git a/tests/test_twin_conformance.py b/tests/test_twin_conformance.py index c54775f..5b4d26d 100644 --- a/tests/test_twin_conformance.py +++ b/tests/test_twin_conformance.py @@ -51,6 +51,7 @@ SelectionPolicy, resolve_selection_policy, evaluate_selection_policy, EvalContext, evaluate_filter, predict_peptide_occurrences, + match_self_peptides, self_matches_in_windows, Affinity, Column, apply_filter, apply_sort, ) from topiary.io_isovar import _check_isovar @@ -61,6 +62,18 @@ from topiary.cli.sv_interest import main as sv_interest_main +SELF_MATCH_TWINS = ( + match_self_peptides, + lambda reference, peptides, **options: reference.match_candidates(peptides, **options), +) +SELF_WINDOW_TWINS = ( + lambda reference, peptides, **options: match_self_peptides(reference, peptides, **options), + lambda reference, peptides, **options: self_matches_in_windows( + {str(index): peptide for index, peptide in enumerate(peptides)}, reference, + peptide_lengths=[len(peptides[0])], **options), +) + + @dataclass(frozen=True) class Twin: """Two callables that must behave identically on the args they share. diff --git a/topiary/__init__.py b/topiary/__init__.py index 4f41bd9..a269121 100644 --- a/topiary/__init__.py +++ b/topiary/__init__.py @@ -141,6 +141,7 @@ unique_fragments, ) from .self_proteome import SelfProteome +from .self_evidence import match_self_peptides, self_matches_in_windows from .io import Metadata, read_csv, read_tsv, to_csv, to_tsv from .io_protein_fragment import read_fragments, write_fragments, iter_fragments from .io_lens import detect_lens_version, read_lens @@ -204,7 +205,7 @@ encode_amino_acids, ) -__version__ = "5.91.0" +__version__ = "5.92.0" __all__ = [ "normalize_isovar_rna_support", @@ -266,6 +267,8 @@ "PredictorSetupError", "mhcflurry_composite_version", "SelfProteome", + "match_self_peptides", + "self_matches_in_windows", "Affinity", "BinOp", "BoolOp", diff --git a/topiary/self_evidence.py b/topiary/self_evidence.py new file mode 100644 index 0000000..83a9c56 --- /dev/null +++ b/topiary/self_evidence.py @@ -0,0 +1,287 @@ +"""Bounded self-sequence searches and supplied presentation evidence.""" + +import hashlib +import json +from collections import defaultdict + +import mhcgnomes +import numpy as np +import pandas as pd + +from .amino_acids import AMINO_ACIDS, encode_amino_acids +from .ranking import mhc_dependence +from .result import TopiaryResult +from .serialization import normalize_python_types + + +def match_self_peptides(reference, peptides, *, max_mismatches=0, alleles=None, + excluded_gene_ids=(), observations=None, predictions=None): + """Return every bounded self match, its origins and supplied evidence. + + Parameters + ---------- + reference : SelfProteome + The explicit corpus searched. Use an unfiltered corpus plus + ``excluded_gene_ids`` to retain both CTA and non-CTA origins of shared + peptides. A filtered corpus cannot recover omitted origins. + peptides : iterable of str + Query sequences, retaining order and duplicates via ``query_index``. + max_mismatches : int + Nonnegative same-length Hamming radius; zero means exact matching. + Indels and TCR recognition are not evaluated. Reference comparisons + are chunked, without a query-by-whole-proteome distance matrix. + alleles : iterable of str or None, optional + One allele per query, parsed with mhcgnomes. Missing alleles remain + unknown; no genotype is borrowed from predictions or other queries. + excluded_gene_ids : iterable of str + Caller-resolved reference exclusions, for example from oncoref. These + mark origins out of scope but never discard them. Topiary owns no CTA + membership or tissue-reference catalog. + observations : DataFrame or iterable of mappings, optional + Supplied records with unique ``evidence_id``, ``peptide``, ``source``, + ``evidence_kind`` (observed/predicted), and ``allele_assignment`` + (confirmed/predicted/unknown). Confirmed or predicted assignments need + an ``allele``; an unresolved ``allele_set`` is not confirmed restriction. + Optional gene_id restricts attribution to that origin. Tissue, assay, + source version and other JSON-compatible facts survive unchanged. + predictions : DataFrame or iterable of mappings, optional + Long-form records with peptide, allele, kind, value, + prediction_method_name and predictor_version. Unknown provenance stays + unknown. All same-peptide records survive, including conflicts; only + finite, per-allele measurements at the query allele count toward + same-allele coverage. No prediction or model selection is run here. + + Returns + ------- + TopiaryResult + One row per query/matched-sequence/reference occurrence. No-hit and + unassessed queries retain a row with a missing self peptide. Search + status, reason and completeness are separate from observation and + prediction coverage. ``self_in_scope`` marks non-excluded origins. + Structured observation/prediction cells and complete search/evidence + identity round-trip through ordinary Topiary CSV/TSV. + + Observation flags describe the supplied records, not biological + absence. Unreported observations remain unknown. Peptide observations + without gene attribution do not establish which matching gene produced + them. Presentation, similarity and binding do not establish TCR + recognition. ``no_match_in_scope`` is bounded negative search evidence, + never a statement of no risk. Empty input returns an empty schema. + """ + if isinstance(peptides, str): + raise TypeError("peptides must be an iterable of sequences, not one string") + peptides = list(peptides) + if any(not isinstance(peptide, str) or not peptide for peptide in peptides): + raise ValueError("Query peptides must be nonempty strings") + if type(max_mismatches) is not int or max_mismatches < 0: + raise ValueError("max_mismatches must be a nonnegative integer") + if isinstance(alleles, str): + raise TypeError("alleles must supply one allele per query") + alleles = [None] * len(peptides) if alleles is None else list(alleles) + if len(alleles) != len(peptides): + raise ValueError("alleles must supply one allele per query") + alleles = [_allele(value) for value in alleles] + if isinstance(excluded_gene_ids, str): + raise TypeError("excluded_gene_ids must be an iterable of gene IDs") + excluded = sorted(set(excluded_gene_ids)) + if any(not isinstance(gene, str) or not gene for gene in excluded): + raise ValueError("Excluded gene IDs must be nonempty strings") + observation_records = _evidence_records(observations, observation=True) + prediction_records = _evidence_records(predictions, observation=False) + by_observation, by_prediction = defaultdict(list), defaultdict(list) + for row in observation_records: + by_observation[row["peptide"]].append(row) + for row in prediction_records: + by_prediction[row["peptide"]].append(row) + supported = set(AMINO_ACIDS) + excluded = set(excluded) + reference_coverage = {} + for length in sorted({len(peptide) for peptide in peptides}): + array = reference._reference_arrays.get(length, []) + unsupported = sum(int((array[start:start + 65536] >= len(AMINO_ACIDS)).any(axis=1).sum()) + for start in range(0, len(array), 65536)) + reference_coverage[str(length)] = dict(n_peptides=len(array), n_unsupported=unsupported) + searches, rows = {}, [] + for query_index, (peptide, allele) in enumerate(zip(peptides, alleles)): + if peptide not in searches: + searches[peptide] = _search(reference, peptide, max_mismatches, supported, + reference_coverage[str(len(peptide))]["n_unsupported"]) + matches, complete, reason = searches[peptide] + status = "matched" if matches else "no_match_in_scope" if complete else "unassessed" + origins = [(match, distance, origin) for match, distance in matches + for origin in reference._provenance.get(match, [(None, None, None)])] + for match, distance, (gene, transcript, offset) in origins or [(None, None, (None, None, None))]: + observed = [record for record in by_observation[match] + if record.get("gene_id") is None or record["gene_id"] == gene] + predicted = by_prediction[match] + same_allele = [record for record in predicted if allele is not None + and record["allele"] == allele and record["prediction_mhc_dependence"] == "single_allele" + and isinstance(record["value"], (int, float)) and not isinstance(record["value"], bool)] + measured = [record for record in observed if record["evidence_kind"] == "observed"] + identity = None if match is None else _digest([reference.reference_version, match, gene, transcript, offset]) + rows.append(dict( + query_index=query_index, peptide=peptide, peptide_length=len(peptide), allele=allele, + self_match_id=identity, self_peptide=match, self_mismatches=distance, + self_gene_id=gene, self_transcript_id=transcript, self_reference_offset=offset, + self_in_scope=None if gene is None else gene not in excluded, + self_reference_version=reference.reference_version, + self_search_status=status, self_search_complete=complete, self_search_reason=reason, + self_observation_status="not_supplied" if observations is None else "reported" if observed else "not_reported", + self_observed=bool(measured) if observed else None, + self_observed_confirmed_same_allele=(any(record["allele_assignment"] == "confirmed" + and record.get("allele") == allele for record in measured) if observed and allele else None), + self_observed_predicted_same_allele=(any(record["allele_assignment"] == "predicted" + and record.get("allele") == allele for record in measured) if observed and allele else None), + self_prediction_status=("not_supplied" if predictions is None else "query_allele_unknown" + if allele is None else "reported_for_allele" if same_allele else "not_reported_for_allele"), + self_observations=observed, self_predictions=predicted, self_same_allele_predictions=same_allele, + )) + columns = ["query_index", "peptide", "peptide_length", "allele", "self_match_id", "self_peptide", "self_mismatches", + "self_gene_id", "self_transcript_id", "self_reference_offset", "self_in_scope", + "self_reference_version", "self_search_status", "self_search_complete", "self_search_reason", + "self_observation_status", "self_observed", "self_observed_confirmed_same_allele", + "self_observed_predicted_same_allele", "self_prediction_status", "self_observations", + "self_predictions", "self_same_allele_predictions"] + scope = dict(reference_version=reference.reference_version, search="same_length_hamming", + reference_coverage=reference_coverage, + max_mismatches=max_mismatches, include_indels=False, excluded_gene_ids=sorted(excluded), + observation_input="not_supplied" if observations is None else "supplied", + observation_sha256=_digest(observation_records), prediction_input="not_supplied" + if predictions is None else "supplied", prediction_sha256=_digest(prediction_records)) + return TopiaryResult(pd.DataFrame(rows, columns=columns), extra={"self_search": scope}) + + +def self_matches_in_windows(windows, reference, *, peptide_lengths, alleles=None, + excluded_gene_ids=(), observations=None, predictions=None): + """Find exact reference peptides at every requested position in windows. + + Parameters + ---------- + windows : mapping of str to str + Window ID to complete proposed sequence. IDs and sequences must be + nonempty strings. Repeated peptides and every source offset survive. + reference : SelfProteome + Corpus passed to :func:`match_self_peptides` without further filtering. + peptide_lengths : iterable of int + Explicit positive lengths to enumerate. Missing reference lengths + remain unassessed; they are never silently removed from the request. + alleles : mapping of str to iterable of str, optional + Explicit alleles for each window. An absent/empty entry produces + unknown allele context; another window's genotype is never borrowed. + excluded_gene_ids, observations, predictions + Passed to :func:`match_self_peptides` with its evidence semantics. + + Returns + ------- + TopiaryResult + Exact-match evidence with window_id, window_sequence and zero-based + peptide_offset. The peptide columns denote window occurrences, not + administered molecules for a cleavage model. Metadata also records + requested window/length combinations too short to contain a peptide. + Empty windows input returns an empty schema. No predictor runs and no + window is chosen, trimmed or declared safe by this function. + """ + windows = dict(windows) + lengths = list(peptide_lengths) + if not lengths or any(type(length) is not int or length <= 0 for length in lengths): + raise ValueError("peptide_lengths must contain positive integers") + lengths = sorted(set(lengths)) + alleles = {} if alleles is None else dict(alleles) + if set(alleles) - set(windows): + raise ValueError("Alleles name an unknown window") + queries, coverage = [], [] + for name, sequence in windows.items(): + if not isinstance(name, str) or not name or not isinstance(sequence, str) or not sequence: + raise ValueError("Window IDs and sequences must be nonempty strings") + named = alleles.get(name, ()) + if isinstance(named, str): + raise TypeError("Window alleles must be an iterable, not one string") + genotype = list(dict.fromkeys(_allele(allele) for allele in named)) or [None] + for length in lengths: + count = max(0, len(sequence) - length + 1) + coverage.append(dict(window_id=name, peptide_length=length, n_occurrences=count, + reason=None if count else "window_shorter_than_peptide_length")) + for offset in range(count): + for allele in genotype: + queries.append(dict(window_id=name, window_sequence=sequence, peptide_offset=offset, + peptide=sequence[offset:offset + length], allele=allele)) + result = match_self_peptides(reference, [row["peptide"] for row in queries], + alleles=[row["allele"] for row in queries], excluded_gene_ids=excluded_gene_ids, + observations=observations, predictions=predictions) + frame = result.df.copy() + for column in ("window_id", "window_sequence", "peptide_offset"): + frame[column] = [queries[index][column] for index in frame.query_index] + return TopiaryResult(frame, extra={**result.extra, "self_window_coverage": coverage}) + + +def _allele(value): + return None if value is None or (pd.api.types.is_scalar(value) and pd.isna(value)) or value == "" else mhcgnomes.parse(str(value)).to_string() + + +def _digest(value): + return hashlib.sha256(json.dumps(value, sort_keys=True, separators=(",", ":"), allow_nan=False).encode()).hexdigest() + + +def _evidence_records(data, *, observation): + if data is None: + return [] + frame = pd.DataFrame(data).copy() + if not frame.columns.is_unique: + raise ValueError("Evidence columns must be unique") + if frame.empty: + return [] + required = {"peptide", "evidence_id", "source", "evidence_kind", "allele_assignment"} if observation else { + "peptide", "allele", "kind", "value", "prediction_method_name", "predictor_version"} + if required - set(frame): + raise ValueError(f"Missing self evidence columns: {sorted(required - set(frame))}") + records = [] + for raw in frame.to_dict("records"): + row = {key: None if pd.api.types.is_scalar(value) and pd.isna(value) else normalize_python_types(value) + for key, value in raw.items()} + if not isinstance(row["peptide"], str) or not row["peptide"]: + raise ValueError("Evidence peptide must be a nonempty string") + row["allele"] = _allele(row.get("allele")) + if row.get("allele_set"): + from .ranking import format_allele_set, split_allele_set + row["allele_set"] = format_allele_set(_allele(name) for name in split_allele_set(row["allele_set"])) + if observation: + if any(not isinstance(row[key], str) or not row[key] for key in ("evidence_id", "source")): + raise ValueError("Observation evidence_id and source must be nonempty strings") + if row["evidence_kind"] not in {"observed", "predicted"} or row["allele_assignment"] not in {"confirmed", "predicted", "unknown"}: + raise ValueError("Unknown observation evidence kind or allele assignment") + if row["allele_assignment"] != "unknown" and not row["allele"]: + raise ValueError("An assigned observation allele must be supplied explicitly") + else: + dependence = row.get("prediction_mhc_dependence") or mhc_dependence(row["kind"], rows=pd.DataFrame([row])) + if dependence not in {"single_allele", "haplotype", "none"}: + raise ValueError("Unknown prediction MHC dependence") + row["prediction_mhc_dependence"] = dependence + _digest(row) # Reject non-JSON facts rather than silently stringifying them. + records.append(row) + if observation and len({row["evidence_id"] for row in records}) != len(records): + raise ValueError("Observation evidence_id must be unique") + return records + + +def _search(reference, peptide, radius, supported, unsupported): + if not set(peptide) <= supported: + return [], False, "query_sequence_unsupported" + length = len(peptide) + candidates = reference._reference_peptides.get(length, []) + if not candidates: + return [], False, "reference_length_unavailable" + complete, reason = not unsupported, "unsupported_reference_sequences" if unsupported else None + if radius == 0: + # Exact matching requires no sequence-model assumption or distance matrix. + return ([(peptide, 0)] if peptide in reference._reference_set(length) else []), complete, reason + query = encode_amino_acids([peptide], length)[0] + array = reference._reference_arrays[length] + matches = [] + for start in range(0, len(array), 65536): + distances = (array[start:start + 65536] != query).sum(axis=1) + for index in np.flatnonzero(distances <= radius): + candidate = candidates[start + int(index)] + if not set(candidate) <= supported: + continue + matches.append((candidate, int(distances[index]))) + return matches, complete, reason diff --git a/topiary/self_proteome.py b/topiary/self_proteome.py index a37b6e3..7101aff 100644 --- a/topiary/self_proteome.py +++ b/topiary/self_proteome.py @@ -34,10 +34,10 @@ edit_distance=1 beats a same-length substitution match at edit_distance≥2. -Binding-aware axes (``self_mimic_*``, ``self_strongest_nearby_*``, -``self_nearest_candidates``) are not implemented yet — they require -MHC prediction on candidate peptides, which is architecturally -separate from the sequence-only ``nearest()`` method. Tracked under +``match_candidates()`` retains all candidates within a same-length Hamming +radius, their origins and supplied observation/prediction evidence. Binding-ranked +axes (``self_mimic_*``, ``self_strongest_nearby_*``) still require an explicit +model-selection and ranking contract, separate from sequence matching. See `#412 `_. Algorithm @@ -179,6 +179,36 @@ def reference_version(self) -> str: # --- lookup --- + def match_candidates(self, peptides, *, max_mismatches=0, alleles=None, + excluded_gene_ids=(), observations=None, predictions=None): + """Return all bounded self matches with origins and supplied evidence. + + Parameters + ---------- + peptides : iterable of str + Query sequences, preserving order and repeats. + max_mismatches : int + Nonnegative same-length Hamming radius. Defaults to exact matches. + alleles : iterable of str or None, optional + One explicit allele per query; missing remains unknown. + excluded_gene_ids : iterable of str + Caller-resolved exclusions; origins survive with an out-of-scope flag. + observations, predictions : DataFrame or iterable of mappings, optional + Supplied evidence with the schema of :func:`topiary.match_self_peptides`. + + Returns + ------- + TopiaryResult + All matches, source origins, supplied evidence and coverage reasons, + through the shared :func:`topiary.match_self_peptides` implementation. + Empty input returns an empty schema. A missing match never establishes + absence of risk. Structured evidence survives ordinary CSV/TSV IO. + """ + from .self_evidence import match_self_peptides + return match_self_peptides( + self, peptides, max_mismatches=max_mismatches, alleles=alleles, + excluded_gene_ids=excluded_gene_ids, observations=observations, predictions=predictions) + def nearest( self, peptides: Iterable[str], @@ -217,9 +247,9 @@ def nearest( hash-set lookup. The **first** indel hit found (deletions checked before insertions, in position order) is returned — not necessarily the best among all indel candidates. This is - fast (~200 lookups per 9-mer) but not exhaustive; the upcoming - ``self_nearest_candidates`` structured column will expose the - full candidate set. + fast (~200 lookups per 9-mer) but not exhaustive. ``match_candidates`` + exposes all candidates within a same-length Hamming radius; it does + not enumerate indels. """ if metric not in ("blosum62", "hamming"): raise ValueError(