From 46a7f2d86713fc2a201a253fd24b178aa548d77d Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung@Romans-MacBook-Air-3.local>
Date: Thu, 25 Jun 2026 12:18:48 -0700
Subject: [PATCH 01/15] feat: add 1-indel search via exhaustive bidirectional
 DFS in Rust

Implements insertion/deletion matching as a new search path in the Rust
engine. Uses the pigeonhole principle for seed selection, bidirectional
DFS with full budget splitting, and two-level deduplication across seeds.
Adds max_indels parameter to Matcher with mutual exclusion against
max_mismatches. max_indels > 1 raises ValueError pending validation.
---
 pepmatch/indel.py               | 114 ++++++++++++++++++
 pepmatch/matcher.py             |  27 ++++-
 pepmatch/rs-engine/src/lib.rs   |   6 +
 pepmatch/rs-engine/src/match.rs | 207 ++++++++++++++++++++++++++++++++
 4 files changed, 352 insertions(+), 2 deletions(-)
 create mode 100644 pepmatch/indel.py

diff --git a/pepmatch/indel.py b/pepmatch/indel.py
new file mode 100644
index 0000000..4ad6833
--- /dev/null
+++ b/pepmatch/indel.py
@@ -0,0 +1,114 @@
+def _is_terminal_deletion(q_idx, query_len, p_idx, protein_len):
+  """Return True if a deletion at this position is at a query or protein boundary.
+
+  Terminal deletions are forbidden because a deletion at the first or last
+  query position produces a match indistinguishable from a shorter peptide,
+  and a deletion at a protein boundary arises from an edge effect rather
+  than a true biological indel event.
+  """
+  if q_idx == 0 or q_idx == query_len - 1:
+    return True
+  if p_idx <= 0 or p_idx >= protein_len - 1:
+    return True
+  return False
+
+
+def _dfs(query, q_idx, protein, p_idx, indels_left, direction):
+  """Exhaustive depth-first search extending one direction from a seed hit.
+
+  Explores three branches at each step: a match branch (both pointers
+  advance), a deletion branch (query pointer advances, protein stays,
+  consumes 0 protein chars), and an insertion branch (protein pointer
+  advances, query stays, consumes 1 protein char). Returns a list of
+  integers, each the number of protein characters consumed by one valid
+  complete alignment path.
+
+  Args:
+    query: the full query peptide string.
+    q_idx: current query position, advances by direction.
+    protein: the full protein sequence string.
+    p_idx: current protein position, advances by direction.
+    indels_left: remaining indel budget for this direction.
+    direction: +1 for right extension, -1 for left extension.
+  """
+  if (direction == 1 and q_idx >= len(query)) or \
+     (direction == -1 and q_idx < 0):
+    return [0]
+
+  all_paths = []
+
+  if 0 <= p_idx < len(protein) and query[q_idx] == protein[p_idx]:
+    for consumed in _dfs(query, q_idx + direction, protein, p_idx + direction, indels_left, direction):
+      all_paths.append(consumed + 1)
+
+  if indels_left > 0:
+    if not _is_terminal_deletion(q_idx, len(query), p_idx, len(protein)):
+      all_paths.extend(
+        _dfs(query, q_idx + direction, protein, p_idx, indels_left - 1, direction)
+      )
+
+    if 0 <= p_idx < len(protein):
+      for consumed in _dfs(query, q_idx, protein, p_idx + direction, indels_left - 1, direction):
+        all_paths.append(consumed + 1)
+
+  return all_paths
+
+
+def extend_bidirectional(query, q_seed_start, p_hit_idx, protein, k, max_indels=1):
+  """Verify and enumerate all indel-consistent matches from a seed hit.
+
+  For each allocation of the indel budget between left and right extensions,
+  runs _dfs in both directions and combines results. Returns all unique
+  (start_0based, matched_sequence) pairs found at this seed hit position.
+
+  Args:
+    query: the full query peptide string.
+    q_seed_start: 0-based start of the seed k-mer within the query.
+    p_hit_idx: 0-based start of the seed hit within protein.
+    protein: the full protein sequence string.
+    k: seed k-mer length.
+    max_indels: maximum total indels allowed across both extensions.
+  """
+  seen = set()
+  results = []
+
+  for r_budget in range(max_indels + 1):
+    l_budget = max_indels - r_budget
+
+    r_paths = _dfs(query, q_seed_start + k, protein, p_hit_idx + k, r_budget, direction=1)
+    l_paths = _dfs(query, q_seed_start - 1, protein, p_hit_idx - 1, l_budget, direction=-1)
+
+    for r_consumed in r_paths:
+      for l_consumed in l_paths:
+        start = p_hit_idx - l_consumed
+        end = p_hit_idx + k + r_consumed
+        matched_seq = protein[start:end]
+        key = (start, matched_seq)
+        if key not in seen:
+          seen.add(key)
+          results.append(key)
+
+  return results
+
+
+def minimal_coverage_seeds(query, k):
+  """Return non-overlapping k-mer seeds plus a final overlapping seed for
+  full query coverage.
+
+  The pigeonhole principle guarantees that for a query of length L with at
+  most d indels, at least one of floor(L / (d+1)) non-overlapping seeds
+  must hit exactly. Appends one final seed anchored at the query end if
+  the last strided seed does not reach it.
+
+  Args:
+    query: the query peptide string.
+    k: seed k-mer length.
+  """
+  seeds = []
+  query_len = len(query)
+  for j in range(0, query_len - k + 1, k):
+    seeds.append((query[j:j + k], j))
+  last_start = query_len - k
+  if seeds and seeds[-1][1] != last_start:
+    seeds.append((query[last_start:], last_start))
+  return seeds
\ No newline at end of file
diff --git a/pepmatch/matcher.py b/pepmatch/matcher.py
index e44a23c..ddd52b6 100755
--- a/pepmatch/matcher.py
+++ b/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,14 @@ 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.')
+
+    self.max_indels = max_indels
+
     if output_format not in VALID_OUTPUT_FORMATS:
       raise ValueError(
         f'Invalid output format. Choose from: {VALID_OUTPUT_FORMATS}'
@@ -128,7 +137,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 +166,18 @@ 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"Preprocessing {self.proteome_name} with k={k}...")
+      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)
+
   def _pepidx_path(self, k):
     return os.path.join(
       self.preprocessed_files_path, f'{self.proteome_name}_{k}mers.pepidx'
diff --git a/pepmatch/rs-engine/src/lib.rs b/pepmatch/rs-engine/src/lib.rs
index 71124dd..fa0ea0b 100644
--- a/pepmatch/rs-engine/src/lib.rs
+++ b/pepmatch/rs-engine/src/lib.rs
@@ -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(())
 }
diff --git a/pepmatch/rs-engine/src/match.rs b/pepmatch/rs-engine/src/match.rs
index 2173c77..9592cd4 100644
--- a/pepmatch/rs-engine/src/match.rs
+++ b/pepmatch/rs-engine/src/match.rs
@@ -526,6 +526,213 @@ fn metadata_columns(index: &PepIndex) -> MetaColumns {
     (pnum, a, b, c, d, e, g, h, i, j)
 }
 
