Text search#

Index answers a question about a set of reference strings. TextIndex answers a question about a text — a proteome, a genome, a set of transcripts — where the query is short and the reference is long and has no length of its own:

Given a reference text and a large batch of short queries of many different lengths, find every position where a query matches within m substitutions. No gaps, full length, exact answer.

This is peptide-to-proteome mapping. It is not homology search: no scoring matrix in the predicate, no E-value, no local alignment.

Why not Index#

To ask a text question with Index you have to enumerate every length-L window of the text as a separate reference string — and do it again for every distinct query length. The human proteome has 68,389,335 nine-mer windows, so one such index costs several gigabytes, and a query set spanning 45 distinct lengths costs 45 of them.

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

from seqtree import TextIndex

ix = TextIndex.build(records, alphabet="aa", k=4)
res = ix.search_batch(["SIINFEKL", "KTAYIAKQRQISFVKSHFSRQ"], max_subs=2, threads=0)

for query, hits in zip(queries, res):
    for h in hits:
        h.ref_id, h.offset, h.n_subs, h.mismatches

How it answers exactly#

One rule, applied per query. Split the query into b disjoint blocks of width bw = L / b >= k, and probe block j’s leading k-mer at radius c_j — every k-mer within c_j substitutions of it.

Lossless if and only if sum(c_j) >= max_subs - b + 1.

The scheme can only miss an occurrence whose errors satisfy e_j > c_j for every block. The lightest such error vector is e_j = c_j + 1, of total weight sum(c_j) + b, so no occurrence within max_subs can hide once sum(c_j) + b > max_subs. Errors falling outside the blocks only lower the total, so an uncovered tail is free.

Probing one block at radius c costs N(c) = sum(C(k,i)·(A-1)^i, i <= c), whose increments climb steeply in c and are the same for every block — so the cheapest legal scheme takes as many blocks as fit and spreads the budget as thinly as possible:

b = min(max_subs + 1, L // k)
r = max(0, max_subs - b + 1)          # budget units to place
c_j = r // b, plus one for the last (r % b) blocks

At b = max_subs + 1 every c_j is 0 and this is plain pigeonhole — one exact seed per error, at least one of which must survive. At b = 1 it is a single substitution ball over q[0:k]. The useful schemes are in between, and they are where the short queries live: at L = 8, max_subs = 2, k = 4 two disjoint 4-mers fit, so c = (0, 1) probes 94 variants where a one-block ball probes 3,267.

Two details are not incidental. The spare units go on the last blocks, worth a further 1.4×: a candidate already matches its own block’s k-mer, so left-to-right verification meets the unconstrained prefix first and rejects after a residue or two — and the high-budget block is the one contributing nearly all the candidates. And a start that two blocks both reach is returned once, emitted by the lowest-indexed block that could have produced it — a test read straight off the mismatch positions verification has already computed, so deduplication needs no sort and no hash set.

None of this is a heuristic. Completeness is checked against a brute-force scan over a grid of query length, max_subs and k, asserting set equality of (ref_id, offset, n_subs) and of the mismatch detail — not merely that hits were found — and separately that the answer is identical across k, so a dispatch bug cannot pass by agreeing with itself.

Choosing k#

k never changes the answer, only the work: a bigger k makes one probe more specific (each 4-mer bucket on the human proteome holds ~208 positions, each 5-mer bucket ~9) but needs L >= b·k to fit the same number of blocks. Measured on the human proteome (UP000005640, 147,506 records / 69,578,135 residues), one thread, Apple M-series, bench/bench_text_index.py:

L

m

ms/query, k = 4

ms/query, k = 5

8

2

1.468

4.694

9

2

1.334

4.526

10

2

1.298

0.129

11

2

1.138

0.108

12

2

0.097

0.126

15

2

0.086

0.014

10

3

3.595

0.347

12

3

1.567

0.394

15

3

1.514

0.138

Hit counts are identical in every cell — the benchmark asserts it rather than reporting it, so a dispatch that got b or a budget wrong would fail there on real text and not only in a unit test.

The crossover is L = 2k: below it only one block fits, so the query pays a full radius-m ball over its leading k-mer and the bigger k is much worse. Pick the largest k with 2k <= L for the bulk of the query set — k = 4 down to length 8, k = 5 from 10 up. At k = 5 every L >= 10 cell above is at or under 0.4 ms/query. Building is cheap either way (0.55 s at k = 4, 0.59 s at k = 5), so a corpus spanning both can simply hold both indexes — still one build per k, never one per query length.

Threads scale well — human proteome, L = 12, max_subs = 2, 50,000 queries: 0.103 ms/query on one thread, 0.013 on eight (7.8x), 0.008 on all cores (12.6x). The ratio is capped by the serial CSR flatten at the end of a batch rather than by the search, which is why it falls short of linear on a workload the search scheme made this fast.

Results are flat arrays#

A 445,000-query run returns a handful of arrays, not millions of Python objects. TextResult holds parallel arrays in CSR form and materialises hit objects only for the query you ask for:

res.num_hits                 # total across the batch
res[i]                       # TextHit list for query i, built on demand
res.arrays()                 # every array as a zero-copy ArrayView (buffer protocol)
res.to_numpy()               # the same, as numpy views; numpy imported on the call

Hits for query i are [query_begin[i], query_begin[i+1]); mismatches for hit h are [mm_begin[h], mm_begin[h+1]). Order is (n_subs, ref_id, offset) — stable across runs and across thread counts, so a downstream digest is stable.

Beyond a position#

Mismatch detail is a pair, not a position. Each mismatch is (pos, query_aa, text_aa), so a caller ranking by chemistry can tell L→I from L→D without re-fetching the window.

Optional scoring. Pass matrix= and each hit also carries a score: the substitution matrix’s similarity summed over the mismatched positions, so 0 is an exact match and less-negative is more conservative. It only scores hits the Hamming predicate has already accepted — it never changes which hits come back.

Nearest shell, whole. best_only=True walks the distance upward and stops at the first shell that has anything in it, returning all of it. A query resolved at m = 0 never enters the m = 1 pass, and no second index is needed.

Group folding. build(..., group_ids=[...]) labels each record with an arbitrary integer — gene ids, species, cluster labels; seqtree does not know what they mean. Then group_by=True collapses hits onto them:

res.groups(i)   # [(group_id, min_subs, n_hits), ...], sorted by group id

More than one row means the nearest parents disagree, which makes a tie a first-class output rather than something every caller re-derives from the raw hits.

Truncation is reported. A degenerate query can match a repeat region a very large number of times, so max_hits caps it — after the sort, so the best hits survive — and sets res.truncated[i]. A cap that is invisible is a recall bug wearing a performance costume.

The alphabet, which is where a silent wrong answer would live#

Symbols are compared through the alphabet’s codec, not as raw bytes. For alphabet="aa" that codec is the 24 symbols ARNDCQEGHILKMFPSTWYVBZX*, so:

  • B, Z, X and * are real symbols, not wildcards — X matches X and is one substitution from A. The human proteome carries 8,417 X residues and the mouse proteome 4,263, so this is not hypothetical.

  • U (selenocysteine) is not in the codec. It encodes to the invalid marker, so no hit can cross a U in the text, and a query containing one is refused by name:

    ValueError: queries[42] ('CASSUGQYF'): symbol 'U' is not in the alphabet

    The same marker separates records, which is what structurally prevents a hit — or a seed — from spanning a record boundary, rather than a check that could be forgotten.

Every query must be at least k residues long. Seeds are indexed only where a whole k-mer fits inside a record, so a shorter query could match within k - L of a record end and be missed; TextIndex raises rather than answering incompletely.

Persistence#

The file is flat, with the arrays at known aligned offsets, so loading it is a map rather than a parse:

ix.save("human_k4.sti")
ix = TextIndex.load("human_k4.sti")            # mmap=True by default

With mmap=True the pages are shared read-only across processes, so a fan-out of workers costs one copy of the index rather than one per worker — and, since there is nothing to parse, the load does not scale with the file.

text

k

residues

file

× text

save

load

mouse

4

23,131,234

116.6 MB

5.0

0.14 s

0.10 ms (mmap)

human

4

69,578,135

348.2 MB

5.0

0.46 s

0.10 ms (mmap)

human

5

69,578,135

378.1 MB

5.4

0.54 s

0.10 ms (mmap)

Reading the file instead (mmap=False) costs 0.12 s on the human index. Roughly 5× the text is what a positional seed index costs: 79 % of the file is post_ids, one uint32 per in-record k-mer start, and every one of those 69.1 M positions has to be addressable for the answer to be exact. Small texts pay a different bill — the bucket table is directly addressed, so it has a floor of A^k × 4 bytes (1.3 MB at k = 4, 31.8 MB at k = 5) regardless of how little text there is.

A .sti is treated as untrusted input: the header’s record, bucket and posting counts are cross-checked against each other and against the file length before any of them is used to size an allocation or to index into the mapping.

Insertions and deletions#

max_indels allows gaps. The two caps are independent, so max_subs=2, max_indels=1 accepts two substitutions and one gap — not three edits of any kind:

res = ix.search_batch(peptides, max_subs=1, max_indels=1, threads=0)
for h in res[0]:
    print(h.ref_id, h.offset, h.length, h.n_subs, h.n_ins, h.n_dels)
    print(ix.ref_seq(h.ref_id)[h.offset : h.offset + h.length])

The seed table does not change, and neither does the scheme. With b = max_subs + max_indels + 1 disjoint blocks the total edit count cannot reach the block count, so the same pigeonhole argument still puts a zero-error block somewhere — and a zero-error block matches the text exactly, gaps or not, so the existing exact lookup finds it. Nothing about the index, the build, or the probe enumeration is different.

What changes is verification. A seed no longer pins the start: the prefix before the block may have gained or lost up to max_indels residues, so every start in [pos - base - max_indels, pos - base + max_indels] is aligned, by a banded DP over the states (query position, net offset, indels spent).

Two consequences worth knowing before you call it:

  • A match is no longer len(query) residues wide. Read length; it is len(query) + n_dels - n_ins.

  • Every seed block must be exact, so a query must be at least (max_subs + max_indels + 1) * k long. At k = 4, max_subs=1, max_indels=1 needs 12 residues. A shorter query raises, naming the bound — it is not answered partially. Build a second index at a smaller k if your queries are shorter than that.

An alignment must begin and end on an aligned pair. A gap at either edge would stretch the reported interval over a residue that matches nothing, and a leading gap is just the same occurrence starting one residue over — allowing it would report one match two or three times.

mismatches is empty on this path: recovering which columns were substituted needs an alignment traceback the verifier does not keep. The counts are exact; the per-column detail is only available at max_indels = 0. matrix= is likewise ignored.

Cost, measured on a synthetic 5,000,000-residue corpus, k = 4, one thread, against the same query set with max_indels = 0:

query length

max_subs

ms/query, no indels

ms/query, max_indels=1

16

1

0.001

0.028

16

2

0.001

0.043

20

3

0.002

0.065

24

3

0.002

0.067

Gapped search costs 25–37× an ungapped one at the same length and max_subs: the banded DP replaces a byte comparison that early-exits after a residue or two, and it runs once per start in the window. It is still well under 0.1 ms/query. Completeness is pinned the same way the substitution path is — set equality against an independent brute force, over k ∈ {3, 4, 5} × max_subs 0–3 × max_indels 1–2, 6,732 occurrences, zero missing and zero extra.

Limits#

TextIndex answers a Hamming question over a text. Gaps, scored alignment and edit distance stay with Index and seqtree.gapblock — indels change the length and break the pigeonhole argument the exactness rests on. Postings are 32-bit, capping the text at 2\ 32 symbols; k is capped so the direct-addressed table stays under 2\ 27 buckets (k <= 5 for amino acids).