API Reference#

Batch CPU budgets#

threads=1 runs inline. threads=0 selects the detected available CPU budget: Linux affinity, root cgroup v2 quota and SLURM allocation; macOS physical cores; Windows process affinity. Legacy nonpositive values also select automatic mode. Positive counts bound workers; a batch launches at most its query count. A request above the detected allocation raises ValueError; explicit counts are never clamped. The caller performs one share of the work: a budget of B uses B compute threads in total, creating only B-1 additional OS threads when there is enough work. Pools are joined before the call returns; use one batch rather than an outer pool of internally threaded searches. Cgroup layouts outside the detected root quota and Windows CPU-rate job limits require an explicit thread budget.

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.

alignment_max_batch(self, queries: collections.abc.Sequence[str], params: seqtree._core.SearchParams, similarity: collections.abc.Sequence[collections.abc.Sequence[float]], similarity_symbols: str, max_indels: int, *, query_target_offsets: collections.abc.Sequence[int], query_target_ids: collections.abc.Sequence[int], reference_target_offsets: collections.abc.Sequence[int], reference_target_ids: collections.abc.Sequence[int], target_class_ids: collections.abc.Sequence[int], reference_class_offsets: collections.abc.Sequence[int], reference_class_ids: collections.abc.Sequence[int], query_group_ids: collections.abc.Sequence[collections.abc.Sequence[int]], reference_group_ids: collections.abc.Sequence[collections.abc.Sequence[int]], group_offsets: collections.abc.Sequence[collections.abc.Sequence[collections.abc.Sequence[float]]], reference_other_counts: collections.abc.Sequence[int], reference_multiplicities: collections.abc.Sequence[int], coefficients: collections.abc.Sequence[float], threads: int = 0, exclude_exact: bool = False) → seqtree._core.AlignmentMaxResult#: Fused numeric requested-label max of cloglog(offset+c3*(rawNW-maxself)+gapcoef*g)*finite information. Candidate scope is independent of float64 global NW; selection gap=max(abs(similarity)); strict ties use M/X/Y for M and X, M/Y/X for Y, terminal M/X/Y. Group ID -1 is unsupported, not imputed. Target filtering precedes NW. Use bounded query chunks against the same immutable Index. GIL released.

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

edit_histogram_batch(self, queries: collections.abc.Sequence[str], params: seqtree._core.SearchParams, threads: int = 0, exclude_exact: bool = False, max_indels: int = -1) → list[list[int]]#: Exact reference counts by total n_subs+n_ins+n_dels, bins 0..max_total_edits. Requires explicit positive max_total_edits and uncapped seqtm/auto all-hit search. max_indels>=0 also caps total insertions+deletions; -1 preserves directional-only caps. Duplicates count separately; exclude_exact retains bin 0 with count zero. threads<=0 uses the detected available CPU budget. Releases the GIL; one reusable hit buffer per worker, no all-query hit vectors.

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.

ref_seqs(self) → list[str]#: Return all reference strings in ref_id order in one native call, preserving duplicates.

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_batch_consume(self, queries: collections.abc.Sequence[str], params: seqtree._core.SearchParams, consume: object, threads: int = 0, chunk_size: int = 32, max_hit_bytes: int = 67108864) → None#: One population/pool, serialized consume(start, stop, owned_buffer) callbacks. Read-only uint32 rows: query_index, ref_id, score, n_subs, n_ins, n_dels. Every query is complete; block order may vary across threads. Output bytes per block never exceed max_hit_bytes; an oversized query raises, never truncates. Native worker scratch also holds the current block of hits. Consumer-retained buffers remain owned; callers control their retention. Parameters are snapshotted.

search_paired_batch_consume(self, other: seqtree._core.Index, queries: collections.abc.Sequence[str], other_queries: collections.abc.Sequence[str], params: seqtree._core.SearchParams, consume: object, threads: int = 0, chunk_size: int = 32, max_hit_bytes: int = 67108864, skip_empty: bool = False) → None#: One pool with two Searchers; consume(start, stop, owned_buffer, other_owned_buffer). Both six-column buffers cover the same complete query block. Combined output bytes are bounded by max_hit_bytes; block order may vary. Same-alphabet indexes and equal-length query populations required. Parameters are snapshotted. skip_empty=True emits no hits for empty query axes; other axes remain searched.

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

