fix: honor explicit k in indel search, clamp above the pigeonhole optimal (#29)

indel_search ignored the user's k and always used the optimized value, so an explicit k was silently swapped and its k-mer table rebuilt. It now honors an explicit k up to the pigeonhole optimal (reusing an existing table rather than rebuilding) and clamps a larger k down to the optimal with a warning, preserving complete recall. _auto_k is generalized to take the edit count so the mismatch and indel paths share one pigeonhole formula.

Fixes #28.

Author: Roman-Young

pepmatch/matcher.py

@@ -182,7 +182,7 @@ def match(self):
     discontinuous_df = pl.DataFrame()
 
     if self.counts_only:
-      k = self.k if self.k_specified else self._auto_k()
+      k = self.k if self.k_specified else self._auto_k(self.max_mismatches)
       df = self._counts_to_dataframe(self._search_counts(k, self.max_mismatches))
       if self.output_format == 'dataframe':
         return df
@@ -197,7 +197,7 @@ def match(self):
       elif self.best_match:
         linear_df = self.best_match_search()
       else:
-        k = self.k if self.k_specified else self._auto_k()
+        k = self.k if self.k_specified else self._auto_k(self.max_mismatches)
         results = self._search(k, self.max_mismatches)
         linear_df = self._to_dataframe(results)
 
@@ -220,8 +220,20 @@ def match(self):
       output_matches(df, self.output_format, self.output_name)
 
   def indel_search(self):
-    min_len = min(len(seq) for _, seq in self.query)
-    k = max(2, min_len // (self.max_indels + 1))
+    # Honor an explicitly-passed k, but never above the pigeonhole ceiling: a k
+    # larger than the optimal drops below max_indels+1 disjoint seeds and forfeits
+    # complete recall, so clamp it down. (The constructor already rejects
+    # k > shortest query length, so an explicit k reaching here is always <= min_len.)
+    optimized_k = self._auto_k(self.max_indels)
+    if self.k_specified:
+      k = self.k
+      if k > optimized_k:
+        print(f"Requested k={k} exceeds k={optimized_k}, the largest k that "
+              f"guarantees complete recall for these queries (pigeonhole); "
+              f"using k={optimized_k}.")
+        k = optimized_k
+    else:
+      k = optimized_k
     pepidx_path = self._pepidx_path(k)
     if not os.path.isfile(pepidx_path):
       print(f"No preprocessed file found for k={k}, building index now "
@@ -232,14 +244,14 @@ def indel_search(self):
     results = rs_indel_match(pepidx_path, self.query, self.max_indels)
     return self._to_dataframe(results, is_indels=True)
 
-  def _auto_k(self):
-    """Optimal k for a seed-based mismatch search (PEPMatch paper / pigeonhole):
-    a length-L peptide with m mismatches must split into at least m+1 disjoint
-    k-mers so one is guaranteed mismatch-free, i.e. k = floor(L / (m + 1)). Uses
-    the shortest query peptide so the guarantee holds for every peptide; floored
-    at 2 (the preprocessor minimum)."""
+  def _auto_k(self, edits):
+    """Optimal k for a seed-based search (PEPMatch paper / pigeonhole): a length-L
+    peptide with `edits` allowed edits (mismatches or indels) must split into at
+    least edits+1 disjoint k-mers so one is guaranteed edit-free, i.e.
+    k = floor(L / (edits + 1)). Uses the shortest query peptide so the guarantee
+    holds for every peptide; floored at 2 (the preprocessor minimum)."""
     min_len = min(len(seq) for _, seq in self.query)
-    return max(2, min_len // (self.max_mismatches + 1))
+    return max(2, min_len // (edits + 1))
 
   def _pepidx_path(self, k):
     return os.path.join(

pepmatch/tests/test_indel_search.py

@@ -254,3 +254,59 @@ def test_indel_positions_insertion_end_to_end(tmp_path):
   row = df.filter(pl.col('Matched Sequence') == 'NALVEXATRFC')
   assert row.height == 1
   assert row['Indel Positions'].item() == 'i: X[6]'
+
+
+def test_indel_honors_explicit_k_up_to_optimal_and_reuses_table(tmp_path, capsys):
+  # Issue #28: a user who already built a k-mer table (e.g. for a mismatch run) and
+  # passes that same k to indel search should have it reused, not silently swapped
+  # for the optimized k and rebuilt. Query len 10, max_indels=1 -> optimized k = 5;
+  # k=3 <= 5 so it is honored, and the pre-built 3-mer index is used as-is.
+  proteome_path = tmp_path / 'proteome.fasta'
+  proteome_path.write_text('>ProtDel\nMKVNALVETRFCGHI\n')
+  # the user "already has" a 3-mer table from an earlier run
+  Matcher(query=['NALVEATRFC'], proteome_file=str(proteome_path), max_mismatches=0,
+          k=3, preprocessed_files_path=str(tmp_path), output_format='dataframe').match()
+  assert (tmp_path / 'proteome_3mers.pepidx').exists()
+  capsys.readouterr()  # discard the pre-build output
+  df = Matcher(query=['NALVEATRFC'], proteome_file=str(proteome_path), max_indels=1,
+               k=3, preprocessed_files_path=str(tmp_path), output_format='dataframe').match()
+  out = capsys.readouterr().out
+  assert '(k=3,' in out                       # honored the user's k, not the optimal 5
+  assert 'No preprocessed file' not in out    # reused the existing table, no rebuild
+  assert 'NALVETRFC' in df['Matched Sequence'].to_list()  # recall preserved at k=3
+
+
+def test_indel_clamps_k_above_optimal_and_warns(tmp_path, capsys):
+  # An explicit k above the pigeonhole optimal drops below max_indels+1 disjoint
+  # seeds and would forfeit complete recall, so indel search clamps it to the
+  # optimal and warns. Query len 10, max_indels=1 -> optimized k = 5; k=7 (still
+  # <= min_len, so it passes the constructor guard) is clamped down to 5.
+  proteome_path = tmp_path / 'proteome.fasta'
+  proteome_path.write_text('>ProtDel\nMKVNALVETRFCGHI\n')
+  df = Matcher(query=['NALVEATRFC'], proteome_file=str(proteome_path), max_indels=1,
+               k=7, preprocessed_files_path=str(tmp_path), output_format='dataframe').match()
+  out = capsys.readouterr().out
+  assert 'Requested k=7 exceeds k=5' in out    # warned about the clamp
+  assert '(k=5,' in out                        # searched at the clamped optimal
+  assert 'NALVETRFC' in df['Matched Sequence'].to_list()   # recall still complete
+  assert (tmp_path / 'proteome_5mers.pepidx').exists()     # built the optimal index
+  assert not (tmp_path / 'proteome_7mers.pepidx').exists() # never the too-large one
+
+
+def test_indel_unspecified_k_derives_optimal(tmp_path, capsys):
+  # With no k passed, indel search derives the pigeonhole optimal: max(2, 10//2) = 5.
+  proteome_path = tmp_path / 'proteome.fasta'
+  proteome_path.write_text('>ProtDel\nMKVNALVETRFCGHI\n')
+  Matcher(query=['NALVEATRFC'], proteome_file=str(proteome_path), max_indels=1,
+          preprocessed_files_path=str(tmp_path), output_format='dataframe').match()
+  out = capsys.readouterr().out
+  assert '(k=5,' in out
+  assert (tmp_path / 'proteome_5mers.pepidx').exists()
+
+
+def test_indel_k_exceeding_shortest_query_raises():
+  # The constructor rejects an explicit k larger than the shortest query peptide
+  # (shorter peptides would be silently skipped). The guard is global, so it also
+  # governs indel search -- an explicit k reaching indel_search is always <= min_len.
+  with pytest.raises(ValueError, match='cannot exceed the shortest query peptide'):
+    Matcher(query=['NALVEATRFC'], proteome_file='unused.fasta', max_indels=1, k=11)