Computed self_nearest_* columns against a curated non-CTA proteome

[iskandr]

Motivation

For substitutions, wt_peptide gives a principled 1:1 pair for differential binding. For indels, frameshifts, fusions, viral, and synthetic fragments there's no meaningful WT pair — but cross-reactivity / autoimmune risk is still a concern. The principled generalization is nearest-self: the closest peptide to the mutant in a curated healthy-tissue proteome, paired with that peptide's MHC binding at the same allele.

docs/fragments.md already reserves the self_nearest_* column namespace and DSL scope (self_nearest_peptide, self_nearest_edit_distance, self_nearest_gene, self_nearest_value / _score / _percentile_rank, etc.) but currently instructs producers to populate externally. This issue is about computing those columns inside topiary, opt-in.

Proposal

1. Bundled reference proteome. Curated healthy-tissue human proteome with CTAs (cancer-testis antigens) explicitly excluded — CTAs are legitimate tumor-specific targets whose testis expression shouldn't count as "self" for cross-reactivity. Shipped as a downloadable artifact (too big for the wheel) with a loader.
2. Pluggable matching. Default: exact k-mer index + Hamming fallback (fast, no BLAST dependency). Interface for callers to plug in BLAST or BLOSUM-weighted edit distance.
3. Population path. TopiaryPredictor(compute_self_nearest=True, self_reference=…) — runs the nearest-self lookup after prediction and populates the full reserved namespace. Rows with no hit within a configurable max distance stay NaN; caller decides via fillna or filters.
4. Paired MHC binding. For each self_nearest_peptide, run _predict_raw_peptides at the mutant's allele to populate self_nearest_value / _score / _percentile_rank. Same primitive as the predict_wt path.

Out of scope

• Healthy-tissue expression data (self_nearest_tissues) — separate issue; needs GTEx / HPA integration.
• Non-human reference proteomes.

Why this instead of indel wt_peptide remapping

For non-SNV effects there is no coordinate-aligned "WT peer" that means what wt_peptide means for an SNV. Remapping would give you a residue string, but the biological question ("is this cross-reactive with a healthy self peptide at the same MHC?") is better answered by nearest-self across the whole self proteome than by the single residue span that happens to line up pre-shift.

[iskandr]

Expanding the two open design questions from our earlier discussion with concrete implementation notes.

Scale ceiling — when pure NumPy stops being enough

The naive vectorized approach: encode reference 9-mers as an `(M, 9)` int8 array, process queries in chunks, use broadcasting for `(chunk, M, 9)` difference tensors, sum over the last axis for Hamming. Memory is `chunk * M * 9` bytes. For `M = 200k` 9-mers and `chunk = 1000`, that's ~1.7GB per chunk — fine. Below are the rough knee points:

Workload	Shape	Pure NumPy?	Notes
~10k queries × 200k refs	2B ops	✓ seconds	Baseline feasible on laptop
100k queries × 200k refs	20B ops	✓ ~minute	Still tractable if batched
1M queries × 200k refs	200B ops	✗	Need Cython/Numba or ANN
Whole-proteome × whole-proteome	4T ops	✗	ANN index required

Decision rule: benchmark before reaching for deps. A Numba-jitted inner loop typically buys 5–20× over NumPy with zero new packaging cost (Numba is already transitively available in most topiary-adjacent envs). ANN (nmslib / FAISS) is a significant dep and only justified if the benchmark shows we're in the 100M+ query regime — which no topiary consumer is today.

Suggest PR 1 ships pure NumPy + a benchmark script. PR 2 adds Numba if the benchmark shows we need it. ANN stays in the "file an issue if you hit it" bucket.

Reference-proteome versioning — lift the mhcflurry pattern

`CachedPredictor` already solved the "predictions from different model weights must not mix silently" problem with a composite-version string that's auto-composed from the local install. Same pattern applies here:

```python
from topiary.self_nearest import (
SelfNearestComputer, self_reference_composite_version,
)

ref = SelfNearestComputer.from_downloaded()

predictor_version = f"uniprot-{release}+cta-list-{pirlygenes.version}"

```

The composite encodes:

• UniProt release (e.g. `2024_05`) — whatever was fetched
• CTA exclusion list identity — pirlygenes version, or a hash of a user-supplied list
• Peptide-length set that was pre-computed (subset of `{8,9,10,11}`)

Stamped onto every `self_nearest_*` row's provenance (new column: `self_nearest_reference_version`). Two runs agreeing on this string produce interchangeable `self_nearest_peptide` values; downstream code can then use the equality as a reproducibility check.

This mirrors `mhcflurry_composite_version()` exactly — user never enumerates the bundle manually; `SelfNearestComputer.from_downloaded()` discovers the active reference and composes the string automatically.

Shape of PR 1 (minimal viable)

