diff --git a/src/raven_toolbox/reconstruction/kegg/assemble.py b/src/raven_toolbox/reconstruction/kegg/assemble.py index a2b5eb9..8a0025e 100644 --- a/src/raven_toolbox/reconstruction/kegg/assemble.py +++ b/src/raven_toolbox/reconstruction/kegg/assemble.py @@ -35,11 +35,17 @@ def assemble_model_from_ko_genes( model_id: str | None = None, model_name: str | None = None, note: str | None = None, + prune_orthology: bool = False, ) -> tuple[cobra.Model, dict[str, list[str]]]: """Build a draft model from a ``{ko: [gene, ...]}`` assignment. Returns ``(model, gpr_map)`` where ``gpr_map`` is the kept reactions' gene lists, so callers can add gene annotations afterwards. + + ``prune_orthology`` restricts each gene-backed reaction's ``kegg.orthology`` + annotation to the KOs that actually contributed a gene, matching RAVEN + ``getKEGGModelForOrganism``'s HMM branch (the organism-annotation path keeps + the full reference KO list, so this is off by default). """ rxn_to_kos: dict[str, set[str]] = {} for ko, rid in zip(ko_reaction["ko"], ko_reaction["reaction"], strict=True): @@ -53,15 +59,19 @@ def assemble_model_from_ko_genes( } gpr_map: dict[str, list[str]] = {} + matched_kos: dict[str, set[str]] = {} spontaneous_kept: set[str] = set() for rxn in reference_model.reactions: rid = rxn.id # Quality filters first: dropped even if it would have genes. if any(not keep_flag and rid in flagged for keep_flag, flagged in drop_if.values()): continue - genes = sorted({g for ko in rxn_to_kos.get(rid, ()) for g in ko_to_genes.get(ko, ())}) + rxn_kos = rxn_to_kos.get(rid, ()) + genes = sorted({g for ko in rxn_kos for g in ko_to_genes.get(ko, ())}) if genes: gpr_map[rid] = genes + # KOs that actually contributed a gene, for optional annotation pruning. + matched_kos[rid] = {ko for ko in rxn_kos if ko_to_genes.get(ko)} elif rid in spontaneous and keep_spontaneous: spontaneous_kept.add(rid) @@ -76,6 +86,20 @@ def assemble_model_from_ko_genes( ) for rid, genes in gpr_map.items(): model.reactions.get_by_id(rid).gene_reaction_rule = " or ".join(genes) + if prune_orthology: + # Restrict each gene-backed reaction's kegg.orthology to the KOs that + # matched a gene, preserving the reference annotation's order. (The + # order-preserving intersection also avoids the index misalignment in + # RAVEN getKEGGModelForOrganism's rxnMiriams pruning.) + for rid, kos in matched_kos.items(): + annotation = model.reactions.get_by_id(rid).annotation + orthology = annotation.get("kegg.orthology") + if orthology is None: + continue + if isinstance(orthology, list): + annotation["kegg.orthology"] = [k for k in orthology if k in kos] + elif orthology not in kos: + annotation["kegg.orthology"] = [] if note is not None: for rid in keep: model.reactions.get_by_id(rid).notes["note"] = note diff --git a/src/raven_toolbox/reconstruction/kegg/query.py b/src/raven_toolbox/reconstruction/kegg/query.py index bbb3ec6..2de4e04 100644 --- a/src/raven_toolbox/reconstruction/kegg/query.py +++ b/src/raven_toolbox/reconstruction/kegg/query.py @@ -31,7 +31,7 @@ from raven_toolbox.reconstruction.kegg.assemble import assemble_model_from_ko_genes from raven_toolbox.reconstruction.kegg.parse import _resolve_artefact, read_kegg_table -_NOTE = "Included by get_kegg_model_from_sequences (using HMMs)" +_NOTE = "Included by KEGG HMM reconstruction" _MIN_EVALUE = 1e-250 # floor for a reported E-value of 0, to keep logs finite @@ -198,6 +198,11 @@ def get_kegg_model_from_sequences( (:func:`assign_kos`), and assembles the model against ``reference_model`` / ``ko_reaction``. Genes are the query proteome's identifiers. """ + if model_id is None: + # RAVEN always sets model.id (to organismID); default it here so the + # draft never inherits the reference model's id. Pass model_id explicitly + # for byte-identical parity with a specific MATLAB organismID. + model_id = Path(fasta).stem hits = parse_hmmsearch_tblout(run_hmmsearch(fasta, library, threads=threads, hmmsearch=hmmsearch)) ko_to_genes = assign_kos( hits, @@ -216,6 +221,7 @@ def get_kegg_model_from_sequences( keep_general=keep_general, model_id=model_id, note=_NOTE, + prune_orthology=True, ) return model diff --git a/tests/test_reconstruction_kegg_query.py b/tests/test_reconstruction_kegg_query.py index 0dd23f0..65c981f 100644 --- a/tests/test_reconstruction_kegg_query.py +++ b/tests/test_reconstruction_kegg_query.py @@ -127,6 +127,50 @@ def test_get_model_from_sequences(reference_and_tables, monkeypatch): assert model.id == "myorg" r = model.reactions.get_by_id("R90010") assert set(r.gene_reaction_rule.split(" or ")) == {"myGeneA", "myGeneB"} - assert r.notes["note"].endswith("(using HMMs)") + assert r.notes["note"] == "Included by KEGG HMM reconstruction" + # Its single KO matched, so the kegg.orthology annotation is unchanged. + assert r.annotation["kegg.orthology"] == ["K90001"] # R90200/R90300 had no matched KOs and are not spontaneous -> absent. assert "R90200" not in model.reactions + + +def test_model_id_defaults_to_fasta_stem(reference_and_tables, monkeypatch): + """RAVEN always sets model.id; with no model_id we default it to the FASTA stem + rather than inheriting the reference model's id.""" + model_ref, tables = reference_and_tables + monkeypatch.setattr( + "raven_toolbox.reconstruction.kegg.query.run_hmmsearch", + lambda *a, **k: "myGeneA - K90001 - 1e-120 400 0\n", + ) + model = get_kegg_model_from_sequences( + "/some/path/eco.faa", model_ref, tables["ko_reaction"], "ignored.hmm", + rxn_flags=tables["rxn_flags"], + ) + assert model.id == "eco" + + +def test_prune_orthology_keeps_only_matched_kos(): + """The FASTA path prunes a kept reaction's kegg.orthology to the KOs that + matched a gene (RAVEN getKEGGModelForOrganism HMM branch), preserving order.""" + import cobra + + from raven_toolbox.reconstruction.kegg.assemble import assemble_model_from_ko_genes + + ref = cobra.Model("ref") + met = cobra.Metabolite("C1", compartment="s") + rxn = cobra.Reaction("R1") + rxn.add_metabolites({met: -1.0}) + rxn.annotation["kegg.orthology"] = ["K1", "K2", "K3"] + ref.add_reactions([rxn]) + ko_reaction = pd.DataFrame( + [("K1", "R1"), ("K2", "R1"), ("K3", "R1")], columns=["ko", "reaction"] + ) + ko_to_genes = {"K2": ["g"]} # only K2 matched a gene + + pruned, _ = assemble_model_from_ko_genes( + ref, ko_reaction, ko_to_genes, prune_orthology=True + ) + assert pruned.reactions.get_by_id("R1").annotation["kegg.orthology"] == ["K2"] + # Default (organism path) keeps the full reference KO list. + full, _ = assemble_model_from_ko_genes(ref, ko_reaction, ko_to_genes) + assert full.reactions.get_by_id("R1").annotation["kegg.orthology"] == ["K1", "K2", "K3"]