+fn minimal_coverage_seeds(query: &[u8], k: usize) -> Vec<(&[u8], usize)> {
+    let query_len = query.len();
+    let mut seeds: Vec<(&[u8], usize)> = Vec::new();
+    let mut j = 0;
+    while j + k <= query_len {
+        seeds.push((&query[j..j + k], j));
+        j += k;
+    }
+    if query_len >= k {
+        let last_start = query_len - k;
+        if seeds.last().map(|(_, s)| *s) != Some(last_start) {
+            seeds.push((&query[last_start..], last_start));
+        }
+    }
+    seeds
+}
+
+fn is_terminal_deletion(q_idx: isize, query_len: usize, p_idx: isize, protein_len: usize) -> bool {
+    if q_idx == 0 || q_idx == query_len as isize - 1 {
+        return true;
+    }
+    if p_idx <= 0 || p_idx >= protein_len as isize - 1 {
+        return true;
+    }
+    false
+}
+
+fn dfs(
+    query: &[u8],
+    q_idx: isize,
+    protein: &[u8],
+    p_idx: isize,
+    indels_left: usize,
+    direction: isize,
+) -> Vec<usize> {
+    if (direction == 1 && q_idx >= query.len() as isize)
+        || (direction == -1 && q_idx < 0)
+    {
+        return vec![0];
+    }
+
+    let mut all_paths: Vec<usize> = Vec::new();
+    let query_len = query.len();
+    let protein_len = protein.len();
+
+    // Match branch: both pointers advance, consumes 1 protein char.
+    if p_idx >= 0
+        && (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) {
+            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) {
+            all_paths.extend(dfs(
+                query, q_idx + direction, protein, p_idx, indels_left - 1, direction,
+            ));
+        }
+        // Insertion branch: protein pointer advances, query stays, 1 protein char consumed.
+        if p_idx >= 0 && (p_idx as usize) < protein_len {
+            for consumed in dfs(
+                query, q_idx, protein, p_idx + direction, indels_left - 1, direction,
+            ) {
+                all_paths.push(consumed + 1);
+            }
+        }
+    }
+
+    all_paths
+}
+
+fn extend_bidirectional(
+    query: &[u8],
+    q_seed_start: usize,
+    p_hit_idx: usize,
+    protein: &[u8],
+    k: usize,
+    max_indels: usize,
+) -> Vec<(usize, Vec<u8>)> {
+    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));
+                }
+            }
+        }
+    }
+
+    results
+}
+
+fn indel_search_peptide(
+    query_id: &str,
+    peptide: &str,
+    indels_allowed: usize,
+    index: &PepIndex,
+) -> Vec<HitRecord> {
+    let pep_bytes = peptide.as_bytes();
+    let peptide_len = pep_bytes.len();
+    let k = index.k;
+
+    if peptide_len < k {
+        return vec![miss_record(query_id, peptide)];
+    }
+
+    let seeds = minimal_coverage_seeds(pep_bytes, k);
+    let mut seen: HashSet<(usize, usize, Vec<u8>)> = HashSet::new();
+    let mut records: Vec<HitRecord> = Vec::new();
+
+    for (seed_bytes, q_seed_start) in &seeds {
+        if let Some(hit_positions) = index.lookup(seed_bytes) {
+            for encoded_hit in hit_positions {
+                let prot_num = (encoded_hit / PROTEIN_INDEX_MULTIPLIER) as usize;
+                if prot_num == 0 || prot_num > index.num_proteins {
+                    continue;
+                }
+                let local_p_idx = (encoded_hit % PROTEIN_INDEX_MULTIPLIER) as usize;
+
+                let prot_base = index.protein_offset(prot_num) as usize;
+                let prot_end = if prot_num < index.num_proteins {
+                    index.protein_offset(prot_num + 1) as usize
+                } else {
+                    index.seq_len
+                };
+                let protein = &index.mmap[index.seq_offset + prot_base..index.seq_offset + prot_end];
+
+                for (true_start, matched) in
+                    extend_bidirectional(pep_bytes, *q_seed_start, local_p_idx, protein, k, indels_allowed)
+                {
+                    let dedup_key = (prot_num, true_start, matched.clone());
+                    if seen.contains(&dedup_key) {
+                        continue;
+                    }
+                    seen.insert(dedup_key);
+
+                    let matched_str = String::from_utf8_lossy(&matched).into_owned();
+                    let n_indels = matched.len().abs_diff(peptide_len);
+
+                    records.push(HitRecord {
+                        query_id: query_id.to_string(),
+                        query_seq: peptide.to_string(),
+                        matched: Some(matched_str),
+                        protein_num: Some(prot_num as u32),
+                        mismatches: Some(n_indels as i64),
+                        mutated: String::from("[]"),
+                        idx_start: Some((true_start + 1) as i64),
+                        idx_end: Some((true_start + matched.len()) as i64),
+                    });
+                }
+            }
+        }
+    }
+
+    if records.is_empty() {
+        vec![miss_record(query_id, peptide)]
+    } else {
+        records
+    }
+}
+
+pub(crate) fn run_indel(
+    pepidx_path: &str,
+    peptides: Vec<(String, String)>,
+    indels_allowed: usize,
+) -> Columns {
+    let index = PepIndex::open(pepidx_path);
+
+    let records: Vec<HitRecord> = peptides
+        .par_iter()
+        .flat_map_iter(|(query_id, peptide)| {
+            indel_search_peptide(query_id, peptide, indels_allowed, &index).into_iter()
+        })
+        .collect();
+
+    unzip_records(records)
+}
+
 pub(crate) fn run_metadata(pepidx_path: &str) -> MetaColumns {
     let index = PepIndex::open(pepidx_path);
     metadata_columns(&index)

From 27e48f144311f42e9c67822a50e4f92195908f2f Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung@Romans-MacBook-Air-3.local>
Date: Thu, 25 Jun 2026 12:19:02 -0700
Subject: [PATCH 02/15] test: add indel search test covering 1-deletion and
 1-insertion cases

Tests verified by hand against Q8V336 (NCAP_DUGBA) in the test proteome.
---
 pepmatch/tests/data/indel_expected.csv |  3 +++
 pepmatch/tests/data/indel_query.fasta  |  4 +++
 pepmatch/tests/test_indel_search.py    | 36 ++++++++++++++++++++++++++
 3 files changed, 43 insertions(+)
 create mode 100644 pepmatch/tests/data/indel_expected.csv
 create mode 100644 pepmatch/tests/data/indel_query.fasta
 create mode 100644 pepmatch/tests/test_indel_search.py

diff --git a/pepmatch/tests/data/indel_expected.csv b/pepmatch/tests/data/indel_expected.csv
new file mode 100644
index 0000000..5c02269
--- /dev/null
+++ b/pepmatch/tests/data/indel_expected.csv
@@ -0,0 +1,3 @@
+Query Sequence,Matched Sequence,Protein ID
+NALVEARFC,NALVEATRFC,Q8V336.1
+EQLVPSYQQA,EQLVSYQQA,Q8V336.1
diff --git a/pepmatch/tests/data/indel_query.fasta b/pepmatch/tests/data/indel_query.fasta
new file mode 100644
index 0000000..4611623
--- /dev/null
+++ b/pepmatch/tests/data/indel_query.fasta
@@ -0,0 +1,4 @@
+>1
+NALVEARFC
+>2
+EQLVPSYQQA
diff --git a/pepmatch/tests/test_indel_search.py b/pepmatch/tests/test_indel_search.py
new file mode 100644
index 0000000..5316b3e
--- /dev/null
+++ b/pepmatch/tests/test_indel_search.py
@@ -0,0 +1,36 @@
+import pytest
+import polars as pl
+import polars.testing as plt
+from pathlib import Path
+from pepmatch import Matcher
+
+@pytest.fixture
+def proteome_path() -> Path:
+  return Path(__file__).parent / 'data' / 'proteome.fasta'
+
+@pytest.fixture
+def query_path() -> Path:
+  return Path(__file__).parent / 'data' / 'indel_query.fasta'
+
+@pytest.fixture
+def expected_path() -> Path:
+  return Path(__file__).parent / 'data' / 'indel_expected.csv'
+
+def test_indel_search(proteome_path, query_path, expected_path):
+  df = Matcher(
+    query=query_path,
+    proteome_file=proteome_path,
+    max_indels=1,
+    output_format='dataframe'
+  ).match()
+
+  expected_df = pl.read_csv(expected_path)
+  cols_to_compare = ['Query Sequence', 'Matched Sequence', 'Protein ID']
+  sort_key = ['Query Sequence', 'Matched Sequence']
+  df_sorted = df.sort(sort_key)
+  expected_df_sorted = expected_df.sort(sort_key)
+
+  plt.assert_frame_equal(
+    df_sorted.select(cols_to_compare),
+    expected_df_sorted.select(cols_to_compare)
+  )

From eb5e4fb861b14cd5b76e5ec6ee8f2a5bb910c183 Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung@Romans-MacBook-Air-3.local>
Date: Thu, 25 Jun 2026 12:56:55 -0700
Subject: [PATCH 03/15] test: add terminal deletion guard test; exclude pepidx
 from git

QNALVEATRFC is constructed so the only match in the test proteome requires
deleting Q at query position 0 (terminal). Verified that removing the guard
produces a spurious NALVEATRFC hit; the guard correctly blocks it.
Also adds *.pepidx to .gitignore to prevent index files from being staged.
---
 .gitignore                          |  1 +
 pepmatch/tests/test_indel_search.py | 17 +++++++++++++++++
 2 files changed, 18 insertions(+)

diff --git a/.gitignore b/.gitignore
index 006f268..75d39b1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,6 @@
 .~*
 *.db
+*.pepidx
 *.pickle
 *.pkl
 *.pyc
diff --git a/pepmatch/tests/test_indel_search.py b/pepmatch/tests/test_indel_search.py
index 5316b3e..0eafa02 100644
--- a/pepmatch/tests/test_indel_search.py
+++ b/pepmatch/tests/test_indel_search.py
@@ -4,6 +4,9 @@
 from pathlib import Path
 from pepmatch import Matcher
 
+# QNALVEATRFC is designed so the only match in the test proteome would require
+# deleting Q at query position 0 (a terminal deletion). The guard must block it.
+
 @pytest.fixture
 def proteome_path() -> Path:
   return Path(__file__).parent / 'data' / 'proteome.fasta'
@@ -34,3 +37,17 @@ def test_indel_search(proteome_path, query_path, expected_path):
     df_sorted.select(cols_to_compare),
     expected_df_sorted.select(cols_to_compare)
   )
