Skip to content

fix: parse grRules into a tree; fix removeGenes, expandModel, findPotentialErrors #673

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

Merged
edkerk merged 1 commit into from
Jul 17, 2026
Merged

fix: parse grRules into a tree; fix removeGenes, expandModel, findPotentialErrors #673

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
118 changes: 64 additions & 54 deletions manipulation/expandModel.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,36 +15,75 @@
% model structure with separate reactions for iso-enzymes, where the
% reaction ids are renamed as id_EXP_1, id_EXP_2, etc.
% rxnToCheck : cell
% cell array with original reaction identifiers for those that
% contained nested and/or-relationships in grRules.
% cell array with original reaction identifiers for those whose grRule
% was not in disjunctive normal form, and so had to be expanded by
% distributivity, plus any whose grRule could not be parsed.
%
% Examples
% --------
% [newModel, rxnToCheck]=expandModel(model);
%
% Notes
% -----
% grRules strings that involve nested expressions of 'and' and 'or' might
% not be parsed correctly if they are not standardized (if the
% standardizeGrRules function was not first run on the model). For those
% reactions, it is therefore advisable to inspect the reactions in
% rxnToCheck to confirm correct model expansion.
% The isozymes are the AND-clauses of the grRule in disjunctive normal form,
% so nested expressions are expanded correctly: "g1 and (g2 or g3)" yields
% "g1 and g2" and "g1 and g3". Reactions listed in rxnToCheck are expanded
% correctly too; they are reported only because a rule needing
% distributivity is often a sign that the GPR was not what its author
% intended. Use standardizeGrRules/findPotentialErrors to inspect them.

%Check how many reactions we will create (the number of or:s in the GPRs).
%This way, we can preallocate all fields and save much computation time

