diff --git a/src/geckopy/ec_model/make_ec_model.py b/src/geckopy/ec_model/make_ec_model.py index 8a4f3db..69f7d35 100644 --- a/src/geckopy/ec_model/make_ec_model.py +++ b/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. 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 d88db2e..d3c975c 100644 --- a/src/geckopy/gather_kcats/merge_dlkcat_and_fuzzy_kcats.py +++ b/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 diff --git a/src/geckopy/limit_proteins/calculate_f_factor.py b/src/geckopy/limit_proteins/calculate_f_factor.py index e363f44..8e456d3 100644 --- a/src/geckopy/limit_proteins/calculate_f_factor.py +++ b/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( diff --git a/src/geckopy/utilities/load_conventional_gem.py b/src/geckopy/utilities/load_conventional_gem.py index 812d8b7..f89c0b2 100644 --- a/src/geckopy/utilities/load_conventional_gem.py +++ b/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) diff --git a/tests/test_calculate_f_factor.py b/tests/test_calculate_f_factor.py index 8c0d322..cae1d40 100644 --- a/tests/test_calculate_f_factor.py +++ b/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)) # --------------------------------------------------------------------------- # diff --git a/tests/test_load_conventional_gem.py b/tests/test_load_conventional_gem.py index 6a83b73..de02c10 100644 --- a/tests/test_load_conventional_gem.py +++ b/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 diff --git a/tests/test_load_ec_model.py b/tests/test_load_ec_model.py index ad9d17a..01febe0 100644 --- a/tests/test_load_ec_model.py +++ b/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 diff --git a/tests/test_make_ec_model.py b/tests/test_make_ec_model.py index 88a2601..70a4a44 100644 --- a/tests/test_make_ec_model.py +++ b/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"]