+
+
+def test_terminal_deletion_blocked(proteome_path):
+  df = Matcher(
+    query=['QNALVEATRFC'],
+    proteome_file=proteome_path,
+    max_indels=1,
+    output_format='dataframe'
+  ).match()
+  assert df['Matched Sequence'].is_null().all(), (
+    'Terminal deletion should be blocked but got match(es): '
+    + str(df.filter(pl.col('Matched Sequence').is_not_null())
+            .select(['Query Sequence', 'Matched Sequence']))
+  )

From 09a6757ffa046879a46ba4a3a53d959133c768b7 Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung9981@gmail.com>
Date: Thu, 25 Jun 2026 15:17:26 -0700
Subject: [PATCH 04/15] fix: restrict terminal deletion guard to protein
 boundaries only

Deletions at query position 0 or query_len-1 are now permitted when the
protein position is not at a sequence boundary, allowing valid biological
indel matches that were previously blocked.
---
 pepmatch/indel.py                   | 15 ++++++---------
 pepmatch/rs-engine/src/match.rs     |  9 +++------
 pepmatch/tests/test_indel_search.py | 11 ++++++-----
 3 files changed, 15 insertions(+), 20 deletions(-)

diff --git a/pepmatch/indel.py b/pepmatch/indel.py
index 4ad6833..71476ee 100644
--- a/pepmatch/indel.py
+++ b/pepmatch/indel.py
@@ -1,13 +1,10 @@
-def _is_terminal_deletion(q_idx, query_len, p_idx, protein_len):
-  """Return True if a deletion at this position is at a query or protein boundary.
+def _is_terminal_deletion(p_idx, protein_len):
+  """Return True if a deletion at this protein position is at a sequence boundary.
 
