Score exact peptide occurrences with retained source context

CHANGELOG.md

@@ -7,6 +7,20 @@ between published tags; older pre-5.0 notes are retained below.
 
 For current interfaces, see the [consumer guide](docs/consumer-guide.md).
 
+## 5.91.0
+
+- Score explicitly identified peptide occurrences with their original flanks,
+  coordinates and source evidence, without scanning additional windows. Share
+  identical inference inputs while preserving every observation; retain model,
+  allele-free/per-allele/haplotype scope and explicit comparator context (#367).
+- Reuse occurrence prediction in additive table rescoring, batching compatible
+  peptides while retaining original measurements and selection behavior.
+- Preserve the canonical prediction schema for empty named-peptide output,
+  including explicitly reported all-missing cache batches (#450).
+- Match fragment and explicit-occurrence WT scores by declared MHC scope, so
+  a different haplotype presenter does not erase the comparator score. Export
+  `prediction_mhc_scope` for consumers sharing that identity rule (#453).
+
 ## 5.90.0
 
 - Evaluate saved policies at explicit occurrence/allele identity before selecting

docs/api.md

@@ -91,11 +91,24 @@ return fresh read-only int8 arrays, so callers never share mutable matrix state.
 |--------|-------|----------|
 | `predict_from_named_sequences(dict)` | `{name: sequence}` | Sliding-window scan |
 | `predict_from_named_peptides(dict)` | `{name: peptide}` | Score as-is |
+| `predict_from_peptide_occurrences(table)` | `prediction_id`, peptide, optional flanks/coordinates/annotations; unique by ID, peptide, offset and sample | Score exact occurrences with their context; no scanning. |
 | `predict_from_sequences(list)` | `[sequence, ...]` | Sliding-window scan |
 | `predict_from_fragments(fragments)` | `[ProteinFragment]` | Universal path — any origin; fragment-level metadata and `target_intervals` threaded through. |
 | `predict_from_variants(variants)` | VariantCollection | Variant pipeline (builds `ProteinFragment`s internally and delegates). |
 | `predict_from_mutation_effects(effects)` | EffectCollection | Same as `predict_from_variants` but starting from pre-computed effects. |
 
+`predict_peptide_occurrences(table, model, use_flanks=True, predict_wt=False)`
+exposes the same occurrence prediction for one configured mhctools model.
+It is exported from `topiary` and is also used by `rescore_candidates`.
+See [exact occurrence prediction](combined-sources.md#predict-exact-peptide-occurrences)
+for input identity, unknown flanks, comparator and coverage rules.
+
+`prediction_mhc_scope(allele, dependence=..., allele_set=...)` returns a canonical
+MHC identity tuple: one allele, the full genotype, or an allele-free scope.
+Missing required identity returns `None`, which must remain unmatchable.
+Comparator joins use this public rule in both fragment and occurrence prediction;
+a haplotype's selected presenter is not its genotype identity.
+
 ## TopiaryResult
 
 `TopiaryResult` is the semantic result object for Topiary prediction tables.

docs/combined-sources.md

@@ -534,6 +534,70 @@ after reloading. The composed workflow test covers original and additive scoring
 retained hypothesis selection, and count union without running a predictor during
 assembly.
 
+## Predict exact peptide occurrences
+
+Topiary 5.91.0 adds `predict_peptide_occurrences` and
+`TopiaryPredictor.predict_from_peptide_occurrences` for a selected peptide
+universe. Supply one record per `(prediction_id, peptide, peptide_offset)` within
+each sample (when `sample_name` or `candidate_sample` is supplied). An ID may name
+several peptide windows in the same source. Use the same occurrence
+identity Vaxrank and the policy evaluator consume. Different occurrences of one
+peptide remain separate even when they share a gene, coordinate or HLA candidate.
+
+```python
+import pandas as pd
+from topiary import (
+    TopiaryPredictor, TopiaryResult, SelectionPolicy, evaluate_selection_policy,
+)
+
+# model is an explicitly configured mhctools predictor with patient alleles.
+occurrences = pd.DataFrame([
+    dict(prediction_id="gene-a:12", peptide="SIINFEKL", peptide_offset=12,
+         n_flank="AAA", c_flank="GGG", gene="gene-a"),
+    dict(prediction_id="gene-b:20", peptide="SIINFEKL", peptide_offset=20,
+         n_flank="TTT", c_flank="CCC", gene="gene-b"),
+])
+predictor = TopiaryPredictor(models=model)
+predictions = predictor.predict_from_peptide_occurrences(occurrences)
+evaluation = evaluate_selection_policy(
+    TopiaryResult(predictions), SelectionPolicy("binding", "1 / affinity.value"),
+    group_keys=["prediction_id", "peptide", "peptide_offset", "allele"],
+    kind_support=predictor.kind_support,
+)
+evaluation.evidence.to_tsv("occurrence-evidence.tsv")
+```
+
+The model scores only the supplied peptides against its configured alleles.
+Within a shared sample/flank/genotype context, distinct peptides are batched and
+identical inference inputs are scored once, then linked to each original
+occurrence. Source annotations and read support are copied, never summed.
+`prediction_mhc_dependence` records each kind's scope; haplotype predictions carry
+the configured `allele_set`. An explicitly supplied haplotype set must match the
+configured model. Missing kind/allele coverage and ambiguous outputs raise.
+
+Flank-dependent models require both `n_flank` and `c_flank`. An empty string is a
+known terminus; a missing value is unknown. `use_flanks=False` explicitly requests
+peptide-only inference while preserving the original context. The output's
+`prediction_flanks_supplied` records whether flanks were supplied in the model call; it does not imply that every returned prediction kind depends on them.
+Missing source offsets remain missing. Model peptide-length/domain checks still
+apply; no padding, comparator sequence, or additional peptide window is inferred.
+
+With `TopiaryPredictor(predict_wt=True)` (or `predict_wt=True` on the standalone
+function), supplied `wt_peptide` values are scored using their own
+`wt_n_flank` / `wt_c_flank`. Flank-dependent comparator prediction requires those
+flanks or explicit `use_flanks=False`; primary-peptide context is not borrowed.
+Absent comparators keep missing scores. Explicit filter/sort settings operate on
+occurrences; `only_novel_epitopes` does not reinterpret this already selected
+peptide universe. Cache failures are strict by default; an explicit
+`cache_miss_handler` / `on_miss` retains the existing reported-partial-batch
+contract, with occurrence IDs in the failure records.
+
+Prediction inputs contain context and annotations, not historical measurement
+columns. Keep imported measurements in their original table and use the additive
+operation below. Context-aware live prediction requires mhctools' flank-aware
+API; sparse historical tables with unknown model versions remain the separate
+prediction-source work in [#368](https://github.com/openvax/topiary/issues/368).
+
 ## Add prediction features on demand
 
 ```python
@@ -564,6 +628,9 @@ models, versions, measurement semantics, UTC creation time, selection and discov
 Discovery source and prediction producer are separate facts.
 
 Configured mhctools models supply `predict_dataframe` and `kind_support()`.
+This path reuses `predict_peptide_occurrences` to batch compatible requests and
+preserve source identity; it appends primary-peptide features and leaves all
+historical comparator measurements intact.
 Supplied flanks are forwarded by default; flank-dependent models require both
 flanks or explicit `use_flanks=False`. Haplotype models require a matching
 stated `allele_set`. Existing candidate identities are preserved; ORF-only rows

scripts/check_vaxrank_candidates.py

@@ -22,6 +22,47 @@
 from vaxrank.vaccine_peptide import VaccinePeptide
 
 
+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
+
+    sequence = "AAASIINFEKLTSIINFEKLTTGILGFVFTL"
+    windows = pd.DataFrame(dict(
+        prediction_id=["one-protein"] * 3,
+        peptide=["SIINFEKL", "SIINFEKL", "GILGFVFTL"],
+        peptide_offset=[3, 12, 22], n_flank=["AAA", "T", "TT"], c_flank=["T", "TT", ""],
+        source_sequence=[sequence] * 3,
+    ))
+    observed = []
+    for use_flanks, expected in ((True, [803., 801., 12.]), (False, [800., 800., 10.])):
+        frame = TopiaryPredictor(models=Model()).predict_from_peptide_occurrences(windows, use_flanks=use_flanks)
+        assert frame.value.tolist() == expected
+        assert frame.prediction_id.tolist() == windows.prediction_id.tolist()
+        rows = [dict(
+            peptide=row.peptide, source=row.source_sequence, offset=row.peptide_offset,
+            prediction_id=row.prediction_id, n_flank=row.n_flank, c_flank=row.c_flank,
+            mutant=Prediction(kind=row.kind, peptide=row.peptide, allele=row.allele,
+                              value=row.value, score=row.score,
+                              predictor_name=row.prediction_method_name, predictor_version=row.predictor_version),
+            source_class="mutation", overlaps_targetable=True, patient_alleles=[row.allele],
+        ) for row in frame.itertuples()]
+        epitopes = candidate_epitopes_from_rows(rows)
+        policy = SelectionPolicy("explicit-windows", "1 / affinity.value")
+        cfg = EpitopeConfig(score_expr=policy.score_by, min_epitope_score=0.)
+        scored = attach_per_allele_scores(epitopes, cfg, topiary_df=frame)
+        keys = prediction_group_columns(frame)
+        evaluated = evaluate_selection_policy(TopiaryResult(frame), policy, group_keys=keys)
+        scores = evaluated.occurrences.set_index(keys).score.to_dict()
+        assert len(scored) == 3
+        for epitope, value in zip(scored, expected):
+            assert epitope.prediction_id == "one-protein"
+            score = epitope.per_allele_scores["HLA-A*02:01"]
+            assert score == pytest.approx(1 / value)
+            assert score == scores[(*epitope.prediction_group_key, "HLA-A*02:01")]
+        observed.append([epitope.per_allele_scores for epitope in scored])
+    assert observed[0] != observed[1]
+
+
 @pytest.mark.parametrize("antigen_kind", ["mutation", "fusion", "splice", "CTA", "ERV", "viral"])
 @pytest.mark.parametrize("policy", ["original", "rescored", "rna_overlay"])
 def test_candidate_features_reach_vaxrank_scoring_and_vaccine_construction(antigen_kind, policy, tmp_path):

tests/test_candidate_tables.py

@@ -123,7 +123,7 @@ def test_rescoring_only_selected_subset_preserves_original_ranking_and_adds_dsl_
     model = Model()
     enriched = rescore_candidates(output, model, prefix="fresh", select="source_label == 'one'")
     feature = "fresh__testmodel__pMHC_affinity__value"
-    assert [call[0] for call in model.calls] == [["SIINFEKL"], ["GILGFVFTL"]]
+    assert [call[0] for call in model.calls] == [["SIINFEKL", "GILGFVFTL"]]
     assert enriched.df.loc[enriched.df.source_label.eq("two"), feature].isna().all()
     pd.testing.assert_frame_equal(enriched.df[before.columns], before)
     pd.testing.assert_frame_equal(output.df, before)
@@ -149,7 +149,7 @@ def test_contexts_are_not_reused_and_flanks_reach_the_predictor():
     output = combined(n_flank="OTHER")
     model = Model()
     scored = rescore_candidates(output, model, prefix="new")
-    assert len(model.calls) == 4
+    assert len(model.calls) == 2
     assert {c[1]["n_flanks"][0] for c in model.calls} == {"AAA", "OTHER"}
     assert scored.df.new__testmodel__pMHC_affinity__value.tolist() == [803, 13, 805, 15]
 
@@ -158,7 +158,7 @@ def test_deduplicated_calls_keep_distinct_observations_and_evidence():
     output = combined()
     model = Model()
     scored = rescore_candidates(output, model, prefix="new")
-    assert len(model.calls) == 2
+    assert len(model.calls) == 1
     assert len(scored.df) == len(output.df)
     assert scored.df.source_observation_id.nunique() == 4
 
@@ -327,7 +327,7 @@ def test_selective_rescoring_does_not_leak_across_genotypes():
     feature = "new__testmodel__pMHC_affinity__value"
     assert scored.df[feature].iloc[:2].notna().all()
     assert scored.df[feature].iloc[2:].isna().all()
-    assert len(model.calls) == 2
+    assert len(model.calls) == 1
 
 
 def test_conflicting_predictions_need_source_or_run_identity_before_dsl_evaluation():
@@ -394,15 +394,17 @@ def kind_support(self):
 
     def predict_dataframe(self, peptides, **kwargs):
         base = super().predict_dataframe(peptides, **kwargs)
-        frames = []
+        records = []
         for kind in self.emitted:
             frame = base.copy()
             frame["kind"] = kind
             if kind == "antigen_processing":
                 frame["allele"] = None
                 frame["value_unit"] = None
-            frames.append(frame)
-        return pd.concat(frames, ignore_index=True)
+            records.extend(frame.to_dict("records"))
+        # Infer one dtype across both allele-bearing and allele-free rows,
+        # including pandas' inferred-string mode (#454).
+        return pd.DataFrame(records)
 
 
 @pytest.mark.parametrize("alleles", [["HLA-B*07:02"] * 2, ["HLA-A*02:01", "HLA-B*07:02"]])

tests/test_consumer_workflows.py

@@ -277,6 +277,85 @@ def test_peptide_level_measurement_conflicts_are_rejected_before_projection():
 )
 
 
+@pytest.mark.parametrize("extension", ["csv", "tsv"])
+def test_explicit_occurrences_from_mixed_sources_rescore_and_replay(tmp_path, monkeypatch, extension):
+    from dataclasses import replace
+    from topiary import (
+        ProteinFragment, SelectionPolicy, combine_sources, evaluate_selection_policy,
+        predict_peptide_occurrences, read_csv, read_tsv, rescore_candidates, replay_selection_policy,
+        select_policy_representatives,
+    )
+    from .test_candidate_tables import Model
+
+    fragments = [
+        ProteinFragment(fragment_id="first", sequence="MAAASIINFEKLAAA", gene="G1",
+                        reference_sequence="MAAAGIINFEKLAAA"),
+        ProteinFragment(fragment_id="second", sequence="MSIINFEKLGGG", gene="G2"),
+    ]
+    inputs = pd.DataFrame([dict(
+        prediction_id=fragment.fragment_id, fragment_id=fragment.fragment_id,
+        peptide=fragment.sequence[offset:offset + 8], peptide_offset=offset,
+        n_flank=fragment.sequence[max(0, offset - 3):offset], c_flank=fragment.sequence[offset + 8:offset + 11],
+        source_sequence=fragment.sequence, gene=fragment.gene, n_rna_alt=7,
+        wt_peptide=fragment.effective_baseline[offset:offset + 8] if fragment.effective_baseline else None,
+        wt_n_flank=fragment.effective_baseline[max(0, offset - 3):offset] if fragment.effective_baseline else None,
+        wt_c_flank=fragment.effective_baseline[offset + 8:offset + 11] if fragment.effective_baseline else None,
+    ) for fragment, offset in zip(fragments, [4, 1])])
+    original = predict_peptide_occurrences(inputs, Model(), predict_wt=True)
+    assert original.peptide.tolist() == ["SIINFEKL", "SIINFEKL"]
+    assert original.wt_value.iloc[0] == 13. and pd.isna(original.wt_value.iloc[1])
+    # Controlled original scores prefer the first occurrence; contextual fresh
+    # scores prefer the second. Neither operation changes the input peptides.
+    original["value"] = original["affinity"] = [10., 20.]
+    combined = combine_sources({"fragment": original, "lens": read_lens(LENS),
+                                "pvacseq": read_pvacseq(PVACSEQ)}, sample_name="p")
+    before = combined.long_df.copy(deep=True)
+    policy = SelectionPolicy("original", "1 / affinity.value", filter_by="source_label == 'fragment'", duplicates="best")
+    def chosen(evaluation, field):
+        representatives = select_policy_representatives(evaluation)
+        return [evaluation.evidence.df.iloc[rows[0]][field] for rows in representatives.evidence_rows]
+    assert chosen(evaluate_selection_policy(combined, policy), "fragment_id") == ["first"]
+
+    model = Model()
+    enriched = rescore_candidates(combined, model, prefix="fresh", select="source_label == 'fragment'")
+    assert [peptides for peptides, _ in model.calls] == [["SIINFEKL"], ["SIINFEKL"]]
+    pd.testing.assert_frame_equal(enriched.long_df[before.columns], before)
+    assert chosen(evaluate_selection_policy(enriched, policy), "fragment_id") == ["first"]
+    fresh = replace(policy, name="fresh", score_by="1 / fresh__testmodel__pMHC_affinity__value")
+    evaluated = evaluate_selection_policy(enriched, fresh)
+    second = select_policy_representatives(evaluated)
+    assert chosen(evaluated, "fragment_id") == ["second"]
+    assert chosen(evaluated, "n_rna_alt") == [7]
+    assert set(evaluated.evidence.df.source_label) == {"fragment", "lens", "pvacseq"}
+    assert evaluated.evidence.df.wt_peptide.dropna().str.len().gt(0).all()
+
+    # Import-only records can also request exact peptide inference explicitly,
+    # without treating their unknown flanks as known termini.
+    imported = combined.df.loc[combined.df.source_label.ne("fragment")]
+    columns = [column for column in ("source_observation_id", "peptide", "peptide_offset", "n_flank", "c_flank", "wt_peptide")
+               if column in imported]
+    requests = imported[columns].drop_duplicates("source_observation_id").rename(columns={"source_observation_id": "prediction_id"})
+    imported_model = Model()
+    with pytest.raises(ValueError, match="both flanks"):
+        predict_peptide_occurrences(requests, imported_model)
+    assert not imported_model.calls
+    exact = predict_peptide_occurrences(requests, imported_model, use_flanks=False)
+    assert set(exact.peptide) == set(requests.peptide)
+    assert exact.prediction_id.tolist() == requests.prediction_id.tolist()
+    assert sum(len(peptides) for peptides, _ in imported_model.calls) == requests.peptide.nunique()
+    pd.testing.assert_series_equal(exact.wt_peptide, requests.wt_peptide.reset_index(drop=True))
+
+    def forbidden(*args, **kwargs):
+        raise AssertionError("Saved evidence replay must not invoke inference")
+    monkeypatch.setattr(Model, "predict_dataframe", forbidden)
+    path = tmp_path / ("occurrence-policy." + extension)
+    getattr(evaluated.evidence, "to_" + extension)(path)
+    loaded = (read_csv if extension == "csv" else read_tsv)(path)
+    replay = replay_selection_policy(loaded)
+    pd.testing.assert_frame_equal(select_policy_representatives(replay), second, check_exact=True)
+    pd.testing.assert_frame_equal(replay.occurrences, evaluated.occurrences, check_exact=True)
+
+
 @pytest.mark.parametrize("wide", [False, True])
 def test_source_tables_combine_rank_and_export_without_predicting(monkeypatch, tmp_path, wide, pandas_string_inference):
     from topiary import (
@@ -334,7 +413,8 @@ def forbidden(*args, **kwargs):
     model = Model()
     enriched = rescore_candidates(restored, model, prefix="fresh", select="source_label == 'direct'",
                                   use_flanks=False)
-    assert len(model.calls) == len(direct)
+    assert len(model.calls) == 1
+    assert set(model.calls[0][0]) == set(direct.peptide)
     original_again = rank_candidates(enriched, expression, ascending=True, duplicates="best")
     assert original_again.candidate_id.tolist() == pooled.candidate_id.tolist()
     np.testing.assert_allclose(original_again.candidate_score, pooled.candidate_score, equal_nan=True)
@@ -3281,7 +3361,7 @@ def forbidden(*args, **kwargs):
                                   result.df.iloc[2].rna_observations[0]])["value"] == 3
     model = Model()
     scored = rescore_candidates(result, model, prefix="fresh")
-    assert len(model.calls) == 2
+    assert len(model.calls) == 1
     assert scored.df.candidate_id.nunique() == 2
     saved = scored.to_wide() if wide else scored
     path = tmp_path / ("reconciled." + extension)
@@ -3327,7 +3407,7 @@ def forbidden(*args, **kwargs):
     assert matches.iloc[0].candidate_id == ranked.iloc[0].candidate_id
     model = Model()
     rescored = rescore_candidates(combined, model, prefix="fresh", select="source_label == 'reported'")
-    assert len(model.calls) == 2
+    assert len(model.calls) == 1
     saved = rescored.to_wide() if wide else rescored
     path = tmp_path / ("exacto-combined." + extension)
     getattr(saved, "to_" + extension)(path)