API Reference#

Grouped by the job you are doing. For a starting point: Index searches a set of sequences, TextIndex searches a long text, and seqtree.pairwise aligns two sequences.

Searching a set of sequences#

The core loop: build an Index from your references, describe the search with SearchParams, get back Hit objects carrying a ref_id you map to your own payload.

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, ref_id: int, query: str, params: seqtree._core.SearchParams) → seqtree._core.Alignment#: Compute a global alignment between query and a reference, on demand.

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, queries: collections.abc.Sequence[str], params: seqtree._core.SearchParams, threads: int = 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).

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

ref_seq(self, ref_id: int) → str#: Return the reference sequence string for a reference id.

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

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

search_batch(self, queries: collections.abc.Sequence[str], params: seqtree._core.SearchParams, threads: int = 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, query: str, params: seqtree._core.SearchParams, k: int = 1) → list#: Return up to k best (lowest-score) hits for one query.

class seqtree.SearchParams(*args, **kwargs)#

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#

it is the only engine that enforces the per-type caps and reports an edit breakdown. Ask for 'seqtrie' by name when a score budget is the entire specification.

Type:: 'seqtm', 'seqtrie', or 'auto' (the default). 'auto' always resolves to seqtm

property gap_extend#: Cost of each additional residue in a gap run. Equal to gap_open means linear gaps; there is no separate mode flag.

property gap_open#: Cost of opening a gap. A run of length L costs gap_open + (L - 1) * gap_extend. The default of 1 is only right for unit cost – with any real matrix use 2 * matrix.scale().

property matrix#

a SubstitutionMatrix, or the name of a builtin ('blosum62', 'blosum45', 'blosum80', 'pam250', 'pam100', 'structural', 'identity'). '' (the default) means unit cost, where a score is a plain edit count. Setting a matrix makes gap_open matter: pass 2 * matrix.scale(), or gaps come out ~14x cheaper than substitutions.

Type:: Substitution matrix

property max_dels#: Maximum deletions (reference residues absent from the query). Enforced exactly by seqtm; seqtrie ignores it.

property max_ins#: Maximum insertions (residues in the query absent from the reference). Enforced exactly by seqtm; seqtrie ignores it. Indels widen the search frontier far more than substitutions do.

property max_penalty#

seqtrie prunes on penalty alone, so without a finite budget it walks the entire index.

Type:: Cap on the accumulated substitution/gap penalty. 0 means unset. Required when engine='seqtrie' is paired with a matrix

property max_subs#: Maximum substitutions. Enforced exactly by seqtm; seqtrie ignores it.

property max_total_edits#: Cap on the total edit count, independent of the per-type caps rather than clamped by them. 0 means no separate total cap – the sum of the three per-type caps applies. This is the only edit limit seqtrie honours.

property mode#: 'all' (the default) returns every hit inside the scope; 'top' sorts by (score, ref_id) and keeps the best max_hits.

property pos_matrix#: PositionalMatrix giving per-position penalties, or None. Applies only on the seqtm Hamming path and only when its width equals the query length. Setting it forces seqtm regardless of engine.

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#: Deletions in the best path. Exact under seqtm, always 0 under seqtrie.

property n_ins#: Insertions in the best path. Exact under seqtm, always 0 under seqtrie.

property n_subs#: Substitutions in the best path to this reference. Exact under seqtm, always 0 under seqtrie, which cannot see edit types.

property ref_id#: Index of the matched reference in the list passed to Index.build().

property score#: Accumulated penalty, always >= 0, with 0 an exact match. Under unit cost this is the edit distance; with a matrix it is the summed substitution and gap penalty. Lower is better.

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#: The query with - inserted at deletion columns. Same length as aligned_ref and ops.

property aligned_ref#: The reference with - inserted at insertion columns.

property ops#

M match, S substitution, I insertion, D deletion.

Type:: One character per alignment column

property score#: Alignment penalty, >= 0 and 0 for an exact match – the same convention as Hit.score, and the opposite of seqtree.pairwise.score(), which returns a signed similarity.

Searching a long text#

For a short query against a proteome or genome, where enumerating every window as its own reference is not affordable. One index answers every query length – see Text search.