-  Terminal deletions are forbidden because a deletion at the first or last
-  query position produces a match indistinguishable from a shorter peptide,
-  and a deletion at a protein boundary arises from an edge effect rather
-  than a true biological indel event.
+  Deletions at p_idx == 0 or p_idx == protein_len - 1 are forbidden because
+  they arise from a database boundary edge effect rather than a true biological
+  indel event.
   """
-  if q_idx == 0 or q_idx == query_len - 1:
-    return True
   if p_idx <= 0 or p_idx >= protein_len - 1:
     return True
   return False
@@ -42,7 +39,7 @@ def _dfs(query, q_idx, protein, p_idx, indels_left, direction):
       all_paths.append(consumed + 1)
 
   if indels_left > 0:
-    if not _is_terminal_deletion(q_idx, len(query), p_idx, len(protein)):
+    if not _is_terminal_deletion(p_idx, len(protein)):
       all_paths.extend(
         _dfs(query, q_idx + direction, protein, p_idx, indels_left - 1, direction)
       )
diff --git a/pepmatch/rs-engine/src/match.rs b/pepmatch/rs-engine/src/match.rs
index 9592cd4..c4bd3e7 100644
--- a/pepmatch/rs-engine/src/match.rs
+++ b/pepmatch/rs-engine/src/match.rs
@@ -543,10 +543,8 @@ fn minimal_coverage_seeds(query: &[u8], k: usize) -> Vec<(&[u8], usize)> {
     seeds
 }
 
-fn is_terminal_deletion(q_idx: isize, query_len: usize, p_idx: isize, protein_len: usize) -> bool {
-    if q_idx == 0 || q_idx == query_len as isize - 1 {
-        return true;
-    }
+fn is_terminal_deletion(p_idx: isize, protein_len: usize) -> bool {
+    // Protein boundary edge effect — no biological indel can be inferred here.
     if p_idx <= 0 || p_idx >= protein_len as isize - 1 {
         return true;
     }
@@ -568,7 +566,6 @@ fn dfs(
     }
 
     let mut all_paths: Vec<usize> = Vec::new();
-    let query_len = query.len();
     let protein_len = protein.len();
 
     // Match branch: both pointers advance, consumes 1 protein char.
@@ -583,7 +580,7 @@ fn dfs(
 
     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 !is_terminal_deletion(p_idx, protein_len) {
             all_paths.extend(dfs(
                 query, q_idx + direction, protein, p_idx, indels_left - 1, direction,
             ));
diff --git a/pepmatch/tests/test_indel_search.py b/pepmatch/tests/test_indel_search.py
index 0eafa02..d00cb7b 100644
--- a/pepmatch/tests/test_indel_search.py
+++ b/pepmatch/tests/test_indel_search.py
@@ -39,15 +39,16 @@ def test_indel_search(proteome_path, query_path, expected_path):
   )
 
 
-def test_terminal_deletion_blocked(proteome_path):
+def test_query_terminal_deletion_found(proteome_path):
+  # Deletion of Q at query position 0 hits NALVEATRFC at protein position 3
+  # in Q8V336 — not a protein boundary, so the match is valid.
   df = Matcher(
     query=['QNALVEATRFC'],
     proteome_file=proteome_path,
     max_indels=1,
     output_format='dataframe'
   ).match()
-  assert df['Matched Sequence'].is_null().all(), (
-    'Terminal deletion should be blocked but got match(es): '
-    + str(df.filter(pl.col('Matched Sequence').is_not_null())
-            .select(['Query Sequence', 'Matched Sequence']))
+  matches = df.filter(pl.col('Matched Sequence').is_not_null())
+  assert 'NALVEATRFC' in matches['Matched Sequence'].to_list(), (
+    'Expected query-terminal deletion match NALVEATRFC not found'
   )

From 1ca61d3d487d4f9819a2da07b7ee7b6d2aff1f9f Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung9981@gmail.com>
Date: Wed, 1 Jul 2026 11:48:44 -0700
Subject: [PATCH 05/15] refactor: move indel.py to tests/ as reference
 implementation

Not imported anywhere in the package; relocating per PR #24 review
(item 2) since its role is validating the Rust indel implementation,
not shipping with it.
---
 pepmatch/rs-engine/Cargo.lock | 2 +-
 pepmatch/{ => tests}/indel.py | 0
 2 files changed, 1 insertion(+), 1 deletion(-)
 rename pepmatch/{ => tests}/indel.py (100%)

diff --git a/pepmatch/rs-engine/Cargo.lock b/pepmatch/rs-engine/Cargo.lock
index 01e4263..33a2a50 100644
--- a/pepmatch/rs-engine/Cargo.lock
+++ b/pepmatch/rs-engine/Cargo.lock
@@ -92,7 +92,7 @@ checksum = "42f5e15c9953c5e4ccceeb2e7382a716482c34515315f7b03532b8b4e8393d2d"
 
 [[package]]
 name = "pepmatch_rs_engine"
-version = "1.16.3"
+version = "1.16.4"
 dependencies = [
  "memmap2",
  "pyo3",
diff --git a/pepmatch/indel.py b/pepmatch/tests/indel.py
similarity index 100%
rename from pepmatch/indel.py
rename to pepmatch/tests/indel.py

From e0b45ec75690daf642924b83f482fcc72ec4eb5f Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung9981@gmail.com>
Date: Wed, 1 Jul 2026 11:56:04 -0700
Subject: [PATCH 06/15] feat: add dedicated Indels column to output schema

Previously the indel count rode in the Mismatches slot. Split it into
its own column (right after Mismatches in FINAL_COLUMNS) so indel rows
report Mismatches=0/Indels=N and non-indel rows report Indels=0,
without disturbing existing Mismatches semantics for exact/mismatch/
discontinuous search. Per PR #24 review item 3.
---
 pepmatch/matcher.py | 18 +++++++++++++-----
 1 file changed, 13 insertions(+), 5 deletions(-)

diff --git a/pepmatch/matcher.py b/pepmatch/matcher.py
index ddd52b6..79ab4c2 100755
--- a/pepmatch/matcher.py
+++ b/pepmatch/matcher.py
@@ -176,7 +176,7 @@ def indel_search(self):
     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)
+    return self._to_dataframe(results, is_indels=True)
 
   def _pepidx_path(self, k):
     return os.path.join(
@@ -352,28 +352,36 @@ def _metadata_table(self) -> pl.DataFrame:
 
   FINAL_COLUMNS = [
     'Query ID','Query Sequence','Matched Sequence','Protein ID','Protein Name','Species',
-    'Taxon ID','Gene','Mismatches','Mutated Positions','Index start','Index end',
+    'Taxon ID','Gene','Mismatches','Indels','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
 
     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'):
+      for c in ('Mismatches','Indels','Index start','Index end','Protein Existence Level','Gene Priority'):
         schema[c] = pl.Int64
       schema['SwissProt Reviewed'] = pl.Boolean
       return pl.DataFrame(schema=schema)
 
+    # rs_indel_match packs the indel count into the same slot rs_match/rs_discontinuous
+    # use for mismatches; split it out here so Mismatches and Indels are never both set.
+    if is_indels:
+      indels, mismatches = mm, [0 if v is not None else None for v in mm]
+    else:
+      mismatches, indels = mm, [0 if v is not None else None for v in mm]
+
     base = pl.DataFrame({
       'Query ID': qid,
       'Query Sequence': qseq,
       'Matched Sequence': matched,
       'protein_num': pl.Series(pnum, dtype=pl.UInt32),
-      'Mismatches': pl.Series(mm, dtype=pl.Int64),
+      'Mismatches': pl.Series(mismatches, dtype=pl.Int64),
+      'Indels': pl.Series(indels, dtype=pl.Int64),
       'Mutated Positions': mutated,
       'Index start': pl.Series(istart, dtype=pl.Int64),
       'Index end': pl.Series(iend, dtype=pl.Int64),

From eb8533c332d5c388e53935fb5a524b767ccba1b3 Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung9981@gmail.com>
Date: Wed, 1 Jul 2026 12:47:32 -0700
Subject: [PATCH 07/15] fix: off-by-one in terminal deletion guard blocked
 valid matches

is_terminal_deletion treated p_idx == 0 and p_idx == protein_len - 1 as
boundary artifacts, but both are real, existing residues, not
out-of-bounds. This required 2 residues of protein context around a
deletion's gap when 1 is sufficient (a real adjacent residue already
confirms the query residue is genuinely absent, not just missing
because the recorded sequence happens to end there). Only strictly
out-of-bounds references (p_idx < 0 or p_idx >= protein_len) are
actually ambiguous.
---
 pepmatch/rs-engine/src/match.rs     |  5 +++--
 pepmatch/tests/test_indel_search.py | 18 ++++++++++++++++++
 2 files changed, 21 insertions(+), 2 deletions(-)

diff --git a/pepmatch/rs-engine/src/match.rs b/pepmatch/rs-engine/src/match.rs
index c4bd3e7..4295178 100644
--- a/pepmatch/rs-engine/src/match.rs
+++ b/pepmatch/rs-engine/src/match.rs
@@ -544,8 +544,9 @@ fn minimal_coverage_seeds(query: &[u8], k: usize) -> Vec<(&[u8], usize)> {
 }
 
 fn is_terminal_deletion(p_idx: isize, protein_len: usize) -> bool {
-    // Protein boundary edge effect — no biological indel can be inferred here.
-    if p_idx <= 0 || p_idx >= protein_len as isize - 1 {
+    // p_idx == 0 and p_idx == protein_len - 1 are real, existing residues — only
+    // an out-of-bounds reference is a genuine protein boundary edge effect.
+    if p_idx < 0 || p_idx >= protein_len as isize {
         return true;
     }
     false
diff --git a/pepmatch/tests/test_indel_search.py b/pepmatch/tests/test_indel_search.py
index d00cb7b..962bfb8 100644
--- a/pepmatch/tests/test_indel_search.py
+++ b/pepmatch/tests/test_indel_search.py
@@ -39,6 +39,24 @@ def test_indel_search(proteome_path, query_path, expected_path):
   )
 
 
+def test_terminal_deletion_allowed_with_single_residue_buffer(tmp_path):
+  # Regression test: p_idx == 0 and p_idx == protein_len - 1 are real,
+  # existing residues, not out-of-bounds — a deletion's gap only needs 1
+  # residue of protein context on either side to be a genuine indel, not 2.
+  proteome_path = tmp_path / 'proteome.fasta'
+  proteome_path.write_text('>Front\nXBCDE\n>Back\nABCDX\n')
+  df = Matcher(
+    query=['ABCDE'],
+    proteome_file=str(proteome_path),
+    max_indels=1,
+    preprocessed_files_path=str(tmp_path),
+    output_format='dataframe'
+  ).match()
+  hits = set(zip(df['Protein ID'].to_list(), df['Matched Sequence'].to_list()))
+  assert ('Front.1', 'BCDE') in hits
+  assert ('Back.1', 'ABCD') in hits
+
+
 def test_query_terminal_deletion_found(proteome_path):
   # Deletion of Q at query position 0 hits NALVEATRFC at protein position 3
   # in Q8V336 — not a protein boundary, so the match is valid.

From 4a78675194d01bdd9416a408016f2cce1bf8c4dc Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung9981@gmail.com>
Date: Wed, 1 Jul 2026 12:52:34 -0700
Subject: [PATCH 08/15] test: strengthen indel test suite per PR #24 review
 item 4

- Replace indel.py's seed/DFS port with a genuinely independent
  brute-force oracle (direct window enumeration, no seeding/recursion)
  so it can catch bugs the DFS and its Python port would share
- Add ValueError path tests for max_indels>1 and indel+mismatch combos
- Add peptide_len < k edge case test
- Add multi-hit test covering two different proteins in one query
- Add randomized property test comparing Rust output against the
  brute-force oracle across 25 seeded random peptide/proteome pairs
---
 pepmatch/tests/indel.py               | 143 +++++++-------------------
 pepmatch/tests/test_indel_property.py |  74 +++++++++++++
 pepmatch/tests/test_indel_search.py   |  47 +++++++++
 3 files changed, 160 insertions(+), 104 deletions(-)
 create mode 100644 pepmatch/tests/test_indel_property.py

diff --git a/pepmatch/tests/indel.py b/pepmatch/tests/indel.py
index 71476ee..50252b6 100644
--- a/pepmatch/tests/indel.py
+++ b/pepmatch/tests/indel.py
@@ -1,111 +1,46 @@
-def _is_terminal_deletion(p_idx, protein_len):
-  """Return True if a deletion at this protein position is at a sequence boundary.
-
-  Deletions at p_idx == 0 or p_idx == protein_len - 1 are forbidden because
-  they arise from a database boundary edge effect rather than a true biological
-  indel event.
-  """
-  if p_idx <= 0 or p_idx >= protein_len - 1:
-    return True
-  return False
-
-
-def _dfs(query, q_idx, protein, p_idx, indels_left, direction):
-  """Exhaustive depth-first search extending one direction from a seed hit.
-
-  Explores three branches at each step: a match branch (both pointers
-  advance), a deletion branch (query pointer advances, protein stays,
-  consumes 0 protein chars), and an insertion branch (protein pointer
-  advances, query stays, consumes 1 protein char). Returns a list of
-  integers, each the number of protein characters consumed by one valid
-  complete alignment path.
-
-  Args:
-    query: the full query peptide string.
-    q_idx: current query position, advances by direction.
-    protein: the full protein sequence string.
-    p_idx: current protein position, advances by direction.
-    indels_left: remaining indel budget for this direction.
-    direction: +1 for right extension, -1 for left extension.
+def _terminal_deletion_blocked(gap_left, gap_right, protein_len):
+  """A deletion at query position d found at protein offset `start` has its gap
+  sitting between protein indices (start+d-1) and (start+d). Both neighbors are
+  real, existing residues as long as they fall within [0, protein_len) — only
+  an out-of-bounds reference means the gap is indistinguishable from the
+  sequence simply ending there rather than a genuine biological indel.
   """
