API Reference#

Index#

class seqtree.Index#

Immutable search index over a set of reference sequences. Build once, then query concurrently; reference id is the position in refs.

align(self: seqtree._core.Index, ref_id: SupportsInt | SupportsIndex, query: str, params: seqtree._core.SearchParams) → seqtree._core.Alignment#: Compute a global alignment between query and a reference, on demand.

static build(refs: collections.abc.Sequence[str], alphabet: str = 'aa') → seqtree._core.Index#: Build an index. alphabet is ‘aa’, ‘nt’, or ‘iupac’. Raises ValueError on a symbol outside the alphabet.

collisions_batch(self: seqtree._core.Index, queries: collections.abc.Sequence[str], params: seqtree._core.SearchParams, threads: SupportsInt | SupportsIndex = 0) → list[int]#: Per-query count of seqtm collisions: how often a reference was re-reached via a different edit path during branch-and-bound (0 for seqtrie / substitution-only).

static load(path: str) → seqtree._core.Index#: Load an index previously written with save(); raises on a corrupt/old file.

ref_seq(self: seqtree._core.Index, ref_id: SupportsInt | SupportsIndex) → str#: Return the reference sequence string for a reference id.

save(self: seqtree._core.Index, path: str) → None#: Serialize the index to a flat binary file for fast reload.

search(self: seqtree._core.Index, query: str, params: seqtree._core.SearchParams) → list#: Return all hits for one query within the scope/budget in params.

search_batch(self: seqtree._core.Index, queries: collections.abc.Sequence[str], params: seqtree._core.SearchParams, threads: SupportsInt | SupportsIndex = 0) → list#: Search many queries in parallel (releases the GIL). threads=0 uses all cores. Returns one hit list per query, in input order.

search_top(self: seqtree._core.Index, query: str, params: seqtree._core.SearchParams, k: SupportsInt | SupportsIndex = 1) → list#: Return up to k best (lowest-score) hits for one query.

SearchParams#

class seqtree.SearchParams#

Search scope and budget. Scope: max_subs/max_ins/max_dels (exact, seqtm) and max_total_edits. Budget: max_penalty with an optional matrix (identity/BLOSUM62/PAM250/PAM100/structural) and gap costs. engine is ‘auto’|’seqtrie’|’seqtm’, mode is ‘all’|’top’.

property engine#

property gap_extend#

property gap_open#

property matrix#

property max_dels#

property max_ins#

property max_penalty#

property max_subs#

property max_total_edits#

property mode#

property pos_matrix#

Hit#

class seqtree.Hit#

A search result. Payload-agnostic: map ref_id back to your own payload downstream. score is a non-negative penalty (0 == exact). n_subs/n_ins/n_dels are exact for the seqtm engine and 0 for seqtrie. Iterable as (ref_id, score, n_subs, n_ins, n_dels).

property n_dels#

property n_ins#

property n_subs#

property ref_id#

property score#

Alignment#

class seqtree.Alignment#

Global alignment of a query to a reference. ops has one char per column: ‘M’ match, ‘S’ substitution, ‘I’ insertion, ‘D’ deletion.

property aligned_query#

property aligned_ref#

property ops#

property score#

Scoring#

class seqtree.SubstitutionMatrix#

Non-negative substitution penalties (penalty(a,a)==0). Build a named builtin (blosum62/pam250/pam100/structural, or unit for identity) or a custom one from a similarity grid whose row/column order matches amino_acids() (or alphabet_symbols(alphabet)).

static blosum45() → seqtree._core.SubstitutionMatrix#

static blosum62() → seqtree._core.SubstitutionMatrix#

static blosum80() → seqtree._core.SubstitutionMatrix#

static from_similarity(grid: collections.abc.Sequence[collections.abc.Sequence[SupportsInt | SupportsIndex]]) → seqtree._core.SubstitutionMatrix#: Build from a square similarity grid (higher == more similar), converted to non-negative penalties via the Gram / squared-distance transform s[a,a] + s[b,b] - 2*s[a,b] (clamped at 0). Row/column order must match the target alphabet’s symbol order (see amino_acids()).

static pam100() → seqtree._core.SubstitutionMatrix#

static pam250() → seqtree._core.SubstitutionMatrix#

penalty(self: seqtree._core.SubstitutionMatrix, a: str, b: str) → int#: Gram-distance substitution penalty between two amino acids: 0 when identical, larger when more dissimilar (s(a,a)+s(b,b)-2 s(a,b)). Characters use the amino_acids() order.

scale(self: seqtree._core.SubstitutionMatrix) → int#: Median penalty over all mismatched symbol pairs – this matrix’s natural unit. Gap costs must be on this scale: BLOSUM62 has scale() == 14, so the default gap_open of 1 makes gaps ~14x cheaper than substitutions and the aligner gaps rather than substitutes. Use gap_open = 1-2 * m.scale().

similarity(self: seqtree._core.SubstitutionMatrix, a: str, b: str) → int#: Raw log-odds similarity (signed). penalty() is the non-negative Gram transform of this; the transform is lossy, so both are kept.

size(self: seqtree._core.SubstitutionMatrix) → int#

static structural() → seqtree._core.SubstitutionMatrix#

static unit(size: SupportsInt | SupportsIndex) → seqtree._core.SubstitutionMatrix#

class seqtree.PositionalMatrix#

Per-position penalties pen(pos, a, b) over a fixed frame width. Build from a base SubstitutionMatrix and per-position integer weights: weight 0 masks the position (free, not counted as a substitution – e.g. an anchor); >1 up-weights it (e.g. a TCR hotspot). Used on the seqtm Hamming path when width == query length.

static from_tables(size: SupportsInt | SupportsIndex, width: SupportsInt | SupportsIndex, data: collections.abc.Sequence[SupportsInt | SupportsIndex], masked: collections.abc.Sequence[SupportsInt | SupportsIndex] = []) → seqtree._core.PositionalMatrix#: Full per-position PSSM. data is row-major [width][size][size]; masked is an optional length-width flag array (non-zero == free position). Use this to give different regions different matrices, e.g. a germline-flank matrix and an N-region core matrix in one frame.

static from_weights(base: seqtree._core.SubstitutionMatrix, weights: collections.abc.Sequence[SupportsInt | SupportsIndex]) → seqtree._core.PositionalMatrix#: pen[pos][a][b] = weights[pos] * base.penalty(a, b); weight 0 masks the position. len(weights) is the frame width. NOTE: penalty(a, a) == 0 for every base matrix, so a weight scales MISMATCH cost only – it is a mismatch-tolerance profile, not an information/match weighting.