class seqtree.TextIndex#

Exact k-mismatch (Hamming) search over a CONCATENATED reference text – a proteome, a genome, a transcript set. Unlike Index, which builds a trie over reference strings and so needs one index per query length, k here belongs to the index: ONE build answers every length and every max_subs. Full length, no gaps, no score in the predicate; the answer is exact, not a heuristic.

property alphabet#: ‘aa’, ‘nt’ or ‘iupac’.

build(refs: collections.abc.Sequence[str], alphabet: str = 'aa', k: int = 4, group_ids: collections.abc.Sequence[int] = []) → seqtree._core.TextIndex#: Build from whole records (NOT windows). k is the seed width; the table is direct-addressed, so alphabet_size**k buckets are allocated (24**4 = 331,776 for amino acids) and larger k is refused. group_ids optionally labels each record – gene ids, species, clusters – so hits can be folded onto them; seqtree does not know what a group means. Raises ValueError on a symbol outside the alphabet, naming the record.

property k#: Seed width this index was built with.

load(path: str, mmap: bool = True) → seqtree._core.TextIndex#: Load an index written by save(). With mmap the file is mapped rather than read, so several processes share one copy of the pages.

property num_refs#: Number of records in the text.

property num_residues#: Total residues, excluding the inter-record separators.

property num_unknown#: Text residues outside the alphabet (36 U in the human proteome, 33 in mouse). They are kept as holes – no hit may cross one – and reported here rather than dropped silently. A query containing one is refused.

ref_seq(self, ref_id: int) → str#: Return record ref_id as a string.

save(self, path: str) → None#: Write a flat, mmap-able index file.

search_batch(self, queries: collections.abc.Sequence[str], max_subs: int = 0, max_indels: int = 0, exclude_exact: bool = False, best_only: bool = False, group_by: bool = False, max_hits: int = 0, matrix: seqtree._core.SubstitutionMatrix | None = None, threads: int = 0) → seqtree._core.TextResult#

Find every position matching each query within max_subs substitutions and max_indels insertions/deletions (releases the GIL; threads=0 uses all cores). The two caps are independent: max_subs=2, max_indels=1 accepts two substitutions AND one gap, not three edits of any kind.

max_indels=0 is the pure Hamming predicate: every hit is exactly len(query) residues wide and carries per-column mismatches. Above 0 the match length varies, so read TextHit.length rather than assuming the query length, and mismatches is empty – the counts are reported, the individual substituted columns are not. Indel search also needs every seed block to be exact, so a query must be at least (max_subs + max_indels + 1) * k long; a shorter one raises rather than being answered incompletely.

best_only walks the distance upward and stops at the first shell with any hit, returning ALL of it. exclude_exact drops 0-mismatch hits. max_hits caps a query after sorting, so the best hits survive, and sets truncated. matrix only SCORES hits the predicate already accepted – it never changes which are returned, and it is ignored when max_indels > 0.

class seqtree.TextResult#

Results of a TextIndex batch, held as flat parallel arrays rather than one object per hit – a 445k-query run comes back as a handful of arrays. len(res) is the query count and res[i] builds the TextHit list for query i on demand, so iterating pairs up with the query list. For the whole batch use arrays() (zero-copy views) or to_numpy(). Hits are ordered (n_subs, ref_id, offset), stable across runs and thread counts.

arrays(self) → dict#: Every underlying array as a zero-copy ArrayView, keyed by name: query_begin, ref_id, offset, n_subs, score, mm_begin, mm_pos, mm_query_aa, mm_text_aa, truncated, group_begin, group_id, group_min_subs, group_n_hits.

groups(self, query: int) → list[tuple[int, int, int]]#: (group_id, min_subs, n_hits) per group reached by query i, sorted by group id. More than one row means the nearest parents disagree – the tie is a first-class output rather than something each caller re-derives.

property num_hits#: Total hits across every query.

to_numpy(self) → dict#: The same arrays as numpy views, sharing memory with this result. Requires numpy; it is imported on the call, never at import time.

property truncated#

1 if max_hits capped it. A cap that is invisible is a recall bug wearing a performance costume, so it is always reported.

Type:: One flag per query

