KEGG database, solver tolerance and zero-total f factor as in GECKO

src/geckopy/ec_model/make_ec_model.py

@@ -29,7 +29,7 @@
 import cobra
 from raven_toolbox.utils.sort import sort_identifiers
 
-from ..databases import UniprotDB, load_uniprot_tsv
+from ..databases import UniprotDB, load_kegg_tsv, load_uniprot_tsv
 from ..solver_support import repair_scipy_problem, switch_to_scipy
 from .ec_model import EcModel
 from .pipeline import (
@@ -162,8 +162,9 @@ def make_ec_model(
         UniProt accession carried in the KEGG row, or to the KEGG
         gene id when that accession is empty), and the EC-number
         lookup (``fill_eccodes_from_database``) consults KEGG when
-        the UniProt result is empty or ends with ``-``. If None,
-        KEGG is not used.
+        the UniProt result is empty or ends with ``-``. If None, the
+        function loads ``adapter.params.path / "data" / "kegg.tsv"``
+        when that file exists, and uses no KEGG data otherwise.
 
     Returns
     -------
@@ -246,6 +247,11 @@ def make_ec_model(
         uniprot_path = adapter.params.path / "data" / "uniprot.tsv"
         uniprot_db = load_uniprot_tsv(uniprot_path)
 
+    if kegg_db is None:
+        kegg_path = adapter.params.path / "data" / "kegg.tsv"
+        if kegg_path.is_file():
+            kegg_db = load_kegg_tsv(kegg_path, auto_download=False)
+
     if gecko_light:
         # Light: stage 7 (per-enzyme data) runs first so the coupling
         # builder can index ec.genes when emitting the per-isozyme rows.

src/geckopy/gather_kcats/merge_dlkcat_and_fuzzy_kcats.py

@@ -51,6 +51,7 @@ def merge_dlkcat_and_fuzzy_kcats(
         DataFrame produced by ``fuzzy_kcat_matching``.
     top_origin_limit
         Origin upper bound for database_top. Must be in ``[1, 6]``.
+        Default 4.
     bottom_origin_limit
         Origin upper bound for database_bottom. Must be in ``[1, 6]``.
     wildcard_limit

src/geckopy/limit_proteins/calculate_f_factor.py

@@ -73,9 +73,9 @@ def calculate_f_factor(
     Returns
     -------
     float
-        The f-factor in [0, 1]. Returns 0.0 if the total proteome
-        abundance is zero, or if none of the proteomics UniProt IDs
-        match a model enzyme.
+        The f-factor in [0, 1]. Returns NaN, as MATLAB does, if the total
+        proteome abundance is zero (the factor is 0/0), and 0.0 if none of
+        the proteomics UniProt IDs match a model enzyme.
     """
     if enzymes is None:
         enzymes = model.ec.enzymes
@@ -90,9 +90,9 @@ def calculate_f_factor(
     if total_prot == 0.0:
         logger.warning(
             "calculate_f_factor: total proteome abundance is zero; "
-            "returning f = 0."
+            "returning f = NaN."
         )
-        return 0.0
+        return float("nan")
 
     enzyme_set = set(enzymes)
     in_model = np.array(

src/geckopy/utilities/load_conventional_gem.py

@@ -22,6 +22,7 @@
 
 import cobra
 from raven_toolbox.io import read_yaml_model
+from raven_toolbox.utils.solver import set_raven_tolerance
 
 from ..solver_support import fast_solver_while_building, switch_to_scipy
 
@@ -48,8 +49,8 @@ def load_conventional_gem(adapter: "ModelAdapter") -> cobra.Model:
     Returns
     -------
     cobra.Model
-        The loaded model. Use it as the input to
-        ``make_ec_model``.
+        The loaded model, with the solver tolerance set to RAVEN's 1e-9.
+        Use it as the input to ``make_ec_model``.
 
     Raises
     ------
@@ -73,4 +74,4 @@ def load_conventional_gem(adapter: "ModelAdapter") -> cobra.Model:
             model = cobra.io.read_sbml_model(str(path))
     if on_fast_solver:
         switch_to_scipy(model)
-    return model
+    return set_raven_tolerance(model)

tests/test_calculate_f_factor.py

@@ -70,17 +70,17 @@ def test_subset_match_yields_correct_ratio():
     assert calculate_f_factor(model, pd) == pytest.approx(0.25)
 
 
-def test_empty_proteome_yields_zero():
+def test_empty_proteome_yields_nan():
     model = _ec_model(["P1"])
     pd = _prot_data([])
-    assert calculate_f_factor(model, pd) == 0.0
+    assert np.isnan(calculate_f_factor(model, pd))
 
 
-def test_zero_total_proteome_yields_zero():
-    """All abundances are zero -> f = 0 (defensive)."""
+def test_zero_total_proteome_yields_nan():
+    """All abundances are zero -> the factor is 0/0, NaN as in MATLAB."""
     model = _ec_model(["P1"])
     pd = _prot_data([("P1", 0.0), ("P2", 0.0)])
-    assert calculate_f_factor(model, pd) == 0.0
+    assert np.isnan(calculate_f_factor(model, pd))
 
 
 # --------------------------------------------------------------------------- #

tests/test_load_conventional_gem.py

@@ -95,3 +95,11 @@ def test_missing_file_raises(tmp_path):
     adapter = _adapter_with_conv_gem(tmp_path, "absent.yml")
     with pytest.raises(FileNotFoundError, match="conv_gem"):
         load_conventional_gem(adapter)
+
+
+def test_loaded_gem_has_raven_solver_tolerance(tmp_path):
+    """The model carries RAVEN's solver tolerance whichever reader was used."""
+    json_path = tmp_path / "tiny.json"
+    cobra.io.save_json_model(cobra.Model("tiny"), str(json_path))
+    adapter = _adapter_with_conv_gem(tmp_path, "tiny.json")
+    assert load_conventional_gem(adapter).tolerance == 1e-9

tests/test_load_ec_model.py

@@ -655,3 +655,10 @@ def test_enzymes_omap_without_entries_loads(tmp_path):
     assert model.ec.rxns == ["R1_EXP_1", "R2", "R3"]
     assert model.ec.notes == ["", 'say "hi" now', ""]
     np.testing.assert_array_equal(model.ec.rxn_enz_mat.toarray(), [[0.0], [1.0], [0.0]])
+
+
+def test_loaded_ec_model_has_raven_solver_tolerance(tmp_path):
+    adapter = _adapter(tmp_path)
+    _write_canonical(tmp_path / "models" / "ecModel.yml")
+    model = load_ec_model("ecModel.yml", adapter=adapter)
+    assert model.tolerance == 1e-9

tests/test_make_ec_model.py

@@ -332,3 +332,48 @@ def test_make_ec_model_unmatched_genes_logged_and_annotated(tmp_path, caplog):
     warned = ec_model.reactions.get_by_id("r1")
     assert "geckopy_warning" in warned.notes
     assert "g_missing" in warned.notes["geckopy_warning"]
+
+
+def _toy_with_missing_gene(tmp_path):
+    (tmp_path / "model_adapter.toml").write_text(
+        'conv_gem = "dummy.xml"\n'
+        'org_name = "test"\n'
+        'enzyme_comp = "c"\n'
+    )
+    adapter = ModelAdapter.from_folder(tmp_path)
+    model = cobra.Model("m")
+    model.compartments = {"c": "c"}
+    m_a = cobra.Metabolite("A", compartment="c")
+    m_b = cobra.Metabolite("B", compartment="c")
+    model.add_metabolites([m_a, m_b])
+    r = cobra.Reaction("r1")
+    r.add_metabolites({m_a: -1, m_b: 1})
+    r.lower_bound, r.upper_bound = 0.0, 1000.0
+    r.gene_reaction_rule = "g_known and g_kegg"
+    model.add_reactions([r])
+    db = UniprotDB(
+        ids=["P1"],
+        genes=["g_known"],
+        eccodes=[""],
+        mw=np.array([10.0]),
+        sequences=["S"],
+    )
+    return model, adapter, db
+
+
+def test_make_ec_model_reads_kegg_tsv_of_the_data_folder(tmp_path):
+    """A gene UniProt does not know is matched through data/kegg.tsv when
+    that file exists, without passing kegg_db."""
+    model, adapter, db = _toy_with_missing_gene(tmp_path)
+    (tmp_path / "data").mkdir()
+    (tmp_path / "data" / "kegg.tsv").write_text(
+        'P2,g_kegg,g_kegg,1.1.1.1,20000,"",MSEQ\n'
+    )
+    ec_model = make_ec_model(model, adapter, uniprot_db=db)
+    assert sorted(ec_model.ec.enzymes) == ["P1", "P2"]
+
+
+def test_make_ec_model_without_kegg_tsv_leaves_the_gene_unmatched(tmp_path):
+    model, adapter, db = _toy_with_missing_gene(tmp_path)
+    ec_model = make_ec_model(model, adapter, uniprot_db=db)
+    assert sorted(ec_model.ec.enzymes) == ["P1"]