feat: support max_indels=2 for homogeneous insertions or deletions

.github/workflows/tests.yml

@@ -30,6 +30,8 @@ jobs:
             test-file: pepmatch/tests/test_indel_search.py
           - name: Indel Property
             test-file: pepmatch/tests/test_indel_property.py
+          - name: 2-Indel Property
+            test-file: pepmatch/tests/test_indel2_property.py
 
     steps:
     - name: Checkout Code

README.md

@@ -8,13 +8,13 @@
 
 **Author:** Daniel Marrama
 
-`PEPMatch` is a high-performance peptide search tool for finding short peptide sequences within a reference proteome. It handles three kinds of search out of the box: **exact matches**, **mismatch (substitution)** searches, and **single-indel (insertion or deletion)** matching. Powered by a Rust engine with Python bindings, it delivers sub-second search times across entire proteomes while maintaining a simple Python API.
+`PEPMatch` is a high-performance peptide search tool for finding short peptide sequences within a reference proteome. It handles three kinds of search out of the box: **exact matches**, **mismatch (substitution)** searches, and **indel (insertion or deletion)** matching. Powered by a Rust engine with Python bindings, it delivers sub-second search times across entire proteomes while maintaining a simple Python API.
 
 ### Key Features
 
 * **Blazing Fast**: Rust-powered search engine with automatic multi-core parallelization via Rayon. Search thousands of peptides against the entire human proteome in seconds.
 * **Unified Index Format**: Single `.pepidx` binary format stores sequences, metadata, and k-mer index in one memory-mapped file. Preprocess once, search repeatedly.
-* **Versatile Searching**: Exact matches, mismatch (substitution) searches, single-indel (insertion/deletion) matching, best match mode, and discontinuous epitope support.
+* **Versatile Searching**: Exact matches, mismatch (substitution) searches, indel (insertion/deletion) matching, best match mode, and discontinuous epitope support.
 * **Counts-Only Mode**: Get aggregate hit counts per peptide with `O(unique queries)` memory instead of `O(hits)` — built for massive query sets and dense reference matching.
 * **Simple API**: Two classes — `Preprocessor` and `Matcher` — handle everything.
 * **Flexible I/O**: Accepts queries from FASTA files, text files, or Python lists. Outputs to CSV, TSV, XLSX, JSON, or Polars DataFrame.