class seqtree.TextHit#

One match of a query against the text. offset is the start within record ref_id; n_subs is the Hamming distance; score is the substitution score (0 unless a matrix was passed). mismatches lists (pos, query_aa, text_aa) – the PAIR, so a caller ranking by chemistry can tell L->I from L->D without re-fetching the window. Iterable as (ref_id, offset, n_subs, score).

property length#: Residues the match spans in the text, i.e. ref_seq(ref_id)[offset : offset + length]. Equal to the query length whenever there are no indels; len(query) + n_dels - n_ins in general.

property mismatches#: [(pos, query_symbol, text_symbol), ...], one per mismatch, pos 0-based within the query. Empty for an exact hit.

property n_dels#

text residues with nothing opposite them in the query. Always 0 unless the search passed max_indels > 0.

Type:: Deletions

property n_ins#

query residues with nothing opposite them in the text. Always 0 unless the search passed max_indels > 0.

Type:: Insertions

property n_subs#: Substituted positions, always <= the max_subs that was searched. 0 with n_ins and n_dels also 0 is an exact occurrence.

property offset#: Start of the occurrence within that record, 0-based. The matched substring is ix.ref_seq(ref_id)[offset:offset + len(query)].

property ref_id#: Index of the record this occurrence falls in, in the list passed to TextIndex.build().

property score#: Summed similarity over the mismatched positions under the matrix= passed to the search, or 0 when none was. Scoring never changes which hits return.

class seqtree.ArrayView#

A read-only 1-D view over one of a TextResult’s arrays. Exposes the buffer protocol, so numpy.asarray(v) and memoryview(v) wrap it without copying, and it keeps its TextResult alive.

property format#: The struct format code of the elements – 'i' for int32, 'I' for uint32, 'H' for uint16, 'b' for int8. The same code memoryview(v).format reports.

Scoring#

How a mismatch is priced. Penalties are non-negative and zero on a match, so a score is a distance: lower is better, everywhere in seqtree.

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)).

blosum45() → seqtree._core.SubstitutionMatrix#: BLOSUM45. More permissive than BLOSUM62 – built from more divergent alignments, so distant substitutions cost less. Amino acids only.

blosum62() → seqtree._core.SubstitutionMatrix#: BLOSUM62, the general-purpose default. Amino acids only; scale() is 14.

blosum80() → seqtree._core.SubstitutionMatrix#: BLOSUM80. Stricter than BLOSUM62 – built from closer alignments, so it separates near-identical sequences more sharply. Amino acids only.

from_similarity(grid: collections.abc.Sequence[collections.abc.Sequence[int]]) → 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()).

pam100() → seqtree._core.SubstitutionMatrix#: PAM100, the same model at a shorter evolutionary distance – stricter than PAM250. Amino acids only.

pam250() → seqtree._core.SubstitutionMatrix#: PAM250, an evolutionary model tuned for distant homology. Amino acids only.

penalty(self, 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) → 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, 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) → int#: Alphabet size this matrix scores over – 24 for amino acids, 4 for nucleotides. A matrix must match the alphabet of the index it is used with.

structural() → seqtree._core.SubstitutionMatrix#: A structure-derived matrix, grouping residues by side-chain shape and charge rather than by observed substitution frequency. Amino acids only.

unit(size: int) → seqtree._core.SubstitutionMatrix#: Identity: every mismatch costs 1, every match 0, so a score is a plain edit count. size is the alphabet size (24 for aa, 4 for nt); use len(alphabet_symbols(alphabet)). The only matrix valid for nucleotides.

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.

from_tables(size: int, width: int, data: collections.abc.Sequence[int], masked: collections.abc.Sequence[int] = []) → 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.

from_weights(base: seqtree._core.SubstitutionMatrix, weights: collections.abc.Sequence[int]) → 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, pos: int) → bool#: True if pos is free: mismatches there cost nothing and are not counted as substitutions. This is how an anchor position is excluded from the ball.

penalty(self, pos: int, a: int, b: int) → int#: Penalty for substituting symbol a with b at position pos. Symbols are single-character strings. Always >= 0, and 0 when a == b or pos is masked.

