Skip to content

Drop each complex without an ortholog under the drop policy #25

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
edkerk wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Drop each complex without an ortholog under the drop policy #25

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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`.
Expand Down
86 changes: 59 additions & 27 deletions src/reconstruction/homology/getModelFromHomology.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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});
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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=[];
Expand All @@ -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)
Expand Down
26 changes: 23 additions & 3 deletions tests/function_tests/tReconstruction.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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', {});
Expand Down
Loading