masked(self: seqtree._core.PositionalMatrix, pos: SupportsInt | SupportsIndex) → bool#

penalty(self: seqtree._core.PositionalMatrix, pos: SupportsInt | SupportsIndex, a: SupportsInt | SupportsIndex, b: SupportsInt | SupportsIndex) → int#

size(self: seqtree._core.PositionalMatrix) → int#

width(self: seqtree._core.PositionalMatrix) → int#

class seqtree.KmerIndex#

Seed-and-extend k-mer index for homology. Build from per-peptide k-mer lists (anchor-masked upstream) + optional allele tags; seed_and_gather fuzzy-matches query k-mers and merges posting lists into ranked candidates entirely in C++ (GIL released).

static build(kmers_per_peptide: collections.abc.Sequence[collections.abc.Sequence[str]], alphabet: str = 'aa', allele_ids: collections.abc.Sequence[SupportsInt | SupportsIndex] = []) → seqtree._core.KmerIndex#

static load(path: str) → seqtree._core.KmerIndex#

num_kmers(self: seqtree._core.KmerIndex) → int#

num_peptides(self: seqtree._core.KmerIndex) → int#

save(self: seqtree._core.KmerIndex, path: str) → None#

seed_and_gather(self: seqtree._core.KmerIndex, query_kmers: collections.abc.Sequence[collections.abc.Sequence[str]], params: seqtree._core.SearchParams, min_shared: SupportsInt | SupportsIndex = 1, allele_filter: SupportsInt | SupportsIndex = -1, threads: SupportsInt | SupportsIndex = 0) → list#: For each query (its k-mer list) return ranked Candidates with >= min_shared shared k-mers; allele_filter >= 0 restricts to that allele tag.

class seqtree.Candidate#

A seed-and-gather hit: peptide_id, shared_kmers (distinct query k-mers that hit it), best_score. Iterable as (peptide_id, shared_kmers, best_score).

property best_score#

property peptide_id#

property shared_kmers#

Functions#

seqtree.pairwise_batch(a: collections.abc.Sequence[str], b: collections.abc.Sequence[str], params: seqtree._core.SearchParams, alphabet: str = 'aa', threads: SupportsInt | SupportsIndex = 0) → list#: Batch-vs-batch search. Indexes the larger set internally and streams the smaller; results are a-major (one hit list per a[i]) with Hit.ref_id pointing into b.

seqtree.alphabet_symbols(alphabet: str = 'aa') → str#: Symbols in code order for an alphabet; custom matrices must follow this order.

seqtree.amino_acids() → str#: The amino-acid symbol order used by the built-in matrices and custom AA matrices.

seqtree.load_control(name='human_trb_aa', size=None, cache_dir=None, alphabet='aa', seed=0)[source]#

Build (or load from cache) an Index over a background control set.

The control is the E-value null. It is filtered to productive clonotypes (sanitize()) and, when subsampled, drawn uniformly over unique clonotypes rather than taken from the abundance-sorted head of the table.

The built index is cached under cache_dir. The cache is content-addressed: its filename carries a fingerprint of everything that determines the sequences – the bundled asset’s own bytes, or the download’s source and seed – so an upgrade that changes the control simply misses the old cache instead of silently serving it. Caches superseded that way are deleted.

Safe to call from many processes at once on a cold cache; see Index.save().

Parameters:

name – control identifier (e.g. "human_trb_aa").
size – number of unique clonotypes. None uses the full bundled subset; a value larger than the bundled subset triggers a HuggingFace download.
cache_dir – where to store the serialized index (default ~/.cache/seqtree).
alphabet – sequence alphabet for the index.
seed – reservoir-sampling seed for the download path, so a given (name, size, seed) is reproducible. The bundled path takes a prefix of a pre-shuffled asset and ignores it.

Returns:

An immutable Index of unique control clonotypes.

seqtree.evalues(target, control, queries, params, threads=0, exclude_exact=False)[source]#

Compute control-calibrated E-values for each query.

Parameters:

target – Index to score hits in (e.g. VDJdb), unique clonotypes.
control – background Index (e.g. healthy-donor control).
queries – list of query strings.
params – SearchParams defining the scope/budget (the ball).
threads – worker threads for the batch searches (0 = all cores).
exclude_exact – drop distance-0 (exact / self) hits from both target and control counts. Set this when queries may themselves be members of the target or control (e.g. a VDJdb-vs-VDJdb scan) so the trivial self-match is not counted.

Returns:

One dict per query with n_target, n_control, E, p_any, p_enrichment, rule_of_three.

seqtree.threshold_for_evalue(target, control, queries, params, e_target, threads=0, exclude_exact=False)[source]#

Per-query score cutoff achieving e_target against control.

One control search at params.max_penalty supplies every cutoff, so calling this and then filtering a target search at the same ceiling costs two scans, not two per query.

Parameters:

target – Index the E-value is expressed against (only its size is used).
control – background Index.
queries – list of query strings.
params – SearchParams; max_penalty is the ceiling theta_max and must be positive.
threads – worker threads for the batch search (0 = all cores).
e_target – Desired E-value.
exclude_exact – Drop distance-0 hits, matching evalues().

Returns:

One integer cutoff per query; -1 where e_target is unreachable at this control size.

Raises:

ValueError – If params.max_penalty is not positive.

seqtree.thetas_from_scores(control_scores, n_target, m_control, e_target, theta_max, *, exclude_exact=False)[source]#

Invert E = (N/M) * n_control for the score cutoff, one cutoff per query.

Scores are integers, so the inversion is exact rather than a root-find: sort a query’s control-hit scores and the answer is the value just below the (k+1)-th smallest, where k = floor(e_target * M / N) is the largest control count the target E allows.

Parameters:

control_scores – Per query, the scores of its control hits found at theta_max. Hits above theta_max are irrelevant and may be omitted.
n_target – N, size of the target index.
m_control – M, size of the control index.
e_target – Desired E-value, e.g. 0.05.
theta_max – Score ceiling the control was searched at. Returned cutoffs never exceed it.
exclude_exact – Drop distance-0 hits, matching evalues().

Returns:

One integer cutoff per query, or -1 where no cutoff achieves e_target. That happens when e_target < 3N/M: with only M control sequences, even an empty ball certifies no better than the rule-of-three bound. Enlarge the control.