• `topiary/self_nearest.py`: `SelfNearestComputer(reference_by_length, metric="hamming")` class + pure-NumPy nearest-neighbor search.
• `topiary.downloads` (new): small fetcher for UniProt human reference + CTA-excluded variants. Analogous to mhcflurry's `mhcflurry-downloads fetch`.
• `TopiaryPredictor(compute_self_nearest=True, self_reference=…)` integration: after prediction, populate the `self_nearest_*` reserved namespace.
• MHC-binding pairing via `_predict_raw_peptides` — same primitive as `predict_wt=True` (TopiaryPredictor(predict_wt=True) — first-class WT predictions pass #123) and `CachedPredictor` fallback.

Benchmark output in the PR description so the Numba decision is data-driven.

Non-goals for PR 1 (separate follow-ups):

• BLOSUM62 metric (add as opt-in `metric="blosum62"`).
• Healthy-tissue expression (`self_nearest_tissues`) — needs GTEx / HPA integration.
• Non-human proteomes.

[iskandr]

Expanding scope materially after review. The single "nearest by edit distance" scalar doesn't answer the cross-reactivity question users actually have; we need to add indels to the candidate space, expose the full candidate set, and stamp it with MHC-binding properties.

Candidate generation with 1aa indels

For a query peptide of length L, candidates are drawn from reference peptides at lengths {L−1, L, L+1}. Per query, enumerate perturbations:

• Substitutions: L × 19 variants
• Deletions: L × 1 variants (length L−1)
• Insertions: (L+1) × 20 variants (length L+1)

For L=9: ≤ 380 strings per query. Check each against a hash set of reference k-mers partitioned by length. The algorithm flips from "all-pairs distance scan" to "enumerate query perturbations, intersect with reference k-mer sets" — O(380) lookups per query instead of O(M) distance computations. Much cheaper, and indels come for free because they're just perturbations at L±1.

Three "nearest" scalars (not one)

The candidate set typically has dozens of hits. Each answers a different question. Three reserved scalar columns, each with the full `self_nearest_*` sub-column set (`_gene`, `_transcript_id`, `_edit_distance`, `_edit_type ∈ {substitution, insertion, deletion}`, `_value` / `_score` / `_percentile_rank` at the query's allele):

Column	Definition	Question answered
`self_nearest_peptide`	Lowest edit distance; tie → strongest binder	"What's the most lexically similar self peptide?"
`self_nearest_by_binding`	Most similar MHC-binding profile	"Which self peptide could fool a T cell?"
`self_strongest_nearby`	Strongest binder within candidate set	"What's the cross-reactivity risk peptide?"

Full candidate set as a structured column

`self_nearest_candidates` — list-of-dicts column, one row per query peptide, each entry carrying {peptide, source_gene, source_transcript_id, edit_distance, edit_type, value, score, percentile_rank}. The three scalars above are aggregations over this set. Users drill into the full set with `Column("self_nearest_candidates").map(...)` in the DSL, or explode to long form for downstream analysis.

Provenance reverse-index

Pre-compute `kmer → [(gene, transcript_id, offset)]` at reference-load time. Peptides can appear in multiple proteins (paralogs / repeats); store all sources with a priority rule for picking a canonical one, plus a `self_nearest_is_ambiguous` flag when multiple sources exist.

MHC-binding on the candidate set is where the cost lives

Candidate generation is µs/query. The dominant cost is predicting MHC binding for every candidate at the query's allele. Typical shape: Q queries × ~40 candidates × ~3 alleles ≈ 120·Q predictions. This is where `CachedPredictor` compounds — pre-compute self-proteome predictions once, reuse forever. A bundled artifact makes sense here: `uniprot-human-nonCTA_netmhcpan41.parquet` with every (peptide, allele) at common HLA supertypes pre-scored, fetched via a helper analogous to mhcflurry-downloads.

Reserved-namespace extensions

`docs/fragments.md` currently reserves:

• `self_nearest_peptide`, `_peptide_length`, `_edit_distance`, `_gene`, `_gene_id`, `_transcript_id`, `_tissues`, `_value`, `_score`, `_percentile_rank`

Add:

• `self_nearest_edit_type` — "substitution" / "insertion" / "deletion"
• `self_nearest_is_ambiguous` — bool (peptide appears in >1 source protein)
• `self_nearest_by_binding_*` — full parallel set for the binding-similarity axis
• `self_strongest_nearby_*` — full parallel set for the strongest-binder axis
• `self_nearest_candidates` — structured column (list of dicts)
• `self_nearest_reference_version` — composite-version string for reproducibility

Phased implementation

1. PR A — substitutions only, three scalar nearest columns, `self_nearest_candidates` structured column. No bundled artifact yet; compute MHC on the fly with optional `CachedPredictor` backing. Core architecture.
2. PR B — add 1aa insertion/deletion to the perturbation enumerator. Same public surface, wider candidate space.
3. PR C — ship the bundled self-proteome × common-HLA prediction artifact + fetcher. Speed win, no new semantics.

This is substantially more ambitious than the original "find the closest string" framing. It's also the column set that actually answers the cross-reactivity questions users asking for "nearest self" are implicitly asking.

[iskandr]

Reference-proteome sources and species support

Third expansion of scope — the reference corpus design needs to be pluggable, multi-source by default, and species-aware.

Default pipeline (human)

• Proteome source: pyensembl (already a topiary dep) — `EnsemblRelease(release=N, species="human")`. Per-release reproducibility for free via Ensembl's versioning.
• CTA filter: `pirlygenes` (already a dep) — curated cancer-testis gene list, filters testis-restricted + germline-restricted genes out of "self."
• Tissue-expression filter (opt-in, separate axis from CTA): reuse topiary's existing tissue-exclusion plumbing. Filter out peptides from genes expressed in user-specified tissues so "what the patient's heart makes" doesn't count as "self" for a heart-vaccine context.

Alternate / supplemental CTA source

`pirl-unc/tsarina` — curated shared-antigen data (different definition, partial overlap with pirlygenes). Opt-in, not a hard dep. Users who want the tsarina definition pass `cta_source="tsarina"` and install it separately.

Configurability

All three filters (proteome source, CTA set, tissue exclusion) need override hooks. Rough shape:

```python
SelfNearestComputer.from_ensembl(
species="human",
release=93,
cta_source="pirlygenes", # "pirlygenes" | "tsarina" | set[gene_id] | callable
exclude_tissues=["heart_muscle", "lung"], # list | None | callable
)

Or, user has their own proteome (non-Ensembl):

SelfNearestComputer.from_fasta(
"my_reference.fa",
cta_source=my_gene_set, # user-supplied
exclude_tissues=None,
)
```

Every hook accepts either a named built-in (`"pirlygenes"`, `"tsarina"`) or a caller-supplied value (gene-id set, callable, tissue list). No magic — the default pipeline is just "pirlygenes + human Ensembl" passed through the same hook surface users can override.

Multi-species (mouse, dog, any pyensembl-supported species)

• Proteome: pyensembl handles every species it supports — nothing topiary-specific. `species="mouse"` / `"dog"` / `"rat"` all work.
• CTA list: pirlygenes and tsarina are human-focused. For non-human species, users supply their own CTA-equivalent gene set, or set `cta_source=None` (take the whole proteome as self). The design can't fabricate non-human CTA data; it just has to not block users who don't have it.
• Tissue data: same story — human-specific defaults, overridable.

Composite reference version

Extends the `mhcflurry_composite_version()` pattern to encode every dimension that could change predictions:

```
ensembl-{species}-{release}+cta-{cta_source_id}+tissues-{sha}
```

Examples:

• `ensembl-human-93+cta-pirlygenes-3.12.1+tissues-default`
• `ensembl-human-93+cta-tsarina-0.5.0+tissues-none`
• `ensembl-mouse-102+cta-none+tissues-none` (non-human, no CTA filter)

Stamped onto every `self_nearest_*` row's provenance (new column: `self_nearest_reference_version`). Two runs with matching strings produce interchangeable `self_nearest_peptide` values.

Doc updates that follow

• `docs/cached.md` won't need changes (this is its own feature).
• `docs/self_nearest.md` (new): covers all three defaults, the override hooks, and the composite-version pattern. Mirrors `docs/fragments.md` in shape.
• `docs/fragments.md` gets the extended reserved namespace (`self_nearest_by_binding_`, `self_strongest_nearby_`, `self_nearest_candidates`, `self_nearest_reference_version`, `self_nearest_edit_type`, `self_nearest_is_ambiguous`).
• `docs/api.md` gets the `SelfNearestComputer` table.
• `README.md` gets a top-level "Cross-reactivity analysis" section sibling to "Cached predictions."

Bundled-artifact strategy (PR C)

The pre-computed proteome × HLA prediction artifact (deferred in the earlier plan) now needs per-species / per-release variants:

• `ensembl-human-93-cta_pirlygenes_3.12.1_netmhcpan41.parquet`
• `ensembl-mouse-102_netmhciipan43.parquet` (no CTA filter applicable)
• ...

Fetcher helper (analogous to `mhcflurry-downloads`):

```
topiary-self-nearest-downloads fetch --species=human --release=93 --predictor=netmhcpan41
```

Ship human + common-HLA supertypes first. Non-human / non-common-HLA: compute on demand with fallback to `CachedPredictor` so the first run is slow but subsequent runs are cached.

Summary of scope growth on this issue

Original: "nearest peptide in a curated non-CTA proteome."
Current: "cross-reactivity analysis column set with per-species / per-source / per-filter reproducibility, three orthogonal 'nearest' axes (sequence, binding-similar, strongest-binder), full candidate set as a structured column, and pluggable proteome/CTA/tissue sources defaulting to pirlygenes + Ensembl but extensible to tsarina, user-supplied gene sets, or other species."

This is the issue that replaces "is my peptide near any self string" with "is my peptide a realistic cross-reactivity risk, and what's the provenance." Much more useful. Still phased A → B → C on the implementation side; the scope expansion is about what ships, not when.

[iskandr]

Indexing strategies — reducing lookup count

The naive "enumerate 380 perturbations per query, hit a hash set" approach is correct but not the only option. Three concrete candidates worth benchmarking:

1. Pigeonhole seed-and-extend (classical)

For edit distance ≤ k between length-L strings, at least one of (k+1) disjoint seeds of length ⌊L/(k+1)⌋ must match exactly. For k=1, L=9: 2 disjoint seeds of length ~4.

• Pre-compute: `seed_table[(position, seed_string)] → [reference_peptide_ids]`. Size O(N · 2).
• Per-query lookup: 2 exact-match set lookups; union the candidate lists; verify each candidate's full edit distance. Query cost drops from 380 hash lookups to 2 + len(candidate_list), and candidate lists for real reference 4-mer frequencies are small.
• Indels: 1-position shifts the query relative to the reference, so length-L±1 references need seed lookups at shifted positions. 2 seed positions × 3 length-target offsets = 6 lookups per query. Still much better than 380.

2. SIMD-vectorized naive scan

Modern CPUs crush pairwise Hamming. Encode reference as `(M, L) int8`; per query, a single SIMD-backed `(q != refs).sum(axis=1)` is microseconds. Zero indexing infrastructure, no edge cases for indels (separate passes per length). Tradeoff: O(Q · M) linear scan, but the constant factor is tiny.

For our expected scale (~200k human reference 9-mers × ~10k mutant queries = 2B ops), this is probably already fast enough on a laptop.

3. Radix trie of the reference

Walk the query down a compressed trie; fork on mismatches, cap at k forks. Nice cache behavior, no hash overhead. Overkill at our size unless measurement shows otherwise.

Decision

Benchmark both pigeonhole and SIMD first. The winner depends on scale:

Scale	Winner
Small ref (~200k) + moderate queries (~10k)	SIMD-vectorized naive (simpler, seconds to run)
Large ref (>1M) or huge queries (>1M)	Pigeonhole seed-and-extend (O(2) lookups per query beats O(M) scan)

Implementation plan: ship whatever the benchmark says is simpler for the default use case (human Ensembl ~200k refs, ~10k mutant queries). Keep the other as a documented fallback in the class — `SelfNearestComputer(algorithm="naive" | "seed_and_extend")`. Users at a different scale regime can switch without rewriting their code.

Filed as a benchmark prerequisite before PR A lands.

[iskandr]

Benchmark spec — construction + lookup, both algorithms

Decision between SIMD-naive and pigeonhole seed-and-extend should fall out of measurements, not intuition. Four distinct axes to report, not one.

Harness

• `time.perf_counter()` for wall-clock (repeat N times, report median + min/max).
• `tracemalloc` + `psutil.Process().memory_info().rss` for memory.
• `tracemalloc.get_traced_memory()` peak for "working memory during" vs baseline RSS for "steady memory after."
• Commit the benchmark script as `tests/bench_self_nearest.py`; not part of the pytest run, driven by a separate CLI entry point. Future changes can re-run it.

Construction phase

Measure for each algorithm at each reference size:

Metric	What
Wall-clock	Seconds to build the index from a list of reference peptides.
Peak RSS during construction	How much memory the process touches while building.
Steady RSS of the built index	How much the object retains in steady state (post-GC).

Vary:

• Reference size (peptides per length): 10k, 100k, 1M, 10M.
• Peptide length: 8, 9, 10, 11 separately, plus a combined-all-lengths run.
• Algorithm: `naive` (int8 array), `seed_and_extend` (pigeonhole dict).

Output shape: 4 (size) × 5 (length) × 2 (algorithm) = 40 rows per metric.

Lookup phase

Per algorithm, per reference size, per batch size:

Metric	What
Single-query latency (µs)	One query at a time, median over 10k trials.
Batch throughput (queries/sec)	N-query batch to completion, warm cache.
Working memory during batch	RSS delta from baseline during a batch.
Candidate-list size (seed_and_extend)	Mean / p95 candidates verified per query (sanity check on pigeonhole).

Vary:

• Batch size: 1, 100, 1k, 10k, 100k queries.
• Reference size: 10k, 100k, 1M, 10M (same as construction).
• Algorithm.

Output: 5 (batch) × 4 (ref) × 2 (algo) = 40 rows per metric.

Sanity-check axes

Cheap to add; catches bugs early:

• Accuracy: both algorithms must return identical candidate sets for every query at every scale. Cross-check.
• Indel correctness: separate trial batches for substitution-only, insertion-only, deletion-only queries. Seed-and-extend's length-shifted seed logic is subtle; verify against the SIMD ground truth.

Deciding the default

After the runs, look at:

1. Crossover point in lookup throughput: below which ref size does SIMD beat seed-and-extend? Above?
2. Construction amortization: for a one-shot (run predictor once, discard cache) workload, does seed-and-extend's higher construction cost matter, or does it amortize across the lookup batch?
3. Memory ceiling: at what reference size does SIMD's chunked diff tensor push the total working set past a reasonable laptop ceiling (~8GB)? Does seed-and-extend stay bounded?
4. Reality check on default size: human Ensembl non-CTA × 4 lengths ≈ 800k reference peptides total; typical mutant query set is 1k–10k peptides. Which algorithm wins at that specific point becomes the default; the other ships as `algorithm="…"` opt-in.

What this prevents

Picking the "clever" algorithm and paying in construction complexity or memory footprint to save lookup time we don't actually need, or picking the naive one and discovering later that it doesn't scale past the default case for anyone who brings their own proteome. Both are real failure modes worth checking against measurement.

[iskandr]

Class rename: `SelfProteome` + scope knob

Dropping `SelfNearestComputer` — weak suffix-`Computer` naming. The object is the reference corpus (plus its index); the nearest-lookup is an operation on it. Final: `SelfProteome`.

Scope toggle

Three named modes + callable escape hatch, selected at construction:

```python

Default: strip CTAs

ref = SelfProteome.from_ensembl(species="human")

All proteins, no filter

ref = SelfProteome.from_ensembl(species="human", scope="all")

Only proteins expressed in named tissues (cross-reactivity-risk focus)

ref = SelfProteome.from_ensembl(
species="human",
scope="protected_tissues",
tissues=["heart", "lung", "liver", "brain"], # None → pirlygenes default set
tissue_source="hpa", # "hpa" | "gtex"
min_expression=1.0, # TPM threshold (GTEx only)
)

Arbitrary gene-level filter

ref = SelfProteome.from_ensembl(
species="human",
scope=lambda gene: gene.gene_id not in my_exclusion_set,
)
```

Scope parameter matrix

`scope=`	Behavior	Parameters
`"all"`	Whole proteome, no filter	—
`"non_cta"` (default)	Remove CTA genes	`cta_source="pirlygenes"` / `"tsarina"` / set / callable
`"protected_tissues"`	Keep only genes expressed in named tissues	`tissues=[…]`, `tissue_source="hpa"` / `"gtex"`, `min_expression=`
callable	Arbitrary `gene → bool`	—

Tissue-source selection

• `"hpa"` (Human Protein Atlas) — qualitative annotations. Binary inclusion, no threshold needed. pirlygenes exposes curated sets.
• `"gtex"` — quantitative TPM. `min_expression=` threshold required; ignored for HPA.

Both pulled through pirlygenes; user chooses per-construction.

Non-human species

• `"all"` — works for any pyensembl-supported species.
• `"non_cta"` and `"protected_tissues"` — need species-appropriate reference data. pirlygenes and tsarina are human-only today; non-human users must supply their own filter or fall back to `"all"`.
• When species + scope + data are incompatible, raise with a clear message rather than silently returning unfiltered.

Composite version string

```
ensembl-{species}-{release}+scope-{scope}+{scope_specific_params}
```

Examples:

• `ensembl-human-93+scope-non_cta+cta-pirlygenes-3.12.1`
• `ensembl-human-93+scope-protected_tissues+tissues-vital_organs+source-hpa`
• `ensembl-human-93+scope-protected_tissues+tissues-sha256:abc123+source-gtex+min_tpm-1.0`
• `ensembl-human-93+scope-all`
• `ensembl-mouse-102+scope-all`

Stamped on every row as `self_nearest_reference_version`. Two runs produce interchangeable `self_nearest_peptide` values iff the composite strings match.

When `tissues` is a user-supplied list (not a named pirlygenes set), the composite uses `tissues-sha256:{hash}` so custom tissue sets stay reproducible without ballooning the string.

Non-Ensembl sources

`SelfProteome.from_fasta(path)` for users with their own proteome. Same scope knob, but without pyensembl-backed gene/tissue metadata the only practically useful scopes are `"all"` and a callable. Document this explicitly so users don't expect `scope="non_cta"` to magically work on arbitrary FASTA.

Default story

• Human user: `SelfProteome.from_ensembl(species="human", release=N)` picks up `scope="non_cta"` defaulting to pirlygenes. Zero configuration needed.
• Non-human user: must be explicit — `scope="all"` or supply their own filter. Clear error otherwise.
• Power user: callable for `scope=` + arbitrary data for `cta_source=` / `tissues=`.

This keeps the zero-config path working for the 90% case (human cross-reactivity analysis against vital organs) while leaving every axis overridable.

[iskandr]

Vaxrank perspective

Not vaxrank-blocking today — vaxrank has a simpler equivalent already (ReferenceProteome.contains(peptide) → bool occurs_in_reference column on EpitopePrediction). But this proposal is a principled upgrade to that, and I'd adopt it if the ergonomics are right.

The vaxrank/topiary semantic gap

Vaxrank today uses a user-supplied pyensembl Genome to build the reference proteome (vaxrank/reference_proteome.py). This is deliberate — it matches the annotation set the variant calls were made against, and the user is responsible for picking a sensible reference. Vaxrank's current check is a boolean "is this exact peptide anywhere in the reference proteome" — not a nearest-neighbor distance, not a CTA-aware subset.

Your proposal is semantically different:

• A bundled non-CTA healthy-tissue proteome (not user-chosen)
• Nearest-neighbor with edit distance + paired MHC binding (not a simple membership test)

Both have value; they answer different questions.

API shape preferences for vaxrank adoption

1. Don't force the bundled proteome. Expose the same compute_self_nearest=True, self_reference=<proteome> path with self_reference being any callable or indexed object implementing a find_nearest(peptide, allele) -> dict protocol. Vaxrank would plug its pyensembl-backed ReferenceProteome into that slot and avoid pulling in a second multi-megabyte blob during a pipeline run.
2. Downgrade gracefully for callers who only want membership. self_nearest_edit_distance == 0 is a superset of vaxrank's current occurs_in_reference check. If the matching path short-circuits on exact match (as k-mer index naturally does), vaxrank can drop ReferenceProteome.contains() for self_nearest_edit_distance == 0 at zero cost.
3. CTA carve-out as a filter, not a hard exclude. The bundled artifact is a ProteinFragment collection (or similar); the CTA list is a separate artifact the user can AND/OR with their own set. Vaxrank users frequently want testis expression to count as "self" for certain indications — burying the CTA exclusion inside the loader makes it invisible.
4. self_nearest_value / _score / _percentile_rank populated via _predict_raw_peptides — same primitive as TopiaryPredictor(predict_wt=True) — first-class WT predictions pass #123's predict_wt. Keep the two features architecturally parallel so a caller enabling both doesn't see surprising interactions (e.g. wt_* computed from length-matched baseline, self_nearest_* computed from whole-proteome search, coexisting in the same row).

Indel / frameshift story

The "why this instead of indel wt_peptide remapping" section in the issue body is the right call. Vaxrank today uses the overlaps_mutation=False branch to fall back to wt_peptide = peptide, wt_ic50 = ic50 for non-SNV rows — which is semantically "I don't know, pretend it's self." self_nearest_* replaces that heuristic with a real answer. That alone is the argument to prefer this path over coordinate remapping.

Vaxrank adoption priority

Behind #123 + #126. Those give large code-deletion wins; this one gives a capability upgrade. Sequence:

1. Ship TopiaryPredictor(predict_wt=True) — first-class WT predictions pass #123 → vaxrank adopts predict_wt=True, deletes WT plumbing.
2. Ship evaluate_scores(df, node) -> pd.Series — row-aligned DSL convenience #126 → vaxrank adopts evaluate_scores / evaluate_filter, deletes reindex dance.
3. Ship Computed self_nearest_* columns against a curated non-CTA proteome #124 → vaxrank optionally swaps occurs_in_reference boolean for self_nearest_edit_distance == 0, gains nearest-neighbor as a new ranking signal.

Step 3 is independent of the first two and can land any time. But it requires the user-pluggable self_reference path — if it's bundled-proteome-only, vaxrank won't adopt because the "user picks the reference" principle is load-bearing.

[iskandr]

Cross-species plan — what's universal, what's human-only, and the plug-in protocol

Pulling this out of the general design because it's the axis most likely to surprise non-human users.

What's universal (works for any pyensembl-supported species)

• Proteome retrieval: `pyensembl.EnsemblRelease(species=…, release=…)` returns protein sequences for any species pyensembl knows. Mouse, dog, rat, zebrafish, arabidopsis — all work.
• Peptide indexing + nearest-neighbor search: the SIMD-naive (and eventual seed-and-extend) algorithms don't care about species at all. Bytes are bytes.
• `scope="all"`: no filtering, universal.
• `scope=callable`: user-supplied gene → bool filter. Works for any species because the user owns the logic.

What's human-only today

• `cta_source="pirlygenes"` — pirlygenes curates human CTA gene lists.
• `cta_source="tsarina"` — tsarina is also human-focused.
• `tissue_source="hpa"` — Human Protein Atlas, human-only by definition.
• `tissue_source="gtex"` — GTEx, human-only.
• Default protected-tissue presets for `scope="protected_tissues"` — curated via pirlygenes, human-tissue-specific.

None of these have canonical cross-species equivalents in a form topiary can bundle today. Mouse has Tabula Muris, ENCODE mouse tissue atlas, MGI; dog has scattered tissue RNA-seq; these aren't yet in pirlygenes or topiary.

Plug-in protocol for non-human users

Every human-only string can be replaced with a user-supplied equivalent:

```python

Mouse, user-supplied CTA list as a gene-ID set

ref = SelfProteome.from_ensembl(
species="mouse",
release=102,
scope="non_cta",
cta_source={"ENSMUSG...", "ENSMUSG..."}, # user's own set
)

Mouse, user-supplied tissue data via callable

def mouse_tissue_lookup(gene):
"""Return list of tissues gene is expressed in (e.g., from Tabula Muris)."""
return ...

ref = SelfProteome.from_ensembl(
species="mouse",
release=102,
scope="protected_tissues",
tissues=["heart", "lung"],
tissue_source=mouse_tissue_lookup,
)

Arbitrary species, fully custom filter

ref = SelfProteome.from_ensembl(
species="canis_familiaris",
release=104,
scope=lambda gene: my_judgment(gene),
)
```

Protocol:

• `cta_source=` accepts `"pirlygenes"` (human default), `"tsarina"` (human opt-in), a `Set[gene_id]`, or a `Callable[[Gene], bool]`.
• `tissue_source=` accepts `"hpa"` (human default for HPA mode), `"gtex"` (human opt-in), or a `Callable[[Gene], List[str]]` returning the tissues where the gene is expressed.
• `scope=callable` bypasses the named-mode layer entirely.

Error story when species + scope combination has no data

Fail loudly at construction — silent unfiltered results are dangerous for immunogenicity analysis:

```

This raises:

SelfProteome.from_ensembl(species="mouse", scope="non_cta")

ValueError: scope="non_cta" defaults to cta_source="pirlygenes",

which is human-only. For species="mouse", either pass

cta_source= explicitly, or use scope="all".

SelfProteome.from_ensembl(species="mouse", scope="protected_tissues")

ValueError: scope="protected_tissues" defaults to tissue_source="hpa",

which is human-only. For species="mouse", either pass

tissue_source= explicitly, or use scope="all".

```

What the class needs to know per species

For each species we formally support with defaults:

Data	Human (shipped)	Mouse (future)	Other (user-supplied)
Proteome	pyensembl	pyensembl	pyensembl
CTA list	pirlygenes default	TBD (Tabula Muris?)	callable / set
Tissue expression	HPA / GTEx	TBD	callable
Protected-tissue preset	pirlygenes curated	TBD	caller's `tissues=` list

Implementation-wise: a `_SPECIES_DEFAULTS` registry in `self_proteome.py` maps species → {cta_source_default, tissue_source_default, protected_tissues_preset}. Human is populated; other species raise with a clear message when asked to default. Community can PR additional species entries without touching the core algorithm.

Composite version string handles species cleanly

Already in the design:

```
ensembl-{species}-{release}+scope-{scope}+{scope_params}
```

Species is a first-class axis, so mouse runs with different tissue data never collide with human runs in reproducibility checks. User-supplied data produces an sha256-hashed identity component so custom CTA lists or tissue callables stay stable across reruns.

What this means for PR A scope

PR A ships the architecture:

• `SelfProteome` class with `from_fasta` and `from_ensembl` (species is just a pyensembl kwarg)
• `scope="all"` works universally, zero config
• `scope="non_cta"` works for human via pirlygenes; raises clearly for other species
• `scope=callable` works for anything
• `_SPECIES_DEFAULTS` registry with only human populated

`scope="protected_tissues"` and custom `tissue_source=` callables land in PR B along with indels. PR B is where the tissue-data layer gets filled out.

Non-human defaults (mouse CTA list, etc.) are explicitly out of scope for the topiary repo — they belong in pirlygenes or similar species-specific packages. Topiary provides the hooks; bundled data is someone else's problem.

[iskandr]

Binding-awareness and multi-allele — concrete design

Clarifying what part A ships vs. what binding-aware nearest actually needs to do in part B.

Part A (current PR #135) is allele-independent sequence-nearest

```python
df = predictor.predict_from_variants(...) # one row per (peptide, allele)

self_nearest_peptide populated via .merge(on='peptide')

→ identical across every allele of a given peptide.

```

This is correct for a sequence-only baseline signal ("is my mutant near anything in self at all?") but it is not a cross-reactivity risk signal — a self peptide that's close lexically but doesn't bind the patient's MHC isn't actually presentable and can't drive cross-reactivity.

Part B introduces binding-aware columns, one value per (peptide, allele)

Four axes in the reserved namespace, each exposed with its full
`_peptide` / `_peptide_length` / `_edit_distance` / `_gene_id` /
`_transcript_id` / `_value` / `_score` / `_percentile_rank`
sub-columns:

Axis	Scope	Semantic
`self_nearest_*` (shipped in part A)	`on="peptide"` — allele-independent	Sequence-closest reference peptide. Fast baseline.
`self_nearest_by_binding_*`	`on=["peptide", "allele"]`	Within the edit-distance-≤1 set, the candidate whose binding profile is most similar to the query's at this allele. "Which self peptide could fool a T cell trained on the mutant?"
`self_nearest_strength_gated_*`	`on=["peptide", "allele"]`	Within the edit-distance-≤1 set, filter to candidates that bind at least as strongly as the query at this allele, then pick sequence-closest. None if no candidate passes the gate — which correctly means "no presentable cross-reactivity risk at this allele."
`self_strongest_nearby_*`	`on=["peptide", "allele"]`	Within the edit-distance-≤1 set, the candidate with the strongest absolute binding at this allele, regardless of sequence distance. "What's the peak cross-reactivity risk peptide?"

The strength-gated axis ships with a tunable metric

```python
predictor = TopiaryPredictor(
models=...,
alleles=...,
self_proteome=ref,
self_nearest_metric="percentile_rank", # default; lower = stronger
# OR: "affinity" (lower nM = stronger), "score" (higher = stronger)
)
```

Metric is used for `self_nearest_strength_gated_` and
`self_strongest_nearby_`. Default is `percentile_rank` because
it's the most cross-tool-comparable metric and doesn't bias across
predictors with different affinity calibrations.

Multi-allele example

Mutant `SIINFEXLG` with two candidate self peptides `SIINFEKLG` and
`TIINFEELG` (both edit distance 1):

	`self_nearest` (seq-only)	`self_strength_gated`	`self_strongest_nearby`
HLA-A*02:01, mutant rank=1.5	`SIINFEKLG` (dist=1)	`SIINFEKLG` if its A*02:01 rank ≤ 1.5, else `None`	whichever candidate has the lower A*02:01 rank
HLA-B*07:02, mutant rank=40	`SIINFEKLG` (dist=1)	probably `SIINFEKLG` (gate is loose at rank=40)	may be `TIINFEELG` if it's a better B*07:02 binder
HLA-C*07:02, mutant rank=99	`SIINFEKLG` (dist=1)	both candidates pass the 99 gate → sequence-closer wins (`SIINFEKLG`)	whichever is strongest at C*07:02

The sequence-only column is the same across all three alleles (by
design — sequence doesn't depend on allele). The binding-aware
columns diverge.

MHC predictions on the candidate set — where the cost lives

Computing any of the three binding-aware columns requires predicting
MHC binding for every (candidate, allele) pair. Typical shape:
Q queries × ~40 edit-distance-1 candidates × ~3 alleles ≈ 120·Q
predictions.

This is exactly where `CachedPredictor` compounds — pre-compute
self-proteome × common HLA once, reuse across every run (the
"bundled artifact" from PR C). Without the artifact, part B's first
run on a given proteome will be slow; subsequent runs hit the cache.

When `None` is the right answer

`self_nearest_strength_gated_*` returning `None` at an allele is a
feature, not a gap: it correctly means "no self peptide within
edit distance 1 that also presents on this MHC." That's a
low-cross-reactivity-risk signal. Don't `fillna(0)` on that column
— users will filter on `self_nearest_strength_gated_peptide.isna()`
as a "safe at this allele" predicate.

Updated PR plan

• Part A (PR Add SelfProteome — nearest-self lookup architecture (part A of #124) #135, in flight): `self_nearest_*` sequence-only,
  allele-independent. Ships the architecture.
• Part B: the three binding-aware axes above + 1aa indel
  candidates + `self_nearest_candidates` structured column + the
  `self_nearest_metric=` knob.
• Part C: bundled self-proteome × common-HLA prediction artifact
  • fetcher. Turns part B's "first run is slow" into "first run is
    fast for anyone using the defaults."

[iskandr]

Design revision — binding + distance + position layering

Revising the earlier sketch after working through the "what does nearest actually mean immunologically" angle. Three substantive changes:

1. Rename `self_nearest_strength_gated_` → `self_mimic_`

The cross-reactivity-risk signal is the peptide most likely to be mistaken for the mutant at this MHC — molecular mimicry, standard immunology terminology. `self_mimic_*` captures it directly and is 3× shorter.

2. Drop `self_nearest_by_binding_*`

It was "binding-similar to the query" — fuzzy axis (what is "similar binding"? within X nM? rank delta < Y?). Cut.

3. Revised four-axis design

Column family	Merge key	Distance metric	Binding gate	Semantic
`self_nearest_*` (part A, shipped)	`on="peptide"`	Hamming	None	Fast allele-independent baseline
`self_mimic_*` (part B)	`on=["peptide", "allele"]`	TCR-weighted BLOSUM62	binding ≥ query binding at this allele	Primary cross-reactivity signal
`self_strongest_nearby_*` (part B)	`on=["peptide", "allele"]`	n/a (selected on binding)	edit-distance ≤ 1	Peak risk — strongest-binding near-neighbor
`self_nearest_candidates` (part B)	`on="peptide"` (list col)	both Hamming + BLOSUM stored	None	Full set for downstream analysis

Multiple candidates within edit distance 1 is the common case

For 9-mers in human non-CTA Ensembl, single-substitution neighbors are common (10–50 is typical for conserved regions). Reference-order tie-breaking is essentially random noise. Need a real metric.

BLOSUM62-weighted distance as the default tie-breaker

```
distance(query, ref) = Σ_i (BLOSUM62_max − BLOSUM62[query[i], ref[i]])
```

• Identical residues → 0 distance per position.
• Conservative subs (I↔L, V↔I, D↔E, K↔R) → small positive (~1–3).
• Non-conservative (W→A, P→G) → large positive (~6–10).

Vectorizable: precompute a 21×21 BLOSUM62 int8 lookup, then scores = table[q_ints[:, None, :], ref_ints[None, :, :]].sum(axis=2). Same SIMD shape as Hamming, ~2× constant factor.

Hamming stays available as an opt-in `metric="hamming"` because it's faster and sometimes what the caller wants ("count positions that differ at all").

Same substitution at different positions is not equivalent

Position matters because the immunological role of each position in a 9-mer varies:

Pos	Role	TCR-facing?
P1	Secondary anchor / solvent	mostly anchor
P2	Primary anchor	no
P3	variable	partial
P4–P7	TCR contact	yes
P8	variable	partial
P9	Primary anchor (C-term)	no

Class II is looser (9-mer core inside a longer peptide, broader TCR contact).

TCR-facing vs. anchor — which direction is "more similar"?

Inverts the naive intuition:

• Substitution at a TCR-facing position: TCR sees a different side chain → distinguishes → LOW cross-reactivity risk → candidate is less similar immunologically.
• Substitution at an anchor position: TCR sees the same face; if candidate still presents, TCR can't distinguish → HIGH cross-reactivity risk → candidate is more similar immunologically.

So the position weights for the distance metric should weight TCR-facing positions heavily and anchor positions lightly. This is opposite of "weight anchors because they're anchors" — anchors matter for presentation (the binding gate), not for recognition (the distance metric).

Two separate questions, two separate axes:

1. Presentation filter (binding gate): does the candidate present at this allele? Checked via MHC binding score.
2. Recognition distance (TCR-aware): given it presents, does the TCR see it as the same? TCR-face-weighted BLOSUM62.

Position weights

Expose `position_weights=` at `SelfProteome` construction — a list of length L defaulting to uniform weights. Sensible alternative profile for MHC-I 9-mers (near-zero at anchors, full at TCR contacts):

```
[0.5, 0.0, 0.5, 1.0, 1.0, 1.0, 1.0, 0.5, 0.0]
```

Defaulting to uniform because position weighting is allele-dependent and contentious — better to ship a sensible default than a clever-wrong one. Advanced users supply their own profile.

MHC binding: per-allele scope, percentile rank default

Allele scope

Per-(peptide, allele) only. A self peptide presented on allele A can only drive T-cell cross-reactivity against a mutant presented on the same allele A — TCRs are allele-specific. Don't cross alleles; it's immunologically nonsense.

```python

For each (peptide, allele) row in the prediction output:

candidates = edit_distance_≤1_neighbors(peptide)

for c in candidates:

binding_at_this_allele = predictor(c, allele)

gated = [c for c in candidates if passes_gate(c, allele)]

self_mimic[peptide, allele] = closest_by_tcr_distance(gated)

```

Binding metric

Default: `percentile_rank`. Cross-tool comparable; threshold semantics well-known (rank ≤ 2 = weak binder, rank ≤ 0.5 = strong binder); robust to affinity-calibration drift across predictor versions.

Exposed options via `binding_metric=`:

Metric	Semantic	Caveat
`"percentile_rank"` (default)	lower stronger	cross-tool comparable
`"affinity"` (raw nM)	lower stronger	predictor-calibration-specific
`"score"` (normalized)	higher stronger	within-tool comparable
`"presentation_rank"` (mhcflurry)	lower stronger	mhcflurry-specific

Binding gate shape

```python
TopiaryPredictor(
...,
self_proteome=ref,
self_mimic_binding_metric="percentile_rank",
self_mimic_binding_gate="at_least_as_strong_as_query", # or float, or callable
)
```

Preset modes:

• `"at_least_as_strong_as_query"` — candidate's binding ≥ query's at this allele. `None` when no candidate passes; semantically correct "no presentable cross-reactivity risk at this allele" signal.
• `"strong_binder"` — absolute threshold (default rank ≤ 2.0).
• float — custom threshold at the chosen metric.
• callable — arbitrary (candidate_prediction, query_prediction) -> bool predicate.

Returning None is a feature

`self_mimic_peptide is None` at an (peptide, allele) row correctly means "no self peptide within edit distance 1 that also presents on this MHC." Don't `fillna` — users filter on `self_mimic_peptide.isna()` as a "safe at this allele" predicate. Same rule applies to `self_strongest_nearby_*` when no edit-distance-≤1 candidate exists.

Computational implications summary

• BLOSUM distance: same SIMD shape as Hamming, ~2× constant factor. Still within budget.
• Position weights: one broadcast multiply, O(L) overhead. Negligible.
• Binding gate: needs `predict(candidate, allele)` for every candidate × allele. This is the expensive step. The bundled self-proteome × common-HLA prediction artifact (part C) turns this into a cached lookup. Without the bundle, users pay a one-time predictor cost per unique candidate; `CachedPredictor` amortizes it across runs.

Part B scope (revised)

1. 1aa insertion / deletion candidates (edit-distance-1 perturbation enumeration).
2. BLOSUM62 distance metric + `metric="hamming"` opt-in.
3. `position_weights=` kwarg (uniform default).
4. `self_mimic_*` columns — TCR-weighted BLOSUM62 + per-allele binding gate.
5. `self_strongest_nearby_*` columns — edit-distance-1 set + max binding at allele.
6. `self_nearest_candidates` structured column — full set with both distance metrics and binding scores.
7. `binding_metric=` / `self_mimic_binding_gate=` construction knobs.

Part C remains bundled-artifact + seed-and-extend indexing.