-  if (direction == 1 and q_idx >= len(query)) or \
-     (direction == -1 and q_idx < 0):
-    return [0]
-
-  all_paths = []
-
-  if 0 <= p_idx < len(protein) and query[q_idx] == protein[p_idx]:
-    for consumed in _dfs(query, q_idx + direction, protein, p_idx + direction, indels_left, direction):
-      all_paths.append(consumed + 1)
-
-  if indels_left > 0:
-    if not _is_terminal_deletion(p_idx, len(protein)):
-      all_paths.extend(
-        _dfs(query, q_idx + direction, protein, p_idx, indels_left - 1, direction)
-      )
-
-    if 0 <= p_idx < len(protein):
-      for consumed in _dfs(query, q_idx, protein, p_idx + direction, indels_left - 1, direction):
-        all_paths.append(consumed + 1)
-
-  return all_paths
+  return gap_left < 0 or gap_right >= protein_len
 
 
-def extend_bidirectional(query, q_seed_start, p_hit_idx, protein, k, max_indels=1):
-  """Verify and enumerate all indel-consistent matches from a seed hit.
+def brute_force_search(query, protein):
+  """Enumerate every (start, matched_sequence) pair where query aligns to a
+  window of protein using at most one single-residue insertion or deletion.
 
-  For each allocation of the indel budget between left and right extensions,
-  runs _dfs in both directions and combines results. Returns all unique
-  (start_0based, matched_sequence) pairs found at this seed hit position.
-
-  Args:
-    query: the full query peptide string.
-    q_seed_start: 0-based start of the seed k-mer within the query.
-    p_hit_idx: 0-based start of the seed hit within protein.
-    protein: the full protein sequence string.
-    k: seed k-mer length.
-    max_indels: maximum total indels allowed across both extensions.
+  Checks each candidate directly against the definition of a valid 1-indel
+  match (no seeding, no recursive extension) so it can serve as an
+  independent oracle for the Rust/seed-based implementation.
   """
-  seen = set()
-  results = []
-
-  for r_budget in range(max_indels + 1):
-    l_budget = max_indels - r_budget
-
-    r_paths = _dfs(query, q_seed_start + k, protein, p_hit_idx + k, r_budget, direction=1)
-    l_paths = _dfs(query, q_seed_start - 1, protein, p_hit_idx - 1, l_budget, direction=-1)
-
-    for r_consumed in r_paths:
-      for l_consumed in l_paths:
-        start = p_hit_idx - l_consumed
-        end = p_hit_idx + k + r_consumed
-        matched_seq = protein[start:end]
-        key = (start, matched_seq)
-        if key not in seen:
-          seen.add(key)
-          results.append(key)
+  query_len = len(query)
+  protein_len = len(protein)
+  results = set()
+
+  for start in range(protein_len - query_len + 1):
+    if protein[start:start + query_len] == query:
+      results.add((start, protein[start:start + query_len]))
+
+  window_len = query_len - 1
+  for d in range(query_len):
+    shortened = query[:d] + query[d + 1:]
+    for start in range(protein_len - window_len + 1):
+      if protein[start:start + window_len] == shortened:
+        if _terminal_deletion_blocked(start + d - 1, start + d, protein_len):
+          continue
+        results.add((start, shortened))
+
+  # i in [1, query_len) keeps the inserted residue interior to the alignment;
+  # an insertion at i=0 or i=query_len would just pad an exact match with an
+  # arbitrary flanking residue and isn't a real edit.
+  window_len = query_len + 1
+  for i in range(1, query_len):
+    for start in range(protein_len - window_len + 1):
+      window = protein[start:start + window_len]
+      if window[:i] == query[:i] and window[i + 1:] == query[i:]:
+        results.add((start, window))
 
   return results
-
-
-def minimal_coverage_seeds(query, k):
-  """Return non-overlapping k-mer seeds plus a final overlapping seed for
-  full query coverage.
-
-  The pigeonhole principle guarantees that for a query of length L with at
-  most d indels, at least one of floor(L / (d+1)) non-overlapping seeds
-  must hit exactly. Appends one final seed anchored at the query end if
-  the last strided seed does not reach it.
-
-  Args:
-    query: the query peptide string.
-    k: seed k-mer length.
-  """
-  seeds = []
-  query_len = len(query)
-  for j in range(0, query_len - k + 1, k):
-    seeds.append((query[j:j + k], j))
-  last_start = query_len - k
-  if seeds and seeds[-1][1] != last_start:
-    seeds.append((query[last_start:], last_start))
-  return seeds
\ No newline at end of file
diff --git a/pepmatch/tests/test_indel_property.py b/pepmatch/tests/test_indel_property.py
new file mode 100644
index 0000000..831a40b
--- /dev/null
+++ b/pepmatch/tests/test_indel_property.py
@@ -0,0 +1,74 @@
+import random
+import pytest
+from pepmatch import Matcher
+from indel import brute_force_search
+
+AMINO_ACIDS = 'ACDEFGHIKLMNPQRSTVWY'
+
+
+def _random_sequence(rng, length):
+  return ''.join(rng.choice(AMINO_ACIDS) for _ in range(length))
+
+
+def _with_deletion(rng, seq):
+  d = rng.randrange(len(seq))
+  return seq[:d] + seq[d + 1:]
+
+
+def _with_insertion(rng, seq):
+  i = rng.randrange(1, len(seq))
+  return seq[:i] + rng.choice(AMINO_ACIDS) + seq[i:]
+
+
+def _random_query(rng, proteins):
+  protein_id = rng.choice(list(proteins))
+  protein = proteins[protein_id]
+  qlen = rng.randint(8, 14)
+  start = rng.randrange(0, len(protein) - qlen + 1)
+  base = protein[start:start + qlen]
+  mode = rng.choice(['exact', 'deletion', 'insertion', 'random'])
+  if mode == 'deletion':
+    return _with_deletion(rng, base)
+  if mode == 'insertion':
+    return _with_insertion(rng, base)
+  if mode == 'random':
+    return _random_sequence(rng, qlen)
+  return base
+
+
+@pytest.mark.parametrize('seed', range(25))
+def test_indel_search_matches_brute_force_oracle(tmp_path, seed):
+  rng = random.Random(seed)
+
+  proteins = {f'P{i}': _random_sequence(rng, rng.randint(20, 40)) for i in range(4)}
+  proteome_path = tmp_path / 'proteome.fasta'
+  proteome_path.write_text(
+    ''.join(f'>{pid}\n{seq}\n' for pid, seq in proteins.items())
+  )
+
+  queries = [(f'q{i}', _random_query(rng, proteins)) for i in range(8)]
+  query_path = tmp_path / 'queries.fasta'
+  query_path.write_text(
+    ''.join(f'>{qid}\n{seq}\n' for qid, seq in queries)
+  )
+
+  df = Matcher(
+    query=str(query_path),
+    proteome_file=str(proteome_path),
+    max_indels=1,
+    preprocessed_files_path=str(tmp_path),
+    output_format='dataframe'
+  ).match()
+
+  actual = set()
+  for row in df.iter_rows(named=True):
+    if row['Matched Sequence'] is not None:
+      actual.add((row['Query Sequence'], row['Protein ID'], row['Matched Sequence']))
+
+  expected = set()
+  for _, qseq in queries:
+    for pid, pseq in proteins.items():
+      for _, matched in brute_force_search(qseq, pseq):
+        expected.add((qseq, f'{pid}.1', matched))
+
+  assert actual == expected
diff --git a/pepmatch/tests/test_indel_search.py b/pepmatch/tests/test_indel_search.py
index 962bfb8..de97d85 100644
--- a/pepmatch/tests/test_indel_search.py
+++ b/pepmatch/tests/test_indel_search.py
@@ -70,3 +70,50 @@ def test_query_terminal_deletion_found(proteome_path):
   assert 'NALVEATRFC' in matches['Matched Sequence'].to_list(), (
     'Expected query-terminal deletion match NALVEATRFC not found'
   )
+
+
+def test_max_indels_greater_than_one_raises():
+  with pytest.raises(ValueError, match='max_indels > 1 is not yet supported'):
+    Matcher(query=['NALVEATRFC'], proteome_file='unused.fasta', max_indels=2)
+
+
+def test_indels_and_mismatches_mutually_exclusive_raises():
+  with pytest.raises(ValueError, match='mutually exclusive'):
+    Matcher(
+      query=['NALVEATRFC'], proteome_file='unused.fasta',
+      max_indels=1, max_mismatches=1
+    )
+
+
+def test_indel_peptide_shorter_than_k(proteome_path):
+  # A single-residue query forces k = max(2, min_len // 2) = 2 while
+  # peptide_len=1 < k; must miss cleanly rather than erroring.
+  df = Matcher(
+    query=['A'],
+    proteome_file=proteome_path,
+    max_indels=1,
+    output_format='dataframe'
+  ).match()
+  assert df['Matched Sequence'].is_null().all()
+
+
+def test_indel_multi_hit_different_proteins(tmp_path):
+  # NALVEATRFC (10 residues) with 1 indel matches two DIFFERENT proteins:
+  # ProtDel is missing the second A (a deletion), ProtIns has an extra X
+  # inserted after E (an insertion). Both hand-verified against the DFS's
+  # seed ("NALVE") + bidirectional extension.
+  proteome_path = tmp_path / 'proteome.fasta'
+  proteome_path.write_text(
+    '>ProtDel\nMKVNALVETRFCGHI\n'
+    '>ProtIns\nMKVNALVEXATRFCGHI\n'
+  )
+  df = Matcher(
+    query=['NALVEATRFC'],
+    proteome_file=str(proteome_path),
+    max_indels=1,
+    preprocessed_files_path=str(tmp_path),
+    output_format='dataframe'
+  ).match()
+  hits = set(zip(df['Protein ID'].to_list(), df['Matched Sequence'].to_list()))
+  assert ('ProtDel.1', 'NALVETRFC') in hits
+  assert ('ProtIns.1', 'NALVEXATRFC') in hits

From 5edd3307fd8f63858d7b6e14cdd9d167414d7136 Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung9981@gmail.com>
Date: Wed, 1 Jul 2026 12:59:56 -0700
Subject: [PATCH 09/15] fix: make indel_search's index-build message explicit

Silently kicking off a full-proteome preprocessing pass with a terse
message is a surprise for large proteomes. Per PR #24 review item 5.
---
 pepmatch/matcher.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pepmatch/matcher.py b/pepmatch/matcher.py
index 79ab4c2..be7a26c 100755
--- a/pepmatch/matcher.py
+++ b/pepmatch/matcher.py
@@ -171,7 +171,8 @@ def indel_search(self):
     k = max(2, min_len // (self.max_indels + 1))
     pepidx_path = self._pepidx_path(k)
     if not os.path.isfile(pepidx_path):
-      print(f"Preprocessing {self.proteome_name} with k={k}...")
+      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})...")

From 190125b8dc424e32a95c575133ab0cc2509698c0 Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung9981@gmail.com>
Date: Wed, 1 Jul 2026 13:05:00 -0700
Subject: [PATCH 10/15] docs: add --max_indels CLI flag and README section for
 indel mode

Per PR #24 review item 6.
---
 README.md         | 12 ++++++++++++
 pepmatch/shell.py |  2 ++
 2 files changed, 14 insertions(+)

diff --git a/README.md b/README.md
index b44a12f..5a27f30 100644
--- a/README.md
+++ b/README.md
@@ -93,6 +93,17 @@ 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 need to preprocess or specify `k` yourself. Currently limited to `max_indels=1`, and mutually exclusive with `max_mismatches`.
+```python
+df = Matcher(
+  query='neoepitopes.fasta',
+  proteome_file='human.fasta',
+  max_indels=1
+).match()
+```
+
 #### 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 +174,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.
diff --git a/pepmatch/shell.py b/pepmatch/shell.py
index 9cd57c7..ef96748 100644
--- a/pepmatch/shell.py
+++ b/pepmatch/shell.py
@@ -21,6 +21,7 @@ def run_matcher():
   parser.add_argument('-q', '--query', required=True)
   parser.add_argument('-p', '--proteome_file', required=True)
   parser.add_argument('-m', '--max_mismatches', type=int, default=0)
+  parser.add_argument('-i', '--max_indels', type=int, default=0)
   parser.add_argument('-k', '--kmer_size', type=int, default=5)
   parser.add_argument('-P', '--preprocessed_files_path', default='.')
   parser.add_argument('-b', '--best_match', action='store_true', default=False)
@@ -33,6 +34,7 @@ def run_matcher():
     query=args.query,
     proteome_file=args.proteome_file,
     max_mismatches=args.max_mismatches,
+    max_indels=args.max_indels,
     k=args.kmer_size,
     preprocessed_files_path=args.preprocessed_files_path,
     best_match=args.best_match,

From 7bd7092ab361806d899c16551246ece0f1335a07 Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung9981@gmail.com>
Date: Wed, 1 Jul 2026 13:10:16 -0700
Subject: [PATCH 11/15] docs: match indel section's caveat placement to house
 style

Trailing constraint note (after the code block) rather than folded
into the intro sentence, mirroring how Counts-Only Mode documents its
best_match incompatibility.
---
 README.md | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 5a27f30..26ce09e 100644
--- a/README.md
+++ b/README.md
@@ -95,7 +95,7 @@ df = Matcher(
 
 #### Indel Searching
 
-Search allowing insertions and deletions (indels) instead of substitution mismatches. The k-mer size is chosen automatically from your query lengths — no need to preprocess or specify `k` yourself. Currently limited to `max_indels=1`, and mutually exclusive with `max_mismatches`.
+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',
@@ -103,6 +103,7 @@ df = Matcher(
   max_indels=1
 ).match()
 ```
+Currently limited to `max_indels=1`; mutually exclusive with `max_mismatches`.
 
 #### Best Match
 

From 63cbb92cc096e467f797346b7d7312ca8af29353 Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung9981@gmail.com>
Date: Wed, 1 Jul 2026 15:15:52 -0700
Subject: [PATCH 12/15] fix: raise ValueError for max_indels combined with
 best_match or counts_only

Both combinations previously produced silent wrong behavior: best_match
was ignored (indel_search() results were never passed through
_best_match_filter), and counts_only silently returned exact-match
counts (max_mismatches=0) rather than indel counts. Neither is
implemented yet, so both now fail loudly like the existing
max_indels/max_mismatches mutual-exclusivity check.
---
 pepmatch/matcher.py                 |  6 ++++++
 pepmatch/tests/test_indel_search.py | 16 ++++++++++++++++
 2 files changed, 22 insertions(+)

diff --git a/pepmatch/matcher.py b/pepmatch/matcher.py
index be7a26c..e641744 100755
--- a/pepmatch/matcher.py
+++ b/pepmatch/matcher.py
@@ -60,6 +60,12 @@ def __init__(
     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:
diff --git a/pepmatch/tests/test_indel_search.py b/pepmatch/tests/test_indel_search.py
index de97d85..0cd03f9 100644
--- a/pepmatch/tests/test_indel_search.py
+++ b/pepmatch/tests/test_indel_search.py
@@ -85,6 +85,22 @@ def test_indels_and_mismatches_mutually_exclusive_raises():
     )
 
 
+def test_max_indels_and_best_match_raises():
+  with pytest.raises(ValueError, match='not yet supported together'):
+    Matcher(
+      query=['NALVEATRFC'], proteome_file='unused.fasta',
+      max_indels=1, best_match=True
+    )
+
+
+def test_max_indels_and_counts_only_raises():
+  with pytest.raises(ValueError, match='not yet supported together'):
+    Matcher(
+      query=['NALVEATRFC'], proteome_file='unused.fasta',
+      max_indels=1, counts_only=True
+    )
+
+
 def test_indel_peptide_shorter_than_k(proteome_path):
   # A single-residue query forces k = max(2, min_len // 2) = 2 while
   # peptide_len=1 < k; must miss cleanly rather than erroring.

From 231bcb2540c96323499c1d8510cfec35ba45a195 Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung9981@gmail.com>
Date: Thu, 2 Jul 2026 10:47:40 -0700
Subject: [PATCH 13/15] docs: clarify indel property test's mutation-helper
 naming

_with_deletion/_with_insertion are named for the edit applied to the
query string, which is the opposite of the Indels label the search
engine reports for the resulting match (removing a query residue
means the protein has extra content there, an insertion match, and
vice versa). Noted during audit to prevent misreading.
---
 pepmatch/tests/test_indel_property.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/pepmatch/tests/test_indel_property.py b/pepmatch/tests/test_indel_property.py
index 831a40b..2b4e1e0 100644
--- a/pepmatch/tests/test_indel_property.py
+++ b/pepmatch/tests/test_indel_property.py
@@ -10,6 +10,9 @@ def _random_sequence(rng, length):
   return ''.join(rng.choice(AMINO_ACIDS) for _ in range(length))
 
 
+# Named for the edit applied to the query string, not the resulting Indels label —
+# e.g. _with_deletion removes a query residue, so the protein has "extra" content
+# there relative to the query, which the search reports as an insertion match.
 def _with_deletion(rng, seq):
   d = rng.randrange(len(seq))
   return seq[:d] + seq[d + 1:]

From b8f349db726669d63f40a0dfd21d77f630429114 Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung9981@gmail.com>
Date: Thu, 2 Jul 2026 15:09:00 -0700
Subject: [PATCH 14/15] fix: block query-terminal deletions and insertions
 unconditionally

A deletion at the query's own first or last residue has no query-side
context to confirm the residue is genuinely absent, as opposed to the
alignment simply starting/ending one residue short -- so it's now
blocked regardless of protein-side buffer, reverting the earlier
relaxation. Applied identically to match.rs's is_terminal_deletion and
the brute-force oracle, so both stay in lockstep.

Added test_terminal_insertion_blocked and
test_insertion_at_second_to_last_position_found to cover the analogous
insertion case: a padded leading/trailing insertion is unverifiable
and blocked, while an interior insertion one position before the
query's own boundary (e.g. LPDGVWEESS) has real query-side context on
both sides and remains valid.

Renamed indel.py to indel_brute_force.py for clarity ahead of future
2-indel brute-force work.
---
 pepmatch/rs-engine/src/match.rs               | 10 ++-
 .../tests/{indel.py => indel_brute_force.py}  | 14 ++--
 pepmatch/tests/test_indel_property.py         |  4 +-
 pepmatch/tests/test_indel_search.py           | 65 +++++++++++++++----
 4 files changed, 68 insertions(+), 25 deletions(-)
 rename pepmatch/tests/{indel.py => indel_brute_force.py} (73%)

diff --git a/pepmatch/rs-engine/src/match.rs b/pepmatch/rs-engine/src/match.rs
index 4295178..8c830db 100644
--- a/pepmatch/rs-engine/src/match.rs
+++ b/pepmatch/rs-engine/src/match.rs
@@ -543,7 +543,12 @@ fn minimal_coverage_seeds(query: &[u8], k: usize) -> Vec<(&[u8], usize)> {
     seeds
 }
 
-fn is_terminal_deletion(p_idx: isize, protein_len: usize) -> bool {
+fn is_terminal_deletion(q_idx: isize, query_len: usize, p_idx: isize, protein_len: usize) -> bool {
+    // A deletion at the query's own boundary has no query-side context to confirm
+    // it's genuinely absent, not just the alignment starting/ending short.
+    if q_idx == 0 || q_idx == query_len as isize - 1 {
+        return true;
+    }
     // p_idx == 0 and p_idx == protein_len - 1 are real, existing residues — only
     // an out-of-bounds reference is a genuine protein boundary edge effect.
     if p_idx < 0 || p_idx >= protein_len as isize {
@@ -567,6 +572,7 @@ fn dfs(
     }
 
     let mut all_paths: Vec<usize> = Vec::new();
+    let query_len = query.len();
     let protein_len = protein.len();
 
     // Match branch: both pointers advance, consumes 1 protein char.
@@ -581,7 +587,7 @@ fn dfs(
 
     if indels_left > 0 {
         // Deletion branch: query pointer advances, protein stays, 0 protein chars consumed.
-        if !is_terminal_deletion(p_idx, protein_len) {
+        if !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,
             ));
diff --git a/pepmatch/tests/indel.py b/pepmatch/tests/indel_brute_force.py
similarity index 73%
rename from pepmatch/tests/indel.py
rename to pepmatch/tests/indel_brute_force.py
index 50252b6..2f4c0d9 100644
--- a/pepmatch/tests/indel.py
+++ b/pepmatch/tests/indel_brute_force.py
@@ -1,10 +1,10 @@
-def _terminal_deletion_blocked(gap_left, gap_right, protein_len):
-  """A deletion at query position d found at protein offset `start` has its gap
-  sitting between protein indices (start+d-1) and (start+d). Both neighbors are
-  real, existing residues as long as they fall within [0, protein_len) — only
-  an out-of-bounds reference means the gap is indistinguishable from the
-  sequence simply ending there rather than a genuine biological indel.
+def _terminal_deletion_blocked(d, query_len, gap_left, gap_right, protein_len):
+  """Blocks a deletion at the query's own boundary (no query-side context to
+  confirm it's genuinely absent) or where its gap references an out-of-bounds
+  protein index (indistinguishable from the sequence simply ending there).
   """
