E-values for TextIndex substring hits

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

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

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.

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

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",
 ]

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,
+    )

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

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"}