From a037baa2bbac392f3989c78e24d66a0cdae531d3 Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven <eduardk@chalmers.se>
Date: Thu, 5 Nov 2026 12:44:26 +0100
Subject: [PATCH] Drop the complexes without an ortholog under the drop policy

The policy drop removes each complex that has a subunit without an ortholog and removes the reaction when no complex remains. It removed the reaction as soon as one subunit had no ortholog.
---
 CHANGELOG.md                                  |  8 ++
 .../homology/getModelFromHomology.m           | 86 +++++++++++++------
 tests/function_tests/tReconstruction.m        | 26 +++++-
 3 files changed, 90 insertions(+), 30 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index f20651e1..d5200eb8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,14 @@
 
 ## Unreleased
 
+### Breaking changes
+
+* `getModelFromHomology` with `complexPolicy` `drop` removes each complex that has a subunit without an ortholog, and removes the reaction only when no complex remains (was: the reaction was removed as soon as one subunit had no ortholog). For the rule `(A and B) or (C and D)` without an ortholog of `B`, the draft has the orthologs of `C` and `D`.
+
+### Fixes
+
+* `getModelFromHomology` reads the rule as a tree, in which `and` binds tighter than `or`, when it removes the placeholder genes of isozymes without an ortholog (was: it removed a placeholder that stood beside an `or` token, also when the placeholder belongs to a complex). A placeholder in an `or` group is removed when the group has a member that is not a placeholder, in nested groups as well. For the rule `(A or B and C) and D` without an ortholog of `B`, where `nA`, `nC` and `nD` are the orthologs of `A`, `C` and `D`, the draft rule is `(nA or (OLD_T_B and nC)) and nD` under `flag`, `(nA or nC) and nD` under `keep` and `nA and nD` under `drop`; it was `nA and nC and nD` under all three. The result differs from before only for rules with unbracketed `and` and `or` or with nested `or` groups.
+
 ### Data
 
 * The KEGG assets are the `kegg120` release of `toolbox-assets` (was: `kegg118`), and `fetchToolboxAsset` reads the manifest `manifest-v5`.
diff --git a/src/reconstruction/homology/getModelFromHomology.m b/src/reconstruction/homology/getModelFromHomology.m
index 41101910..4dba5084 100755
--- a/src/reconstruction/homology/getModelFromHomology.m
+++ b/src/reconstruction/homology/getModelFromHomology.m
@@ -61,8 +61,16 @@
 % complexPolicy : char
 %     what happens to an "and"-linked subunit without an ortholog (case
 %     insensitive): "flag" keeps the reaction with the subunit renamed
-%     OLD_MODELID_geneName, for review; "keep" drops that subunit;
-%     "drop" drops the whole reaction (default "flag").
+%     OLD_MODELID_geneName, for review; "keep" drops that subunit; "drop"
+%     drops each complex that has such a subunit and drops the reaction when
+%     no complex remains (default "flag"). With every policy, a group of
+%     isozymes ("or") counts as a missing subunit only when none of its members
+%     has an ortholog, and "and" binds tighter than "or". Under "drop", a
+%     complex is one "and"-linked alternative of the rule, as in its
+%     disjunctive normal form. The rule is not rewritten to its normal form,
+%     and a complex that contains another complex is kept. For the rule
+%     (A and B) or (C and D) without an ortholog of B, the draft has the
+%     orthologs of C and D.
 % keepGeneFree : logical
 %     also transfer the template reactions that have no grRule (spontaneous,
 %     transport, exchange, artificial). By default they are dropped, because
@@ -611,30 +619,19 @@
 %complexes
 draftModel=mergeModels(models,'metNames');
 
-%Remove OLD_ genes that appear only in OR branches (they were added for
-%AND-complex genes with no orthologs; in OR branches they are redundant).
-%Uses token-level removal to avoid partial-match hazards.
+%Remove OLD_ genes that are a member of an isozyme group (an "or" group) with
+%a member that is not an OLD_ gene: an isozyme without an ortholog is redundant
+%when another member of its group is present. When every member of a group is an
+%OLD_ gene, the last one stays for the group. The rule is parsed so that "and"
+%binds tighter than "or".
 for oldI=1:numel(draftModel.grRules)
     if ~contains(draftModel.grRules{oldI},'OLD_')
         continue
     end
     tkns=gprTokenize(draftModel.grRules{oldI});
