From 06679928976dde1bd7b8fe1dfa2f49932b4d5056 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Semid=C3=A1n=20Robaina=20Est=C3=A9vez?= Date: Fri, 19 Jun 2026 13:08:01 +0100 Subject: [PATCH] Add --best_hmm_wins to fix paralog cross-hits dropping syntenic blocks (#96) When the HMM panel contains paralogous models, the same peptide can score above threshold under more than one HMM. Pynteny kept every hit, so the duplicate rows interleaved (by contig/gene_pos sort) between true syntenic neighbours and broke the rolling-window matcher's position-diff check, silently dropping a genuine operon from synteny_matched.tsv. - Add an opt-in `--best_hmm_wins` flag (default off, preserving current behaviour). When enabled, SyntenyHMMfilter.get_all_HMM_hits deduplicates by peptide (full_label), keeping the highest-scoring HMM, before the rolling-window scan. bitscore is carried only on this path and dropped right after dedup, so the default code path is unchanged. - Thread the flag CLI (`-b/--best_hmm_wins`) -> synteny_search -> filter_FASTA_by_synteny_structure -> SyntenyHMMfilter, and through the api.Search class. - Upgrade the duplicate-hits warning to report the count and the top offending HMM combinations, and to point at --best_hmm_wins. - Docs: search.md gains a "note on paralogous HMMs" section and the new flag in the options table. - Tests: new tests/test_filter.py reproduces the silent drop on cross-hits and verifies --best_hmm_wins recovers the operon. Closes #96 Co-Authored-By: Claude Opus 4.8 --- docs/subcommands/search.md | 9 ++++ src/pynteny/api.py | 5 +++ src/pynteny/cli.py | 14 +++++++ src/pynteny/filter.py | 65 +++++++++++++++++++++++++--- src/pynteny/subcommands.py | 1 + tests/test_filter.py | 86 ++++++++++++++++++++++++++++++++++++++ 6 files changed, 175 insertions(+), 5 deletions(-) create mode 100644 tests/test_filter.py diff --git a/docs/subcommands/search.md b/docs/subcommands/search.md index d53ccb1..45cadfd 100644 --- a/docs/subcommands/search.md +++ b/docs/subcommands/search.md @@ -20,6 +20,7 @@ usage: pynteny search [-h] [args] |`-a`|`--hmmsearch_args`|`None`|list of comma-separated additional arguments to hmmsearch for each input hmm. A single argument may be provided, in which case the same additional argument is employed in all hmms.| |`-g`|`--gene_ids`||use gene symbols in synteny structure instead of HMM names. If set, a path to the hmm database metadata file must be provided in argument '--hmm_meta'| |`-u`|`--unordered`||whether the HMMs should be arranged in the exact same order displayed in the synteny_structure or in any order. If ordered, the filters will filter collinear rather than syntenic structures. If more than two HMMs are employed, the largest maximum distance among any pair is considered to run the search.| +|`-b`|`--best_hmm_wins`||if the same peptide is hit by more than one HMM (e.g. paralogous models that cross-hit), keep only the highest-scoring HMM for that peptide. Use this when your HMM panel contains close paralogs to avoid silently dropping genuine syntenic blocks. Disabled by default.| |`-r`|`--reuse`||reuse hmmsearch result table in following synteny searches. Do not delete hmmer_outputs subdirectory for this option to work.| |`-m`|`--hmm_meta`|`None`|path to hmm database metadata file| |`-l`|`--log`|`None`|path to log file. Log not written by default.| @@ -36,3 +37,11 @@ where $n_{ab}$ corresponds to the maximum number of genes between $HMM_a$ and $H Several HMMs can be assigned to the same ORF, in which case the search is performed for all of them. In this case, HMM names must be separated by "|" and grouped within parentheses, as shown above. If the PGAP database is employed (see `pynteny download` below), synteny blocks can also be specified by gene symbols, such as $$\lt leuD \space 0 \space \lt leuC \space 1 \space \lt leuA.$$ In that case, the program will try to match gene symbols to HMM names in the PGAP database before running the search. + +### A note on paralogous HMMs + +By default, Pynteny's synteny matcher does **not** disambiguate paralogous HMMs. If your panel contains close paralogs (e.g. nitrogenase α/β subunits `nifD`/`nifK`, MoFe-cofactor biosynthesis `nifE`/`nifN`, or light-independent protochlorophyllide reductase `bchL`/`nifH`), the same peptide can score above threshold under more than one HMM. These cross-hits add duplicate rows that can interleave with true syntenic neighbours and cause a genuine syntenic block to be **silently dropped** (a warning is printed). To handle paralog-aware panels, prefer one of: + +- a disjoint, orthogonal HMM set; +- `|`-grouped HMMs declared in the synteny structure (e.g. `>(nifD|nifK)`); or +- the `--best_hmm_wins` flag, which assigns each peptide to its single highest-scoring HMM before matching. diff --git a/src/pynteny/api.py b/src/pynteny/api.py index 63294db..6066a2f 100644 --- a/src/pynteny/api.py +++ b/src/pynteny/api.py @@ -71,6 +71,7 @@ def __init__( synteny_struc: str, gene_ids: bool = False, unordered: bool = False, + best_hmm_wins: bool = False, reuse: bool = False, hmm_dir: Path = None, hmm_meta: Path = None, @@ -100,6 +101,9 @@ def __init__( exact same order displayed in the synteny_structure or in any order If ordered, the filters would filter collinear rather than syntenic structures. Defaults to False. + best_hmm_wins (bool, optional): if True, when the same peptide is hit + by more than one HMM (paralog cross-hits), keep only the + highest-scoring HMM for that peptide. Defaults to False. reuse (bool, optional): if True then HMMER3 won't be run again for HMMs already searched in the same output directory. Defaults to False. hmm_dir (Path, optional): path to directory containing input HMM files. @@ -132,6 +136,7 @@ def __init__( logfile=Path(logfile) if logfile is not None else logfile, processes=processes, unordered=unordered, + best_hmm_wins=best_hmm_wins, reuse=reuse, ) diff --git a/src/pynteny/cli.py b/src/pynteny/cli.py index 75648d1..09544fe 100644 --- a/src/pynteny/cli.py +++ b/src/pynteny/cli.py @@ -287,6 +287,20 @@ def search() -> argparse.ArgumentParser: "pair is considered to run the search." ), ) + optional.add_argument( + "-b", + "--best_hmm_wins", + dest="best_hmm_wins", + default=False, + action="store_true", + help=( + "if the same peptide is hit by more than one HMM (e.g. paralogous \n" + "models that cross-hit), keep only the highest-scoring HMM for that \n" + "peptide. Use this when your HMM panel contains close paralogs \n" + "(e.g. nitrogenase nifD/nifK, nifE/nifN) to avoid silently dropping \n" + "genuine syntenic blocks. Disabled by default." + ), + ) optional.add_argument( "-r", "--reuse", diff --git a/src/pynteny/filter.py b/src/pynteny/filter.py index 5f5f5eb..2d6a10a 100644 --- a/src/pynteny/filter.py +++ b/src/pynteny/filter.py @@ -145,7 +145,11 @@ class SyntenyHMMfilter: """Tools to search for synteny structures among sets of hmm models""" def __init__( - self, hmm_hits: dict, synteny_structure: str, unordered: bool = True + self, + hmm_hits: dict, + synteny_structure: str, + unordered: bool = True, + best_hmm_wins: bool = False, ) -> None: """Search for contigs that satisfy the given gene synteny structure. @@ -166,8 +170,14 @@ def __init__( exact same order displayed in the synteny_structure or in any order. If ordered, the filters would filter collinear rather than syntenic structures. Defaults to True. + best_hmm_wins (bool, optional): if True, when the same peptide is hit + by more than one HMM (e.g. paralogous models that cross-hit), + keep only the highest-scoring HMM for that peptide. This prevents + duplicate rows from breaking the rolling-window synteny matcher. + Defaults to False. """ self._unordered = unordered + self._best_hmm_wins = best_hmm_wins self._hmm_hits = hmm_hits self._hmms = list(hmm_hits.keys()) self._synteny_structure = synteny_structure @@ -263,8 +273,15 @@ def get_all_HMM_hits(self) -> pd.DataFrame: if not labels: logger.error(f"No records found in database matching HMM: {hmm}") sys.exit(1) - hit_labels[hmm] = labelparser.parse_from_list(labels) - hit_labels[hmm]["hmm"] = hmm + parsed = labelparser.parse_from_list(labels) + parsed["hmm"] = hmm + # Carry the HMMER bitscore along (only when needed) so paralog + # cross-hits on the same peptide can be disambiguated by score. + # It is dropped again right after deduplication so the default code + # path keeps exactly the same columns as before. + if self._best_hmm_wins and "bitscore" in hits.columns: + parsed["bitscore"] = hits.bitscore.values + hit_labels[hmm] = parsed # Create single dataframe with new column corresponding to HMM and all hits # Remove contigs with less hits than the number of hmms in synteny structure (drop if enabling partial searching) all_hit_labels = ( @@ -274,6 +291,20 @@ def get_all_HMM_hits(self) -> pd.DataFrame: .sort_values(["contig", "gene_pos"], ascending=True) ) all_hit_labels = all_hit_labels.reset_index(drop=True) + if self._best_hmm_wins: + # Keep only the highest-scoring HMM per peptide. Without this, + # paralogous HMMs that cross-hit the same peptide create duplicate + # rows that interleave with true syntenic neighbours and break the + # rolling-window matcher (silent false negatives). + if "bitscore" not in all_hit_labels.columns: + all_hit_labels["bitscore"] = float("nan") + all_hit_labels = ( + all_hit_labels.sort_values("bitscore", ascending=False) + .drop_duplicates(subset=["full_label"], keep="first") + .sort_values(["contig", "gene_pos"], ascending=True) + .reset_index(drop=True) + .drop(columns="bitscore") + ) if self._contains_hmm_groups: all_hit_labels = self._merge_hits_by_HMM_group(all_hit_labels) all_hit_labels = all_hit_labels.reset_index(drop=True) @@ -291,8 +322,22 @@ def filter_hits_by_synteny_structure(self) -> dict: ) all_hit_labels = self.get_all_HMM_hits() if all_hit_labels.full_label.duplicated().any(): + n_dup = int(all_hit_labels.full_label.duplicated().sum()) + dup_mask = all_hit_labels.full_label.duplicated(keep=False) + affected_combinations = ( + all_hit_labels[dup_mask] + .groupby("full_label") + .hmm.apply(lambda s: tuple(sorted(set(s)))) + .value_counts() + .head(5) + .to_dict() + ) logger.warning( - "At least two different HMMs produced identical sequence hits" + f"{n_dup} peptide(s) hit by >=2 HMMs (paralog cross-hits). " + f"Top HMM combinations: {affected_combinations}. " + "Syntenic patterns adjacent to these peptides may be silently " + "dropped. Re-run with --best_hmm_wins to keep only each " + "peptide's highest-scoring HMM." ) if all_hit_labels.full_label.duplicated().sum() == all_hit_labels.shape[0] / 2: logger.error( @@ -505,6 +550,7 @@ def filter_FASTA_by_synteny_structure( input_fasta: Path, input_hmms: list[Path], unordered: bool = False, + best_hmm_wins: bool = False, hmm_meta: Path = None, hmmer_output_dir: Path = None, reuse_hmmer_results: bool = True, @@ -532,6 +578,10 @@ def filter_FASTA_by_synteny_structure( displayed in the synteny structure string or not, i.e., whether to search for only synteny (colocation) or collinearity as well (same order). Defaults to False. + best_hmm_wins (bool, optional): if True, when the same peptide is hit by + more than one HMM (paralog cross-hits), keep only the highest-scoring + HMM for that peptide before matching the synteny structure. Defaults + to False. hmm_meta (Path, optional): path to PGAP's metadata file. Defaults to None. hmmer_output_dir (Path, optional): output directory to store HMMER3 output files. Defaults to None. reuse_hmmer_results (bool, optional): if True then HMMER3 won't be run again for HMMs already @@ -601,7 +651,12 @@ def filter_FASTA_by_synteny_structure( sys.exit(1) logger.info("Filtering results by synteny structure") - syntenyfilter = SyntenyHMMfilter(hmm_hits, synteny_structure, unordered=unordered) + syntenyfilter = SyntenyHMMfilter( + hmm_hits, + synteny_structure, + unordered=unordered, + best_hmm_wins=best_hmm_wins, + ) hits_by_contig = syntenyfilter.filter_hits_by_synteny_structure() if hmm_meta is not None: return SyntenyHits.from_hits_dict(hits_by_contig).add_HMM_meta_info_to_hits( diff --git a/src/pynteny/subcommands.py b/src/pynteny/subcommands.py index ee58461..41ccdf8 100644 --- a/src/pynteny/subcommands.py +++ b/src/pynteny/subcommands.py @@ -152,6 +152,7 @@ def synteny_search(args: Union[CommandArgs, ArgumentParser]) -> SyntenyHits: synteny_hits = filter_FASTA_by_synteny_structure( synteny_structure=args.synteny_struc, unordered=args.unordered, + best_hmm_wins=getattr(args, "best_hmm_wins", False), input_fasta=args.data, input_hmms=input_hmms, hmm_meta=args.hmm_meta, diff --git a/tests/test_filter.py b/tests/test_filter.py new file mode 100644 index 0000000..6878e48 --- /dev/null +++ b/tests/test_filter.py @@ -0,0 +1,86 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" +Unit tests for the synteny filter, including the paralog cross-hit case +(see issue #96). +""" + +import unittest + +import pandas as pd + +from pynteny.filter import SyntenyHMMfilter + + +def _paralog_crosshit_hits() -> dict: + """Build a minimal hmm_hits dict reproducing a paralog cross-hit scenario. + + A textbook syntenic operon on a single contig, all on the positive strand, + at gene positions 33 -> 34 -> 35. Each peptide has a single fixed label + (independent of which HMM hits it). The paralogous nifD/nifK models also + cross-hit each other's peptides, so peptides p34 and p35 each appear in two + HMM hit tables. The highest-scoring HMM per peptide is exactly the canonical + nifH -> nifD -> nifK operon. + """ + p33 = "p33__ctg_33_1_100_pos" # true nifH peptide + p34 = "p34__ctg_34_1_100_pos" # true nifD peptide + p35 = "p35__ctg_35_1_100_pos" # true nifK peptide + nifH = pd.DataFrame({"id": [p33], "bitscore": [452.7]}) + nifD = pd.DataFrame( + {"id": [p34, p35], "bitscore": [856.5, 95.0]} # p34 real, p35 cross-hit + ) + nifK = pd.DataFrame( + {"id": [p34, p35], "bitscore": [88.0, 861.6]} # p34 cross-hit, p35 real + ) + return {"nifH": nifH, "nifD": nifD, "nifK": nifK} + + +SYNTENY_STRUC = ">nifH 0 >nifD 0 >nifK" + + +def _matched_labels(matched_hits: dict) -> set: + """Flatten the {contig: {hmm: [labels]}} structure into a set of labels.""" + return { + label + for hmm_to_labels in matched_hits.values() + for labels in hmm_to_labels.values() + for label in labels + } + + +class TestParalogCrossHits(unittest.TestCase): + def test_default_drops_operon_on_crosshit(self): + """Without best_hmm_wins, cross-hits silently drop the operon (issue #96).""" + synteny_filter = SyntenyHMMfilter( + _paralog_crosshit_hits(), SYNTENY_STRUC, unordered=False + ) + labels = _matched_labels(synteny_filter.filter_hits_by_synteny_structure()) + self.assertNotIn( + "p33__ctg_33_1_100_pos", + labels, + "Expected the operon to be dropped under cross-hits with default behaviour", + ) + + def test_best_hmm_wins_recovers_operon(self): + """With best_hmm_wins, each peptide is assigned its best HMM and the + canonical operon is recovered.""" + synteny_filter = SyntenyHMMfilter( + _paralog_crosshit_hits(), + SYNTENY_STRUC, + unordered=False, + best_hmm_wins=True, + ) + labels = _matched_labels(synteny_filter.filter_hits_by_synteny_structure()) + self.assertTrue( + { + "p33__ctg_33_1_100_pos", + "p34__ctg_34_1_100_pos", + "p35__ctg_35_1_100_pos", + }.issubset(labels), + "Expected best_hmm_wins to recover the full syntenic operon", + ) + + +if __name__ == "__main__": + unittest.main()