diff --git a/pepmatch/matcher.py b/pepmatch/matcher.py index e103303..5a4e7ae 100755 --- a/pepmatch/matcher.py +++ b/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( diff --git a/pepmatch/tests/test_indel_search.py b/pepmatch/tests/test_indel_search.py index bc18a6e..d1d674f 100644 --- a/pepmatch/tests/test_indel_search.py +++ b/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)