For edit-count summaries, Index.edit_histogram_batch(queries, params, threads=0, exclude_exact=False, max_indels=-1) returns one list of exact counts per query, indexed by total substitutions, insertions and deletions. Set params.max_total_edits explicitly to a positive value; bins run from zero through that cap. Reference duplicates count separately. exclude_exact=True retains bin zero with count zero. Nonpositive threads uses the detected available CPU budget, matching the other batch APIs. This requires seqtm or auto and all-hit mode; the search uses one reusable native hit buffer per worker rather than materializing all hits in Python. Counts are not cumulative; sum the bins up to a radius to obtain its ball count. The optional final max_indels argument limits total insertions plus deletions; -1 preserves the directional-only limits. Bounded unit-cost profiles reuse bit-vector rejection and exact constrained dynamic programming; other profiles use the generic search worker. Existing eligible count paths are retained when the bound is redundant. Explicit seqtrie remains unsupported.

For full hits with bounded output, use Index.search_batch_consume(queries, params, consume, threads=0, chunk_size=32, max_hit_bytes=64*1024**2). It searches the complete population with one native pool. Serialized callbacks, which may run on native worker threads, receive (start, stop, owned_buffer) with read-only uint32 columns query_index, ref_id, score, n_subs, n_ins, n_dels. Every query is complete; block order may vary across threads, while hit order within each query is unchanged. Each output buffer obeys the explicit byte ceiling. An oversized single query raises rather than truncating its hits. Workers also retain the current block’s native hit scratch. Consume and release buffers to bound total memory; buffers explicitly retained by a callback remain owned. Parameters and matrices are snapshotted, callbacks are serialized, and workers join before errors propagate.

For coordinated searches over two same-alphabet indexes, use Index.search_paired_batch_consume(other, queries, other_queries, params, consume, threads=0, chunk_size=32, max_hit_bytes=64*1024**2, skip_empty=False). Both equal-length query populations run in one pool with two reusable Searchers per worker. Serialized callbacks receive (start, stop, owned_buffer, other_owned_buffer) for the same complete query block. The byte ceiling applies to the sum of both output buffers; ownership, snapshots and error cleanup follow the single-index consumer contract. This API provides no pairing or score fusion; the caller owns those interpretations. skip_empty=True emits no hits for an empty query axis while retaining searches on its partner; the default preserves ordinary empty-string search semantics.

Floating alignment target reductions#

Index.alignment_max_batch filters unit-cost search candidates by requested numeric target IDs, aligns each relevant reference sequence once, and reduces annotation scores without exporting the query-by-reference pair table. It releases the GIL and uses the same batch CPU budget as search. Build the index with unique sequences. Bound output memory with contiguous query chunks against that index. This is a numeric reduction API; callers supply their own group mappings and counts.

The positional arguments are queries, params, similarity, similarity_symbols, max_indels. similarity is a square finite float64 grid in the explicit similarity_symbols order; asymmetric grids are allowed. Every input residue must be represented. Candidate params must describe uncapped unit-cost all-hit seqtm/auto search. max_indels independently caps total insertion and deletion residues. Candidate existence does not constrain the selected alignment. Eligible short unit-cost profiles reuse the exact bounded-count predicate to collect full-index candidate IDs; other profiles retain the existing searcher. Global candidate counts are preserved before target filtering.

Global Needleman–Wunsch maximizes U - lambda*g, where U sums nongap similarities, g counts gap residues including terminal gaps, and lambda = max(abs(similarity)) > 0. Strict predecessor ties use M/X/Y for the match and insertion states, M/Y/X for the deletion state, and M/X/Y at the terminal cell. The normalized feature is U - max(H_query, H_reference), where H is the sum of diagonal self similarities. No integer scaling or Gram transform is used.

Required keyword arrays describe the annotations:

query_target_offsets (Q+1) and query_target_ids (M) are requested-target CSR; repeated targets within a query are rejected.
reference_target_offsets (R+1) and reference_target_ids (L) map each reference to annotation lanes. target_class_ids maps each target to a class.
reference_class_offsets (R+1) and reference_class_ids enumerate each reference’s unique class roles, independently of target annotations.
query_group_ids (C by Q), reference_group_ids (C by L), and group_offsets (C small rectangular matrices) describe additive components. Each matrix is indexed by query and reference group ID; -1 is unsupported.
reference_other_counts (L) supplies distinct neighboring references with another class, and reference_multiplicities (L) supplies positive observation counts. These counts are supplied and checked, never estimated by this reducer.
coefficients contains [intercept, normalized_coefficient, gap_coefficient].