-    again=true;
-    while again
-        again=false;
-        for kk=1:numel(tkns)
-            if startsWith(tkns{kk},'OLD_')
-                if kk>1 && strcmp(tkns{kk-1},'or')
-                    tkns(kk-1:kk)=[];
-                    again=true; break
-                elseif kk<numel(tkns) && strcmp(tkns{kk+1},'or')
-                    tkns(kk:kk+1)=[];
-                    again=true; break
-                end
-            end
-        end
-    end
-    draftModel.grRules{oldI}=gprBuildRule(tkns);
+    [node,~]=gprParseOr(tkns,1);
+    node=gprDropOldIsozymes(node);
+    draftModel.grRules{oldI}=gprRender(node,'');
 end
 
 %Apply the complex policy to the OLD_ genes left in "and" relations, and
@@ -646,7 +643,8 @@
         continue
     end
     if contains(rule,'OLD_') && strcmp(complexPolicy,'drop')
-        toRemove(j)=true;
+        draftModel.grRules{j}=gprSimplify(rule,'OLD_',true);
+        toRemove(j)=isempty(draftModel.grRules{j});
     elseif contains(rule,'OLD_') && strcmp(complexPolicy,'keep')
         draftModel.grRules{j}=gprSimplify(rule,'OLD_');
         toRemove(j)=isempty(draftModel.grRules{j});
@@ -717,6 +715,32 @@
 end
 
 
+function node=gprDropOldIsozymes(node)
+% Remove the OLD_ genes that are direct members of an "or" group together with
+% a member that is not an OLD_ gene. A group of OLD_ genes alone keeps its last
+% member.
+if strcmp(node.op,'gene')
+    return
+end
+for k=1:numel(node.kids)
+    node.kids{k}=gprDropOldIsozymes(node.kids{k});
+end
+if strcmp(node.op,'or')
+    isOld=cellfun(@(kid) strcmp(kid.op,'gene') && startsWith(kid.gene,'OLD_'),node.kids);
+    if any(isOld)
+        if all(isOld)
+            node.kids=node.kids(end);
+        else
+            node.kids=node.kids(~isOld);
+        end
+    end
+    if numel(node.kids)==1
+        node=node.kids{1};
+    end
+end
+end
+
+
 function tokens=gprReplaceToken(tokens,oldGene,newGenes)
 % Replace each occurrence of token oldGene with the gene(s) in newGenes.
 % newGenes is a cell array; multiple genes are wrapped in parentheses and
@@ -746,18 +770,22 @@
 end
 
 
-function rule=gprSimplify(rule,removePrefix)
+function rule=gprSimplify(rule,removePrefix,dropComplex)
 % Rewrite a grRule as a tree: genes starting with removePrefix (if any) are
 % dropped from their group, groups of the same operator are flattened,
 % repeated genes or groups are removed, and brackets appear only around a
-% group of the other operator. A rule with no gene left becomes ''.
+% group of the other operator. With dropComplex, an "and" group that loses
+% a gene is removed as a whole. A rule with no gene left becomes ''.
+if nargin<3
+    dropComplex=false;
+end
 tokens=gprTokenize(rule);
 if isempty(tokens)
     rule='';
     return
 end
 [node,~]=gprParseOr(tokens,1);
-node=gprPrune(node,removePrefix);
+node=gprPrune(node,removePrefix,dropComplex);
 if isempty(node)
     rule='';
 else
@@ -809,7 +837,7 @@
 end
 
 
-function node=gprPrune(node,removePrefix)
+function node=gprPrune(node,removePrefix,dropComplex)
 if strcmp(node.op,'gene')
     if ~isempty(removePrefix) && startsWith(node.gene,removePrefix)
         node=[];
@@ -819,8 +847,12 @@
 kids={};
 seen={};
 for k=1:numel(node.kids)
-    kid=gprPrune(node.kids{k},removePrefix);
+    kid=gprPrune(node.kids{k},removePrefix,dropComplex);
     if isempty(kid)
+        if dropComplex && strcmp(node.op,'and')
+            node=[];
+            return
+        end
         continue
     end
     if strcmp(kid.op,node.op)