@@ -103,9 +103,9 @@ df = Matcher(
   max_indels=1
 ).match()
 ```
-Each indel hit is annotated in an **`Indel Positions`** column with the edit type, residue, and 1-based query position — e.g. `d: A[6]` (deletion of `A` at query position 6) or `i: X[6]` (insertion of `X` before query position 6); an exact match reports `[]`. In a repeat the exact position is ambiguous, so the inclusive range of all valid positions is reported, e.g. `d: A[2,4]`.
+Each indel hit is annotated in an **`Indel Positions`** column with the edit type, residue, and 1-based query position — e.g. `d: A[6]` (deletion of `A` at query position 6) or `i: X[6]` (insertion of `X` before query position 6); an exact match reports `[]`. In a repeat the exact position is ambiguous, so the inclusive range of all valid positions is reported, e.g. `d: A[2,4]`. When two edits fall at one site they are reported as a single chunk (`d: YA[2]`); when they fall apart they are reported separately (`d: C[3], d: E[5]`).
 
-Currently limited to `max_indels=1`; mutually exclusive with `max_mismatches`.
+Supports `max_indels=1` and `max_indels=2`; mutually exclusive with `max_mismatches`. At `max_indels=2` the search is **homogeneous** — a match may use two insertions or two deletions, but never one of each (a mixed pair is a substitution, which mismatch search already covers). The two edits need not be adjacent.
 
 #### Best Match
 
@@ -177,7 +177,7 @@ pepmatch-match -q peptides.fasta -p human.fasta -m 0 -k 5
 * `-q`, `--query` (Required): Path to the query file.
 * `-p`, `--proteome_file` (Required): Path to the proteome FASTA file.
 * `-m`, `--max_mismatches`: Maximum mismatches allowed (default: 0).
-* `-i`, `--max_indels`: Maximum indels allowed (default: 0). Currently limited to 1, and mutually exclusive with `-m`.
+* `-i`, `--max_indels`: Maximum indels allowed (default: 0). Supports 1 or 2, and mutually exclusive with `-m`. At 2, a match uses two insertions or two deletions, never one of each.
 * `-k`, `--kmer_size`: K-mer size (default: 5).
 * `-P`, `--preprocessed_files_path`: Directory containing preprocessed files.
 * `-b`, `--best_match`: Enable best match mode.

pepmatch/matcher.py

@@ -1,5 +1,6 @@
 import os
 import polars as pl
+from itertools import combinations
 from pathlib import Path
 from Bio import SeqIO
 from ._rs import rs_preprocess, rs_match, rs_discontinuous, rs_metadata, rs_match_counts, rs_indel_match
@@ -23,41 +24,88 @@ def output_matches(df: pl.DataFrame, output_format: str, output_name: str) -> No
     df.write_json(path)
 
 
-def _indel_edit(query, matched):
-  """Return (kind, residues, low, high) for a single-indel match, or None for an
-  exact match: kind 'd'/'i', the deleted query or inserted protein residue(s), and
-  [low, high] the inclusive 1-based range of valid positions (a run of equivalent
-  positions in a repeat)."""
-  if matched == query:
-    return None
-  L = len(query)
-  if len(matched) < L:
-    # interior positions only — query-terminal deletions are barred
-    positions = [i for i in range(2, L) if query[:i - 1] + query[i:] == matched]
-    if not positions:
-      return None
-    return ('d', query[positions[0] - 1], positions[0], positions[-1])
-  # interior matched residues only — boundary insertions are barred
-  positions, residue = [], None
-  for k in range(1, len(matched) - 1):
-    if matched[:k] + matched[k + 1:] == query:
-      positions.append(k + 1)   # 1-based query position the inserted residue precedes
-      residue = matched[k]
-  if not positions:
-    return None
-  return ('i', residue, positions[0], positions[-1])
+def _indel_placements(query, matched):
+  """Every valid way the edits could be placed, as a tuple of (1-based query
+  position, residue) per edit. A repeat makes several placements equivalent."""
+  L, M = len(query), len(matched)
+  n = abs(M - L)
+  if n == 0:
+    return []
+  out = []
+  if M < L:
+    # deletions: interior query positions only — query-terminal deletions are barred
+    for combo in combinations(range(1, L - 1), n):
+      if ''.join(query[i] for i in range(L) if i not in combo) == matched:
+        out.append(tuple((i + 1, query[i]) for i in combo))
+  else:
+    # insertions: the residue comes from the protein; its query position is the
+    # number of query residues preceding it. Boundary insertions are barred.
+    for combo in combinations(range(M), n):
+      if ''.join(matched[i] for i in range(M) if i not in combo) != query:
+        continue
+      placement = []
+      for k in combo:
+        p = sum(1 for x in range(k) if x not in combo) + 1
+        if p == 1 or p == L + 1:
+          placement = None
+          break
+        placement.append((p, matched[k]))
+      if placement:
+        out.append(tuple(placement))
+  return out
+
+
+def _indel_edits(query, matched):
+  """Entries for the Indel Positions column as (kind, residues, low, high)."""
+  placements = _indel_placements(query, matched)
+  if not placements:
+    return []
+  kind = 'd' if len(matched) < len(query) else 'i'
+
+  if len(placements[0]) == 1:
+    positions = sorted(p for ((p, _),) in placements)
+    residue = placements[0][0][1]
+    return [(kind, residue, positions[0], positions[-1])]
+
+  # Two edits at one site form a chunk: adjacent query positions for a deletion,
+  # the same query gap for an insertion.
+  chunks = sorted({
+    (p1, r1 + r2) for (p1, r1), (p2, r2) in placements
+    if (p2 == p1 + 1 if kind == 'd' else p2 == p1)
+  })
+  if chunks:
+    entries = []
+    for start, residues in chunks:
+      # Only range contiguous starts that remove the SAME residues — in a periodic
+      # repeat (ABABAB) the removed pair changes along the range, and ranging them
+      # would misreport which residues went missing.
+      if entries and entries[-1][1] == residues and entries[-1][3] + 1 == start:
+        k, r, low, _ = entries[-1]
+        entries[-1] = (k, r, low, start)
+      else:
+        entries.append((kind, residues, start, start))
+    return entries
+
+  # No chunk: two independent edits. Their positions decouple, and each edit's own
+  # positions form a run of one repeated residue, so a per-edit range is exact.
+  entries = []
+  for slot in (0, 1):
+    positions = sorted({pl[slot][0] for pl in placements})
+    entries.append((kind, placements[0][slot][1], positions[0], positions[-1]))
+  return entries
 
 
 def format_indel_positions(query, matched):
-  """Render the Indel Positions column, e.g. 'd: A[6]', 'i: X[2,4]', or '[]' for
-  an exact match. Positions are 1-based; a range [low,high] collapses to [n] when
-  the position is unambiguous."""
-  edit = _indel_edit(query, matched)
-  if edit is None:
+  """Render the Indel Positions column, e.g. 'd: A[6]', 'd: AA[2,4]',
+  'd: C[3], d: E[5]', or '[]' for an exact match. Positions are 1-based; a range
+  collapses to a single number when the placement is unambiguous."""
+  edits = _indel_edits(query, matched)
+  if not edits:
     return '[]'
-  kind, residues, low, high = edit
-  span = f'[{low}]' if low == high else f'[{low},{high}]'
-  return f'{kind}: {residues}{span}'
+  return ', '.join(
+    f'{kind}: {residues}[{low}]' if low == high else f'{kind}: {residues}[{low},{high}]'
+    for kind, residues, low, high in edits
+  )
 
 
 class Matcher:
@@ -89,8 +137,10 @@ def __init__(
     if k != 0 and k < 2:
       raise ValueError('k must be >= 2.')
 
-    if max_indels > 1:
-      raise ValueError('max_indels > 1 is not yet supported. Only max_indels=1 has been validated.')
+    # max_indels=2 searches homogeneous edits only: two insertions OR two deletions,
+    # never one of each (a mixed pair is a substitution, which mismatch search covers).
+    if max_indels > 2:
+      raise ValueError('max_indels > 2 is not yet supported. Only max_indels<=2 has been validated.')
 
     if max_indels > 0 and max_mismatches > 0:
       raise ValueError('max_indels and max_mismatches are mutually exclusive.')

pepmatch/rs-engine/src/match.rs

@@ -557,13 +557,18 @@ fn is_terminal_deletion(q_idx: isize, query_len: usize, p_idx: isize, protein_le
     false
 }
 
+// `allow_del` / `allow_ins` mask the edit branches so a single alignment uses deletions
+// OR insertions, never both. extend_bidirectional fixes the mask for a seed's two
+// extensions, which is what stops a left-deletion mixing with a right-insertion.
 fn dfs(
     query: &[u8],
     q_idx: isize,
     protein: &[u8],
     p_idx: isize,
     indels_left: usize,
     direction: isize,
+    allow_del: bool,
+    allow_ins: bool,
 ) -> Vec<usize> {
     if (direction == 1 && q_idx >= query.len() as isize)
         || (direction == -1 && q_idx < 0)
@@ -580,22 +585,27 @@ fn dfs(
         && (p_idx as usize) < protein_len
         && query[q_idx as usize] == protein[p_idx as usize]
     {
-        for consumed in dfs(query, q_idx + direction, protein, p_idx + direction, indels_left, direction) {
+        for consumed in dfs(
+            query, q_idx + direction, protein, p_idx + direction, indels_left, direction,
+            allow_del, allow_ins,
+        ) {
             all_paths.push(consumed + 1);
         }
     }
 
     if indels_left > 0 {
         // Deletion branch: query pointer advances, protein stays, 0 protein chars consumed.
-        if !is_terminal_deletion(q_idx, query_len, p_idx, protein_len) {
+        if allow_del && !is_terminal_deletion(q_idx, query_len, p_idx, protein_len) {
             all_paths.extend(dfs(
                 query, q_idx + direction, protein, p_idx, indels_left - 1, direction,
+                allow_del, allow_ins,
             ));
         }
         // Insertion branch: protein pointer advances, query stays, 1 protein char consumed.
-        if p_idx >= 0 && (p_idx as usize) < protein_len {
+        if allow_ins && p_idx >= 0 && (p_idx as usize) < protein_len {
             for consumed in dfs(
                 query, q_idx, protein, p_idx + direction, indels_left - 1, direction,
+                allow_del, allow_ins,
             ) {
                 all_paths.push(consumed + 1);
             }
@@ -616,33 +626,39 @@ fn extend_bidirectional(
     let mut seen: HashSet<(usize, Vec<u8>)> = HashSet::new();
     let mut results: Vec<(usize, Vec<u8>)> = Vec::new();
 
-    for r_budget in 0..=max_indels {
-        let l_budget = max_indels - r_budget;
-
-        let r_paths = dfs(
-            query, (q_seed_start + k) as isize, protein, (p_hit_idx + k) as isize,
-            r_budget, 1,
-        );
-        let l_paths = dfs(
-            query, q_seed_start as isize - 1, protein, p_hit_idx as isize - 1,
-            l_budget, -1,
-        );
-
-        for &r_consumed in &r_paths {
-            for &l_consumed in &l_paths {
-                if l_consumed > p_hit_idx {
-                    continue;
-                }
-                let start = p_hit_idx - l_consumed;
-                let end = p_hit_idx + k + r_consumed;
-                if end > protein.len() {
-                    continue;
-                }
-                let matched = protein[start..end].to_vec();
-                let key = (start, matched.clone());
-                if !seen.contains(&key) {
-                    seen.insert(key);
-                    results.push((start, matched));
+    // Two homogeneous passes: deletions-only, then insertions-only. Fixing the edit type
+    // across BOTH extensions of the seed is what keeps an alignment homogeneous — masking
+    // per-side would still let a left-deletion pair with a right-insertion. `seen`
+    // collapses the exact match, which both passes find.
+    for &(allow_del, allow_ins) in &[(true, false), (false, true)] {
+        for r_budget in 0..=max_indels {
+            let l_budget = max_indels - r_budget;
+
+            let r_paths = dfs(
+                query, (q_seed_start + k) as isize, protein, (p_hit_idx + k) as isize,
+                r_budget, 1, allow_del, allow_ins,
+            );
+            let l_paths = dfs(
+                query, q_seed_start as isize - 1, protein, p_hit_idx as isize - 1,
+                l_budget, -1, allow_del, allow_ins,
+            );
+
+            for &r_consumed in &r_paths {
+                for &l_consumed in &l_paths {
+                    if l_consumed > p_hit_idx {
+                        continue;
+                    }
+                    let start = p_hit_idx - l_consumed;
+                    let end = p_hit_idx + k + r_consumed;
+                    if end > protein.len() {
+                        continue;
+                    }
+                    let matched = protein[start..end].to_vec();
+                    let key = (start, matched.clone());
+                    if !seen.contains(&key) {
+                        seen.insert(key);
+                        results.push((start, matched));
+                    }
                 }
             }
         }

pepmatch/tests/indel2_brute_force.py

@@ -0,0 +1,75 @@
+def _out_of_bounds_deletion(p_idx, protein_len):
+  """Blocks a deletion whose protein index is truly out of bounds — a database edge
+  effect rather than an indel. p_idx 0 and protein_len-1 are real residues, allowed."""
+  return p_idx < 0 or p_idx >= protein_len
+
+
+def _align_deletions(query, q_idx, window, w_idx, dels_left, protein, w_offset):
+  """query[q_idx:] aligns to window[w_idx:] using only deletions — no insertions, no
+  mismatches. Keeping the two edit kinds in separate passes is what makes a match
+  homogeneous: an alignment can never mix a deletion with an insertion."""
+  if q_idx == len(query):
+    return w_idx == len(window)
+
+  if w_idx < len(window) and query[q_idx] == window[w_idx]:
+    if _align_deletions(query, q_idx + 1, window, w_idx + 1, dels_left, protein, w_offset):
+      return True
+
+  if dels_left > 0:
+    # A deletion at the query's own first/last residue has no query-side context to
+    # confirm the residue is genuinely absent, so it is barred.
+    query_terminal = q_idx == 0 or q_idx == len(query) - 1
+    if not query_terminal and not _out_of_bounds_deletion(w_offset + w_idx, len(protein)):
+      if _align_deletions(query, q_idx + 1, window, w_idx, dels_left - 1, protein, w_offset):
+        return True
+
+  return False
+
+
+def _align_insertions(query, q_idx, window, w_idx, ins_left):
+  """query[q_idx:] aligns to window[w_idx:] using only insertions."""
+  if q_idx == len(query):
+    return w_idx == len(window)
+
+  if w_idx < len(window) and query[q_idx] == window[w_idx]:
+    if _align_insertions(query, q_idx + 1, window, w_idx + 1, ins_left):
+      return True
+
+  # q_idx > 0 bars an insertion before the first query residue; one after the last is
+  # unreachable because the base case returns as soon as the query is consumed. Both
+  # would only pad an exact match with an arbitrary flanking residue.
+  if ins_left > 0 and q_idx > 0 and w_idx < len(window):
+    if _align_insertions(query, q_idx, window, w_idx + 1, ins_left - 1):
+      return True
+
+  return False
+
+
+def brute_force_search(query, protein, max_indels=2):
+  """Enumerate every (start, matched_sequence) where query aligns to a window of
+  protein using at most max_indels homogeneous indels — insertions only or deletions
+  only, never both.
+
+  Checks each candidate window directly against the definition of a valid match (no
+  seeding, no bidirectional extension), so it can serve as an independent oracle for
+  the Rust engine.
+  """
+  query_len = len(query)
+  protein_len = len(protein)
+  results = set()
+
+  for start in range(protein_len):
+    for window_len in range(max(0, query_len - max_indels), query_len + max_indels + 1):
+      end = start + window_len
+      if end > protein_len:
+        continue
+      window = protein[start:end]
+      if window_len <= query_len:
+        if _align_deletions(
+          query, 0, window, 0, query_len - window_len, protein, start
+        ):
+          results.add((start, window))
+      elif _align_insertions(query, 0, window, 0, window_len - query_len):
+        results.add((start, window))
+
+  return results