Example

>>> # E(4) = (100/100)*1 = 1.0 is allowed; E(5) = 2.0 is not.
>>> thetas_from_scores([[2, 5, 9]], n_target=100, m_control=100, e_target=1.0,
...                    theta_max=10)
[4]

Pairwise alignment (Needleman-Wunsch / Smith-Waterman)#

Needleman-Wunsch and Smith-Waterman: ordinary protein alignment, without BioPython.

Everything else in seqtree minimises a non-negative penalty – that is what a search ball and an E-value need. This module does the opposite: it maximises a raw log-odds similarity, the way BLAST and BioPython do, because that is what an ordinary pairwise alignment means and what downstream code expects to get back.

The two views live on the same SubstitutionMatrix:

mat.penalty(a, b)     >= 0, zero on the diagonal   -- search, E-values, gap blocks
mat.similarity(a, b)  signed log-odds              -- the aligners here

The penalty is the Gram transform of the similarity, pen = s(a,a) + s(b,b) - 2·s(a,b), which is lossy: it forces the diagonal to zero and destroys s(a,a). So the raw grid is kept rather than reconstructed, and a similarity score cannot be recovered from a penalty.

Conventions, all verified against BioPython (tests/python/test_pairwise.py runs 6,720 comparisons across three matrices, ten gap settings and both modes; zero disagreements):

a gap run of length L costs gap_open + (L-1)·gap_extend – gap_open is the cost of the first gap column, not a surcharge on top of it;
gap_open == gap_extend gives linear gaps. There is no separate mode for it;
mode="global" charges end gaps like any other (true Needleman-Wunsch, not semi-global);
mode="local" never lets the score fall below zero and takes the best cell anywhere (Smith-Waterman).

Gap costs are positive magnitudes and are subtracted. BLAST’s protein defaults are gap_open=11, gap_extend=1; BioPython’s PairwiseAligner("blastp") preset uses 12, 1.

Example

>>> import seqtree
>>> from seqtree.pairwise import align, score
>>> mat = seqtree.SubstitutionMatrix.blosum62()
>>> score("CASSLGQAYEQYF", "CASSPGQAYEQF", mat)          # global, BLAST defaults
45
>>> score("CASSLGQAYEQYF", "CASSPGQAYEQF", mat, gap_open=12)   # BioPython's 'blastp' preset
44
>>> aln = align("WWWAAAWWW", "KKKAAAKKK", mat, mode="local")   # Smith-Waterman
>>> aln.score, aln.aligned_query, aln.aligned_ref
(12, 'AAA', 'AAA')

seqtree.pairwise.score(query, ref, matrix, mode='global', gap_open=11, gap_extend=1, alphabet='aa')[source]#

Optimal alignment score of query against ref.

Parameters:

query (str) – First sequence.
ref (str) – Second sequence.
matrix (SubstitutionMatrix) – Scoring matrix; its similarity view is used, not its penalty.
mode (str) – "global" for Needleman-Wunsch, "local" for Smith-Waterman. "nw" and "sw" are accepted too.
gap_open (int) – Cost of the first column of a gap. Positive; it is subtracted.
gap_extend (int) – Cost of each further column. Equal to gap_open means linear gaps.
alphabet (str) – "aa", "nt" or "iupac".

Returns:

The score, signed. Higher is more similar – the opposite sense to the rest of seqtree.

Raises:

ValueError – On a negative gap cost, an unknown mode, or a symbol outside the alphabet.

Return type:

int

Example

>>> m = SubstitutionMatrix.blosum62()
>>> score("AAA", "AAA", m)
12
>>> score("AAA", "AAAAA", m)          # a length-2 gap: 11 + 1*1 = 12
0

seqtree.pairwise.align(query, ref, matrix, mode='global', gap_open=11, gap_extend=1, alphabet='aa')[source]#

As score(), but also returns the aligned strings and the edit ops.

Returns:

An Alignment. Note Alignment.score here is a similarity (signed, higher is better), whereas the same field from seqtree.Index.align() is a penalty. In local mode the aligned strings are the matched sub-sequences only.

Parameters:

query (str)
ref (str)
matrix (SubstitutionMatrix)
mode (str)
gap_open (int)
gap_extend (int)
alphabet (str)

Return type:

Alignment

Example

>>> m = SubstitutionMatrix.blosum62()
>>> a = align("CASSLGQAYEQYF", "CASSPGQAYEQF", m)
>>> a.aligned_query, a.aligned_ref
('CASSLGQAYEQYF', 'CASSPGQAYEQ-F')

seqtree.pairwise.score_matrix(queries, refs, matrix, mode='global', gap_open=11, gap_extend=1, alphabet='aa', threads=0)[source]#

Every query against every reference, in C++ with the GIL released.

Returns:

A ScoreMatrix of shape (len(queries), len(refs)) holding signed similarity scores. numpy.asarray wraps it without copying.

Parameters:

queries (Sequence[str])
refs (Sequence[str])
matrix (SubstitutionMatrix)
mode (str)
gap_open (int)
gap_extend (int)
alphabet (str)
threads (int)

Return type:

ScoreMatrix

seqtree.pairwise.dist_matrix(queries, refs, matrix, mode='global', gap_open=11, gap_extend=1, alphabet='aa', threads=0)[source]#

Alignment distances: d(a, b) = s(a,a) + s(b,b) - 2·s(a,b).

The Gram transform, applied at the sequence level to the alignment scores rather than per residue. Non-negative, symmetric, zero on the diagonal – so it is a distance, and it is what a prototype-distance embedding actually wants. This is the quantity users of BioPython hand-roll, and it is computed here without a Python loop: the self-scores are taken once per sequence, not once per pair.

Returns:

A ScoreMatrix of shape (len(queries), len(refs)).

Parameters:

queries (Sequence[str])
refs (Sequence[str])
matrix (SubstitutionMatrix)
mode (str)
gap_open (int)
gap_extend (int)
alphabet (str)
threads (int)

Return type:

ScoreMatrix

Example

>>> m = SubstitutionMatrix.blosum62()
>>> d = dist_matrix(["CASSLGQAYEQYF"], ["CASSLGQAYEQYF", "CASSPGQAYEQF"], m)
>>> d[0, 0], d[0, 1] > 0
(0, True)

Plain edit distances (Hamming / Levenshtein)#

Plain string edit distances: Hamming and Levenshtein, in C++, without a dependency.