diff --git a/tests/function_tests/tReconstruction.m b/tests/function_tests/tReconstruction.m
index 53eb1742..494b5aaf 100644
--- a/tests/function_tests/tReconstruction.m
+++ b/tests/function_tests/tReconstruction.m
@@ -87,13 +87,19 @@ function getModelFromHomologyComplexPolicyAndOptions(testCase)
 
         function getModelFromHomologyComplexPolicyOnNestedRules(testCase)
             % An isozyme group whose members all lack an ortholog leaves one flagged
-            % subunit; "drop" removes every reaction that keeps a flagged gene.
+            % subunit; "drop" removes each complex that keeps a flagged gene and the
+            % reaction when no complex remains.
             % Columns: rule, mapped template genes, policy, genes in the rule ({} = reaction removed).
             cases = {'(A or B) and C', {'C'}, 'flag', {'OLD_T_B', 'nC'}; ...
                      '(A or B) and C', {'C'}, 'keep', {'nC'}; ...
                      '(A or B) and C', {'C'}, 'drop', {}; ...
-                     'A or (B and C)', {'A', 'C'}, 'drop', {}; ...
-                     '(A and B) or C', {'C'}, 'drop', {}};
+                     'A or (B and C)', {'A', 'C'}, 'drop', {'nA'}; ...
+                     '(A and B) or C', {'C'}, 'drop', {'nC'}; ...
+                     'A and B', {'A'}, 'drop', {}; ...
+                     '(A and B) or (C and D)', {'A', 'C', 'D'}, 'drop', {'nC', 'nD'}; ...
+                     '(A or B and C) and D', {'A', 'C', 'D'}, 'drop', {'nA', 'nD'}; ...
+                     '(A and B or C) and D', {'A', 'C', 'D'}, 'drop', {'nC', 'nD'}; ...
+                     'D and (A or B and C)', {'A', 'C', 'D'}, 'drop', {'nA', 'nD'}};
             for k = 1:size(cases, 1)
                 t = tReconstruction.homologyTemplate('T', {'R1'}, cases(k, 1));
                 pairs = cell(numel(cases{k, 2}), 2);
@@ -108,6 +114,20 @@ function getModelFromHomologyComplexPolicyOnNestedRules(testCase)
                     testCase.verifyFalse(contains(draft.grRules{strcmp(draft.rxns, 'R1')}, ' or '), cases{k, 1});
                 end
             end
+            % "and" binds tighter than "or": B belongs to the complex B and C, not to the isozymes of A
+            t = tReconstruction.homologyTemplate('T', {'R1'}, {'(A or B and C) and D'});
+            draft = tReconstruction.homologyDraft({t}, {{'A', 'nA'; 'C', 'nC'; 'D', 'nD'}}, 'flag', {});
+            testCase.verifyEqual(draft.grRules{1}, '(nA or (OLD_T_B and nC)) and nD');
+            % complete complexes that share a subunit with an incomplete one stay, and a
+            % complex that contains another is not merged into it
+            t = tReconstruction.homologyTemplate('T', {'R1'}, {'(A and B) or (A and C) or (D and E) or (A and B and F)'});
+            draft = tReconstruction.homologyDraft({t}, {{'A', 'nA'; 'B', 'nB'; 'C', 'nC'; 'D', 'nD'; 'F', 'nF'}}, 'drop', {});
+            testCase.verifyEqual(draft.grRules{1}, '(nA and nB) or (nA and nC) or (nA and nB and nF)');
+            t = tReconstruction.homologyTemplate('T', {'R1'}, {'(A and B) or B'});
+            draft = tReconstruction.homologyDraft({t}, {{'A', 'nA'; 'B', 'nB'}}, 'drop', {});
+            testCase.verifyEqual(draft.grRules{1}, '(nA and nB) or nB');
+            draft = tReconstruction.homologyDraft({t}, {{'B', 'nB'}}, 'drop', {});
+            testCase.verifyEqual(draft.grRules{1}, 'nB');
             % the flagged subunit stays next to its isozyme in the other branch
             t = tReconstruction.homologyTemplate('T', {'R1'}, {'A or (B and C)'});
             draft = tReconstruction.homologyDraft({t}, {{'A', 'nA'; 'C', 'nC'}}, 'flag', {});