From a037baa2bbac392f3989c78e24d66a0cdae531d3 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 5 Nov 2026 12:44:26 +0100 Subject: [PATCH 1/2] 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 Date: Sat, 7 Nov 2026 20:31:05 +0100 Subject: [PATCH 2/2] Treat a gene taken by an earlier model as a subunit without an ortholog In the preferred-order step, a gene that maps to a new gene an earlier model already maps loses its mappings and follows complexPolicy. For the rule (a and b) or x with a taken, flag gives (OLD_MODELID_a and nb) or nx. A reaction with no gene that keeps a mapping is removed, and a model without reactions is left out of the merge. --- CHANGELOG.md | 1 + .../homology/getModelFromHomology.m | 42 ++++++++++----- tests/function_tests/tReconstruction.m | 53 +++++++++++++++++++ 3 files changed, 84 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d5200eb8..6892c2d4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ ### Breaking changes +* `getModelFromHomology` with `preferredOrder` treats a gene of a later model that maps to a new gene an earlier model already maps as a subunit without an ortholog, and `complexPolicy` decides what happens to it (was: the gene and every part of the rule that needs it were removed). For the rule `(a and b) or x` with `a` taken by an earlier model, `flag` gives `(OLD_MODELID_a and nb) or nx`, `keep` gives `nb or nx` and `drop` gives `nx`. A reaction that has genes, none of which keeps a mapping, is not included. * `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 diff --git a/src/reconstruction/homology/getModelFromHomology.m b/src/reconstruction/homology/getModelFromHomology.m index 4dba5084..2768b5fb 100755 --- a/src/reconstruction/homology/getModelFromHomology.m +++ b/src/reconstruction/homology/getModelFromHomology.m @@ -24,7 +24,12 @@ % the order in which reactions are added from the models, as a cell array % of model IDs. If not supplied, reactions are included from all models, % otherwise one gene results in reactions from one model only -% (default {}). +% (default {}). A gene of a later model that maps to a new gene which an +% earlier model already maps loses its mappings and counts as a subunit +% without an ortholog, so complexPolicy decides what happens to it. For +% the rule (a and b) or x with a taken by an earlier model, "flag" gives +% (OLD_MODELID_a and nb) or nx, "keep" gives nb or nx and "drop" gives nx. +% A reaction that has genes, none of which keeps a mapping, is not included. % strictness : double % integer that specifies which reactions are included (default 1): % @@ -503,16 +508,24 @@ [~, genesToDelete]=find(finalMappings{i}(usedGenes(a),:)); %IMPORTANT! IS it really correct to remove all genes that map? genesToDelete=unique(genesToDelete); %Maybe not needed, but for clarity if nothing else - %Remove all the genes that were already found and add the other - %ones to allUsedGenes - if ~isempty(genesToDelete) - models{useOrderIndexes(i)}=removeGenes(models{useOrderIndexes(i)},allGenes{i+1}(genesToDelete),true,true,false); - end + %The genes that were already found lose their mappings and are + %treated as genes without an ortholog. Add the other ones to + %allUsedGenes allUsedGenes(usedGenes)=true; %Remove the deleted genes from finalMappings and allGenes. finalMappings{i}(:,genesToDelete)=[]; allGenes{i+1}(genesToDelete)=[]; + + %Remove the reactions that have no gene with a mapping left + if ~isempty(genesToDelete) + a=ismember(models{useOrderIndexes(i)}.genes,allGenes{i+1}); + [rxnsToKeep,~] = find(models{useOrderIndexes(i)}.rxnGeneMat(:,a)); + rxnsToRemove = repmat(1,numel(models{useOrderIndexes(i)}.rxns),1); + rxnsToRemove(rxnsToKeep) = 0; + rxnsToRemove(ismember(models{useOrderIndexes(i)}.rxns,geneFreeRxns{useOrderIndexes(i)})) = 0; + models{useOrderIndexes(i)}=removeReactions(models{useOrderIndexes(i)},find(rxnsToRemove),true,true,true); + end end end @@ -606,7 +619,7 @@ else models{useOrderIndexes(i)}.genes=allGenes{1}(unique(newGenes)); end - if isfield(models{useOrderIndexes(i)},'geneComps') + if isfield(models{useOrderIndexes(i)},'geneComps') && ~isempty(models{useOrderIndexes(i)}.geneComps) geneComps=models{useOrderIndexes(i)}.geneComps(1); models{useOrderIndexes(i)}.geneComps=zeros(numel(models{useOrderIndexes(i)}.genes),1); %Assume that all genes are in the same compartment, and this @@ -616,7 +629,12 @@ end %Now merge the models. All information should be correct except for "or" -%complexes +%complexes. A model that has lost all its reactions adds nothing to the draft +modelIds=cellfun(@(m) m.id,models,'UniformOutput',false); +hasRxns=cellfun(@(m) ~isempty(m.rxns),models); +if any(hasRxns) + models=models(hasRxns); +end draftModel=mergeModels(models,'metNames'); %Remove OLD_ genes that are a member of an isozyme group (an "or" group) with @@ -659,11 +677,11 @@ %Change name of the resulting model draftModel.id=getModelFor; name='Generated by getModelFromHomology using '; -for i=1:numel(models) - if i