size(self) → int#: Alphabet size this matrix scores over.

width(self) → int#: Frame width in positions. A query must be exactly this long for the positional path to apply.

Dense matrices#

Every query against every reference in one GIL-released call, for when nothing can be pruned. Returned by seqtree.pairwise.score_matrix(), seqtree.distance.hamming_matrix(), and seqtree.gapblock.score_matrix().

class seqtree.ScoreMatrix#

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, i: int) → list[int]#: Row i as a list of penalties, one per reference.

property shape#: (n_rows, n_cols) – queries by references, the same shape numpy.asarray(sm).shape reports.

Seed-and-extend#

Candidate generation at million scale: match query k-mers, merge the posting lists, rank. Used by seqtree.pmhc.

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).

build(kmers_per_peptide: collections.abc.Sequence[collections.abc.Sequence[str]], alphabet: str = 'aa', allele_ids: collections.abc.Sequence[int] = []) → seqtree._core.KmerIndex#: Build the index. kmers_per_peptide[i] is the k-mer list for peptide i – produce it upstream (seqtree.layout.kmers masks anchors first), since which k-mers represent a peptide is a domain decision. allele_ids is an optional parallel tag per peptide, used by seed_and_gather’s allele_filter; pass [] for none. The k-mers need not all be the same width.

load(path: str) → seqtree._core.KmerIndex#: Read an index written by save(). Raises on a truncated file or a format version this build does not know.

num_kmers(self) → int#: Number of distinct k-mers across all peptides. The postings lists total more than this whenever a k-mer occurs in several peptides.

num_peptides(self) → int#: Number of peptides indexed – the length of kmers_per_peptide. Also len(index).

save(self, path: str) → None#: Write the index to path. Writes a temporary and renames it into place, so a concurrent reader never sees a partial file.