+  if d == 0 or d == query_len - 1:
+    return True
   return gap_left < 0 or gap_right >= protein_len
 
 
@@ -29,7 +29,7 @@ def brute_force_search(query, protein):
     shortened = query[:d] + query[d + 1:]
     for start in range(protein_len - window_len + 1):
       if protein[start:start + window_len] == shortened:
-        if _terminal_deletion_blocked(start + d - 1, start + d, protein_len):
+        if _terminal_deletion_blocked(d, query_len, start + d - 1, start + d, protein_len):
           continue
         results.add((start, shortened))
 
diff --git a/pepmatch/tests/test_indel_property.py b/pepmatch/tests/test_indel_property.py
index 2b4e1e0..59a96f6 100644
--- a/pepmatch/tests/test_indel_property.py
+++ b/pepmatch/tests/test_indel_property.py
@@ -1,7 +1,7 @@
 import random
 import pytest
 from pepmatch import Matcher
-from indel import brute_force_search
+from indel_brute_force import brute_force_search
 
 AMINO_ACIDS = 'ACDEFGHIKLMNPQRSTVWY'
 
@@ -39,7 +39,7 @@ def _random_query(rng, proteins):
   return base
 
 
-@pytest.mark.parametrize('seed', range(25))
+@pytest.mark.parametrize('seed', range(15))
 def test_indel_search_matches_brute_force_oracle(tmp_path, seed):
   rng = random.Random(seed)
 