numOrs = count(model.grRules, ' or ');
toAdd = sum(numOrs);
%Work out the isozymes of every reaction up front. The number of copies is
%the number of DNF clauses, which is not the number of ' or ' substrings:
%"(g1 or g2) and (g3 or g4)" has two or:s but four isozymes.
prevNumRxns = length(model.rxns);
clauses = cell(prevNumRxns,1);
nCopies = zeros(prevNumRxns,1);
rxnToCheck={};
for i=1:prevNumRxns
rule = model.grRules{i};
if isempty(rule)
continue
end
try
c = grRuleToDNF(rule);
catch ME
if ~strcmp(ME.identifier,'RAVEN:badGrRule')
rethrow(ME)
end
%An unparseable rule is left untouched rather than guessed at.
rxnToCheck{end+1,1}=model.rxns{i}; %#ok<AGROW>
continue
end
if numel(c) <= 1
continue
end
%DEVIATION from raven-toolbox expand.py, which has no such guard: the
%cross-product is exponential in the number of ORed complexes, and
%MATLAB grows the model fields eagerly, so a pathological rule would
%exhaust memory instead of raising. Fail with the culprit named.
if numel(c) > 10000
error('RAVEN:grRuleTooComplex',['Reaction ' model.rxns{i} ' has a grRule ' ...
'with ' num2str(numel(c)) ' isozymes after expansion. This is almost ' ...
'certainly a malformed rule; check it with findPotentialErrors.']);
end
clauses{i} = c;
nCopies(i) = numel(c)-1;
if ~isDnfGrRule(rule)
rxnToCheck{end+1,1}=model.rxns{i}; %#ok<AGROW>
end
end
toAdd = sum(nCopies);
if toAdd > 0
%Calculate indices to expand
%For example, if a reaction with index x has 2 or:s, meaning it has 3
%reactions after the split, we should add two copies of this reaction
%For fields that should just be copied to the new reactions, we just keep
%track of that there are two copies, i.e., we add x x to this vector.
%That is exactly what repelem does for us.
cpyIndices = repelem(1:prevNumRxns, numOrs);
%(:) forces a column. Indexing a 1x1 field takes the orientation of the
%index rather than of the field, so a row here made every concatenation
%below fail on a single-reaction model.
cpyIndices = reshape(repelem((1:prevNumRxns)', nCopies), [], 1);

%Copy all fields that should just be copied
model.S=[model.S model.S(:,cpyIndices)];
Expand Down Expand Up @@ -92,54 +131,25 @@
%Loop throught those reactions and fill in the expanded data
nextIndex = prevNumRxns + 1;
for i=1:prevNumRxns
if (numOrs(i) > 0)
%Check that it doesn't contain nested 'and' and 'or' relations and
%print a warning if it does
if ~isempty(strfind(model.grRules{i},' and '))
rxnToCheck{end+1,1}=model.rxns{i};
end

%Get rid of all '(' and ')' since I'm not looking at complex stuff
%anyways
geneString=model.grRules{i};
geneString=strrep(geneString,'(','');
geneString=strrep(geneString,')','');
geneString=strrep(geneString,' or ',';');
if (nCopies(i) > 0)
geneSets = clauses{i};

%Split the string into gene names
geneNames=regexp(geneString,';','split');

%Update the reaction to only use the first gene
model.grRules{i}=['(' geneNames{1} ')'];
%Find the gene in the gene list If ' and ' relationship, first
%split the genes
%The first isozyme replaces the original reaction
model.grRules{i}=strjoin(geneSets{1},' and ');
model.rxnGeneMat(i,:)=0;
if ~isempty(strfind(geneNames(1),' and '))
andGenes=regexp(geneNames{1},' and ','split');
model.rxnGeneMat(i,ismember(model.genes,andGenes)) = 1;
else
[~, index]=ismember(geneNames(1),model.genes);
model.rxnGeneMat(i,index)=1;
end
model.rxnGeneMat(i,ismember(model.genes,geneSets{1})) = 1;

%Insert the reactions at the end of the model and without
%allocating space. This is not nice, but ok for now
for j=2:numel(geneNames)
%The rest are appended, named after the original reaction, which
%is why it is only renamed once the loop is done
for j=2:numel(geneSets)
ind = nextIndex+j-2;
model.rxns{ind}=[model.rxns{i} '_EXP_' num2str(j)];

model.grRules{ind}=['(' geneNames{j} ')'];

if ~isempty(strfind(geneNames(j),' and '))
andGenes=regexp(geneNames{j},' and ','split');
model.rxnGeneMat(ind,ismember(model.genes,andGenes)) = 1;
else
model.rxnGeneMat(ind,ismember(model.genes,geneNames(j))) = 1;
end
model.grRules{ind}=strjoin(geneSets{j},' and ');
model.rxnGeneMat(ind,ismember(model.genes,geneSets{j})) = 1;
end
model.rxns{i}=[model.rxns{i}, '_EXP_1'];
nextIndex = nextIndex + numOrs(i);
end
nextIndex = nextIndex + nCopies(i);
end
end
newModel=model;
else
Expand Down
100 changes: 25 additions & 75 deletions manipulation/findPotentialErrors.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,15 +25,18 @@
%
% Notes
% -----
% Detection uses a parse-tree walk rather than substring search, so gene
% IDs that contain substrings such as "and" or "or" do not cause false
% positives or negatives.
% Detection walks the parse tree from parseGrRule rather than searching the
% raw string, so gene IDs that contain substrings such as "and" or "or" do
% not cause false positives or negatives.
%
% A rule is non-DNF when the top-level operator is AND and one of its
% operands contains OR, or equivalently when any OR subexpression appears
% inside a bracket group that is joined to another group by AND.
% A rule is non-DNF when an AND operator has an OR anywhere beneath it.
% Bracketing alone never makes a rule non-DNF: "((G1 and G2) or G3)" and
% "(G1 or G2) or (G3 and G4)" are both fine.
% Use standardizeGrRules to attempt automatic repair.
%
% Rules that cannot be parsed at all are reported as issues too, with the
% parser's message as the reason.
%
% Examples
% --------
% issues = findPotentialErrors(model);
Expand All @@ -53,36 +56,26 @@
if isempty(rule)
continue
end
% Normalize operator keywords to single-character symbols so that gene
% IDs containing " and " or " or " as substrings are not misread.
ruleN=lower(rule);
ruleN=regexprep(ruleN,' and ',' & ');
ruleN=regexprep(ruleN,' or ',' | ');

% Rules with only one operator type are always DNF.
hasAnd=contains(ruleN,' & ');
hasOr =contains(ruleN,' | ');
if ~hasAnd || ~hasOr
continue
end

% Split at top-level | (OR) to get each conjunction term.
terms=splitTopLevel(ruleN,'|');
nonDnf=false;
for j=1:numel(terms)
term=strtrim(terms{j});
% Strip a single pair of enclosing parens if present.
if ~isempty(term) && term(1)=='(' && matchClose(term,1)==numel(term)
term=term(2:end-1);
end
% If this AND-term still contains an OR operator the rule is non-DNF.
if contains(term,' | ')
nonDnf=true;
break
try
isDnf=isDnfGrRule(rule);
catch ME
if ~strcmp(ME.identifier,'RAVEN:badGrRule')
rethrow(ME)
end
% A rule that cannot be parsed is itself worth reporting: it is
% certainly not usable by the isoenzyme/complex reasoning, and
% silently skipping it would hide the very thing this function
% exists to surface.
issues(end+1,1)=struct( ...
'index', i, ...
'rxn', model.rxns{i}, ...
'grRule', rule, ...
'reason', ['Cannot be parsed: ' ME.message]); %#ok<AGROW>
continue
end

if nonDnf
if ~isDnf
issues(end+1,1)=struct( ...
'index', i, ...
'rxn', model.rxns{i}, ...
Expand All @@ -91,46 +84,3 @@
end
end
end


function parts=splitTopLevel(str,sep)
% Split str on sep only at bracket depth zero.
parts={};
depth=0;
start=1;
n=numel(str);
sn=numel(sep);
k=1;
while k<=n
c=str(k);
if c=='('
depth=depth+1;
elseif c==')'
depth=depth-1;
elseif depth==0 && k+sn-1<=n && strcmp(str(k:k+sn-1),sep)
parts{end+1}=str(start:k-1); %#ok<AGROW>
start=k+sn;
k=k+sn-1;
end
k=k+1;
end
parts{end+1}=str(start:end);
end


function pos=matchClose(str,openPos)
% Return the position of the closing ")" that matches str(openPos).
depth=0;
pos=-1;
for k=openPos:numel(str)
if str(k)=='('
depth=depth+1;
elseif str(k)==')'
depth=depth-1;
if depth==0
pos=k;
return
end
end
end
end
94 changes: 67 additions & 27 deletions manipulation/removeGenes.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,10 +75,9 @@
for j = 1:numel(geneRxns)
index = geneRxns(j);
grRule = reducedModel.grRules{index};
ruleGenes = reducedModel.genes(logical(rxnGeneMat(index,:)));
if ~ismember(index,toCheck) && canCarryFlux(index) && ~isempty(grRule)
%Check if rxn can carry flux without this gene:
canCarryFlux(index) = canRxnCarryFlux(ruleGenes,grRule,genes{i});
canCarryFlux(index) = canRxnCarryFlux(grRule,genes{i});
%Adapt gene rule & gene matrix:
grRule = removeGeneFromRule(grRule,genes{i});
reducedModel.grRules{index} = grRule;
Expand All @@ -105,33 +104,74 @@
end
end

function canIt = canRxnCarryFlux(ruleGenes,geneRule,geneToRemove)
%This function converts a gene rule to a logical statement, and then
%asseses if the rule is true (i.e. rxn can still carry flux) or not (cannot
%carry flux).
geneRule = [' ', geneRule, ' '];
for i = 1:length(ruleGenes)
if strcmp(ruleGenes{i},geneToRemove)
geneRule = strrep(geneRule,[' ' ruleGenes{i} ' '],' false ');
geneRule = strrep(geneRule,['(' ruleGenes{i} ' '],'(false ');
geneRule = strrep(geneRule,[' ' ruleGenes{i} ')'],' false)');
else
geneRule = strrep(geneRule,[' ' ruleGenes{i} ' '],' true ');
geneRule = strrep(geneRule,['(' ruleGenes{i} ' '],'(true ');
geneRule = strrep(geneRule,[' ' ruleGenes{i} ')'],' true)');
end
function canIt = canRxnCarryFlux(geneRule,geneToRemove)
%Evaluate the rule with geneToRemove absent and every other gene present.
%A complex (AND) needs all of its subunits, isozymes (OR) need only one.
tree = parseGrRule(geneRule);
if isempty(tree)
canIt = true;
return
end
canIt = evalWithout(tree,geneToRemove);
end

function tf = evalWithout(node,geneToRemove)
switch node.type
case 'gene'
tf = ~strcmp(node.id,geneToRemove);
case 'and'
tf = all(cellfun(@(c) evalWithout(c,geneToRemove),node.children));
case 'or'
tf = any(cellfun(@(c) evalWithout(c,geneToRemove),node.children));
otherwise
error('RAVEN:badGrRule',['Unexpected grRule node type: ' node.type]);
end
geneRule = strtrim(geneRule);
geneRule = strrep(geneRule,'and','&&');
geneRule = strrep(geneRule,'or','||');
canIt = eval(geneRule);
end

function geneRule = removeGeneFromRule(geneRule,geneToRemove)
%This function receives a standard gene rule and it returns it without the
%chosen gene.
geneSets = strsplit(geneRule,' or ');
hasGene = ~cellfun(@isempty,strfind(geneSets,geneToRemove));
geneSets = geneSets(~hasGene);
geneRule = strjoin(geneSets,' or ');
%This function receives a gene rule and it returns it without the chosen
%gene. A complex that loses a subunit is dropped whole; an isozyme is
%dropped from its OR without disturbing the alternatives.
tree = parseGrRule(geneRule);
geneRule = grRuleToString(pruneGene(tree,geneToRemove));
end

function node = pruneGene(node,geneToRemove)
%Remove geneToRemove from the tree, returning [] when nothing is left that
%can catalyse the reaction.
if isempty(node)
node = [];
return
end
switch node.type
case 'gene'
if strcmp(node.id,geneToRemove)
node = [];
end
case 'and'
%A complex missing any subunit cannot form at all.
for k = 1:numel(node.children)
if isempty(pruneGene(node.children{k},geneToRemove))
node = [];
return
end
end
case 'or'
kept = {};
for k = 1:numel(node.children)
child = pruneGene(node.children{k},geneToRemove);
if ~isempty(child)
kept{end+1} = child; %#ok<AGROW>
end
end
if isempty(kept)
node = [];
elseif isscalar(kept)
node = kept{1};
else
node.children = kept;
end
otherwise
error('RAVEN:badGrRule',['Unexpected grRule node type: ' node.type]);
end
end
Loading