The final predictor is the intercept plus group offsets plus the fitted normalized and gap terms. Its complementary-log-log transform is multiplied by -log10((k+1)/N), where N is the number of unique reference/class roles and k is the supplied other-class count for that annotation. With exclude_exact=True, matching query/reference sequences are excluded only from candidate hits. Reference/class roles and supplied other-class counts remain fixed. The native module reports this contract through _core.exact_exclusion_hits_only; version strings alone do not establish support. Inconsistent model counters raise instead of clipping. These supplied-count semantics do not make the output a calibrated probability.

The result’s arrays() method returns read-only owned zero-copy buffers. Per requested target: scores, raw_similarity, normalized_similarity, pfit, informativeness, gap_counts, best_ref_ids, best_lane_ids, counts, modeled_counts, aa_counts, and exact_counts. Counts distinguish observation multiplicity, supported multiplicity, distinct reference sequences, and raw exact multiplicity before exclusion. IDs are -1 if unavailable. Per query: candidate_counts counts all full-index candidates before filtering, nw_pairs counts actual alignments, and query_supported records group support. Supported no-hit scores are zero; unsupported scores are NaN while counts remain. Tied maxima choose the smallest annotation lane ID. Invalid CSR, IDs, duplicate references or roles, nonfinite coefficients, and overflow raise explicitly.

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.

record_distances(self) → list[list[int]]#: Sorted best substitution distances per distinct record for each query. Reduces offsets in one native call; rejects truncated results.

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.

The same question for substring hits in a text, where counts are per distinct record rather than per offset and the ball is a mismatch radius rather than a penalty. The module docstring of seqtree.text_evalue explains both choices.

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

Control-calibrated E-values for each query’s substring neighbourhood.

Parameters:

target – TextIndex to score hits in, e.g. a proteome.
control – background TextIndex, matched to target in composition and record-length distribution. See the module note – this choice is the assumption the calibration rests on.
queries – query strings, each at least target.k long.
max_subs – substitution radius defining the ball, searched identically in both texts.
threads – worker threads for the batch searches (0 = all cores).
exclude_exact – drop 0-mismatch hits from both counts. Set this when a query may occur verbatim in either text. A record with any exact occurrence is excluded, even if it also contains inexact occurrences.

Returns:

One dict per query with n_target, n_control, E, p_any, p_enrichment and rule_of_three. Counts are distinct records, not offsets.

Raises:

ValueError – If the control text is empty, or if a search was truncated, which would make a count a floor rather than a count.

Example

>>> import seqtree
>>> proteome = seqtree.TextIndex.build(["MKVLAAGIVGLNLGGK"], "aa", k=4)
>>> shuffled = seqtree.TextIndex.build(["GKLNVGALMGIVKGLA"], "aa", k=4)
>>> seqtree.text_evalues(proteome, shuffled, ["AAGIVGL"])[0]["n_target"]
1

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

Per-query mismatch cutoff achieving e_target against control.

A fixed radius is not calibrated: at max_subs=2 a low-complexity query picks up thousands of chance occurrences while a compositionally unusual one picks up none, so the radius that buys a given E-value has to be found per query. One control search at max_subs supplies every cutoff.

Parameters:

target – TextIndex the E-value is expressed against. Only num_refs is used.
control – background TextIndex.
queries – query strings.
max_subs – radius the control is searched at, and the ceiling for the returned cutoffs.
e_target – desired E-value, e.g. 0.05.
threads – worker threads for the batch search (0 = all cores).
exclude_exact – drop 0-mismatch hits, matching text_evalues().

Returns:

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

Raises:

ValueError – If max_subs is negative, or a search was truncated.

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, position_weights_by_length=None, query_group_ids=None, reference_group_ids=None, group_distances=None, query_position_weights=None, reference_position_weights=None, feature_offsets=None)[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).
position_weights_by_length – Optional mapping of full sequence lengths to nonnegative integer weights. Weights index the longer sequence; both axes must be covered. None or an empty mapping preserves unweighted defaults.
query_position_weights – Optional nonnegative integer rows, one per query and residue. Supply with reference_position_weights; incompatible with length profiles. Substitutions use max(query_weight, reference_weight) at aligned coordinates. Gap costs remain unchanged.
reference_position_weights – Per-reference rows with the same contract.
query_group_ids – Optional integer IDs, one per query; supply all group options together.
reference_group_ids – Optional integer IDs, one per reference.
group_distances – Nonempty rectangular nonnegative integer penalty table. Added to each pair score in the same units; no symmetry/zero diagonal imposed.
feature_offsets – Optional PairedFeatureOffsets descriptor. Adds paired-code costs averaged at shared depth, ties-even rounded, then scaled.
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++. Evaluated only for unequal input lengths and valid block starts; it must be a pure nonnegative function. 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.ball_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_extend=1, gap_prior=None, alphabet='aa', threads=0, exclude_exact=False, position_weights_by_length=None, query_group_ids=None, reference_group_ids=None, group_distances=None, query_position_weights=None, reference_position_weights=None, feature_offsets=None)[source]#

Exhaustive bounded ball rows, reference order, unset edit fields; no Q*N matrix.

One integer cutoff per query. Negative cutoffs yield empty rows. Optional integer positional weights index the longer full sequence, symmetrically, at every length. Optional additive group distances match score_matrix(); exact exclusion remains full encoded identity, irrespective of group distance. This reduces output memory; it is not a trie-based candidate filter.

seqtree.gapblock.group_audit_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_extend=1, gap_prior=None, alphabet='aa', threads=0, exclude_exact=False, position_weights_by_length=None, query_group_ids=None, reference_group_ids=None, group_distances=None, query_position_weights=None, reference_position_weights=None, feature_offsets=None, query_partition_ids=None, reference_partition_ids=None)[source]#

Stream component-ball count, full-score count and bounded nearest in one batch.

Geometry/offset options match ball_batch. Component excludes the nonnegative group offset; full score includes it. For each query threshold, return three lists: component count, full count and minimum full score among component hits. Empty component balls have nearest=-1. Duplicate references count separately; full encoded exact puncturing is identical. Optional nonnegative partition IDs restrict each query to references with the same ID. Both arrays are required; unmatched partitions return zero counts and nearest=-1. Partition restriction precedes scoring and is independent of additive group/feature distances. O(Q) output; no sparse hit objects.

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

Exhaustive single-gap-block top-k penalties, one native batch, GIL released.

Return one list of Hit(ref_id, score) per query, ordered by (score, ref_id). At most k reference ROWS survive; repeated reference sequences remain distinct rows. Hit edit fields are unset (zero), not an alignment decomposition. Empty references yield empty query rows.

Arguments match score_matrix(), with positive integer k added. Exact exclusion compares case-insensitive alphabet-encoded sequences BEFORE selection; a nonidentity zero penalty remains eligible.

Shared storage holds encoded inputs and the existing optional prior cube. Each worker holds O(k + max_sequence_length) scratch; results take O(Q*k), never a dense Q*N matrix. This remains exhaustive O(Q*N*sequence_length) scoring; no radius predicate, approximation or statistical model is added.

seqtree.gapblock.count_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_extend=1, gap_prior=None, alphabet='aa', threads=0, exclude_exact=False, position_weights_by_length=None, query_group_ids=None, reference_group_ids=None, group_distances=None, linear_mass=False, query_position_weights=None, reference_position_weights=None, feature_offsets=None, shared_thresholds=False)[source]#

Count all reference scores <= each query’s requested integer thresholds.

Optional positional weights and additive group distances match score_matrix(). Exact exclusion remains full encoded identity, irrespective of group distance. thresholds has one row per query; rows may be empty, unsorted or contain duplicates. With shared_thresholds=True, supply one flat cutoff list; it is validated and sorted once for the entire batch. Return integer counts in their original order, including ALL score ties. Negative cutoffs count zero. Same scorer/exclusion as topk_batch(); one native batch with O(Q*S) output, no dense score matrix. Worker scratch is O(max_sequence_length + maximum_threshold_row_length). With linear_mass=True, return (counts, masses) from the same traversal; each mass is sum(max(0, threshold - score)) over nonexcluded references. Boundary ties still count but contribute zero mass. Scores are nonnegative; integer mass is uint64, with overflow bounds checked before workers. Negative and zero thresholds have zero mass. No normalization or statistical model is applied; default counts and paired APIs remain unchanged.

