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/README.md b/README.md index b44a12f..26ce09e 100644 --- a/README.md +++ b/README.md @@ -93,6 +93,18 @@ df = Matcher( ).match() ``` +#### Indel Searching + +Search allowing insertions and deletions (indels) instead of substitution mismatches. The k-mer size is chosen automatically from your query lengths — no manual preprocessing needed. +```python +df = Matcher( + query='neoepitopes.fasta', + proteome_file='human.fasta', + max_indels=1 +).match() +``` +Currently limited to `max_indels=1`; mutually exclusive with `max_mismatches`. + #### Best Match Automatically finds the optimal match for each peptide by trying different k-mer sizes and mismatch thresholds. No manual preprocessing required. @@ -163,6 +175,7 @@ pepmatch-match -q peptides.fasta -p human.fasta -m 0 -k 5 * `-q`, `--query` (Required): Path to the query file. * `-p`, `--proteome_file` (Required): Path to the proteome FASTA file. * `-m`, `--max_mismatches`: Maximum mismatches allowed (default: 0). +* `-i`, `--max_indels`: Maximum indels allowed (default: 0). Currently limited to 1, and mutually exclusive with `-m`. * `-k`, `--kmer_size`: K-mer size (default: 5). * `-P`, `--preprocessed_files_path`: Directory containing preprocessed files. * `-b`, `--best_match`: Enable best match mode. diff --git a/pepmatch/matcher.py b/pepmatch/matcher.py index e44a23c..bc2db4f 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,20 @@ def __init__( if k != 0 and k < 2: raise ValueError('k must be >= 2.') + if max_indels > 1: + raise ValueError('max_indels > 1 is not yet supported. Only max_indels=1 has been validated.') + + if max_indels > 0 and max_mismatches > 0: + raise ValueError('max_indels and max_mismatches are mutually exclusive.') + + if max_indels > 0 and best_match: + raise ValueError('max_indels and best_match are not yet supported together.') + + if max_indels > 0 and counts_only: + raise ValueError('max_indels and counts_only are not yet supported together.') + + self.max_indels = max_indels + if output_format not in VALID_OUTPUT_FORMATS: raise ValueError( f'Invalid output format. Choose from: {VALID_OUTPUT_FORMATS}' @@ -73,6 +88,8 @@ def __init__( self.query = self._parse_query(query) self.discontinuous_epitopes = self._find_discontinuous_epitopes() + if self.max_indels > 0 and self.discontinuous_epitopes: + raise ValueError('max_indels and discontinuous epitopes are not supported together.') self.query = self._clean_query() if not self.query and not self.discontinuous_epitopes: raise ValueError('Query is empty.') @@ -128,7 +145,9 @@ def match(self): return output_matches(df, self.output_format, self.output_name) if self.query: - if self.best_match and self.k_specified: + if self.max_indels > 0: + linear_df = self.indel_search() + elif self.best_match and self.k_specified: results = self._search(self.k, self.max_mismatches) linear_df = self._best_match_filter(self._to_dataframe(results)) elif self.best_match: @@ -155,6 +174,19 @@ def match(self): else: output_matches(df, self.output_format, self.output_name) + def indel_search(self): + min_len = min(len(seq) for _, seq in self.query) + k = max(2, min_len // (self.max_indels + 1)) + pepidx_path = self._pepidx_path(k) + if not os.path.isfile(pepidx_path): + print(f"No preprocessed file found for k={k}, building index now " + f"(this may take a moment)...") + rs_preprocess(self.proteome_file, k, pepidx_path) + print(f"Searching {len(self.query)} peptides against {self.proteome_name} " + f"(k={k}, max_indels={self.max_indels})...") + results = rs_indel_match(pepidx_path, self.query, self.max_indels) + return self._to_dataframe(results, is_indels=True) + def _pepidx_path(self, k): return os.path.join( self.preprocessed_files_path, f'{self.proteome_name}_{k}mers.pepidx' @@ -327,20 +359,31 @@ def _metadata_table(self) -> pl.DataFrame: pl.col('SwissProt Reviewed').cast(pl.Int64).cast(pl.Boolean), ]) - FINAL_COLUMNS = [ - 'Query ID','Query Sequence','Matched Sequence','Protein ID','Protein Name','Species', - 'Taxon ID','Gene','Mismatches','Mutated Positions','Index start','Index end', - 'Protein Existence Level','Gene Priority','SwissProt Reviewed', - ] + def _final_columns(self, is_indels): + # One edit-count column per mode, in a fixed position: Indels for indel search, + # Mismatches for every other mode. Never both — an all-zero twin column would + # imply a search that didn't run. + edit_col = 'Indels' if is_indels else 'Mismatches' + return [ + 'Query ID','Query Sequence','Matched Sequence','Protein ID','Protein Name','Species', + 'Taxon ID','Gene', edit_col, 'Mutated Positions','Index start','Index end', + 'Protein Existence Level','Gene Priority','SwissProt Reviewed', + ] - def _to_dataframe(self, cols): + def _to_dataframe(self, cols, is_indels=False): """Build the result DataFrame from columnar Rust output, reconstructing protein metadata via a single join (instead of cloning it into every hit row).""" qid, qseq, matched, pnum, mm, mutated, istart, iend = cols + # rs_indel_match packs the indel count into the same slot rs_match/rs_discontinuous + # use for mismatches, so the values are identical in shape — only the column name + # differs by mode (Indels vs Mismatches), and only one is ever emitted. + edit_col = 'Indels' if is_indels else 'Mismatches' + final_columns = self._final_columns(is_indels) + if not qid: - schema = {c: pl.Utf8 for c in self.FINAL_COLUMNS} - for c in ('Mismatches','Index start','Index end','Protein Existence Level','Gene Priority'): + schema = {c: pl.Utf8 for c in final_columns} + for c in (edit_col, 'Index start', 'Index end', 'Protein Existence Level', 'Gene Priority'): schema[c] = pl.Int64 schema['SwissProt Reviewed'] = pl.Boolean return pl.DataFrame(schema=schema) @@ -350,7 +393,7 @@ def _to_dataframe(self, cols): 'Query Sequence': qseq, 'Matched Sequence': matched, 'protein_num': pl.Series(pnum, dtype=pl.UInt32), - 'Mismatches': pl.Series(mm, dtype=pl.Int64), + edit_col: pl.Series(mm, dtype=pl.Int64), 'Mutated Positions': mutated, 'Index start': pl.Series(istart, dtype=pl.Int64), 'Index end': pl.Series(iend, dtype=pl.Int64), @@ -366,4 +409,4 @@ def _to_dataframe(self, cols): .alias('Protein ID') ) - return df.drop('Sequence Version').select(self.FINAL_COLUMNS) + return df.drop('Sequence Version').select(final_columns) 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/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..8c830db 100644 --- a/pepmatch/rs-engine/src/match.rs +++ b/pepmatch/rs-engine/src/match.rs @@ -526,6 +526,217 @@ 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 { + // 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 { + return true; + } + false +} + +fn dfs( + query: &[u8], + q_idx: isize, + protein: &[u8], + p_idx: isize, + indels_left: usize, + direction: isize, +) -> Vec { + if (direction == 1 && q_idx >= query.len() as isize) + || (direction == -1 && q_idx < 0) + { + return vec![0]; + } + + let mut all_paths: Vec = 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)> { + let mut seen: HashSet<(usize, Vec)> = HashSet::new(); + let mut results: Vec<(usize, Vec)> = 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 { + 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)> = HashSet::new(); + let mut records: Vec = 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 = 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) 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, 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/indel_brute_force.py b/pepmatch/tests/indel_brute_force.py new file mode 100644 index 0000000..2f4c0d9 --- /dev/null +++ b/pepmatch/tests/indel_brute_force.py @@ -0,0 +1,46 @@ +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 + + +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. + + 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. + """ + 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(d, query_len, 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 diff --git a/pepmatch/tests/test_indel_property.py b/pepmatch/tests/test_indel_property.py new file mode 100644 index 0000000..59a96f6 --- /dev/null +++ b/pepmatch/tests/test_indel_property.py @@ -0,0 +1,77 @@ +import random +import pytest +from pepmatch import Matcher +from indel_brute_force import brute_force_search + +AMINO_ACIDS = 'ACDEFGHIKLMNPQRSTVWY' + + +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:] + + +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(15)) +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 new file mode 100644 index 0000000..20771c3 --- /dev/null +++ b/pepmatch/tests/test_indel_search.py @@ -0,0 +1,207 @@ +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) + ) + + +def test_terminal_deletion_allowed_with_single_residue_buffer(tmp_path): + # 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\nACDEFG\n>Back\nXYABCE\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', 'ACDE') in hits + assert ('Back.1', 'ABCE') in hits + + +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, + max_indels=1, + output_format='dataframe' + ).match() + matches = df.filter(pl.col('Matched Sequence').is_not_null()) + 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(): + 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_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_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. + 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 + + +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