These are the unweighted distances – unit costs, no substitution matrix, no gap model, no alphabet. That is the whole point: when all you need is “how many edits apart are these two strings”, you should not have to build an SubstitutionMatrix or reach for python-Levenshtein / rapidfuzz. seqtree still needs nothing at runtime.

hamming() – number of differing positions; defined only for equal-length sequences, and raises ValueError otherwise;
levenshtein() – the classic insertion / deletion / substitution edit distance, each edit costing 1;
hamming_matrix() / levenshtein_matrix() – every a against every b in one GIL-released, multi-threaded C++ call, returned as a zero-copy ScoreMatrix.

Comparison is case-sensitive, byte for byte – unlike the search engines, which fold case. For a weighted alignment (a substitution matrix, affine gaps, local mode), use seqtree.pairwise instead.

Example

>>> from seqtree.distance import hamming, levenshtein
>>> hamming("CASSLGQYF", "CASSPGQYF")
1
>>> levenshtein("kitten", "sitting")
3

seqtree.distance.hamming(a, b)[source]#

Number of positions at which a and b differ.

Parameters:

a (str) – First sequence.
b (str) – Second sequence, of the same length as a.

Returns:

The count of differing positions (0 when identical).

Raises:

ValueError – If a and b have different lengths – Hamming distance is undefined for unequal lengths; use levenshtein() for that.

Return type:

int

Example

>>> hamming("AAAA", "AAAA")
0
>>> hamming("AAAA", "ATAT")
2

seqtree.distance.levenshtein(a, b)[source]#

Edit distance: fewest single-character insert / delete / substitute steps from a to b.

Parameters:

a (str) – First sequence.
b (str) – Second sequence; may be any length.

Returns:

The edit distance (0 when identical, max(len(a), len(b)) at most).

Return type:

int

Example

>>> levenshtein("flaw", "lawn")
2
>>> levenshtein("CASSLGQAYEQYF", "CASSPGQAYEQF")
2

seqtree.distance.hamming_matrix(a, b, threads=0)[source]#

Hamming distance of every a against every b, in parallel C++.

Parameters:

a (Sequence[str]) – Query sequences (the rows).
b (Sequence[str]) – Reference sequences (the columns).
threads (int) – Worker threads; 0 uses all cores.

Returns:

A ScoreMatrix of shape (len(a), len(b)) of int32 distances. numpy.asarray wraps it without copying.

Raises:

ValueError – If any (a[i], b[k]) pair has mismatched lengths.

Return type:

ScoreMatrix

Example

>>> import numpy as np
>>> d = np.asarray(hamming_matrix(["AAAA", "AAAT"], ["AAAA", "TTTT"]))
>>> d.tolist()
[[0, 4], [1, 3]]

seqtree.distance.levenshtein_matrix(a, b, threads=0)[source]#

Levenshtein distance of every a against every b, in parallel C++.

Unlike hamming_matrix(), sequences may differ in length freely.

Parameters:

a (Sequence[str]) – Query sequences (the rows).
b (Sequence[str]) – Reference sequences (the columns).
threads (int) – Worker threads; 0 uses all cores.

Returns:

A ScoreMatrix of shape (len(a), len(b)) of int32 distances, zero-copy through numpy.asarray.

Return type:

ScoreMatrix

Gap-block alignment#

Single-gap-block alignment for anchored loops (CDR3 / junction).

CDR3 length variation comes from V/J trimming plus N-addition – one contiguous indel event, not scattered indels. So we restrict the alignment to exactly one gap block of length d = abs(len(q) - len(r)), placed anywhere, and pick its position by score.

Two facts make this cheap and exact:

The optimum over all block positions is a prefix/suffix sum, so gapblock_score() is O(min(m, n)) rather than the O(m*n) of a full DP.
Each contiguous-d-deletion variant of the query is one block position. Hamming- matching those variants against the ordinary trie therefore reproduces the same optimum with no new engine – see GapBlockIndex.

The score is a non-negative penalty with s(q, q) == 0, so it defines a ball in the sense of appendix/evalue.tex and flows through seqtree.evalues() unchanged.

Gap costs must be on the matrix’s scale. The Gram transform puts a typical BLOSUM62 mismatch at SubstitutionMatrix.scale() == 14; a gap_open of 1 would make gaps ~14x cheaper than substitutions and every alignment would degenerate to gaps.

A sequence score alone does not locate the gap. Measured against 58 pairs of TCR-pMHC crystal structures sharing peptide and MHC, the minimum-BLOSUM62 block position agreed with the structurally correct one 8.6% of the time – indistinguishable from picking at random (8.6%). A central prior lifts that to 25.9% and cuts loop CA-RMSD from 2.15 A to 1.62 A (oracle: 1.52 A). Pass gap_prior=central_prior(...) unless you have a better one.

A prior is also what makes a column frame possible. Pairwise-optimal gap placement is not transitive: align A to B and B to C independently and the two column assignments do not compose, so a set of unequal-length sequences has no consistent column index – and hence no profile. A rule that maps length to block position supplies one. Only a rule whose block start is constant in d is transitive, i.e. one that pins the block to a fixed frame column c; see frame_prior() and embed_in_frame(). central_prior() is not such a rule – its block start drifts with d and the correspondence it induces between two shorter members splits into two blocks.

seqtree.gapblock.gapblock_score(q, r, matrix=None, gap_open=None, gap_extend=1, gap_prior=None, _pen=None)[source]#

Optimal single-gap-block alignment score of q against r.

Returns (score, block_position). score is a non-negative penalty, zero iff the sequences are identical. block_position indexes the shorter sequence and is inclusive at both ends: 0 is a leading block, min(m, n) a trailing one.

The gap block sits in whichever sequence is shorter, so this is symmetric in q/r. Being a strict restriction of affine alignment, the score is always >= the unrestricted Gotoh optimum; on pairs differing by a pure length change the two agree exactly, and at a calibrated gap_open they agree on 90% of pairs that also carry substitutions.

Parameters:

q (str) – Query sequence.
r (str) – Reference sequence.
matrix (SubstitutionMatrix | None) – Substitution penalties. None means unit cost (1 per mismatch).
gap_open (int | None) – Cost of opening the block. Defaults to 2 * matrix.scale(), or 1 for unit cost. Must be >= 0.
gap_extend (int) – Cost of each additional gap column. Must be >= 0.
gap_prior (Callable[[int, int, int], int] | None) – GapPrior, added to each candidate position. See central_prior(), profile_prior(), frame_prior(). None disables it. It applies only when there is a block to place (d > 0); otherwise s(q, q) would be non-zero and the score would no longer define a ball.
_pen (dict[tuple[str, str], int] | None) – Precomputed penalty lookup, for hot loops.