seqtree.gapblock.paired_topk_batch(queries_a, queries_b, refs_a, refs_b, k=10, matrix=None, gap_open=None, gap_extend=1, gap_prior=None, alphabet='aa', threads=0, exclude_exact=False)[source]#

Top-k linked two-view rows, ordered by (max(score_a,score_b), ref_id).

Return one list of (ref_id, score_a, score_b) tuples per linked query. Both query axes and both reference axes must have matching lengths. Exclusion removes only rows identical in BOTH views; one-view identities remain eligible. Scoring arguments and memory bounds match topk_batch(). No cross-row combinations, payloads or statistical assumptions are introduced.

seqtree.gapblock.paired_count_batch(queries_a, queries_b, refs_a, refs_b, thresholds, matrix=None, gap_open=None, gap_extend=1, gap_prior=None, alphabet='aa', threads=0, exclude_exact=False)[source]#

Count all linked two-view rows with max(lane penalties) <= each cutoff.

Shape, threshold order, ties and negative-cutoff semantics match count_batch(); exact exclusion is the both-view predicate described in paired_topk_batch(). One native batch; no dense score matrix.

seqtree.gapblock.paired_sum_count_batch(queries_a, queries_b, refs_a, refs_b, thresholds, matrix=None, gap_open=None, gap_extend=1, gap_prior=None, alphabet='aa', threads=0, exclude_exact=False, position_weights_by_length_alpha=None, position_weights_by_length_beta=None, query_group_ids_alpha=None, query_group_ids_beta=None, reference_group_ids_alpha=None, reference_group_ids_beta=None, group_distances_alpha=None, group_distances_beta=None, query_position_weights_alpha=None, query_position_weights_beta=None, reference_position_weights_alpha=None, reference_position_weights_beta=None, feature_offsets_alpha=None, feature_offsets_beta=None, shared_thresholds=False)[source]#

Return (counts, linear masses) for linked rows under summed lane distances.

Shared matrix/gap costs/prior and separate full-coordinate lane weights and group tables use the existing kernel. Sum distance is computed in int64; thresholds remain int32. Count every score <= cutoff and sum max(0, cutoff - score); ties count with zero mass. Negative cutoffs return zero. Preserve duplicate reference observations and input threshold order. Exact exclusion removes a row only when both encoded full strings match. One GIL-released native traversal, O(Q*S) outputs and O(max_length + max_threshold_row) worker scratch; no matrix/edges. shared_thresholds=True accepts one flat cutoff list validated and sorted once for all queries. The paired-max API is unchanged.

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]

class seqtree.gapblock.PairedFeatureOffsets(*args, **kwargs)#

Bases: object

Immutable factored paired-feature offsets; validated by each batch.

seqtree.gapblock.grouped_count_batch(queries, refs, thresholds, *, reference_target_offsets, reference_target_ids, target_peptide_ids, reference_exclusive_peptide_ids, query_gene_ids, reference_gene_ids, chunk_size=64, matrix=None, gap_open=None, gap_extend=1, gap_prior=None, alphabet='aa', threads=0, exclude_exact=False, position_weights_by_length=None, query_group_ids=None, reference_group_ids=None, group_distances=None, query_position_weights=None, reference_position_weights=None, feature_offsets=None, query_withheld_peptide_ids=None, reference_peptide_offsets=None, reference_peptide_ids=None, score_quantum=1, alignment_mode='gapblock', alignment_scale=1, trim_left=0, trim_right=0, additional_feature_offsets=None, return_hits=False, query_target_ids=None, query_reference_ids=None, max_hit_bytes=67108864)[source]#

Yield bounded label curves from one immutable prepared scoring population.

return_hits=True instead yields chunk.hits[H,3] read-only uint32 rows (global query index, reference index, raw distance) at thresholds[-1], ordered by query then reference; counts/sizes/exact_counts are empty. Requires one query_target_ids entry per query and scans only that target’s CSR references. Optional query_reference_ids keeps reference indices strictly greater than each query’s supplied index (-1 keeps all). Exact exclusion still compares full sequences. Supports the single non-withheld population API only. max_hit_bytes bounds actual hit output per chunk (default64MiB); exceeding it raises without dropping hits. Row assembly plus final output temporarily uses at most roughly three times that cap, excluding inputs and worker DP scratch.

