diff --git a/README.md b/README.md index 1a966c4..bec0279 100644 --- a/README.md +++ b/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 diff --git a/pepmatch/matcher.py b/pepmatch/matcher.py index bc2db4f..285f9b2 100755 --- a/pepmatch/matcher.py +++ b/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), }) diff --git a/pepmatch/tests/test_indel_search.py b/pepmatch/tests/test_indel_search.py index 20771c3..bc18a6e 100644 --- a/pepmatch/tests/test_indel_search.py +++ b/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 `: []`, 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]'