Returns:

(score, block_position).

Raises:

ValueError – If a gap cost is negative.

Return type:

tuple[int, int]

Example

>>> m = SubstitutionMatrix.blosum62()
>>> gapblock_score("CASSLGQAYEQYF", "CASSLGQAYEQYF", m)
(0, 0)

seqtree.gapblock.score_matrix(queries, refs, matrix=None, gap_open=None, gap_extend=1, gap_prior=None, alphabet='aa', threads=0)[source]#

Gap-block penalty of every query against every reference, in C++ with the GIL released.

The exhaustive counterpart of GapBlockIndex.search(): no budget, no trie, every cell scored. This is the shape a prototype-distance embedding wants – n clonotypes against a few thousand fixed references – and it is where gapblock_score() stops being fast enough, at roughly 0.4 M pairs/s in Python against ~50 M in the kernel.

Parameters:

queries (Sequence[str]) – Query sequences (the matrix rows).
refs (Sequence[str]) – Reference sequences (the matrix columns).
matrix (SubstitutionMatrix | None) – Substitution penalties; None means unit cost.
gap_open (int | None) – Block-opening cost. Defaults to 2 * matrix.scale(). See the module note: leaving this at 1 with a real matrix makes gaps ~14x cheaper than substitutions.
gap_extend (int) – Cost per additional gap column.
gap_prior (Callable[[int, int, int], int] | None) – GapPrior, materialized once into a lookup cube and then read from C++. None lets the score alone choose the block position.
alphabet (str) – "aa", "nt", or "iupac". Symbols outside it raise.
threads (int) – Worker threads; 0 means one per core. Rows are disjoint, so this scales.

Returns:

A ScoreMatrix of shape (len(queries), len(refs)). It holds 4 * len(queries) * len(refs) bytes – 1.2 GB at 100k x 3000 – so chunk the queries if that does not fit. numpy.asarray wraps it without copying.

Raises:

ValueError – If a gap cost is negative, or the prior returns a negative value.

Return type:

ScoreMatrix

Example

>>> m = SubstitutionMatrix.blosum62()
>>> sm = score_matrix(["CASSLGQAYEQYF"], ["CASSLGQAYEQYF", "CASSLGAYEQYF"], m)
>>> sm.shape
(1, 2)
>>> sm[0, 0]
0

seqtree.gapblock.deletion_variants(q, d)[source]#

Every contiguous-d-deletion variant of q, as (block_position, variant).

Variant i is q[:i] + q[i+d:]; it is exactly the query as seen through a gap block opened at position i. d == 0 has no block, so it yields a single identity variant (not len(q) + 1 copies of it).

Example

>>> deletion_variants("CAST", 1)
[(0, 'AST'), (1, 'CST'), (2, 'CAT'), (3, 'CAS')]
>>> deletion_variants("CAST", 0)
[(0, 'CAST')]

Parameters:

q (str)
d (int)

Return type:

list[tuple[int, str]]

seqtree.gapblock.central_prior(lam)[source]#

Penalise blocks whose midpoint sits away from the centre of the longer sequence.

lam * abs(block_midpoint - m/2), where the block spans [i, i + d). lam ~ 1.5 * matrix.scale() reproduces the structurally-fitted optimum. Returns an integer so the total score stays an exact non-negative penalty.

Not a transitive frame rule: its block start (m - d) // 2 moves with d.

Parameters:: lam (int)
Return type:: Callable[[int, int, int], int]

seqtree.gapblock.profile_prior(lam, w)[source]#

Charge lam per unit of positional weight the block deletes.

lam * sum(w(j, m) for j in block). With w(j, m) the probability that position j of a length-m sequence is germline-templated, this reads as lam times the expected number of templated residues the gap had to remove – deleting conserved framework is implausible, deleting a non-templated insert is free.

Parameters:

lam (int) – Cost per unit weight. Must be >= 0.
w (Callable[[int, int], float] | Sequence[float]) – w(j, m) -> float in [0, 1], or a fixed sequence indexed by j when the frame has a fixed width. Prefer the callable: loop length varies.

Returns:

A GapPrior. Zero at d == 0 (the sum is empty) and non-negative. Unlike central_prior() it is also monotone non-decreasing in d: a longer block can only delete more weight.

Raises:

ValueError – If lam is negative.

Return type:

Callable[[int, int, int], int]

seqtree.gapblock.frame_prior(lam, c)[source]#

Pin the block to frame column c: lam * abs(i - c).

The block start does not depend on d, which makes this the only kind of rule under which embedding two sequences into a common frame reproduces their pairwise single-block alignment (see embed_in_frame()). Equivalent to left-anchoring the first c residues and right-anchoring the rest.

A large lam makes the pin hard: exactly one layout survives.

Parameters:

lam (int)
c (int)

Return type:

Callable[[int, int, int], int]

seqtree.gapblock.positions_prior(starts)[source]#

Allow the block to open only at starts; let the score choose among them.

A non-negative start counts from the sequence’s beginning, a negative one from the end of the shorter sequence – so (3, 4, -4, -3) reproduces the fixed gap set that mir.distances.aligner.JunctionAligner hardcodes for every locus. Starts outside [0, shorter] clamp into range, so at least one layout always survives.

This is the “score several candidate placements and keep the best” rule. It is weaker than it looks: measured on human TRB retrieval at a matched false-positive rate, candidates (3, 4, mid) reached precision 0.156 against 0.414 for a single hard-pinned centre. The freer the placement, the more readily an unrelated reference manufactures a low score.

Parameters:: starts (Iterable[int]) – Permitted block starts. Negative values index from the end.
Returns:: A GapPrior returning 0 at a permitted start and UNREACHABLE elsewhere.
Raises:: ValueError – If starts is empty.
Return type:: Callable[[int, int, int], int]

seqtree.gapblock.embed_in_frame(seq, width, c, gap='-')[source]#

Place seq into a width-column frame, gaps blocked at column c.

Columns [0, c) hold the sequence’s own prefix (left-anchored) and columns [c + d, width) its suffix (right-anchored), with d = width - len(seq) gap columns between them. Applying this to every member of a set yields a multiple alignment whose columns are consistent – which is what a position weight matrix needs.