alignment_mode defaults to gapblock. Global/local reuse the existing affine similarity recurrence and require a raw similarity matrix. For trimmed sequences q/r, junction distance is alignment_scale * max(0, sum(diag(q)) + sum(diag(r)) - 2*similarity(q,r)), using raw diagonal sums rather than self-alignment scores. Nonnegative trim_left/trim_right remove residues for scoring only; exact exclusion and offsets retain full input identities. Gap costs are raw similarity magnitudes in these modes. Position weights and gap_prior are rejected. Paired lanes compute distances independently before summation. Scratch rows are reused per worker; no dense query/reference score matrix is produced.

additional_feature_offsets is an optional iterable of PairedFeatureOffsets. Every term is independently ties-even rounded/scaled, then added to the primary feature and group offsets. All descriptors are snapshotted/validated once. Combined score/offset bounds must fit int32 per lane. In the paired API, use additional_feature_offsets_alpha/beta; gap_prior_alpha/beta independently override the shared gap_prior. None inherits the shared prior for that lane.

score_quantum is a positive int32. Floor the raw alignment score to a multiple of this quantum before adding group/feature offsets. Default1 leaves every score unchanged. Paired lanes floor independently before their sum. This option applies only to grouped APIs.

Cutoffs are one strictly increasing flat integer list. References are a caller’s deduplicated search-key axis; CSR offsets/IDs assign each key to zero or more target IDs. Per-reference exclusive peptide IDs come from the full label set, with -1 for shared/unselected labels; target_peptide_ids maps each target. Canonical gene IDs define the auxiliary same-gene subset independently of allele costs. Channels are whole, targets, same-gene targets, exclusive peptides. Each chunk owns read-only uint32 buffers: counts[C,1+2*T+P,K], sizes[C,1+2*T+P], exact_counts[C,T], with start/stop query offsets. numpy.asarray is zero-copy; no NumPy import/dependency is required. Exact counts retain raw multiplicity; exclude_exact masks hit curves by full sequence identity; model population sizes remain fixed. Optional query_withheld_peptide_ids is one ID/query (-1 means none) and requires full reference_peptide_offsets/reference_peptide_ids CSR. Each sorted unique row includes all raw peptide labels, including sparse/unselected IDs in the target_peptide_ids namespace. Extra IDs beyond output peptide channels are valid. Selected target parents and static exclusives must agree with full CSR. Withholding removes every membership of that peptide, drops keys with no survivors, and recomputes exclusive surviving peptides before exact exclusion. Raw exact counts include only surviving targets. Full CSR without query IDs is rejected; None preserves the existing traversal. Encoding/geometry validation occurs once; each query/reference is scored once under the existing kernel, with GIL-released bounded workers. Memory is bounded by chunk_size, labels and cutoffs, never Q*labels. Retained chunks own their data.

seqtree.gapblock.grouped_paired_sum_count_batch(queries_a, queries_b, refs_a, refs_b, thresholds, *, reference_target_offsets, reference_target_ids, target_peptide_ids, reference_exclusive_peptide_ids, query_gene_ids, reference_gene_ids, chunk_size=64, matrix=None, gap_open=None, gap_extend=1, gap_prior=None, alphabet='aa', threads=0, exclude_exact=False, position_weights_by_length_alpha=None, position_weights_by_length_beta=None, query_group_ids_alpha=None, query_group_ids_beta=None, reference_group_ids_alpha=None, reference_group_ids_beta=None, group_distances_alpha=None, group_distances_beta=None, query_position_weights_alpha=None, query_position_weights_beta=None, reference_position_weights_alpha=None, reference_position_weights_beta=None, feature_offsets_alpha=None, feature_offsets_beta=None, query_withheld_peptide_ids=None, reference_peptide_offsets=None, reference_peptide_ids=None, score_quantum_alpha=1, score_quantum_beta=1, alignment_mode='gapblock', alignment_scale=1, trim_left=0, trim_right=0, additional_feature_offsets_alpha=None, additional_feature_offsets_beta=None, gap_prior_alpha=None, gap_prior_beta=None)[source]#

Grouped buffers under existing linked sum geometry; see grouped_count_batch.

Composite full-key exclusion requires both lanes equal. Canonical paired gene identities are caller-encoded tuple IDs, with equality requiring both genes. Shared matrix/prior and lane-specific weights/group/feature offsets retain the paired_sum_count_batch geometry. Counts, sizes and raw exact counts use uint32.

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.