Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
59 changes: 54 additions & 5 deletions pepmatch/matcher.py
Original file line number Diff line number Diff line change
Expand Up @@ -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."""
Expand Down Expand Up @@ -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',
]

Expand All @@ -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:
Expand All @@ -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),
})
Expand Down
65 changes: 57 additions & 8 deletions pepmatch/tests/test_indel_search.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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]'
Loading