diff --git a/src/geckopy/databases/brenda_loader.py b/src/geckopy/databases/brenda_loader.py index 4fe7827..0e005d0 100644 --- a/src/geckopy/databases/brenda_loader.py +++ b/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 diff --git a/src/geckopy/gather_kcats/merge_dlkcat_and_fuzzy_kcats.py b/src/geckopy/gather_kcats/merge_dlkcat_and_fuzzy_kcats.py index 7455d5c..d88db2e 100644 --- a/src/geckopy/gather_kcats/merge_dlkcat_and_fuzzy_kcats.py +++ b/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: diff --git a/src/geckopy/gather_kcats/merge_kcats.py b/src/geckopy/gather_kcats/merge_kcats.py index 1443fca..860a89d 100644 --- a/src/geckopy/gather_kcats/merge_kcats.py +++ b/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 diff --git a/tests/test_brenda_loader.py b/tests/test_brenda_loader.py index 7585f24..266e874 100644 --- a/tests/test_brenda_loader.py +++ b/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 # --------------------------------------------------------------------------- # diff --git a/tests/test_merge_dlkcat_and_fuzzy_kcats.py b/tests/test_merge_dlkcat_and_fuzzy_kcats.py index f083bde..cd52cc7 100644 --- a/tests/test_merge_dlkcat_and_fuzzy_kcats.py +++ b/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"