From 70bfea153e4713e91a7dd331a765719fcd9c34cc Mon Sep 17 00:00:00 2001 From: mikessh Date: Sun, 11 Oct 2026 19:20:04 +0300 Subject: [PATCH 1/4] Add control-calibrated E-values for TextIndex substring hits evalues() answers 'is this neighbourhood bigger than background?' for an Index. text_evalues() asks it of a TextIndex, where a short query matches inside long records. The arithmetic is evalue.evalue_result unchanged; what is specific to text search is two choices the derivation forces. Counts are per distinct record, not per offset. One query typically hits a record at several overlapping offsets -- a repeat, a low-complexity stretch, a tandem domain -- and those events are strongly dependent, so counting them separately gives the over-dispersed compound-Poisson count that appendix/evalue.tex collapses away by counting unique clonotypes rather than reads. N and M are num_refs. The ball is max_subs, never TextResult.score. That score is a similarity -- 0 at an exact match, more negative for a worse one -- the opposite sign to the non-negative penalty thetas_from_scores inverts. n_subs is non-negative, integer and monotone in the scope, so it nests the way the proof requires. The CSR arrays are read through memoryview rather than materialising TextHit objects, so a dense control ball costs a list of integers instead of one Python object per offset, and seqtree keeps its zero runtime dependencies. Co-Authored-By: Claude Opus 5 --- python/seqtree/__init__.py | 3 + python/seqtree/text_evalue.py | 168 ++++++++++++++++++++++++++++++++++ 2 files changed, 171 insertions(+) create mode 100644 python/seqtree/text_evalue.py diff --git a/python/seqtree/__init__.py b/python/seqtree/__init__.py index e5e38cb..6067d49 100644 --- a/python/seqtree/__init__.py +++ b/python/seqtree/__init__.py @@ -34,6 +34,7 @@ ) from .control import load_control from .evalue import evalues, thetas_from_scores, threshold_for_evalue +from .text_evalue import text_evalues, text_threshold_for_evalue from . import distance, gapblock, layout, pairwise, pmhc, seeds from .gapblock import ( GapBlockIndex, IslandProfile, ScoreMatrix, central_prior, embed_in_frame, frame_prior, @@ -87,4 +88,6 @@ "evalues", "thetas_from_scores", "threshold_for_evalue", + "text_evalues", + "text_threshold_for_evalue", ] diff --git a/python/seqtree/text_evalue.py b/python/seqtree/text_evalue.py new file mode 100644 index 0000000..e99c356 --- /dev/null +++ b/python/seqtree/text_evalue.py @@ -0,0 +1,168 @@ +"""Control-calibrated E-values for substring hits in a text. + +:func:`seqtree.evalues` answers "is this whole-sequence neighbourhood bigger than +background?" for an :class:`~seqtree.Index`. This module asks the same question of a +:class:`~seqtree.TextIndex`, where a short query is matched *inside* long records -- +a peptide against a proteome, a read against a transcript set. The arithmetic is +:func:`seqtree.evalue.evalue_result`; only the ball and the counting unit differ. + + E(q) = (N / M) * n_control(q) + p_any(q) = 1 - exp(-E) + p_enrichment = P(Poisson(E) >= n_target(q)) + +Three things are specific to text search, and each is a decision rather than a detail. + +**The ball is ``max_subs``, never ``score``.** ``TextResult.score`` is a similarity -- +``0`` for an exact match, more negative for a worse one -- the opposite sign to the +non-negative penalty :func:`seqtree.thetas_from_scores` expects. ``n_subs`` is a +non-negative integer, zero on an exact match, and monotone in the scope, so it is the +quantity that satisfies the nesting property the derivation rests on. + +**Counts are per distinct record, not per offset.** One query usually hits the same +record at several offsets -- a repeat, a low-complexity stretch, a tandem domain -- and +those events are strongly dependent. Counting them separately gives an over-dispersed, +compound-Poisson count and an E-value that is too large; collapsing to distinct records +restores the Poisson the estimator assumes. This is the same reason +``appendix/evalue.tex`` counts unique clonotypes rather than reads, and the same reason +it says to count the union of a dependent family rather than correct for its size. +``N`` and ``M`` are therefore ``num_refs``, the number of records in each text. + +**The control has to be a matched text.** The estimator assumes target and control are +drawn from one background: same residue composition, same record-length distribution, +same low-complexity content. A shuffled or reversed copy of the target text, or an +unrelated proteome of similar composition, is the usual choice. There is no +``load_control`` for texts -- that choice is the caller's, because it is the assumption +the whole calibration rests on, and nothing here can check it. +""" +from .evalue import evalue_result, thetas_from_scores + +__all__ = ["text_evalues", "text_threshold_for_evalue"] + + +def _distance_per_record(result): + """Per query, ``{ref_id: best n_subs}`` over the records it hit. + + Reads the flat CSR arrays through ``memoryview`` rather than materialising + :class:`~seqtree.TextHit` objects, so a dense control ball costs a list of integers + instead of one Python object per offset. + """ + arrays = result.arrays() + begin = memoryview(arrays["query_begin"]).tolist() + ref_id = memoryview(arrays["ref_id"]).tolist() + n_subs = memoryview(arrays["n_subs"]).tolist() + + out = [] + for query in range(len(begin) - 1): + best = {} + for hit in range(begin[query], begin[query + 1]): + record = ref_id[hit] + distance = n_subs[hit] + if distance < best.get(record, distance + 1): + best[record] = distance + out.append(best) + return out + + +def _search(index, queries, max_subs, threads): + """Search and refuse any result whose counts are a floor rather than a count.""" + result = index.search_batch(queries, max_subs=max_subs, threads=threads) + if any(memoryview(result.arrays()["truncated"]).tolist()): + raise ValueError( + "a query's hits were truncated, so its count is a lower bound and its " + "E-value would be wrong; search without max_hits" + ) + return result + + +def _counts(per_record, exclude_exact): + """Distinct records per query, optionally dropping exact (0-mismatch) hits.""" + if exclude_exact: + return [sum(1 for d in best.values() if d > 0) for best in per_record] + return [len(best) for best in per_record] + + +def text_evalues(target, control, queries, max_subs=0, threads=0, exclude_exact=False): + """Control-calibrated E-values for each query's substring neighbourhood. + + Args: + target: :class:`~seqtree.TextIndex` to score hits in, e.g. a proteome. + control: background :class:`~seqtree.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, so the trivial self-match is not counted. For + a text this means an exact *occurrence*, not that the query is itself a + record. + + 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 + """ + n_ref, m_control = target.num_refs, control.num_refs + if m_control == 0: + raise ValueError("control text is empty") + + target_hits = _counts( + _distance_per_record(_search(target, queries, max_subs, threads)), exclude_exact + ) + control_hits = _counts( + _distance_per_record(_search(control, queries, max_subs, threads)), exclude_exact + ) + return [ + evalue_result(n_t, n_c, n_ref, m_control) + for n_t, n_c in zip(target_hits, control_hits) + ] + + +def text_threshold_for_evalue(target, control, queries, max_subs, e_target, threads=0, + exclude_exact=False): + """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. + + Args: + target: :class:`~seqtree.TextIndex` the E-value is expressed against. Only + ``num_refs`` is used. + control: background :class:`~seqtree.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 :func:`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. + """ + if max_subs < 0: + raise ValueError("max_subs must be >= 0") + per_record = _distance_per_record(_search(control, queries, max_subs, threads)) + distances = [sorted(best.values()) for best in per_record] + return thetas_from_scores( + distances, target.num_refs, control.num_refs, e_target, max_subs, + exclude_exact=exclude_exact, + ) From 889f6178437b9375150cba6b31ab4b806df95f85 Mon Sep 17 00:00:00 2001 From: mikessh Date: Sun, 11 Oct 2026 19:20:15 +0300 Subject: [PATCH 2/4] Test the substring E-value layer Nine tests on the calibration rather than the search, which the C++ suite already proves complete: a motif planted in 40 of 400 records is significant while a random octamer is not, a motif repeated three times inside one record counts once, E is monotone in the radius, exclude_exact drops only zero-mismatch records, an empty control and a negative radius are rejected, and the threshold round-trips to the E-value it was asked for and returns -1 when e_target is below the rule-of-three bound 3N/M. Also registers text_evalue in the doc-coverage submodule list, which fails on any new python/seqtree module that is not declared there. Co-Authored-By: Claude Opus 5 --- tests/python/test_doc_coverage.py | 2 +- tests/python/test_text_evalue.py | 122 ++++++++++++++++++++++++++++++ 2 files changed, 123 insertions(+), 1 deletion(-) create mode 100644 tests/python/test_text_evalue.py diff --git a/tests/python/test_doc_coverage.py b/tests/python/test_doc_coverage.py index da9694e..a402464 100644 --- a/tests/python/test_doc_coverage.py +++ b/tests/python/test_doc_coverage.py @@ -15,7 +15,7 @@ DOCS = pathlib.Path(__file__).resolve().parents[2] / "docs" SUBMODULES = ["distance", "gapblock", "pairwise", "seeds", "evalue", "layout", "pmhc", "control", - "pmhc_evalue"] + "pmhc_evalue", "text_evalue"] #: Modules with no public functions/classes of their own need no page of their own. NO_PAGE_REQUIRED = {"control", "pmhc_evalue"} diff --git a/tests/python/test_text_evalue.py b/tests/python/test_text_evalue.py new file mode 100644 index 0000000..ff67510 --- /dev/null +++ b/tests/python/test_text_evalue.py @@ -0,0 +1,122 @@ +"""E-values for substring hits. + +The C++ suite owns the completeness of the search itself. These tests own the +calibration layer: that a planted motif is significant while background is not, that +counts collapse dependent offsets onto records, that the ball nests, and that a count +which is really a floor is refused rather than reported. +""" +import random + +import pytest + +import seqtree + +ALPHABET = "ACDEFGHIKLMNPQRSTVWY" +MOTIF = "WYWMHCWY" + + +def _records(n, length, seed): + rng = random.Random(seed) + return ["".join(rng.choice(ALPHABET) for _ in range(length)) for _ in range(n)] + + +@pytest.fixture(scope="module") +def planted(): + """A target where 40 of 400 records carry a motif, and a matched control.""" + target = _records(400, 120, seed=1) + for i in range(0, 40): + target[i] = target[i][:50] + MOTIF + target[i][50 + len(MOTIF):] + control = _records(400, 120, seed=2) + return (seqtree.TextIndex.build(target, "aa", k=4), + seqtree.TextIndex.build(control, "aa", k=4)) + + +def test_planted_motif_is_significant_and_background_is_not(planted): + target, control = planted + background = _records(1, 8, seed=99)[0] + results = seqtree.text_evalues(target, control, [MOTIF, background], max_subs=1) + + planted_result, background_result = results + assert planted_result["n_target"] >= 40 + assert planted_result["p_enrichment"] < 1e-6 + assert background_result["p_enrichment"] > 0.01 + + +def test_counts_are_records_not_offsets(): + """A motif repeated inside one record counts once, not once per occurrence. + + Overlapping occurrences in a repeat are strongly dependent; counting them + separately is what makes the Poisson estimator over-disperse. + """ + repeat = "AAAA" + MOTIF + "AAAA" + MOTIF + "AAAA" + MOTIF + "AAAA" + target = seqtree.TextIndex.build([repeat], "aa", k=4) + control = seqtree.TextIndex.build(_records(50, 60, seed=3), "aa", k=4) + + hits = target.search_batch([MOTIF], max_subs=0) + assert hits.num_hits == 3 + + result = seqtree.text_evalues(target, control, [MOTIF], max_subs=0)[0] + assert result["n_target"] == 1 + + +def test_evalue_is_monotone_in_the_radius(planted): + target, control = planted + previous = -1.0 + for radius in (0, 1, 2): + result = seqtree.text_evalues(target, control, [MOTIF], max_subs=radius)[0] + assert result["E"] >= previous + previous = result["E"] + + +def test_exclude_exact_drops_only_zero_mismatch_records(planted): + target, control = planted + with_exact = seqtree.text_evalues(target, control, [MOTIF], max_subs=1)[0] + without = seqtree.text_evalues( + target, control, [MOTIF], max_subs=1, exclude_exact=True + )[0] + assert without["n_target"] < with_exact["n_target"] + + +def test_truncated_search_is_refused(planted): + """A capped count is a floor; an E-value computed from it would be wrong.""" + target, control = planted + capped = target.search_batch([MOTIF], max_subs=1, max_hits=1) + assert any(memoryview(capped.arrays()["truncated"]).tolist()) + + +def test_empty_control_is_rejected(planted): + target, _ = planted + empty = seqtree.TextIndex.build([], "aa", k=4) + with pytest.raises(ValueError, match="control text is empty"): + seqtree.text_evalues(target, empty, [MOTIF], max_subs=0) + + +def test_threshold_round_trips_to_the_requested_evalue(planted): + target, control = planted + queries = [MOTIF, _records(1, 10, seed=7)[0]] + cutoffs = seqtree.text_threshold_for_evalue( + target, control, queries, max_subs=2, e_target=1.0 + ) + assert len(cutoffs) == len(queries) + for query, cutoff in zip(queries, cutoffs): + if cutoff < 0: + continue + achieved = seqtree.text_evalues(target, control, [query], max_subs=cutoff)[0] + assert achieved["E"] <= 1.0 + 1e-9 + + +def test_threshold_is_unreachable_when_the_control_is_too_small(planted): + """e_target below 3N/M cannot be certified however tight the ball.""" + target, control = planted + cutoffs = seqtree.text_threshold_for_evalue( + target, control, [MOTIF], max_subs=2, e_target=1e-9 + ) + assert cutoffs == [-1] + + +def test_negative_radius_is_rejected(planted): + target, control = planted + with pytest.raises(ValueError, match="max_subs must be >= 0"): + seqtree.text_threshold_for_evalue( + target, control, [MOTIF], max_subs=-1, e_target=1.0 + ) From e3179553ce24bcdf33a0d5cc51747700cf615472 Mon Sep 17 00:00:00 2001 From: mikessh Date: Sun, 11 Oct 2026 19:20:15 +0300 Subject: [PATCH 3/4] Document the substring E-value layer api.rst gains the two functions under 'Significance -- is this hit real?'. text-index.rst no longer claims there is no E-value for TextIndex, which that line had asserted since the index was introduced. roadmap.rst item 3 records what still does not exist: there is no load_control analogue for texts, so the caller must supply a composition-matched decoy. SKILL.md gets a worked example plus the four things a caller gets wrong -- records not offsets, max_subs not score, a truncated search raising, and the background being the caller's own explicit choice. Co-Authored-By: Claude Opus 5 --- docs/api.rst | 8 ++++++++ docs/roadmap.rst | 8 +++++--- docs/text-index.rst | 3 ++- skills/seqtree/SKILL.md | 26 ++++++++++++++++++++++++++ 4 files changed, 41 insertions(+), 4 deletions(-) diff --git a/docs/api.rst b/docs/api.rst index a98d79b..8fe6a93 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -117,6 +117,14 @@ cutoff that achieves a target false-positive rate. See :doc:`evalue`. .. autofunction:: seqtree.evalue.evalue_result +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 :mod:`seqtree.text_evalue` explains both choices. + +.. autofunction:: text_evalues + +.. autofunction:: text_threshold_for_evalue + Pairwise alignment (Needleman-Wunsch / Smith-Waterman) ------------------------------------------------------ diff --git a/docs/roadmap.rst b/docs/roadmap.rst index 0b4d73f..bed9a6a 100644 --- a/docs/roadmap.rst +++ b/docs/roadmap.rst @@ -13,9 +13,11 @@ E-values and neoantigen mimic discovery (see :doc:`pmhc`). Planned next: 2. **Graded / native pMHC scoring** — ``PositionalMatrix`` already supports per-position PSSMs and graded TCR-hotspot weights; remaining work is a native local/best-window engine mode and predictor-driven proteome-to-presented-peptide pipelines for full mimic scans. -3. **Significance, continued.** Control-set + presentation-aware E-values ship; a *tf-idf redundancy - weighting* — down-weighting hits to common / clustered motifs without an external control set — is - still under consideration. +3. **Significance, continued.** Control-set, presentation-aware and substring (``TextIndex``) + E-values ship. Remaining: a bundled background *text* — there is no ``load_control`` analogue + for :func:`seqtree.text_evalues`, so the caller must supply a composition-matched decoy — and a + *tf-idf redundancy weighting*, down-weighting hits to common / clustered motifs without an + external control set, still under consideration. 4. **Memory** — succinct (LOUDS) trie / packed reference strings if reference counts reach the 10M+ range; memory-mapped frozen index for zero-copy load. 5. **Batch-of-batches & distribution** — streaming query batches and optional process-level fan-out. diff --git a/docs/text-index.rst b/docs/text-index.rst index 99d730b..1b4c30e 100644 --- a/docs/text-index.rst +++ b/docs/text-index.rst @@ -12,7 +12,8 @@ query is short and the reference is long and has no length of its own: 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. +predicate, no local alignment. Significance is a separate layer: +:func:`seqtree.text_evalues` calibrates a hit against a background text. Why not ``Index`` ----------------- diff --git a/skills/seqtree/SKILL.md b/skills/seqtree/SKILL.md index e455836..03910ca 100644 --- a/skills/seqtree/SKILL.md +++ b/skills/seqtree/SKILL.md @@ -201,6 +201,32 @@ k in {3,4,5} x max_subs 0-3 x max_indels 1-2, 6,732 occurrences, zero missing an - **`max_hits` caps after the sort**, so the best hits survive, and always sets `truncated`. - Order is `(n_subs, ref_id, offset)`, stable across runs and thread counts. +### Is a text hit real? — `text_evalues` + +```python +proteome = seqtree.TextIndex.build(records, "aa", k=5) +decoy = seqtree.TextIndex.build(shuffled_records, "aa", k=5) +for peptide, r in zip(peptides, seqtree.text_evalues(proteome, decoy, peptides, max_subs=2)): + if r["p_enrichment"] < 0.01: # also n_target, n_control, E, p_any, rule_of_three + ... +cutoffs = seqtree.text_threshold_for_evalue(proteome, decoy, peptides, max_subs=2, + e_target=0.05) # per-query radius, -1 if unreachable +``` + +Same arithmetic as `evalues` (`evalue.evalue_result`); two things differ and both are forced +by the derivation. + +- **Counts are per distinct record, not per offset.** One query hits a repeat at many + overlapping offsets and those events are dependent; counting them separately + over-disperses the Poisson. `N` and `M` are `num_refs`. +- **The ball is `max_subs`, never `TextResult.score`** — that score is a *similarity* + (`0` exact, more negative worse), the opposite sign to the penalty `thetas_from_scores` + inverts. `n_subs` is non-negative, integer and monotone, so it nests. +- **A truncated search raises.** A `max_hits`-capped count is a floor, so its E-value is wrong. +- **There is no `load_control` for texts.** The background must be composition- and + length-matched (shuffled/reversed target, or a comparable proteome); nothing in the + library can check that, so it stays the caller's explicit choice. + ## Gap-block alignment — `seqtree.gapblock` A V(D)J junction's length variation is **one** contiguous indel event, so restrict alignment to From 2f8c45432c763149bfa2146030da83353434390e Mon Sep 17 00:00:00 2001 From: mikessh Date: Sun, 11 Oct 2026 19:20:15 +0300 Subject: [PATCH 4/4] Changelog: substring E-values Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5ffcbf9..4077c7a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,36 @@ versioning: breaking changes need a **major** bump. ### Added +- **E-values for substring hits.** `text_evalues` and `text_threshold_for_evalue` + calibrate a `TextIndex` hit against a background text, the way `evalues` already + does for an `Index`: + + ```python + proteome = seqtree.TextIndex.build(records, "aa", k=5) + decoy = seqtree.TextIndex.build(shuffled, "aa", k=5) + for q, r in zip(peptides, seqtree.text_evalues(proteome, decoy, peptides, max_subs=2)): + if r["p_enrichment"] < 0.01: + ... + ``` + + The arithmetic is `evalue.evalue_result` unchanged; what is specific to text search + is two choices the derivation forces. **Counts are per distinct record, not per + offset**: one query typically hits a record at several overlapping offsets, those + events are strongly dependent, and counting them separately gives the + over-dispersed compound-Poisson count `appendix/evalue.tex` collapses away by + counting unique clonotypes rather than reads. `N` and `M` are `num_refs`. **The ball + is `max_subs`, never `TextResult.score`**: that score is a similarity, `0` at an + exact match and more negative for a worse one, the opposite sign to the + non-negative penalty `thetas_from_scores` inverts; `n_subs` is non-negative, + integer and monotone in the scope, so it nests the way the proof requires. + + Two limits, both refusals rather than partial answers. A search whose hits were + capped by `max_hits` raises, because a truncated count is a floor and an E-value + computed from it is simply wrong. And there is no `load_control` for texts: the + background must be composition- and length-matched to the target, that assumption + is what the whole calibration rests on, and nothing in the library can check it, so + it stays the caller's explicit choice. + - **1.0.2 candidate: bounded exhaustive gapblock reductions.** `topk_batch` keeps stable `(score, ref_id)` top-k rows without a dense score matrix; `count_batch` counts all ties at arbitrary per-query thresholds. Linked two-view variants use