diff --git a/CHANGELOG.md b/CHANGELOG.md index 7fd2718..6e57433 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,17 @@ ## Unreleased +* Fix: `write_dlkcat_input` and the OpenKineticsPredictor input builder + raised `TypeError: cannot use 'list' as a set element` on any model that + had been saved and reloaded. cobra returns an annotation value as a list + for a model read from a file, so `annotation['smiles']` is `"CCO"` on a + model annotated in memory by `find_met_smiles` and `["CCO"]` on the same + model after a round trip, and the ignore-list check hashed it. The new + `databases.met_smiles` reads either shape; a reloaded model now writes + the same DLKcat input as a freshly annotated one. +* Fix: `parse_okp_output` raised `TypeError: cannot use 'list' as a dict + key` on a reloaded model for the same reason, when indexing metabolites + by their SMILES annotation to map each prediction back to an `ec.rxn`. * Fix: evotune's trust tiers never applied to a fuzzy-matched kcat. `ec.source` reads `"brenda (wc=0, origin=1)"` where MATLAB's reads `"brenda"`, and `classify_kcat_source` compared the whole string against diff --git a/src/geckopy/databases/__init__.py b/src/geckopy/databases/__init__.py index 0a71b50..0f88da5 100644 --- a/src/geckopy/databases/__init__.py +++ b/src/geckopy/databases/__init__.py @@ -13,7 +13,7 @@ from .pax_db_loader import ProtData, load_pax_db from .phyl_dist import PhylDist, download_phyl_dist, load_phyl_dist from .prot_data_loader import load_prot_data -from .pubchem import find_met_smiles +from .pubchem import find_met_smiles, met_smiles from .uniprot_download import download_uniprot from .uniprot_loader import UniprotDB, load_uniprot_tsv @@ -41,4 +41,5 @@ "load_phyl_dist", "load_prot_data", "load_uniprot_tsv", + "met_smiles", ] diff --git a/src/geckopy/databases/pubchem.py b/src/geckopy/databases/pubchem.py index 3a933a4..edcee5a 100644 --- a/src/geckopy/databases/pubchem.py +++ b/src/geckopy/databases/pubchem.py @@ -42,6 +42,24 @@ _FIRST_LINE_PATTERN = re.compile(r"^(\S*)\n") +def met_smiles(metabolite) -> str: + """The metabolite's SMILES, whatever shape its annotation is in. + + cobra stores an annotation value as the object it was assigned, and + returns a *list* for a model read back from a file: a metabolite + annotated in memory by :func:`find_met_smiles` carries ``"CCO"``, + and the same metabolite after a save/load round trip carries + ``["CCO"]``. Readers that assume a string get a list that compares + equal to nothing and is not hashable. + + Returns the first non-empty entry, or ``""`` when there is none. + """ + value = metabolite.annotation.get("smiles", "") + if isinstance(value, (list, tuple, set)): + return next((str(v) for v in value if v), "") + return str(value) if value else "" + + def find_met_smiles( model: "cobra.Model", *, diff --git a/src/geckopy/gather_kcats/open_kinetics_predictor/parse_output.py b/src/geckopy/gather_kcats/open_kinetics_predictor/parse_output.py index fd811ef..55b8945 100644 --- a/src/geckopy/gather_kcats/open_kinetics_predictor/parse_output.py +++ b/src/geckopy/gather_kcats/open_kinetics_predictor/parse_output.py @@ -17,6 +17,8 @@ import pandas as pd +from ...databases.pubchem import met_smiles as read_met_smiles + if TYPE_CHECKING: from ...ec_model.ec_model import EcModel @@ -88,7 +90,7 @@ def parse_okp_output( seq_to_proteins = _index(model.ec.sequence) smiles_to_mets: dict[str, list] = {} for met in model.metabolites: - smiles = met.annotation.get("smiles", "") + smiles = read_met_smiles(met) if smiles: smiles_to_mets.setdefault(smiles, []).append(met) diff --git a/src/geckopy/gather_kcats/write_dlkcat_input.py b/src/geckopy/gather_kcats/write_dlkcat_input.py index b1700b0..ea0c1f8 100644 --- a/src/geckopy/gather_kcats/write_dlkcat_input.py +++ b/src/geckopy/gather_kcats/write_dlkcat_input.py @@ -14,6 +14,8 @@ import pandas as pd from scipy import sparse +from ..databases.pubchem import met_smiles as read_met_smiles + if TYPE_CHECKING: from ..databases.dlkcat_ignore_lists import DLKcatIgnoreLists from ..ec_model.ec_model import EcModel @@ -64,7 +66,8 @@ def write_dlkcat_input( SMILES are dropped from the output. Otherwise the SMILES is written as the literal string ``"None"``. - SMILES are read from each metabolite's ``annotation['smiles']``. + SMILES are read from each metabolite's ``annotation['smiles']``, + which cobra returns as a list for a model loaded from a file. Parameters ---------- @@ -194,7 +197,7 @@ def extract_enzyme_substrate_pairs( metabolites = list(model.metabolites) met_ids = [m.id for m in metabolites] met_names = [m.name for m in metabolites] - met_smiles = [m.annotation.get("smiles", "") for m in metabolites] + met_smiles = [read_met_smiles(m) for m in metabolites] met_normalized = [_normalize(n) for n in met_names] # Build the ignore mask. diff --git a/tests/test_open_kinetics_predictor.py b/tests/test_open_kinetics_predictor.py index 3077627..a34fd14 100644 --- a/tests/test_open_kinetics_predictor.py +++ b/tests/test_open_kinetics_predictor.py @@ -169,6 +169,22 @@ def test_parse_maps_back_to_ec_rxns(): assert set(df["rxn_id"]) == {"r1", "r2"} +def test_parse_maps_back_from_a_reloaded_models_list_valued_smiles(): + """cobra returns annotation values as lists for a model read from a + file, which is what `save_ec_model` + `load_ec_model` produces. The + SMILES index has to key on the string inside the list, not on the + list itself.""" + model = _ec_model() + for met in model.metabolites: + if met.annotation.get("smiles"): + met.annotation["smiles"] = [met.annotation["smiles"]] + + df = parse_okp_output(model, _OKP_RESULT_CSV) + assert len(df) == 3 + assert set(df["rxn_id"]) == {"r1", "r2"} + assert all(s in (["ala"], ["beta"]) for s in df["substrates"]) + + def test_parse_source_provenance_stripped_and_verbatim(): df = parse_okp_output(_ec_model(), _OKP_RESULT_CSV) sources = set(df["source"]) diff --git a/tests/test_pubchem.py b/tests/test_pubchem.py index 168f60a..4ba33d2 100644 --- a/tests/test_pubchem.py +++ b/tests/test_pubchem.py @@ -343,3 +343,28 @@ def test_logs_summary_after_run(monkeypatch, tmp_path, caplog): find_met_smiles(model, cache_path=cache_path) assert "SMILES found for" in caplog.text + + +# --------------------------------------------------------------------------- # +# Reading the annotation back +# --------------------------------------------------------------------------- # + +def test_met_smiles_accepts_the_shape_a_loaded_model_carries(): + """cobra returns a list for a model read from a file, a string for one + annotated in memory. Both mean the same molecule.""" + import cobra + + from geckopy.databases import met_smiles + + met = cobra.Metabolite("etoh") + met.annotation["smiles"] = "CCO" + assert met_smiles(met) == "CCO" + + met.annotation["smiles"] = ["CCO"] + assert met_smiles(met) == "CCO" + + met.annotation["smiles"] = [] + assert met_smiles(met) == "" + + del met.annotation["smiles"] + assert met_smiles(met) == "" diff --git a/tests/test_write_dlkcat_input.py b/tests/test_write_dlkcat_input.py index ea2e773..cf9593f 100644 --- a/tests/test_write_dlkcat_input.py +++ b/tests/test_write_dlkcat_input.py @@ -417,3 +417,40 @@ def test_written_file_has_no_header(tmp_path): first_line = out.read_text(encoding="utf-8").splitlines()[0] # Should be the data row, NOT a header. assert first_line.startswith("r1\t") + + +def test_a_reloaded_model_carries_its_smiles_as_a_list(tmp_path): + """cobra returns annotation values as lists for a model read from a + file, which is what `save_ec_model` + `load_ec_model` produces. The + SMILES has to survive that, and the ignore-list check has to not + choke on it.""" + model = _ec_model( + [("r1", [("A", -1.0, "alpha"), ("B", 1.0, "beta")])], + genes_per_rxn={"r1": ["g1"]}, + gene_sequences={"g1": "MASEQ"}, + met_smiles={"A": "C(C)O"}, + ) + for met in model.metabolites: + if met.annotation.get("smiles"): + met.annotation["smiles"] = [met.annotation["smiles"]] + + df = write_dlkcat_input(model, tmp_path / "DLKcat.tsv", _ignore_lists()) + assert len(df) == 1 + assert df.iloc[0]["smiles"] == "C(C)O" + + +def test_a_list_valued_smiles_is_still_matched_against_the_ignore_list(tmp_path): + model = _ec_model( + [("r1", [("A", -1.0, "alpha"), ("B", -1.0, "beta"), ("C", 1.0, "gamma")])], + genes_per_rxn={"r1": ["g1"]}, + gene_sequences={"g1": "MASEQ"}, + met_smiles={"A": "SMI_A", "B": "O"}, + ) + for met in model.metabolites: + if met.annotation.get("smiles"): + met.annotation["smiles"] = [met.annotation["smiles"]] + + df = write_dlkcat_input( + model, tmp_path / "DLKcat.tsv", _ignore_lists(ignore_smiles=["O"]), + ) + assert set(df["substrate"]) == {"alpha"}