seed_and_gather(self, query_kmers: collections.abc.Sequence[collections.abc.Sequence[str]], params: seqtree._core.SearchParams, min_shared: int = 1, allele_filter: int = -1, threads: int = 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#: Lowest penalty over the k-mer matches that reached this peptide. A tiebreak within a shared_kmers group, not a full-length alignment score.

property peptide_id#: Index of the candidate peptide in the list the index was built from.

property shared_kmers#

candidates are ranked on it first.

Type:: How many distinct query k-mers reached this peptide. The seed evidence

Alphabets and batch helpers#

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.pairwise_batch(a: collections.abc.Sequence[str], b: collections.abc.Sequence[str], params: seqtree._core.SearchParams, alphabet: str = 'aa', threads: int = 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.

Significance – is this hit real?#

A score alone is not evidence: a germline-adjacent query collects neighbours by chance. These count the same ball in a background control and report an E-value, or invert it for the score cutoff that achieves a target false-positive rate. See E-values: is this hit real?.

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.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.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]

seqtree.evalue.evalue_result(n_target, n_control, n_ref, m_control)[source]#

The E-value record for one query, shared by every caller that computes one.

E = (N / M) * n_control, or the rule-of-three upper bound 3N / M when the control ball came back empty – an observed zero does not mean the background rate is zero, it means it is below roughly 3 / M.

Parameters:

n_target – hits in the target set for this query.
n_control – hits in the background control for this query.
n_ref – size of the target set, N.
m_control – size of the control, M. Zero yields an infinite E-value.

Returns:

dict with n_target, n_control, E, p_any, p_enrichment, rule_of_three.

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.

The same module also enumerates a Hamming ball rather than scoring a pair you already hold:

neighbourhood() – every sequence within r substitutions of one centre;
neighbourhood_union() – the same over many centres, each distinct sequence emitted once. Deduplication happens during generation, never as a pass over the sum(19*L_i) multiset;
union_size() – the cardinality alone, for sizing a job before running it.

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, union_size
>>> hamming("CASSLGQYF", "CASSPGQYF")
1
>>> levenshtein("kitten", "sitting")
3
>>> union_size(["CASSLGQYF", "CASSPGQYF"])      # 2 * 172 balls, 20 sequences shared
324

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

seqtree.distance.neighbourhood(seq, r=1, alphabet=None, include_self=True, shell=False)[source]#

Every sequence within r substitutions of seq: its closed Hamming ball.

Substitution only, so every member has len(seq) characters – Hamming distance is undefined across lengths (see hamming()), and an indel is a different question. Over a k-letter alphabet the ball holds sum((k-1)^d * C(L, d) for d in 0..r) sequences; at r = 1 that is 19*L + 1.

Parameters:

seq (str) – The centre.
r (int) – Radius in substitutions.
alphabet (str | None) – Substituting alphabet; None is the 20 standard amino acids.
include_self (bool) – Keep seq itself (the d = 0 member).
shell (bool) – Return only the members at distance exactly r, instead of the closed ball. The cognacy-retention profile is estimated per shell, not per ball.

Returns:

The distinct members, ordered by distance and sorted within each shell, so seq itself comes first when include_self.

Raises:

ValueError – If r is negative.

Return type:

list[str]

Example

>>> len(neighbourhood("CASSLGQYF"))                 # 19 * 9 + 1
172
>>> len(neighbourhood("CASSLGQYF", include_self=False))
171
>>> neighbourhood("A", 1, alphabet="ACGT")
['A', 'C', 'G', 'T']
>>> len(neighbourhood("CASSLGQYF", 2, shell=True))   # 19^2 * C(9, 2)
12996

seqtree.distance.neighbourhood_union(seqs, r=1, alphabet=None, include_self=True, shell=False)[source]#

Every sequence within r substitutions of any of seqs, each listed once.

This is the union of the per-sequence Hamming balls, not their concatenation. Near-duplicate centres – co-specific TCR junctions, an error family around one UMI – have balls that overlap heavily, so the difference is the dominant term rather than a correction. The dedup happens as the ball is walked, so the sum(19*L_i) multiset never exists.

Parameters:

seqs (Iterable[str]) – The centres. Duplicates and mixed lengths are fine; substitution-only means balls of different lengths cannot overlap.
r (int) – Radius in substitutions.
alphabet (str | None) – Substituting alphabet; None is the 20 standard amino acids. Characters of seqs outside it are still substitutable, so a centre containing X has |alphabet| neighbours at that position rather than |alphabet| - 1.
include_self (bool) – Keep the centres themselves. False drops every centre, including one that happens to be a neighbour of another – the r = 0 shell of a union is exactly set(seqs).
shell (bool) – Return only the members at distance exactly r from the nearest centre, instead of the whole closed ball.

Returns:

The distinct members, ordered by distance and sorted within each shell (so the centres come first). Empty input gives [].

Raises:

ValueError – If r is negative.
TypeError – If seqs is a single string – its characters would otherwise become the centres, silently. Use neighbourhood() for one sequence.

Return type:

list[str]

Example

>>> neighbourhood_union(["AA", "AC"], 1, alphabet="AC")
['AA', 'AC', 'CA', 'CC']
>>> a = neighbourhood("CASSLGQYF")
>>> b = neighbourhood("CASSPGQYF")                  # one substitution away
>>> len(a) + len(b), len(neighbourhood_union(["CASSLGQYF", "CASSPGQYF"]))
(344, 324)

seqtree.distance.union_size(seqs, r=1, alphabet=None, include_self=True, shell=False)[source]#

How many distinct sequences neighbourhood_union() would return.

Same walk, but the result list is never built or sorted – use it to size a job (an OLGA P_gen pass over the union, say) before committing to it. The dedup set is still held, so this bounds the output, not the peak memory.

Parameters:

seqs (Iterable[str]) – The centres.
r (int) – Radius in substitutions.
alphabet (str | None) – Substituting alphabet; None is the 20 standard amino acids.
include_self (bool) – Count the centres themselves.
shell (bool) – Count only distance-exactly-r members.

Returns:

len(neighbourhood_union(seqs, r, alphabet, include_self, shell)).

Raises:

ValueError – If r is negative.
TypeError – If seqs is a single string, as for neighbourhood_union().

Return type:

int

Example

>>> union_size(["CASSLGQYF", "CASSPGQYF"])          # 344 with double-counting
324
>>> union_size(["CASSLGQYF", "CASSLGQYF"])          # a repeat adds nothing
172

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.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.