feat: add 1-indel search via bidirectional DFS in Rust (#24)

- Adds max_indels=1 search to the Rust engine: pigeonhole-seeded bidirectional DFS (run_indel -> indel_search_peptide -> extend_bidirectional).
- End-position deletions and terminal insertions are disallowed.
- Output carries an Indels column in indel mode (Mismatches elsewhere, never both); indel is
rejected with mismatches/best_match/counts_only/discontinuous. CLI -i/--max_indels
flag + README section. 41 tests incl. a brute-force property oracle.

Author: Roman-Young

.gitignore

@@ -1,5 +1,6 @@
 .~*
 *.db
+*.pepidx
 *.pickle
 *.pkl
 *.pyc

README.md

@@ -93,6 +93,18 @@ df = Matcher(
 ).match()
 ```
 
+#### Indel Searching
+
+Search allowing insertions and deletions (indels) instead of substitution mismatches. The k-mer size is chosen automatically from your query lengths — no manual preprocessing needed.
+```python
+df = Matcher(
+  query='neoepitopes.fasta',
+  proteome_file='human.fasta',
+  max_indels=1
+).match()
+```
+Currently limited to `max_indels=1`; mutually exclusive with `max_mismatches`.
+
 #### Best Match
 
 Automatically finds the optimal match for each peptide by trying different k-mer sizes and mismatch thresholds. No manual preprocessing required.
@@ -163,6 +175,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`.
 * `-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

@@ -2,7 +2,7 @@
 import polars as pl
 from pathlib import Path
 from Bio import SeqIO
-from ._rs import rs_preprocess, rs_match, rs_discontinuous, rs_metadata, rs_match_counts
+from ._rs import rs_preprocess, rs_match, rs_discontinuous, rs_metadata, rs_match_counts, rs_indel_match
 
 VALID_OUTPUT_FORMATS = ['dataframe', 'csv', 'tsv', 'xlsx', 'json']
 FASTA_EXTENSIONS = {
@@ -31,6 +31,7 @@ def __init__(
     query,
     proteome_file,
     max_mismatches=0,
+    max_indels=0,
     k=0,
     preprocessed_files_path='.',
     best_match=False,
@@ -53,6 +54,20 @@ 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.')
+
+    if max_indels > 0 and max_mismatches > 0:
+      raise ValueError('max_indels and max_mismatches are mutually exclusive.')
+
+    if max_indels > 0 and best_match:
+      raise ValueError('max_indels and best_match are not yet supported together.')
+
+    if max_indels > 0 and counts_only:
+      raise ValueError('max_indels and counts_only are not yet supported together.')
+
+    self.max_indels = max_indels
+
     if output_format not in VALID_OUTPUT_FORMATS:
       raise ValueError(
         f'Invalid output format. Choose from: {VALID_OUTPUT_FORMATS}'
@@ -73,6 +88,8 @@ def __init__(
 
     self.query = self._parse_query(query)
     self.discontinuous_epitopes = self._find_discontinuous_epitopes()
+    if self.max_indels > 0 and self.discontinuous_epitopes:
+      raise ValueError('max_indels and discontinuous epitopes are not supported together.')
     self.query = self._clean_query()
     if not self.query and not self.discontinuous_epitopes:
       raise ValueError('Query is empty.')
@@ -128,7 +145,9 @@ def match(self):
       return output_matches(df, self.output_format, self.output_name)
 
     if self.query:
-      if self.best_match and self.k_specified:
+      if self.max_indels > 0:
+        linear_df = self.indel_search()
+      elif self.best_match and self.k_specified:
         results = self._search(self.k, self.max_mismatches)
         linear_df = self._best_match_filter(self._to_dataframe(results))
       elif self.best_match:
@@ -155,6 +174,19 @@ def match(self):
     else:
       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))
+    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 "
+            f"(this may take a moment)...")
+      rs_preprocess(self.proteome_file, k, pepidx_path)
+    print(f"Searching {len(self.query)} peptides against {self.proteome_name} "
+          f"(k={k}, max_indels={self.max_indels})...")
+    results = rs_indel_match(pepidx_path, self.query, self.max_indels)
+    return self._to_dataframe(results, is_indels=True)
+
   def _pepidx_path(self, k):
     return os.path.join(
       self.preprocessed_files_path, f'{self.proteome_name}_{k}mers.pepidx'
@@ -327,20 +359,31 @@ def _metadata_table(self) -> pl.DataFrame:
       pl.col('SwissProt Reviewed').cast(pl.Int64).cast(pl.Boolean),
     ])
 
-  FINAL_COLUMNS = [
-    'Query ID','Query Sequence','Matched Sequence','Protein ID','Protein Name','Species',
-    'Taxon ID','Gene','Mismatches','Mutated Positions','Index start','Index end',
-    'Protein Existence Level','Gene Priority','SwissProt Reviewed',
-  ]
+  def _final_columns(self, is_indels):
+    # One edit-count column per mode, in a fixed position: Indels for indel search,
+    # Mismatches for every other mode. Never both — an all-zero twin column would
+    # imply a search that didn't run.
+    edit_col = 'Indels' if is_indels else 'Mismatches'
+    return [
+      'Query ID','Query Sequence','Matched Sequence','Protein ID','Protein Name','Species',
+      'Taxon ID','Gene', edit_col, 'Mutated Positions','Index start','Index end',
+      'Protein Existence Level','Gene Priority','SwissProt Reviewed',
+    ]
 
-  def _to_dataframe(self, cols):
+  def _to_dataframe(self, cols, is_indels=False):
     """Build the result DataFrame from columnar Rust output, reconstructing protein
     metadata via a single join (instead of cloning it into every hit row)."""
     qid, qseq, matched, pnum, mm, mutated, istart, iend = cols
 
+    # rs_indel_match packs the indel count into the same slot rs_match/rs_discontinuous
+    # use for mismatches, so the values are identical in shape — only the column name
+    # differs by mode (Indels vs Mismatches), and only one is ever emitted.
+    edit_col = 'Indels' if is_indels else 'Mismatches'
+    final_columns = self._final_columns(is_indels)
+
     if not qid:
-      schema = {c: pl.Utf8 for c in self.FINAL_COLUMNS}
-      for c in ('Mismatches','Index start','Index end','Protein Existence Level','Gene Priority'):
+      schema = {c: pl.Utf8 for c in final_columns}
+      for c in (edit_col, 'Index start', 'Index end', 'Protein Existence Level', 'Gene Priority'):
         schema[c] = pl.Int64
       schema['SwissProt Reviewed'] = pl.Boolean
       return pl.DataFrame(schema=schema)
@@ -350,7 +393,7 @@ def _to_dataframe(self, cols):
       'Query Sequence': qseq,
       'Matched Sequence': matched,
       'protein_num': pl.Series(pnum, dtype=pl.UInt32),
-      'Mismatches': pl.Series(mm, dtype=pl.Int64),
+      edit_col: pl.Series(mm, dtype=pl.Int64),
       'Mutated Positions': mutated,
       'Index start': pl.Series(istart, dtype=pl.Int64),
       'Index end': pl.Series(iend, dtype=pl.Int64),
@@ -366,4 +409,4 @@ def _to_dataframe(self, cols):
           .alias('Protein ID')
       )
 
-    return df.drop('Sequence Version').select(self.FINAL_COLUMNS)
+    return df.drop('Sequence Version').select(final_columns)

pepmatch/rs-engine/Cargo.lock

@@ -34,6 +34,11 @@ fn rs_match_counts(pepidx_path: &str, peptides: Vec<(String, String)>, k: usize,
     matching::run_counts(pepidx_path, peptides, k, max_mismatches)
 }
 
+#[pyfunction]
+fn rs_indel_match(pepidx_path: &str, peptides: Vec<(String, String)>, indels_allowed: usize) -> matching::Columns {
+    matching::run_indel(pepidx_path, peptides, indels_allowed)
+}
+
 #[pymodule]
 fn _rs(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(rs_version, m)?)?;
@@ -42,5 +47,6 @@ fn _rs(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(rs_discontinuous, m)?)?;
     m.add_function(wrap_pyfunction!(rs_metadata, m)?)?;
     m.add_function(wrap_pyfunction!(rs_match_counts, m)?)?;
+    m.add_function(wrap_pyfunction!(rs_indel_match, m)?)?;
     Ok(())
 }