feat: annotate indel hits with edit type and query position (Indel Positions column) (#27)

Author: Roman-Young

README.md

@@ -103,6 +103,8 @@ df = Matcher(
   max_indels=1
 ).match()
 ```
+Each indel hit is annotated in an **`Indel Positions`** column with the edit type, residue, and 1-based query position — e.g. `d: A[6]` (deletion of `A` at query position 6) or `i: X[6]` (insertion of `X` before query position 6); an exact match reports `[]`. In a repeat the exact position is ambiguous, so the inclusive range of all valid positions is reported, e.g. `d: A[2,4]`.
+
 Currently limited to `max_indels=1`; mutually exclusive with `max_mismatches`.
 
 #### Best Match

pepmatch/matcher.py

@@ -22,6 +22,44 @@ def output_matches(df: pl.DataFrame, output_format: str, output_name: str) -> No
   elif output_format == 'json':
     df.write_json(path)
 
+
+def _indel_edit(query, matched):
+  """Return (kind, residues, low, high) for a single-indel match, or None for an
+  exact match: kind 'd'/'i', the deleted query or inserted protein residue(s), and
+  [low, high] the inclusive 1-based range of valid positions (a run of equivalent
+  positions in a repeat)."""
+  if matched == query:
+    return None
+  L = len(query)
+  if len(matched) < L:
+    # interior positions only — query-terminal deletions are barred
+    positions = [i for i in range(2, L) if query[:i - 1] + query[i:] == matched]
+    if not positions:
+      return None
+    return ('d', query[positions[0] - 1], positions[0], positions[-1])
+  # interior matched residues only — boundary insertions are barred
+  positions, residue = [], None
+  for k in range(1, len(matched) - 1):
+    if matched[:k] + matched[k + 1:] == query:
+      positions.append(k + 1)   # 1-based query position the inserted residue precedes
+      residue = matched[k]
+  if not positions:
+    return None
+  return ('i', residue, positions[0], positions[-1])
+
+
+def format_indel_positions(query, matched):
+  """Render the Indel Positions column, e.g. 'd: A[6]', 'i: X[2,4]', or '[]' for
+  an exact match. Positions are 1-based; a range [low,high] collapses to [n] when
+  the position is unambiguous."""
+  edit = _indel_edit(query, matched)
+  if edit is None:
+    return '[]'
+  kind, residues, low, high = edit
+  span = f'[{low}]' if low == high else f'[{low},{high}]'
+  return f'{kind}: {residues}{span}'
+
+
 class Matcher:
   """Searches query peptides against a preprocessed proteome index
   and returns matches as a Polars DataFrame or output file."""
@@ -360,13 +398,15 @@ def _metadata_table(self) -> pl.DataFrame:
     ])
 
   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.
+    # One edit-count column and one edit-position column per mode, in fixed
+    # positions: Indels / Indel Positions for indel search, Mismatches / Mutated
+    # Positions for every other mode. Never both — an all-zero/empty twin column
+    # would imply a search that didn't run.
     edit_col = 'Indels' if is_indels else 'Mismatches'
+    pos_col = 'Indel Positions' if is_indels else 'Mutated Positions'
     return [
       'Query ID','Query Sequence','Matched Sequence','Protein ID','Protein Name','Species',
-      'Taxon ID','Gene', edit_col, 'Mutated Positions','Index start','Index end',
+      'Taxon ID','Gene', edit_col, pos_col,'Index start','Index end',
       'Protein Existence Level','Gene Priority','SwissProt Reviewed',
     ]
 
@@ -379,6 +419,7 @@ def _to_dataframe(self, cols, is_indels=False):
     # 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'
+    pos_col = 'Indel Positions' if is_indels else 'Mutated Positions'
     final_columns = self._final_columns(is_indels)
 
     if not qid:
@@ -388,13 +429,21 @@ def _to_dataframe(self, cols, is_indels=False):
       schema['SwissProt Reviewed'] = pl.Boolean
       return pl.DataFrame(schema=schema)
 
+    # In indel mode the edit position is derivable from (query, matched), so we
+    # compute Indel Positions here rather than in Rust; miss rows (no match) stay null.
+    if is_indels:
+      positions = [format_indel_positions(q, m) if m is not None else None
+                   for q, m in zip(qseq, matched)]
+    else:
+      positions = mutated
+
     base = pl.DataFrame({
       'Query ID': qid,
       'Query Sequence': qseq,
       'Matched Sequence': matched,
       'protein_num': pl.Series(pnum, dtype=pl.UInt32),
       edit_col: pl.Series(mm, dtype=pl.Int64),
-      'Mutated Positions': mutated,
+      pos_col: positions,
       'Index start': pl.Series(istart, dtype=pl.Int64),
       'Index end': pl.Series(iend, dtype=pl.Int64),
     })

pepmatch/tests/test_indel_search.py

@@ -3,6 +3,7 @@
 import polars.testing as plt
 from pathlib import Path
 from pepmatch import Matcher
+from pepmatch.matcher import format_indel_positions
 
 @pytest.fixture
 def proteome_path() -> Path:
@@ -182,26 +183,74 @@ def test_indel_multi_hit_different_proteins(tmp_path):
 
 
 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.
+  # One edit-count and one edit-position column per mode: an indel search reports
+  # Indels + Indel Positions and must NOT carry the mismatch-mode twins.
   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
+  assert 'Indels' in df.columns and 'Indel Positions' in df.columns
+  assert 'Mismatches' not in df.columns and 'Mutated Positions' 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.
+  # The mirror: a non-indel search keeps Mismatches + Mutated Positions and must
+  # NOT gain the indel-mode twins 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
+  assert 'Mismatches' in df.columns and 'Mutated Positions' in df.columns
+  assert 'Indels' not in df.columns and 'Indel Positions' not in df.columns
+
+
+def test_indel_positions_annotation_unit():
+  # Hand-verified annotations, format `<d|i>: <residues>[<pos or range>]`, 1-based.
+  # Deletion residue comes from the query, insertion residue from the protein. In a
+  # repeat the exact position is ambiguous, so a range of all valid positions is
+  # reported (collapsing to a single number when unambiguous).
+  assert format_indel_positions('ABCDEF', 'ABCDEF') == '[]'           # exact
+  assert format_indel_positions('YYADGY', 'YADGY') == 'd: Y[2]'       # only pos 2 (1 is terminal)
+  assert format_indel_positions('NALVEATRFC', 'NALVETRFC') == 'd: A[6]'  # the 2nd A, unambiguous
+  assert format_indel_positions('ABCDEF', 'ABXCDEF') == 'i: X[3]'     # X inserted before C
+  assert format_indel_positions('AAAAA', 'AAAA') == 'd: A[2,4]'       # deletable at 2, 3 or 4
+  assert format_indel_positions('AAAAAA', 'AAAAA') == 'd: A[2,5]'
+  assert format_indel_positions('AAAAAA', 'AAAAAAA') == 'i: A[2,6]'   # insertable across the run
+
+
+def test_indel_positions_deletion_end_to_end(tmp_path):
+  # Query NALVEATRFC vs a protein missing the 2nd A -> matched NALVETRFC,
+  # annotated as a deletion of A at query position 6.
+  proteome_path = tmp_path / 'proteome.fasta'
+  proteome_path.write_text('>P\nMKVNALVETRFCGHI\n')
+  df = Matcher(
+    query=['NALVEATRFC'],
+    proteome_file=str(proteome_path),
+    max_indels=1,
+    preprocessed_files_path=str(tmp_path),
+    output_format='dataframe'
+  ).match()
+  row = df.filter(pl.col('Matched Sequence') == 'NALVETRFC')
+  assert row.height == 1
+  assert row['Indel Positions'].item() == 'd: A[6]'
+
+
+def test_indel_positions_insertion_end_to_end(tmp_path):
+  # Query NALVEATRFC vs a protein with an extra X after E -> matched NALVEXATRFC,
+  # annotated as an insertion of X before query position 6.
+  proteome_path = tmp_path / 'proteome.fasta'
+  proteome_path.write_text('>P\nMKVNALVEXATRFCGHI\n')
+  df = Matcher(
+    query=['NALVEATRFC'],
+    proteome_file=str(proteome_path),
+    max_indels=1,
+    preprocessed_files_path=str(tmp_path),
+    output_format='dataframe'
+  ).match()
+  row = df.filter(pl.col('Matched Sequence') == 'NALVEXATRFC')
+  assert row.height == 1
+  assert row['Indel Positions'].item() == 'i: X[6]'