seqtree 1.0.0 — TextIndex, nanobind bindings, first stable release

First stable release. From 1.0.0 the project follows semantic versioning: breaking changes need a major bump.

TextIndex — exact k-mismatch search over a long text

Index is a trie over reference strings, so asking a text question with it means enumerating every length-L window as its own string — and again for every distinct query length. The human proteome has 68,389,335 nine-mer windows; a query set spanning 45 distinct lengths costs 45 multi-gigabyte builds. That is what ran over 2 h without finishing and filled 225 GB of index cache downstream.

TextIndex keys on a k-mer seed table over the flat text, so k belongs to the index, not the query: one build answers every length and every max_subs.

ix  = seqtree.TextIndex.build(records, alphabet="aa", k=4, group_ids=gene_id_per_record)
res = ix.search_batch(peptides, max_subs=2, best_only=True, group_by=True, threads=0)

Exact, not heuristic, by one search scheme: split the query into b = min(m+1, L//k) disjoint blocks and probe block j's leading k-mer at radius c_j. It is lossless exactly when Σ c_j ≥ m − b + 1. Completeness is pinned by brute-force set equality over L 6–30 × max_subs 0–3 × k ∈ {3,4,5}, and by the answer being identical across k.

On the human proteome (147,506 records, 69,578,135 residues), one thread:

query length	max_subs	ms/query
9	2	1.28
12	3	0.39 (k=5)
15	3	0.14 (k=5)

Results come back as flat CSR arrays with zero-copy numpy views, mismatches as (pos, query_aa, text_aa) pairs, an optional fold onto caller-supplied group ids that makes a tie explicit, and a cap that is always reported.

nanobind bindings

The bindings moved from pybind11 to nanobind: smaller module, lower per-call overhead, and generated type stubs. The package has shipped py.typed since the beginning while providing no stubs at all, so every C++ symbol resolved to Any; _core.pyi now ships in the wheel. Externally visible behaviour is unchanged, including ScoreMatrix's buffer protocol and the exception mapping pinned by message text.

Also in this release

• Windows portability for the mmap paths, exercised by CI on all three platforms.
• Audit pass: one parallel_for replaces seven hand-rolled thread-pool copies (one of which had no exception plumbing at all); 59 previously undocumented members of the C++ binding surface are documented and gated by a test; psutil is out of the [bench] extra.
• Docs: two claims that no longer matched the code are corrected — engine="auto" never routed to seqtrie, and the first BLOSUM62 example used gap_open=8 against the library's own 2 * scale() = 28 rule.

Zero runtime dependencies. Wheels for CPython 3.10–3.13 on Linux x86-64, macOS arm64, and Windows x86-64.

Full detail in CHANGELOG.md · docs at antigenomics.github.io/seqtree

seqtree 0.7.0 — deduplicated Hamming-neighbourhood enumeration

Added

seqtree.distance now enumerates a Hamming ball, not just scores one.

• neighbourhood(seq, r=1, alphabet=None, include_self=True, shell=False) — the closed ball (19·L + 1 at r = 1 over the 20 standard residues).
• neighbourhood_union(seqs, ...) — the union over many centres, each distinct sequence emitted once.
• union_size(seqs, ...) — the cardinality, without building or sorting the result list.

Substitution only, fixed length. Deduplication happens during a multi-source breadth-first walk, so the Σ 19·L_i multiset is never materialised — the point, since near-duplicate centres overlap heavily. For 200 length-14 junctions all within distance 1 of a common centre the per-sequence balls double-count 41.7% (53,400 → 31,122); at spread 2, 4.5%; at spread ≥ 3, nothing.

shell=True returns the sphere at exactly r from the nearest centre. Shells partition the closed ball exactly: shells 0..r are pairwise disjoint, their union is the ball, and every member sits at distance exactly d from its nearest centre — so a per-shell quantity sums back to the per-ball one.

neighbourhood_union and union_size raise TypeError on a single string rather than iterating it: union_size("CASSLGQYF") would otherwise take the 8 distinct characters as centres and answer 20 instead of 172.

seqtree.__version__, read from the installed distribution metadata, so pyproject.toml is the single source.

Fixed

The documentation build stamped every page 0.6.1 while the package was 0.7.0 — docs/conf.py hand-copied the version. Both fields now read seqtree.__version__.

Gates

341 → 346 tests passed, 3 skipped · ctest 1/1 · sphinx -W clean · oracle and perf regression within threshold · sdist + wheel build, twine check PASSED.

v0.6.1

Fixed

• gapblock_matrix's out-of-alphabet error named only the bad symbol, not which sequence it
  came from ("symbol '_' is not in the alphabet") — unhelpful when a single malformed row (e.g.
  a junction_aa containing a legacy out-of-frame marker) is buried in a batch of hundreds of
  thousands. The message now names the side, index, and offending string:
  "queries[42] ('CASSIRS_YEQYF'): symbol '_' is not in the alphabet". No change to which symbols
  are accepted (still the standard 20 + B/Z/X/* for alphabet="aa").

seqtree 0.6.0 — the corrected MJ A-N contact energy

Fixed

• structural scored A–N contacts as if they did not interact, because the source table is
  corrupted there. The Miyazawa–Jernigan A–N contact energy was transcribed as 0.00 — the
  generator's comment excused it as a pair the source left unlisted. It is not unlisted. In
  MJ_Keskin_potentials.csv the lower triangle runs A-A, R-A, R-R, [N-A], N-R, …, and
  the N-A slot reads V,1 (mirrored 1,V), where 1 is a mangled residue symbol. The true
  value is 0.15. It is not a stray duplicate of V-N either — the V row separately lists
  N = 0.12.

  Substituting 0.00 for 0.15 understated A's interaction strength and overstated N's
  (q(A) −0.04100 → −0.03350, q(N) +0.04350 → +0.05100). Since structural is rank-1 by
  construction — every cell is a function of the 20 per-residue strengths — this moved the
  strong→weak ordering that the matrix exists to encode, swapping N and P:
  FWCLYMIVHGATNPRSQDEK → FWCLYMIVHGATPNRSQDEK.

  Four cells of the 24×24 grid change: A-W/W-A 6 → 5, and B-C/C-B 4 → 3 (B is
  mean(N, D), so N's shift carries into it). Scores from structural() change for sequences
  containing A, W, or B; all other matrices are untouched. N and P were near-tied, which is why
  a real correction to the ordering moves so few cells.

• The same comment claimed the source is "near- but not perfectly symmetric". It is perfectly
  symmetric — 0 asymmetric directed pairs across all 400. The symmetrisation in
  structural_grid() is a no-op guard, not a repair, and is now documented as such.

seqtree 0.5.0 — plain edit distances

Added — seqtree.distance

Plain, unweighted Hamming and Levenshtein string distances in C++ — unit costs, no substitution matrix, no gap model, no alphabet. When all you need is how many edits apart two strings are, you no longer build a SubstitutionMatrix or add python-Levenshtein / rapidfuzz. seqtree still needs nothing at runtime.

distance.hamming(a, b)	differing positions; equal length only (raises ValueError otherwise)
distance.levenshtein(a, b)	insertions + deletions + substitutions, each cost 1, O(min(m,n)) memory
distance.hamming_matrix(a, b, threads=0)	dense len(a) × len(b) int32, GIL released, zero-copy numpy
distance.levenshtein_matrix(a, b, threads=0)	same, for mixed-length sequences

Comparison is case-sensitive, byte for byte. For a weighted alignment (a matrix, affine gaps, local mode) use seqtree.pairwise instead. Verified against pure-Python oracles over random data.

Full changelog: see CHANGELOG.md.

v0.4.0 — pairwise alignment without BioPython, and two cache-corruption fixes

Pairwise alignment without BioPython — plus two cache-corruption fixes that were tagged 0.3.1 but never published.

Added

• seqtree.pairwise — Needleman–Wunsch and Smith–Waterman. Ordinary protein alignment on the raw log-odds scale, so reaching for BioPython is no longer necessary.

  pairwise.score(q, r, matrix, mode=...)	optimal score, O(min(m,n)) memory
  pairwise.align(...)	plus the aligned strings and ops
  pairwise.score_matrix(queries, refs, ...)	dense n × K, GIL released, zero-copy numpy
  pairwise.dist_matrix(...)	d = s(a,a) + s(b,b) − 2·s(a,b) — non-negative, zero on the diagonal

  mode="global" is Needleman–Wunsch, mode="local" is Smith–Waterman, and gap_open == gap_extend gives linear gaps — no separate mode. A gap run of length L costs gap_open + (L-1)·gap_extend, and global charges end gaps (true NW, not semi-global).

  It is a drop-in. Verified against Bio.Align.PairwiseAligner as an oracle across three matrices × fifteen gap/mode settings × sixty sequence shapes — zero disagreements — and on real germline V genes. BioPython is a test-only dependency; seqtree still has zero required runtime dependencies and never imports it.

  And faster, because there is no Python in the per-pair loop:

  sequence length	seqtree, 1 thread	seqtree, 16 threads	BioPython	speedup
  15 (a junction)	1.7 M pairs/s	20.1 M pairs/s	0.31 M pairs/s	65×
  90 (a germline V gene)	72 k pairs/s	893 k pairs/s	10 k pairs/s	87×

• SubstitutionMatrix.similarity(a, b) — the raw signed log-odds, alongside the existing non-negative penalty(a, b). The Gram transform pen = s(a,a) + s(b,b) − 2·s(a,b) is lossy (it forces the diagonal to zero), so both views are now stored.

• SubstitutionMatrix.blosum45() and .blosum80(), usable wherever a matrix name is accepted.

Fixed

• A cold cache shared by concurrent processes could hand back a half-written index. Index::save wrote straight into the destination, so a process that checked for the cache while another was still writing it loaded a truncated file and raised truncated or corrupt index — 10 times out of 10 on the 45 MB control. Saves now write a temporary and rename it into place (atomic on the same filesystem, POSIX and Windows alike).

  This is the first-use failure of any multi-process fan-out sharing ~/.cache: pytest-xdist, a Snakemake or Nextflow rule, a multiprocessing pool. Warm caches were never at risk, and CI matrix jobs (separate runners) were never affected.

• The control cache is content-addressed, so a stale cache can no longer be served silently. The old key named neither the alphabet, nor the seed, nor the source data — so two calls that must draw different samples shared one file, and an upgrade that changed the bundled control served the previous release's control from a warm cache. You no longer need to clear ~/.cache/seqtree when upgrading.

• Index.align compared raw characters instead of codec-encoded ones, so a lowercase query against an identical uppercase reference scored 12 under unit cost but 0 under a matrix, and both labelled identical residues as substitutions.

Full detail in CHANGELOG.md.

v0.3.0 — gap-block alignment, calibrated cutoffs, batch scoring, island profiles

Gap-block alignment, calibrated cutoffs, batch scoring, and per-island profiles — plus a corrected background control that changes every E-value.

Added

• gapblock — single contiguous-indel alignment for V(D)J junctions, the position set by a gap prior (central_prior, profile_prior, frame_prior, positions_prior) rather than the score alone. Exactly optimal against unrestricted affine alignment on 98.8% of genuinely related pairs at a calibrated gap_open.
• gapblock.score_matrix / ScoreMatrix — every query × every reference in one GIL-released C++ call: ~532 M pairs/s on 16 cores, numpy.asarray is zero-copy, and seqtree keeps its zero runtime dependencies. The shape a prototype-distance embedding needs, where nothing can be pruned.
• gapblock.IslandProfile — a per-island position weight matrix whose column penalty is measured against the column consensus, so it stays ≥ 0, is 0 on the consensus, and flows through threshold_for_evalue unchanged. Beats min-over-members only at a strict cutoff (48.5% vs 37.6% at repertoire scale; indistinguishable at a loose one).
• threshold_for_evalue / thetas_from_scores — the per-query score cutoff that achieves a target E-value, inverted exactly from the integer scores.

Breaking

• The bundled control is regenerated as a uniform sample of productive clonotypes (it was the abundance head, 25.8× too self-similar). Every E-value moves. After upgrading, delete ~/.cache/seqtree/control_*.sqtree so the cached index is rebuilt from the corrected control.
• engine="auto" always resolves to seqtm; passing a matrix to seqtrie without an explicit max_penalty now raises. GapPrior takes (block_start, block_length, longer_length).

Full detail in CHANGELOG.md.

v0.2.0 — Miyazawa-Jernigan structural matrix

Highlights

structural is now a Miyazawa-Jernigan interaction-strength matrix. Each residue's interaction strength q(a) = mean_b e(a,b) is read off the MJ contact potential, and sim(a,b) = 10·(1 − |q̂(a) − q̂(b)|). Substitutions between residues of like interaction strength are cheap, so the matrix separates strong hydrophobic interactors (F W C L Y M I V) from weak polar/charged ones (S Q D E K) — the strong/weak-interactor axis of TCR-recognition models (Košmrlj et al., PNAS 2008, doi:10.1073/pnas.0808081105; MJ energies from Miyazawa & Jernigan, J Mol Biol 1996). This lets seqtrie align dissimilar-but-chemically-equivalent loops.

Other

• bench/gen_matrices.py is now self-contained (embeds the MJ contact matrix; drops the texshade.sty/kpsewhich build dependency).
• Docs, header comments, and test anchors updated.

Full diff: v0.1.0...v0.2.0

seqtree 0.1.0

First non-beta release.

• Add SubstitutionMatrix.penalty(a, b) Python API: char-based Gram-distance substitution penalty (0 on identity, larger for dissimilar pairs), the same per-position cost the seqtm/seqtrie engines apply. Enables downstream MI-weighted BLOSUM scoring (mhcmatch).

Builds wheels (cp310-cp313; Linux/macOS/Windows) and publishes to PyPI via CI.

seqtree v0.0.3

seqtree v0.0.3 — pMHC epitope search, position-aware scoring, KmerIndex, MHC-allele guessing, and a reproducible benchmark pipeline.

Highlights since v0.0.2:

• pMHC epitope homology layer (anchor-masked k-mers, mimics, allele assignment) and C++ KmerIndex seed-and-gather.
• Position-aware scoring (PositionalMatrix) and local mode.
• Control-set E-values with Elhanati selection factor; MHC-I/II ROC-PR guessing benchmark + non-binder filter.
• PAM50 + custom (Gram-distance) substitution matrices; seqtm collision metric.
• Reproducible benchmark pipeline: deterministic table producers (shell+python) → plot scripts → committed oracle, with CI oracle-diff and time/memory regression checks.

Wheels (cp310–cp313; Linux x86-64, macOS arm64, Windows x86-64) and the sdist are published to PyPI by the Publish workflow.