Example

>>> # c = 4: the V-templated CASS stays left, the J-templated EQYF stays right.
>>> for s in ("CASSLGQGAYEQYF", "CASSLGQAYEQYF", "CASSGQAYEQYF"):
...     print(embed_in_frame(s, 14, 4))
CASSLGQGAYEQYF
CASS-LGQAYEQYF
CASS--GQAYEQYF

Parameters:

seq (str)
width (int)
c (int)
gap (str)

Return type:

str

seqtree.gapblock.gap_cost(d, gap_open, gap_extend)[source]#

Affine cost of one gap block of length d. d == 0 costs nothing.

The d == 0 guard is not cosmetic: gap_open + (d-1)*gap_extend evaluates to gap_open - gap_extend at d == 0, which is negative whenever gap_open < gap_extend. A negative score would break s >= 0, scope monotonicity and admissible trie pruning.

Parameters:

d (int)
gap_open (int)
gap_extend (int)

Return type:

int

class seqtree.gapblock.GapBlockIndex(refs, alphabet='aa', d_max=1)[source]#

Bases: object

Search a reference set under the single-gap-block model, reusing the Hamming engine.

Refs shorter than the query are reached by Hamming-matching the query’s deletion variants against the ordinary index (the Hamming path only ever terminates on refs of the query’s length, so no length partitioning is needed). Refs longer than the query are reached by pre-indexing the refs’ own deletion variants, one auxiliary index per block length.

Building costs O(d_max * total_residues) extra index entries – roughly 14x the base index for CDR3 at d_max=1. Build once, query many.

Profiled over the bundled 250k control at d_max=3: 91% of query time is the first branch (the query’s own deletion variants against the base index) and only 9% the auxiliary indices, despite those holding 9.8M entries. Netting the prior out of each variant’s budget already cuts that branch from ~15 sub-searches per query to 2.5. Deduplicating variants would touch 7-10% of them before pruning and fewer after, and bucketing the auxiliary indices by reference length saves no memory – neither is worth the code.

Parameters:

refs (Iterable[str])
alphabet (str)
d_max (int)

search(query, max_penalty, matrix=None, gap_open=None, gap_extend=1, gap_prior=None)[source]#

All refs within max_penalty under the gap-block score.

Returns (ref_id, score, block_length, block_position) per ref, best score kept, sorted by ascending score. block_position indexes the shorter of the two.

Parameters:

query (str)
max_penalty (int)
matrix (SubstitutionMatrix | None)
gap_open (int | None)
gap_extend (int)
gap_prior (Callable[[int, int, int], int] | None)

Return type:

list[tuple[int, int, int, int]]

class seqtree.gapblock.ScoreMatrix#

Bases: pybind11_object

A read-only (n_queries, n_refs) int32 penalty matrix, row-major. Exposes the buffer protocol, so numpy.asarray(sm) and memoryview(sm) both wrap it without copying. Index it with sm[i, k] or pull one row with sm.row(i).

row(self: seqtree._core.ScoreMatrix, i: SupportsInt | SupportsIndex) → list[int]#: Row i as a list of penalties, one per reference.

property shape#

class seqtree.gapblock.IslandProfile(penalties, width, c)[source]#

Bases: object

A position weight matrix over one island, scored as a non-negative penalty.

Once a set of related sequences has been embedded into a common frame (see embed_in_frame()), each column has a residue distribution and a query can be scored column by column instead of against every member. The column penalty is measured against the column’s own consensus:

pen(j, a) = round(lam * log(p_max_j / p_j(a)))

which is what keeps the score usable. A textbook PWM log-odds score is signed; this one is >= 0 and exactly 0 on the consensus sequence, so it still defines a ball – centred on the consensus rather than on any one member – and still flows through seqtree.thetas_from_scores().

The gap is a column symbol like any other, so a column never gapped in the training members charges heavily for a gap there. There is no separate affine gap term: the island’s own members say where a gap is tolerated.

When this is worth it depends on your cutoff, and there are two regimes. The E-value’s k = floor(e_target * M / N) is the number of control neighbours the cutoff may admit, so the false-positive rate is k / M and it moves with N, the size of the set you annotate.

Building islands within one epitope group puts N at the group size (median 88), so k has median 142 of M = 250,000: FPR ~ 5.7e-4.
Annotating a whole repertoire against known islands puts N at ~20,000. Then e_target = 0.05 gives k = 0, which seqtree.thetas_from_scores() reports as -1: the rule of three certifies no E below 3N/M = 0.236. At that smallest certifiable E, k = 3: FPR ~ 1.2e-5.

Recall on held-out members of 108 calibrated VDJdb islands of >= 10 (human TRB, three splits each, paired bootstrap over islands, 250,000 control negatives):

regime	FPR	min-over-members	IslandProfile	difference [95% CI]
loose reference	1%	99.5 %	99.1 %	-0.40 [-1.09, +0.14]
per-epitope	0.0568%	88.3 %	89.3 %	+0.93 [-0.80, +2.79]
repertoire	0.0012%	37.6 %	48.5 %	+10.90 [+7.69, +14.21]

So: no significant difference when you are building the islands, and a large one when you use them to annotate a repertoire. On islands of >= 50 members the repertoire-regime gap is 9.8 % against 22.6 %.

It does not generalise. Junctions specific to the same epitope that fall in a different island are recovered by neither this nor min-over-members (3.5 % vs 3.7 % at a 1 % FPR, by neither at either operating point). Distinct islands share no motif either representation finds.

Nor is it a compression: 14 columns x 21 symbols x 4 B is 1,176 B against 182 B of member strings. An island needs 84 members before the profile is the smaller of the two, which 3.7 % of real islands reach.

Parameters:

penalties (list[dict[str, int]]) – One dict per frame column, mapping symbol to a non-negative integer penalty.
width (int) – Frame width; equals len(penalties).
c (int) – Frame column where the gap block opens.

Example

>>> members = ["CASSLGQAYEQYF", "CASSLGQGYEQYF", "CASSLGQAYEQYF"]
>>> p = IslandProfile.fit(members)
>>> p.score(p.consensus())
0
>>> p.score("CASSLGQAYEQYF") <= p.score("CASSLGQGYEQYF")
True

classmethod fit(members, c=None, lam=1000, pseudocount=0.5, gap='-')[source]#

Fit a profile to an island’s members.

Parameters:

