KEGG FASTA/HMM reconstruction: align model output with MATLAB RAVEN

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

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
 

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"]