diff --git a/pepmatch/tests/test_indel_search.py b/pepmatch/tests/test_indel_search.py
index 0cd03f9..a1d90e5 100644
--- a/pepmatch/tests/test_indel_search.py
+++ b/pepmatch/tests/test_indel_search.py
@@ -4,9 +4,6 @@
 from pathlib import Path
 from pepmatch import Matcher
 
-# QNALVEATRFC is designed so the only match in the test proteome would require
-# deleting Q at query position 0 (a terminal deletion). The guard must block it.
-
 @pytest.fixture
 def proteome_path() -> Path:
   return Path(__file__).parent / 'data' / 'proteome.fasta'
@@ -40,11 +37,11 @@ def test_indel_search(proteome_path, query_path, expected_path):
 
 
 def test_terminal_deletion_allowed_with_single_residue_buffer(tmp_path):
-  # Regression test: p_idx == 0 and p_idx == protein_len - 1 are real,
-  # existing residues, not out-of-bounds — a deletion's gap only needs 1
-  # residue of protein context on either side to be a genuine indel, not 2.
+  # Regression test: p_idx==0/protein_len-1 are real residues, not out-of-bounds —
+  # 1 residue of protein buffer is enough, not 2. Uses interior query positions
+  # (d=1,3 of 'ABCDE') to isolate this from the query-terminal guard below.
   proteome_path = tmp_path / 'proteome.fasta'
-  proteome_path.write_text('>Front\nXBCDE\n>Back\nABCDX\n')
+  proteome_path.write_text('>Front\nACDEFG\n>Back\nXYABCE\n')
   df = Matcher(
     query=['ABCDE'],
     proteome_file=str(proteome_path),
@@ -53,13 +50,13 @@ def test_terminal_deletion_allowed_with_single_residue_buffer(tmp_path):
     output_format='dataframe'
   ).match()
   hits = set(zip(df['Protein ID'].to_list(), df['Matched Sequence'].to_list()))
