Fix paralog cross-hits silently dropping syntenic blocks (--best_hmm_wins)

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.

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,
         )
 

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",

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(

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,

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()