members (Sequence[str]) – The island. Must be non-empty. The frame width is the longest member.
c (int | None) – Frame column for the gap block. None picks the column minimising summed column entropy – the frame the members themselves prefer. On real islands the mode lands at 6, where crystal structures put the block.
lam (int) – Score resolution. Scores are compared against a control-calibrated cutoff, so any monotone rescaling cancels; lam only controls integer rounding.
pseudocount (float) – Added to every symbol count, so an unseen residue is expensive but finite.
gap (str) – The gap symbol used in the frame.

Returns:

A fitted IslandProfile.

Raises:

ValueError – If members is empty, lam is negative, pseudocount is not positive, or an explicit c exceeds the shortest member’s length.

Return type:

IslandProfile

consensus(gap='-')[source]#

The zero-penalty sequence: each column’s most frequent symbol, gaps stripped.

This is the centre of the ball the profile defines. score(consensus()) == 0.

Parameters:: gap (str)
Return type:: str

score(seq, gap='-')[source]#

Penalty of seq in this island’s frame; 0 on the consensus.

A sequence that does not fit the frame – longer than width, or shorter than c – cannot be embedded and scores UNREACHABLE. That is a rejection, not an error: it must still count as a scored sequence when a cutoff is calibrated against a control.

Parameters:

seq (str)
gap (str)

Return type:

int

score_batch(seqs, gap='-')[source]#

score() for many sequences. Score the whole control this way, then hand the result to seqtree.thetas_from_scores() for a calibrated cutoff.

Parameters:

seqs (Iterable[str])
gap (str)

Return type:

list[int]

Seed E-values#

E-values for shared k-mer seeds in the variable core of an anchored loop.

The germline flanks of a CDR3 junction carry almost no evidence: an exact N-terminal 4-mer is shared by 31.0% of the 250k human TRB control repertoire (CASS alone by 56.5%), and a C-terminal 4-mer by 14.1%. A central 4-mer is shared by 0.080% – about 386x more selective. Same four residues, ~2e4-fold different evidence.

So the significance of a shared seed has to be computed, not assumed. For a seed w and a target set of N sequences, the expected number of chance targets sharing w with the query is

\[E_{\mathrm{seed}}(w) = N \cdot n_C(w) / M\]

where n_C(w) counts control sequences containing w and M = |C|. Under the null the target draw is independent of the query, so conditioning on “the query contains w” is vacuous – do not square the probability. Occurrence-weighted over the bundled control with N = 1e5, the median E_seed of a central k-mer is 20.8 (k=4), 2.0 (k=5), 0.40 (k=6). A typical shared central 4-mer is therefore not significant, but 4.9% of them are; the median crosses E_seed < 1 at k=6.

This cannot be modelled. The residual KL divergence of the empirical central-k-mer distribution from a fitted background is 0.85 / 2.28 / 5.46 bits (independent per-position), 0.49 / 1.77 / 4.79 (Markov-1) and 0.43 / 1.60 / 4.48 (Markov-2) at k = 4 / 5 / 6 – growing with k, because D-gene germline runs (GGG, LAGG, SGGG) correlate. Count directly.

Scope. Seeds buy precision, not recall. Among same-epitope VDJdb pairs that lie in different sequence islands (i.e. beyond the reach of any anchored alignment), only 0.5% share a central 4-mer at all and 0.0% share a 6-mer – a real 4x enrichment over cross-epitope pairs, but negligible coverage. Use this to decide whether a seed you already found is meaningful, and to prune uninformative seeds before gathering. Do not expect it to connect distant relatives.

seqtree.seeds.core_kmers(seq, k, flank=4)[source]#

The distinct k-mers of seq’s variable core, excluding flank residues at each end.

flank=4 is the junction default: it drops the conserved Cys plus the first three germline residues at the 5’ end, and the conserved Phe/Trp plus three at the 3’ end. Sequences with a core shorter than k yield nothing.

Example

>>> sorted(core_kmers("CASSLGQAYEQYF", 4))   # core is 'LGQAY'
['GQAY', 'LGQA']

Parameters:

seq (str)
k (int)
flank (int)

Return type:

set[str]

class seqtree.seeds.SeedIndex(seqs, k=5, flank=4)[source]#

Bases: object

Inverted index over the core k-mers of a sequence set.

Built over a control repertoire it calibrates seed significance (evalue(), significant()). Built over a target set it gathers candidates (gather()). Counts are per sequence, not per occurrence: the event is “sequence x contains w”.

Parameters:

seqs (Iterable[str])
k (int)
flank (int)

classmethod from_index(index, k=5, flank=4)[source]#

Build from an existing Index, e.g. the one from load_control().

Parameters:

index (Index)
k (int)
flank (int)

Return type:

SeedIndex

count(seed)[source]#

Number of indexed sequences whose core contains seed.

Parameters:: seed (str)
Return type:: int

evalue(seed, n_target)[source]#

E_seed = n_target * count(seed) / len(self).

The expected number of sequences in a target set of size n_target that share seed with the query by chance. Below 1, the shared seed is itself evidence. Uses the rule of three (3/M) when the seed is absent from the control, matching seqtree.evalue.

Parameters:

seed (str)
n_target (int)

Return type:

float

seed_evalues(query, n_target)[source]#

E_seed for every core k-mer of query, keyed by seed.

Parameters:

query (str)
n_target (int)

Return type:

dict[str, float]

significant(query, n_target, alpha=1.0)[source]#

Core k-mers of query with E_seed < alpha, rarest first.

alpha=1.0 keeps seeds expected to arise fewer than once by chance.

Parameters:

query (str)
n_target (int)
alpha (float)

Return type:

list[str]

union_evalue(query, n_target, seeds=None)[source]#

E-value for sharing at least one of query’s seeds with a target sequence.

The query carries several overlapping core k-mers, so the per-seed E-values form a family. Rather than correct for that multiplicity, count the union directly: it is a single measurable set, so E = n_target * |union of postings| / M is exact and needs no correction constant. (A Boole bound over the seeds is only ~7% loose here, because overlapping string seeds gather nearly disjoint background sets – but it is unnecessary.)

Parameters:

query (str)
n_target (int)
seeds (Sequence[str] | None)

Return type:

float

gather(query, seeds=None)[source]#

Indices of sequences sharing at least one of query’s core k-mers.

Pass seeds=control.significant(query, ...) to gather on informative seeds only: an N-terminal 4-mer would otherwise pull in nearly half the repertoire.

Parameters:

query (str)
seeds (Sequence[str] | None)