-  assert ('Front.1', 'BCDE') in hits
-  assert ('Back.1', 'ABCD') in hits
+  assert ('Front.1', 'ACDE') in hits
+  assert ('Back.1', 'ABCE') in hits
 
 
-def test_query_terminal_deletion_found(proteome_path):
-  # Deletion of Q at query position 0 hits NALVEATRFC at protein position 3
-  # in Q8V336 — not a protein boundary, so the match is valid.
+def test_query_terminal_deletion_blocked(proteome_path):
+  # Deleting Q at query position 0 is the only way QNALVEATRFC could match
+  # NALVEATRFC in Q8V336 — blocked unconditionally, a query-boundary deletion.
   df = Matcher(
     query=['QNALVEATRFC'],
     proteome_file=proteome_path,
@@ -67,9 +64,49 @@ def test_query_terminal_deletion_found(proteome_path):
     output_format='dataframe'
   ).match()
   matches = df.filter(pl.col('Matched Sequence').is_not_null())
-  assert 'NALVEATRFC' in matches['Matched Sequence'].to_list(), (
-    'Expected query-terminal deletion match NALVEATRFC not found'
+  assert 'NALVEATRFC' not in matches['Matched Sequence'].to_list(), (
+    'Query-terminal deletion match NALVEATRFC should be blocked'
+  )
+
+
+def test_terminal_insertion_blocked(tmp_path):
+  # A leading/trailing insertion pads the query with no context to confirm
+  # it's genuine, not an arbitrary flanking residue — the exact match still hits.
+  proteome_path = tmp_path / 'proteome.fasta'
+  proteome_path.write_text(
+    '>Leading\nXNALVEATRFC\n'
+    '>Trailing\nNALVEATRFCX\n'
   )
+  df = Matcher(
+    query=['NALVEATRFC'],
+    proteome_file=str(proteome_path),
+    max_indels=1,
+    preprocessed_files_path=str(tmp_path),
+    output_format='dataframe'
+  ).match()
+  hits = set(zip(df['Protein ID'].to_list(), df['Matched Sequence'].to_list()))
+  assert ('Leading.1', 'NALVEATRFC') in hits
+  assert ('Trailing.1', 'NALVEATRFC') in hits
+  assert ('Leading.1', 'XNALVEATRFC') not in hits
+  assert ('Trailing.1', 'NALVEATRFCX') not in hits
+
+
+def test_insertion_at_second_to_last_position_found(tmp_path):
+  # The repeated S makes the inserted residue's exact position ambiguous, but
+  # the second-to-last placement is verifiable (query's final S still matches
+  # a real residue) and must be found — unlike a genuinely terminal insertion.
+  proteome_path = tmp_path / 'proteome.fasta'
+  proteome_path.write_text('>P\nLPDGVWEESS\n')
+  df = Matcher(
+    query=['LPDGVWEES'],
+    proteome_file=str(proteome_path),
+    max_indels=1,
+    preprocessed_files_path=str(tmp_path),
+    output_format='dataframe'
+  ).match()
+  hits = set(zip(df['Protein ID'].to_list(), df['Matched Sequence'].to_list()))
+  assert ('P.1', 'LPDGVWEES') in hits
+  assert ('P.1', 'LPDGVWEESS') in hits
 
 
 def test_max_indels_greater_than_one_raises():

From 089fc97cd98fd3855392b75aec456db819fa443c Mon Sep 17 00:00:00 2001
From: Roman Young <romanyoung9981@gmail.com>
Date: Wed, 8 Jul 2026 12:53:17 -0700
Subject: [PATCH 15/15] fix: emit one edit-count column per search mode; reject
 indel + discontinuous

Each mode now carries a single edit-count column instead of always showing
both: Indels for indel search, Mismatches for every other mode. The prior
schema surfaced an all-zero Indels column on mismatch runs (and an all-zero
Mismatches column on indel runs), which reads as a search that did not happen;
dropping the twin also restores the original schema for the mismatch, exact,
discontinuous, and best_match outputs. _to_dataframe now names the single
column by mode via _final_columns(is_indels), replacing the fixed
FINAL_COLUMNS list.

Since a single match() call can concat a linear result frame with a
discontinuous one, and the two must share columns, max_indels combined with
discontinuous epitopes is now rejected in __init__ -- a discontinuous indel
search is undefined anyway (a position-anchored epitope has no contiguous
window to seed or extend).

Adds test_indel_mode_emits_indels_column_only,
test_mismatch_mode_emits_mismatches_column_only, and
test_max_indels_and_discontinuous_raises.
---
 pepmatch/matcher.py                 | 39 ++++++++++++++++-------------
 pepmatch/tests/test_indel_search.py | 35 ++++++++++++++++++++++++++
 2 files changed, 57 insertions(+), 17 deletions(-)

diff --git a/pepmatch/matcher.py b/pepmatch/matcher.py
index e641744..bc2db4f 100755
--- a/pepmatch/matcher.py
+++ b/pepmatch/matcher.py
@@ -88,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.')
@@ -357,38 +359,41 @@ 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','Indels','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, 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','Indels','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)
 
-    # rs_indel_match packs the indel count into the same slot rs_match/rs_discontinuous
-    # use for mismatches; split it out here so Mismatches and Indels are never both set.
-    if is_indels:
-      indels, mismatches = mm, [0 if v is not None else None for v in mm]
-    else:
-      mismatches, indels = mm, [0 if v is not None else None for v in mm]
-
     base = pl.DataFrame({
       'Query ID': qid,
       'Query Sequence': qseq,
       'Matched Sequence': matched,
       'protein_num': pl.Series(pnum, dtype=pl.UInt32),
-      'Mismatches': pl.Series(mismatches, dtype=pl.Int64),
-      'Indels': pl.Series(indels, 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),
@@ -404,4 +409,4 @@ def _to_dataframe(self, cols, is_indels=False):
           .alias('Protein ID')
       )
 
-    return df.drop('Sequence Version').select(self.FINAL_COLUMNS)
+    return df.drop('Sequence Version').select(final_columns)
diff --git a/pepmatch/tests/test_indel_search.py b/pepmatch/tests/test_indel_search.py
index a1d90e5..20771c3 100644
--- a/pepmatch/tests/test_indel_search.py
+++ b/pepmatch/tests/test_indel_search.py
@@ -138,6 +138,15 @@ def test_max_indels_and_counts_only_raises():
     )
 
 
+def test_max_indels_and_discontinuous_raises():
+  # A discontinuous (position-anchored) epitope has no contiguous window for the
+  # indel engine to seed/extend, so the combination is undefined — it must error
+  # up front rather than crash later when the indel and discontinuous result
+  # frames (now carrying different edit-count columns) fail to concat.
+  with pytest.raises(ValueError, match='discontinuous'):
+    Matcher(query=['N1, A3, L5'], proteome_file='unused.fasta', max_indels=1)
+
+
 def test_indel_peptide_shorter_than_k(proteome_path):
   # A single-residue query forces k = max(2, min_len // 2) = 2 while
   # peptide_len=1 < k; must miss cleanly rather than erroring.
@@ -170,3 +179,29 @@ def test_indel_multi_hit_different_proteins(tmp_path):
   hits = set(zip(df['Protein ID'].to_list(), df['Matched Sequence'].to_list()))
   assert ('ProtDel.1', 'NALVETRFC') in hits
   assert ('ProtIns.1', 'NALVEXATRFC') in hits
+
+
+def test_indel_mode_emits_indels_column_only(proteome_path):
+  # One edit-count column per mode: an indel search reports counts in an Indels
+  # column and must NOT carry a separate always-zero Mismatches column.
+  df = Matcher(
+    query=['NALVEATRFC'],
+    proteome_file=proteome_path,
+    max_indels=1,
+    output_format='dataframe'
+  ).match()
+  assert 'Indels' in df.columns
+  assert 'Mismatches' not in df.columns
+
+
+def test_mismatch_mode_emits_mismatches_column_only(proteome_path):
+  # The mirror: a non-indel search keeps its Mismatches column and must NOT gain
+  # an always-zero Indels column suggesting an indel search that never ran.
+  df = Matcher(
+    query=['NALVEATRFC'],
+    proteome_file=proteome_path,
+    max_mismatches=0,
+    output_format='dataframe'
+  ).match()
+  assert 'Mismatches' in df.columns
+  assert 'Indels' not in df.columns