BRENDA: filter implausible SA-derived kcats; default top_origin_limit 6 -> 4

src/geckopy/databases/brenda_loader.py

@@ -33,6 +33,7 @@
 from pathlib import Path
 from typing import IO, Literal
 
+import numpy as np
 import pandas as pd
 
 logger = logging.getLogger(__name__)
@@ -95,7 +96,9 @@ def sa_for(self, aggregation: Literal["max", "median"]) -> pd.DataFrame:
         )
 
 
-def load_brenda_data(folder: str | Path) -> BrendaData:
+def load_brenda_data(
+    folder: str | Path, *, filter_outlier_kcats: bool = True,
+) -> BrendaData:
     """Load BRENDA kcat, SA, and MW dumps from a folder.
 
     Ported from GECKO MATLAB:
@@ -129,6 +132,12 @@ def load_brenda_data(folder: str | Path) -> BrendaData:
     ----------
     folder
         Directory containing all three BRENDA dump files.
+    filter_outlier_kcats
+        If True (default), drop a row whose kcat is more than 4 orders
+        of magnitude below the typical kcat of its EC before it can be
+        matched to any reaction -- see ``_filter_outlier_rows`` for the
+        rationale and mechanism. Set False to reproduce the unfiltered
+        1.0-era behaviour.
 
     Returns
     -------
@@ -170,6 +179,11 @@ def load_brenda_data(folder: str | Path) -> BrendaData:
         _sa_view_from_wide(sa_wide, "sa_median"), mw_raw,
     )
 
+    if filter_outlier_kcats:
+        kcat_max, kcat_median, sa_max, sa_median = _filter_outlier_rows(
+            [kcat_max, kcat_median, sa_max, sa_median],
+        )
+
     return BrendaData(
         kcat_max=kcat_max,
         kcat_median=kcat_median,
@@ -371,3 +385,95 @@ def _join_sa_with_mw(sa: pd.DataFrame, mw: pd.DataFrame) -> pd.DataFrame:
     if not out_rows:
         return pd.DataFrame(columns=_SA_COLUMNS)
     return pd.DataFrame(out_rows, columns=_SA_COLUMNS)
+
+
+def _filter_outlier_rows(
+    tables: list[pd.DataFrame], *, decades: float = 4.0, min_ref: int = 5,
+) -> list[pd.DataFrame]:
+    """Drop a row whose kcat is implausibly low for its EC.
+
+    BRENDA's specific-activity x molecular-weight join can produce a
+    derived kcat that is many orders of magnitude below anything else
+    reported for the same EC -- almost always a unit or transcription
+    error upstream, not a real slow enzyme (a documented case: EC
+    2.3.1.297's only molecular-weight row is for a different organism
+    than its lowest specific-activity row, and the derived kcat this
+    produces, ~4e-7 1/s, is 8 orders of magnitude below the yeast SA
+    measurement for the same EC that has no matching MW row and is
+    silently dropped instead). Once matched to a reaction, a kcat that
+    small makes the enzyme-usage coefficient enormous and can make a
+    model untunable.
+
+    This is a *relative* floor, not a fixed threshold: a row is dropped
+    only if it is more than ``decades`` orders of magnitude below the
+    median kcat of its own EC. "Its own EC" widens from the full
+    4-level code to the 3-, 2-, then 1-level prefix if fewer than
+    ``min_ref`` values are available at the current level, so a
+    sparsely-represented EC is still judged against something. The
+    reference pool for one EC is built from every positive-kcat row of
+    every table passed in (plain kcat rows and SA-derived rows alike),
+    computed once up front from the *unfiltered* tables so that
+    dropping a row in one table cannot shift the reference used to
+    judge another.
+
+    Real enzymes can be genuinely slow, so this only catches a row far
+    outside its own EC's spread, not a merely-low one.
+
+    Parameters
+    ----------
+    tables
+        The four BRENDA views (``kcat_max``, ``kcat_median``,
+        ``sa_max``, ``sa_median``); each must have ``ec_code`` and
+        ``kcat`` columns. Order is preserved in the return value.
+    decades
+        How many orders of magnitude below the reference counts as an
+        outlier.
+    min_ref
+        Minimum number of reference values required before a level is
+        trusted; the search widens to a broader EC prefix otherwise.
+
+    Returns
+    -------
+    list[pandas.DataFrame]
+        The same tables, each with outlier rows removed and its index
+        reset.
+    """
+    pool: dict[str, list[float]] = {}
+    for table in tables:
+        for ec, kcat in zip(table["ec_code"], table["kcat"]):
+            if kcat > 0 and "-" not in ec:
+                pool.setdefault(ec, []).append(np.log10(kcat))
+
+    log_values_by_level: dict[tuple[int, str], list[float]] = {}
+    for ec, values in pool.items():
+        parts = ec.split(".")
+        for level in (4, 3, 2, 1):
+            key = (level, ".".join(parts[:level]))
+            log_values_by_level.setdefault(key, []).extend(values)
+
+    reference_cache: dict[str, float | None] = {}
+
+    def reference_log_kcat(ec: str) -> float | None:
+        if ec not in reference_cache:
+            parts = ec.split(".")
+            reference_cache[ec] = None
+            for level in (4, 3, 2, 1):
+                values = log_values_by_level.get((level, ".".join(parts[:level])), [])
+                if len(values) >= min_ref:
+                    reference_cache[ec] = float(np.median(values))
+                    break
+        return reference_cache[ec]
+
+    filtered = []
+    for table in tables:
+        if table.empty:
+            # An empty boolean mask loses the (already-correct) column
+            # set on an empty frame; nothing to filter anyway.
+            filtered.append(table)
+            continue
+        keep = []
+        for ec, kcat in zip(table["ec_code"], table["kcat"]):
+            ref = reference_log_kcat(ec) if kcat > 0 else None
+            keep.append(not (ref is not None and np.log10(kcat) < ref - decades))
+        filtered.append(table[pd.Series(keep, index=table.index)].reset_index(drop=True))
+    return filtered

src/geckopy/gather_kcats/merge_dlkcat_and_fuzzy_kcats.py

@@ -22,7 +22,7 @@ def merge_dlkcat_and_fuzzy_kcats(
     dlkcat_kcats: pd.DataFrame,
     fuzzy_kcats: pd.DataFrame,
     *,
-    top_origin_limit: int = 6,
+    top_origin_limit: int = 4,
     bottom_origin_limit: int = 6,
     wildcard_limit: int = 3,
 ) -> pd.DataFrame:

src/geckopy/gather_kcats/merge_kcats.py

@@ -94,7 +94,7 @@ def _strip_okp_prefix(normalized_source: str) -> str:
 def merge_kcats(
     *kcat_lists: pd.DataFrame,
     source_priority: Sequence[str],
-    top_origin_limit: int = 6,
+    top_origin_limit: int = 4,
     bottom_origin_limit: int = 6,
     wildcard_limit: int = 3,
     database_sources: Iterable[str] = _DEFAULT_DATABASE_SOURCES,
@@ -126,6 +126,13 @@ def merge_kcats(
         source not listed here is dropped (with a warning).
     top_origin_limit
         Origin upper bound for ``database_top`` fuzzy rows. ``[1, 6]``.
+        Default ``4``, not ``6``: origins 5-6 are specific-activity ×
+        molecular-weight derived kcats, which BRENDA's own data quality
+        does not always support treating as better than a prediction
+        (a bad SA/MW join can otherwise silently override a good
+        prediction with an outlier -- see ``load_brenda_data``'s
+        ``filter_outlier_kcats``, which removes the worst of these
+        before they reach this function, but does not catch every one).
     bottom_origin_limit
         Origin upper bound for ``database_bottom`` fuzzy rows. ``[1, 6]``.
     wildcard_limit

tests/test_brenda_loader.py

@@ -88,6 +88,85 @@ def test_all_empty_files_yield_empty_dataframes(tmp_path):
     assert result.sa_median.empty
 
 
+# --------------------------------------------------------------------------- #
+# Outlier filtering (filter_outlier_kcats)
+# --------------------------------------------------------------------------- #
+
+
+def _kcat_row(ec, organism, kcat):
+    return f"{ec}\tm1\t{organism}\t{kcat}\t{kcat}\t1\t*\n"
+
+
+def test_outlier_row_dropped_by_default(tmp_path):
+    # Five real-looking values for the same EC (the reference pool),
+    # plus one 8-decades-lower outlier for that same EC.
+    kcat = "".join(
+        _kcat_row("1.1.1.1", f"org{i}", v)
+        for i, v in enumerate([10, 12, 8, 15, 11])
+    )
+    kcat += _kcat_row("1.1.1.1", "org_outlier", 1e-6)
+    _write_brenda_files(tmp_path, kcat=kcat)
+
+    filtered = load_brenda_data(tmp_path)
+    assert len(filtered.kcat_max) == 5
+    assert 1e-6 not in filtered.kcat_max["kcat"].values
+
+    unfiltered = load_brenda_data(tmp_path, filter_outlier_kcats=False)
+    assert len(unfiltered.kcat_max) == 6
+    assert 1e-6 in unfiltered.kcat_max["kcat"].values
+
+
+def test_reference_pool_widens_to_ec_prefix_when_sparse(tmp_path):
+    # Only one row for 1.1.1.9 itself (too few for its own reference),
+    # but four more under the shared 1.1.1.- prefix -- enough (5) once
+    # widened. The lone 1.1.1.9 row is a real outlier against that
+    # widened reference and should be dropped.
+    kcat = "".join(
+        _kcat_row(ec, f"org{i}", v)
+        for i, (ec, v) in enumerate([
+            ("1.1.1.1", 10), ("1.1.1.2", 12), ("1.1.1.3", 8), ("1.1.1.4", 15),
+        ])
+    )
+    kcat += _kcat_row("1.1.1.9", "org_outlier", 1e-6)
+    _write_brenda_files(tmp_path, kcat=kcat)
+
+    filtered = load_brenda_data(tmp_path)
+    assert "1.1.1.9" not in filtered.kcat_max["ec_code"].values
+    assert len(filtered.kcat_max) == 4
+
+
+def test_outlier_filter_uses_kcat_and_sa_together_as_one_reference_pool(tmp_path):
+    # The reference pool for one EC draws on both the plain kcat table
+    # and the SA-derived table; four kcat rows plus one SA row is
+    # enough (5) to judge a sixth, wildly-off SA row as an outlier.
+    kcat = "".join(
+        _kcat_row("2.3.1.297", f"org{i}", v)
+        for i, v in enumerate([10, 12, 8, 15])
+    )
+    sa = (
+        "2.3.1.297\t*\torg_good\t2.0\t2.0\t1\t*\n"     # ~ SA*MW = 20, in range
+        "2.3.1.297\t*\torg_bad\t1e-8\t1e-8\t1\t*\n"    # ~ SA*MW = 1e-7, outlier
+    )
+    mw = (
+        "2.3.1.297\t*\torg_good\t10000\t1\t*\n"
+        "2.3.1.297\t*\torg_bad\t10000\t1\t*\n"
+    )
+    _write_brenda_files(tmp_path, kcat=kcat, sa=sa, mw=mw)
+
+    filtered = load_brenda_data(tmp_path)
+    assert len(filtered.sa_max) == 1
+    assert filtered.sa_max.iloc[0]["organism"] == "org_good"
+
+
+def test_empty_tables_unaffected_by_outlier_filter(tmp_path):
+    _write_brenda_files(tmp_path)
+    result = load_brenda_data(tmp_path, filter_outlier_kcats=True)
+    assert list(result.kcat_max.columns) == [
+        "ec_code", "substrate", "organism", "kcat", "n",
+    ]
+    assert result.kcat_max.empty
+
+
 # --------------------------------------------------------------------------- #
 # Missing files
 # --------------------------------------------------------------------------- #

tests/test_merge_dlkcat_and_fuzzy_kcats.py

@@ -160,13 +160,30 @@ def test_fuzzy_prio1_wins_when_reaction_also_in_dlkcat():
     assert row["kcat"] == 5.0
 
 
-def test_fuzzy_prio1_uses_default_top_origin_limit_6():
-    """With default top_origin_limit=6, origin=6 still qualifies for prio1."""
+def test_fuzzy_prio1_uses_default_top_origin_limit_4():
+    """Default top_origin_limit=4: an origin=6 (specific-activity x MW
+    derived) match does NOT qualify for prio1, so DLKcat wins instead.
+    Origins 5-6 are demoted by default because BRENDA's SA/MW join can
+    produce a wildly wrong outlier that would otherwise silently
+    outrank a good prediction -- see brenda_loader's
+    filter_outlier_kcats, which this default is meant to pair with."""
     out = merge_dlkcat_and_fuzzy_kcats(
         _dlkcat_df([_dlkcat_row("r1", 100.0)]),
         _fuzzy_df([_fuzzy_row("r1", 5.0, wildcard_level=0, origin=6)]),
     )
     assert len(out) == 1
+    assert out.iloc[0]["source"] == "DLKcat"
+
+
+def test_fuzzy_prio1_origin_6_qualifies_when_limit_explicitly_widened():
+    """The pre-4-default behaviour is still reachable by passing
+    top_origin_limit=6 explicitly."""
+    out = merge_dlkcat_and_fuzzy_kcats(
+        _dlkcat_df([_dlkcat_row("r1", 100.0)]),
+        _fuzzy_df([_fuzzy_row("r1", 5.0, wildcard_level=0, origin=6)]),
+        top_origin_limit=6,
+    )
+    assert len(out) == 1
     assert out.iloc[0]["source"] == "brenda"