Return type:

set[int]

Layout and anchors#

Pluggable index layout for epitopes/CDR3: anchor specs, masking, k-mer extraction.

Homology that drives TCR cross-reactivity is a shared central, TCR-facing motif (Dolton et al., Cell 2023: one HLA-A*02:01 TCR sees EAAGIGILTV / LLLGIGILVL / NLSALGIFST via x-x-x-A/G-I/L-G-I-x-x-x), not the MHC anchors. We therefore mask anchor positions and search anchor-masked k-mers. Anchor positions are parametrized (presets per class, overridable); MASK/GAP are spare symbols in the amino-acid alphabet so masked positions are first-class characters.

class seqtree.layout.AnchorSpec(cls, anchors=())[source]#

Bases: object

1-based anchor positions; negatives count from the C-terminus (-1 == last).

Parameters:

cls (str)
anchors (tuple)

cls: str#

anchors: tuple = ()#

resolve(length)[source]#

0-based anchor indices for a peptide of the given length.

Parameters:: length (int)
Return type:: set

seqtree.layout.spec_for(cls, override=None)[source]#

Resolve the anchor layout for a presentation class, honouring an explicit override.

Parameters:

cls (str) – One of the keys of DEFAULTS ("mhc1", "mhc2", "cdr3").
override (AnchorSpec | None) – Use this spec instead of the default when given.

Returns:

The AnchorSpec to apply.

Raises:

ValueError – If cls is unknown and no override is supplied.

Return type:

AnchorSpec

seqtree.layout.mask_anchors(pep, spec)[source]#

Replace anchor positions with MASK so they don’t drive homology.

Parameters:

pep (str)
spec (AnchorSpec)

Return type:

str

seqtree.layout.kmers(pep, k, spec=None)[source]#

Anchor-masked contiguous k-mers. Peptides shorter than k yield the whole (masked) peptide as a single token so they remain searchable.

Parameters:

pep (str)
k (int)
spec (AnchorSpec | None)

Return type:

list[str]

seqtree.layout.presentation_features(pep, cls, register='anchored')[source]#

Short binding-motif signatures (anchor / pocket residues) for the reverse problem – peptide -> presenting allele. These KEEP the anchors and drop the TCR-facing positions (the opposite of kmers()), so peptides binding the same allele share a signature.

class I: N-pocket P1-P3 + C-pocket P(Ω-1),PΩ -> one 5-residue signature. class II: core anchors P1,P4,P6,P9 of the 9-mer core. register='anchored' (the register trick) picks the single best register by _core_anchor_score so the signature is consistent across peptides of an allele; register='all' keeps every register (register-agnostic, noisier).

Parameters:

pep (str)
cls (str)
register (str)

Return type:

list[str]

seqtree.layout.weight_profile(length, spec, mode='tcr_facing', hotspot=(), hotspot_weight=1)[source]#

Per-position weights for a PositionalMatrix over a length-length frame.

mode=’tcr_facing’: anchors -> 0 (free), others -> 1 (optionally up-weight a central hotspot, e.g. class-I P4-P7). mode=’presentation’: the inverse – anchors -> 1, TCR-facing positions -> 0 – for allele assignment.

Parameters:

length (int)
spec (AnchorSpec)
mode (str)
hotspot (tuple)
hotspot_weight (int)

Return type:

list[int]

Epitope (pMHC) search#

Peptide-MHC epitope homology search and neoantigen molecular-mimicry discovery.

Homology is a shared TCR-facing k-mer (anchors masked); see seqtree.layout. Built on the C++ KmerIndex (seed-and-gather). MHC restriction is a payload filter; significance (E-values) lives in seqtree.pmhc_evalue.

class seqtree.pmhc.EpitopeHit(epitope: str, mhc: str, cls: str, shared_kmers: int, score: int, gene: str = '', species: str = '')[source]#

Bases: object

Parameters:

epitope (str)
mhc (str)
cls (str)
shared_kmers (int)
score (int)
gene (str)
species (str)

epitope: str#

mhc: str#

cls: str#

shared_kmers: int#

score: int#

gene: str = ''#

species: str = ''#

class seqtree.pmhc.PMHCStore(k=4, anchor_overrides=None)[source]#

Bases: object

Searchable epitope store, partitioned by MHC class.

classmethod from_records(records, k=4, anchor_overrides=None)[source]#: records: iterable of dicts with epitope, mhc (or mhc_a), mhc_class[, gene, species].

classmethod from_pmhc(path, classes=('mhc1', 'mhc2'), species=None, k=4, anchor_overrides=None)[source]#: Stream the isalgo/pmhc_data TSV(.gz).

size(cls)[source]#: Number of indexed epitopes for a presentation class ("mhc1" / "mhc2").

search_homologs(query, cls, mhc=None, max_subs=1, matrix='', min_shared=1, exclude_self=True, threads=0)[source]#: TCR-facing homologs of query in class cls, optionally restricted to mhc.

assign_allele(query, cls, top=5)[source]#

Rank alleles by how typical the query’s ANCHOR signature is among each allele’s presented peptides (a lightweight presentation prior, not NetMHCpan). Returns [(allele, score, n_match, n_allele)] sorted by score desc.

score is the log-odds of presentation by the allele vs the marginal background, so it doubles as a non-binder filter: score <= 0 means the anchors are not enriched for that allele (a non-binder for it), and if every allele scores <= 0 the peptide binds nothing in the panel. Class II is promiscuous, so expect several alleles with score > 0 (multi-label).

seqtree.pmhc.spec_has_anchors(spec)[source]#: Does this layout pin any anchor positions? Class II specs do not.

seqtree.pmhc.build_kmer_index(peptides, cls='mhc1', k=4, anchor_overrides=None)[source]#: Build a bare KmerIndex over a peptide set (for find_mimics / custom corpora).

seqtree.pmhc.find_mimics(neoantigen, self_set, bacterial_sets=None, control=None, cls='mhc1', k=4, max_subs=1, matrix='', min_shared=1, anchor_overrides=None, threads=0)[source]#

Discover TCR-facing molecular mimics of a neoantigen, with presentation-aware E-values.

All inputs are assumed presented by the same (compatible) MHC – cross-allele groove similarity is out of scope. self_set is the host presented peptidome (same allele); bacterial_sets maps organism -> presented peptides; control is the per-allele presented background (defaults to self_set). Returns one entry per source with the homolog hits, the expected count E, and the enrichment p-value.