From 14e8e2cc846b6fff80f61c4409324d4a52b8f2ab Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 17 Jul 2026 11:04:51 +0200 Subject: [PATCH] Parse grRules into a tree; fix removeGenes, expandModel, findPotentialErrors Three functions each reasoned about grRules by searching the raw string, and each got it wrong in a different way. They now share one parse tree. New in utils/: parseGrRule (tokenise + recursive descent), isDnfGrRule, grRuleToDNF, grRuleToString. Operators are matched as whole words, so a gene called RAND1 is a gene, and "10" is never confused with "100". raven-toolbox gets this for free from cobra's Python AST (utils/gpr.py, manipulation/expand.py); MATLAB has no equivalent, so the tree is built here. isDnfGrRule mirrors _is_dnf_node/_contains_or and grRuleToDNF mirrors _node_to_dnf node for node. removeGenes (R2): removeGeneFromRule used an unanchored strfind over the ' or '-split rule, so removing gene "10" from "10 or 100" also matched "100" and left the reaction with no grRule at all while gene 100 stayed in the model -- the association was destroyed, silently, and standardizeGrRules then zeroed the rxnGeneMat row. Numeric Entrez/HGNC ids make this the norm rather than the exception in the tINIT/HPA human workflow, and findGeneDeletions and runSimpleOptKnock mispredict essentiality off the result. Both subfunctions now walk the tree, comparing gene ids with strcmp. This also retires the last eval-based grRule evaluation in the codebase. expandModel (EX1): stripped every parenthesis and split on ' or ', which destroys precedence. "g1 and (g2 or g3)" expanded to "g1 and g2" and "g3" -- the second isozyme silently lost its mandatory g1 subunit, so a g1 knockout was predicted viable when it is lethal. Isozymes are now the DNF clauses. The reaction count comes from the clauses too: "(g1 or g2) and (g3 or g4)" has two or:s but four isozymes, and the old preallocation sized for three. rxnToCheck now lists only rules that genuinely needed distributivity, not every rule containing ' and '. findPotentialErrors (S2): the one-level split flagged valid DNF such as "((G1 and G2) or G3)" and "(G1 or G2) or (G3 and G4)". That was worse than cosmetic: standardizeGrRules skips flagged rules, so a false positive meant the rule was never repaired and kept the brackets the function exists to remove. Unparseable rules are now reported as issues rather than skipped. Also fix a pre-existing crash in expandModel: cpyIndices was a row, and indexing a 1x1 field takes the orientation of the index, so every concatenation failed on a single-reaction model. Tests: nothing under testing/ exercised grRules at all. Adds parser coverage to tUtils and regression tests to tManipulation for each defect above, using the reproductions as the anchors. --- manipulation/expandModel.m | 118 +++++++++-------- manipulation/findPotentialErrors.m | 100 ++++---------- manipulation/removeGenes.m | 94 +++++++++---- testing/function_tests/tManipulation.m | 101 ++++++++++++++ testing/function_tests/tUtils.m | 70 ++++++++++ utils/grRuleToDNF.m | 76 +++++++++++ utils/grRuleToString.m | 57 ++++++++ utils/isDnfGrRule.m | 66 ++++++++++ utils/parseGrRule.m | 176 +++++++++++++++++++++++++ 9 files changed, 702 insertions(+), 156 deletions(-) create mode 100644 utils/grRuleToDNF.m create mode 100644 utils/grRuleToString.m create mode 100644 utils/isDnfGrRule.m create mode 100644 utils/parseGrRule.m diff --git a/manipulation/expandModel.m b/manipulation/expandModel.m index 288d2844..1c03dd66 100755 --- a/manipulation/expandModel.m +++ b/manipulation/expandModel.m @@ -15,8 +15,9 @@ % 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 % -------- @@ -24,19 +25,54 @@ % % 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 + 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 + 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 @@ -44,7 +80,10 @@ %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)]; @@ -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 diff --git a/manipulation/findPotentialErrors.m b/manipulation/findPotentialErrors.m index 1b15ff04..48960d9a 100644 --- a/manipulation/findPotentialErrors.m +++ b/manipulation/findPotentialErrors.m @@ -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); @@ -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 + continue end - if nonDnf + if ~isDnf issues(end+1,1)=struct( ... 'index', i, ... 'rxn', model.rxns{i}, ... @@ -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 - 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 diff --git a/manipulation/removeGenes.m b/manipulation/removeGenes.m index be7ded01..e2a50f32 100755 --- a/manipulation/removeGenes.m +++ b/manipulation/removeGenes.m @@ -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; @@ -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 + 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 diff --git a/testing/function_tests/tManipulation.m b/testing/function_tests/tManipulation.m index cebf6a1e..f0aa4483 100644 --- a/testing/function_tests/tManipulation.m +++ b/testing/function_tests/tManipulation.m @@ -251,5 +251,106 @@ function standardizeGrRulesReturnsRules(testCase) testCase.verifyNumElements(grRules, numel(testCase.model.rxns)); end + function findPotentialErrorsFlagsOnlyNonDnf(testCase) + m = struct(); + m.rxns = {'R1';'R2';'R3';'R4'}; + m.grRules = {'((G1 and G2) or G3)' % DNF, just bracketed + '(G1 or G2) or (G3 and G4)' % DNF + 'G1 or G2' % DNF + '(G1 or G2) and (G3 or G4)'};% genuinely non-DNF + issues = findPotentialErrors(m); + testCase.verifyEqual(vertcat(issues.index), 4); + end + + function findPotentialErrorsReportsUnparseable(testCase) + m = struct(); + m.rxns = {'R1'}; + m.grRules = {'(G1 and G2'}; + issues = findPotentialErrors(m); + testCase.verifyNumElements(issues, 1); + testCase.verifySubstring(issues(1).reason, 'Cannot be parsed'); + end + + function standardizeGrRulesRepairsBracketedDnf(testCase) + % A rule that is DNF but redundantly bracketed must be repaired, + % not skipped: standardizeGrRules leaves flagged rules alone, so a + % false positive from findPotentialErrors silently prevents repair. + m = struct(); + m.rxns = {'R1'}; + m.grRules = {'((G1 and G2) or G3)'}; + m.genes = {'G1';'G2';'G3'}; + m.rxnGeneMat = sparse([1 1 1]); + [grRules,~,indexes2check] = standardizeGrRules(m, true); + testCase.verifyEmpty(indexes2check); + testCase.verifyEqual(grRules{1}, '(G1 and G2) or G3'); + end + + function removeGenesMatchesWholeGeneIds(testCase) + % Gene "10" is a prefix of "100". Removing it must not take "100" + % with it, which an unanchored substring search would. + m = testCase.gprTestModel('10 or 100', {'10';'100'}, [1 1]); + r = removeGenes(m, {'10'}); + testCase.verifyEqual(r.grRules{1}, '100'); + testCase.verifyEqual(r.ub(1), 1000); + end + + function removeGenesDropsWholeComplex(testCase) + % A complex missing a subunit cannot form; the other isozyme lives. + m = testCase.gprTestModel('(G1 and G2) or G3', {'G1';'G2';'G3'}, [1 1 1]); + r = removeGenes(m, {'G1'}); + testCase.verifyEqual(r.grRules{1}, 'G3'); + end + + function removeGenesBlocksWhenNoEnzymeLeft(testCase) + m = testCase.gprTestModel('G1 and G2', {'G1';'G2'}, [1 1]); + r = removeGenes(m, {'G1'}); + testCase.verifyEmpty(r.grRules{1}); + testCase.verifyEqual(r.lb(1), 0); + testCase.verifyEqual(r.ub(1), 0); + end + + function expandModelKeepsMandatorySubunit(testCase) + % "g1 and (g2 or g3)" is two isozymes, both needing g1. Stripping + % brackets and splitting on ' or ' loses g1 from the second. + m = testCase.gprTestModel('g1 and (g2 or g3)', {'g1';'g2';'g3'}, [1 1 1]); + e = expandModel(m); + testCase.verifyEqual(sort(e.grRules), {'g1 and g2';'g1 and g3'}); + end + + function expandModelDistributesBothSides(testCase) + % Two or:s, but four isozymes. + m = testCase.gprTestModel('(g1 or g2) and (g3 or g4)', ... + {'g1';'g2';'g3';'g4'}, [1 1 1 1]); + [e, rxnToCheck] = expandModel(m); + testCase.verifyEqual(sort(e.grRules), ... + {'g1 and g3';'g1 and g4';'g2 and g3';'g2 and g4'}); + testCase.verifyEqual(rxnToCheck, {'R1'}); + end + + function expandModelLeavesDnfAlone(testCase) + % An OR of complexes expands without needing distributivity, so it + % must not be reported as needing a check. + m = testCase.gprTestModel('(g1 and g2) or (g3 and g4)', ... + {'g1';'g2';'g3';'g4'}, [1 1 1 1]); + [e, rxnToCheck] = expandModel(m); + testCase.verifyEqual(sort(e.grRules), {'g1 and g2';'g3 and g4'}); + testCase.verifyEmpty(rxnToCheck); + end + + end + + methods (Access = private) + function m = gprTestModel(~, grRule, genes, rxnGeneRow) + % Smallest model that removeGenes/expandModel will operate on. + m = struct(); + m.rxns = {'R1'}; m.rxnNames = {'R1'}; + m.mets = {'A';'B'}; m.metNames = {'A';'B'}; m.metComps = [1;1]; + m.comps = {'c'}; m.compNames = {'c'}; + m.S = sparse([-1;1]); m.lb = 0; m.ub = 1000; m.rev = 0; m.c = 0; + m.b = [0;0]; + m.genes = genes; + m.grRules = {grRule}; + m.rxnGeneMat = sparse(rxnGeneRow); + end end end diff --git a/testing/function_tests/tUtils.m b/testing/function_tests/tUtils.m index 928bd9d1..877c842f 100644 --- a/testing/function_tests/tUtils.m +++ b/testing/function_tests/tUtils.m @@ -87,5 +87,75 @@ function runRAVENtestsExists(testCase) testCase.verifyEqual(exist('runRAVENtests','file'), 2); end + function parseGrRuleBuildsTree(testCase) + t = parseGrRule('(G1 and G2) or G3'); + testCase.verifyEqual(t.type, 'or'); + testCase.verifyEqual(numel(t.children), 2); + testCase.verifyEqual(t.children{1}.type, 'and'); + testCase.verifyEqual(t.children{2}.id, 'G3'); + end + + function parseGrRuleEmptyIsEmpty(testCase) + testCase.verifyEmpty(parseGrRule('')); + testCase.verifyEmpty(parseGrRule(' ')); + end + + function parseGrRuleTakesWholeWordOperators(testCase) + % Gene IDs containing "and"/"or" are genes, not operators. + t = parseGrRule('RAND1 or ORF2'); + testCase.verifyEqual(t.type, 'or'); + testCase.verifyEqual(t.children{1}.id, 'RAND1'); + testCase.verifyEqual(t.children{2}.id, 'ORF2'); + end + + function parseGrRuleAcceptsSymbolOperators(testCase) + testCase.verifyEqual(parseGrRule('G1 & G2').type, 'and'); + testCase.verifyEqual(parseGrRule('G1 | G2').type, 'or'); + testCase.verifyEqual(parseGrRule('G1 AND G2').type, 'and'); + end + + function parseGrRuleRejectsMalformed(testCase) + testCase.verifyError(@() parseGrRule('(G1 and G2'), 'RAVEN:badGrRule'); + testCase.verifyError(@() parseGrRule('G1 and'), 'RAVEN:badGrRule'); + testCase.verifyError(@() parseGrRule('G1 G2'), 'RAVEN:badGrRule'); + end + + function isDnfGrRuleClassifies(testCase) + testCase.verifyTrue(isDnfGrRule('(G1 and G2) or G3')); + testCase.verifyTrue(isDnfGrRule('G1')); + testCase.verifyTrue(isDnfGrRule('')); + testCase.verifyFalse(isDnfGrRule('(G1 or G2) and G3')); + testCase.verifyFalse(isDnfGrRule('G1 and (G2 or G3)')); + end + + function isDnfGrRuleIgnoresRedundantBrackets(testCase) + % Bracketing alone never makes a rule non-DNF. Both of these were + % false positives for the one-level split this replaced. + testCase.verifyTrue(isDnfGrRule('((G1 and G2) or G3)')); + testCase.verifyTrue(isDnfGrRule('(G1 or G2) or (G3 and G4)')); + end + + function grRuleToDNFDistributes(testCase) + testCase.verifyEqual(grRuleToDNF('g1 and (g2 or g3)'), ... + {{'g1','g2'}, {'g1','g3'}}); + % Two or:s, but four isozymes: the count of ' or ' is not the + % count of clauses. + testCase.verifyEqual(numel(grRuleToDNF('(g1 or g2) and (g3 or g4)')), 4); + end + + function grRuleToDNFSimpleCases(testCase) + testCase.verifyEqual(grRuleToDNF(''), {}); + testCase.verifyEqual(grRuleToDNF('G1'), {{'G1'}}); + testCase.verifyEqual(grRuleToDNF('(G1 and G2) or G3'), {{'G1','G2'},{'G3'}}); + end + + function grRuleToStringRoundTrips(testCase) + rules = {'(G1 and G2) or G3', 'G1', 'G1 and G2', 'G1 or G2'}; + for k = 1:numel(rules) + testCase.verifyEqual(grRuleToString(parseGrRule(rules{k})), rules{k}); + end + testCase.verifyEqual(grRuleToString(parseGrRule('')), ''); + end + end end diff --git a/utils/grRuleToDNF.m b/utils/grRuleToDNF.m new file mode 100644 index 00000000..d5b5360c --- /dev/null +++ b/utils/grRuleToDNF.m @@ -0,0 +1,76 @@ +function clauses=grRuleToDNF(rule) +% grRuleToDNF Convert a grRule to disjunctive normal form. +% +% Returns the rule as a list of AND-clauses, applying distributivity where +% needed: "G1 and (G2 or G3)" becomes {{G1,G2},{G1,G3}}. Each clause is one +% isozyme; the genes within it are the subunits of that complex. +% +% Parameters +% ---------- +% rule : char or string or struct +% a grRule, or a tree from parseGrRule. +% +% Returns +% ------- +% clauses : cell +% one cell array of gene IDs per AND-clause. An empty rule yields {}. A +% rule with no OR anywhere yields a single clause. +% +% Notes +% ----- +% Gene order within a clause follows the order the genes appear in the rule. +% Duplicates are not removed: "G1 and G1" stays a two-element clause, since +% collapsing it would silently rewrite the user's rule. +% +% The number of clauses is the product of the OR-arities across a complex, +% so a pathological rule can expand combinatorially. Callers that build one +% object per clause should check numel(clauses) first. +% +% Examples +% -------- +% grRuleToDNF('(G1 and G2) or G3') % {{'G1','G2'}, {'G3'}} +% grRuleToDNF('G1 and (G2 or G3)') % {{'G1','G2'}, {'G1','G3'}} +% +% See also: parseGrRule, isDnfGrRule, expandModel + +if ~isstruct(rule) + rule=parseGrRule(rule); +end +if isempty(rule) + clauses={}; + return +end +clauses=nodeToDNF(rule); +end + + +function clauses=nodeToDNF(node) +% Mirrors _node_to_dnf in raven-toolbox manipulation/expand.py. +switch node.type + case 'gene' + clauses={{node.id}}; + case 'or' + % OR: concatenate the disjuncts' clauses. + clauses={}; + for i=1:numel(node.children) + clauses=[clauses nodeToDNF(node.children{i})]; %#ok + end + case 'and' + % AND: cross-product of the children's clauses. + clauses={{}}; + for i=1:numel(node.children) + childDNF=nodeToDNF(node.children{i}); + newClauses=cell(1,numel(clauses)*numel(childDNF)); + k=0; + for a=1:numel(clauses) + for b=1:numel(childDNF) + k=k+1; + newClauses{k}=[clauses{a} childDNF{b}]; + end + end + clauses=newClauses; + end + otherwise + error('RAVEN:badGrRule',['Unexpected grRule node type: ' node.type]); +end +end diff --git a/utils/grRuleToString.m b/utils/grRuleToString.m new file mode 100644 index 00000000..19706872 --- /dev/null +++ b/utils/grRuleToString.m @@ -0,0 +1,57 @@ +function rule=grRuleToString(tree) +% grRuleToString Render a grRule syntax tree back to a string. +% +% Writes the standard RAVEN format: " and " / " or " as operators, with +% parentheses added only where precedence requires them, so that +% parseGrRule(grRuleToString(tree)) round-trips to the same structure. +% +% Parameters +% ---------- +% tree : struct +% a tree from parseGrRule. An empty tree yields ''. +% +% Returns +% ------- +% rule : char +% the rendered grRule. +% +% Examples +% -------- +% grRuleToString(parseGrRule('(G1 and G2) or G3')) % '(G1 and G2) or G3' +% +% See also: parseGrRule, standardizeGrRules + +if isempty(tree) + rule=''; + return +end +rule=nodeToString(tree); +end + + +function str=nodeToString(node) +switch node.type + case 'gene' + str=node.id; + case 'and' + % An OR child inside an AND must be bracketed to keep precedence. + parts=cellfun(@(c) bracketIf(c,'or'),node.children,'UniformOutput',false); + str=strjoin(parts,' and '); + case 'or' + % An AND child inside an OR needs no brackets for correctness, but + % RAVEN writes complexes bracketed, matching standardizeGrRules and + % the joinOR helper in removeLowScoreGenes. + parts=cellfun(@(c) bracketIf(c,'and'),node.children,'UniformOutput',false); + str=strjoin(parts,' or '); + otherwise + error('RAVEN:badGrRule',['Unexpected grRule node type: ' node.type]); +end +end + + +function str=bracketIf(node,type) +str=nodeToString(node); +if strcmp(node.type,type) + str=['(' str ')']; +end +end diff --git a/utils/isDnfGrRule.m b/utils/isDnfGrRule.m new file mode 100644 index 00000000..f6e7258b --- /dev/null +++ b/utils/isDnfGrRule.m @@ -0,0 +1,66 @@ +function tf=isDnfGrRule(rule) +% isDnfGrRule Test whether a grRule is in disjunctive normal form. +% +% DNF ("OR of AND-complexes") means no AND operator has an OR anywhere +% beneath it, i.e. the rule is a single gene, a pure AND-complex, or an OR +% of those. Standard RAVEN format requires DNF: subunits that must all be +% present are ANDed into a complex, and the reaction is catalysed by at +% least one such complex. +% +% Parameters +% ---------- +% rule : char or string or struct +% a grRule, or a tree from parseGrRule. An empty rule is trivially DNF. +% +% Returns +% ------- +% tf : logical +% true if the rule is in disjunctive normal form. +% +% Examples +% -------- +% isDnfGrRule('(G1 and G2) or G3') % true +% isDnfGrRule('(G1 or G2) and G3') % false +% +% See also: parseGrRule, grRuleToDNF, findPotentialErrors + +if ~isstruct(rule) + rule=parseGrRule(rule); +end +if isempty(rule) + tf=true; + return +end +tf=isDnfNode(rule); +end + + +function tf=isDnfNode(node) +% Mirrors _is_dnf_node in raven-toolbox utils/gpr.py. +switch node.type + case 'gene' + tf=true; + case 'and' + % An AND-complex may not contain an OR anywhere beneath it. + tf=~any(cellfun(@containsOr,node.children)); + case 'or' + % Every disjunct must itself be DNF. + tf=all(cellfun(@isDnfNode,node.children)); + otherwise + % Unknown node type: don't flag it as a problem. + tf=true; +end +end + + +function tf=containsOr(node) +% Mirrors _contains_or in raven-toolbox utils/gpr.py. +switch node.type + case 'or' + tf=true; + case 'and' + tf=any(cellfun(@containsOr,node.children)); + otherwise + tf=false; +end +end diff --git a/utils/parseGrRule.m b/utils/parseGrRule.m new file mode 100644 index 00000000..596a5694 --- /dev/null +++ b/utils/parseGrRule.m @@ -0,0 +1,176 @@ +function tree=parseGrRule(rule) +% parseGrRule Parse a grRule into a syntax tree. +% +% Parses a gene-protein-reaction rule into a tree of AND/OR/gene nodes, so +% that callers can reason about its structure instead of searching the raw +% string. Gene IDs are matched as whole tokens, so IDs that contain "and" or +% "or" as a substring (RAND1), or that are prefixes of one another (10 and +% 100), are never confused with operators or with each other. +% +% Parameters +% ---------- +% rule : char or string +% a grRule, e.g. "(G1 and G2) or G3". Both the word forms ("and", "or", +% any case) and the symbol forms ("&", "&&", "|", "||") are accepted. +% An empty rule yields an empty tree. +% +% Returns +% ------- +% tree : struct +% a scalar struct with fields: +% +% - type : char — 'gene', 'and' or 'or'. +% - id : char — the gene ID; '' for 'and'/'or' nodes. +% - children : cell — child nodes; {} for 'gene' nodes. +% +% An empty rule returns a 0x0 struct array, for which isempty is true. +% +% Notes +% ----- +% The grammar is the usual precedence, with AND binding tighter than OR: +% +% expr := term ( OR term )* +% term := factor ( AND factor )* +% factor := gene | '(' expr ')' +% +% Malformed rules raise RAVEN:badGrRule rather than being silently +% misparsed. Callers that lint user input (findPotentialErrors) catch this +% and report it instead of failing. +% +% Examples +% -------- +% tree = parseGrRule('(G1 and G2) or G3'); +% tree.type % 'or' +% numel(tree.children) % 2 +% +% See also: isDnfGrRule, grRuleToDNF, grRuleToString + +if nargin<1 || isempty(rule) + tree=emptyTree(); + return +end +rule=char(rule); +if all(isspace(rule)) + tree=emptyTree(); + return +end + +toks=tokenizeGrRule(rule); +if isempty(toks) + tree=emptyTree(); + return +end + +[tree,pos]=parseExpr(toks,1,rule); +if pos<=numel(toks) + error('RAVEN:badGrRule',['Unexpected "' tokenText(toks{pos}) '" in grRule: ' rule]); +end +end + + +function tree=emptyTree() +tree=struct('type',{},'id',{},'children',{}); +end + + +function toks=tokenizeGrRule(rule) +% Split into '(' , ')' , 'AND' , 'OR' and gene tokens. Operators are matched +% as whole words, which is what keeps gene IDs containing "and"/"or" intact. +toks={}; +i=1; +n=numel(rule); +while i<=n + c=rule(i); + if isspace(c) + i=i+1; + continue + end + if c=='(' || c==')' + toks{end+1}=c; %#ok + i=i+1; + continue + end + % Read a word: everything up to whitespace or a bracket. + j=i; + while j<=n && ~isspace(rule(j)) && rule(j)~='(' && rule(j)~=')' + j=j+1; + end + word=rule(i:j-1); + switch lower(word) + case {'and','&','&&'} + toks{end+1}='AND'; %#ok + case {'or','|','||'} + toks{end+1}='OR'; %#ok + otherwise + toks{end+1}=struct('gene',word); %#ok + end + i=j; +end +end + + +function [node,pos]=parseExpr(toks,pos,rule) +% expr := term ( OR term )* +[node,pos]=parseTerm(toks,pos,rule); +if pos>numel(toks) || ~isOp(toks{pos},'OR') + return +end +children={node}; +while pos<=numel(toks) && isOp(toks{pos},'OR') + [child,pos]=parseTerm(toks,pos+1,rule); + children{end+1}=child; %#ok +end +node=struct('type','or','id','','children',{children}); +end + + +function [node,pos]=parseTerm(toks,pos,rule) +% term := factor ( AND factor )* +[node,pos]=parseFactor(toks,pos,rule); +if pos>numel(toks) || ~isOp(toks{pos},'AND') + return +end +children={node}; +while pos<=numel(toks) && isOp(toks{pos},'AND') + [child,pos]=parseFactor(toks,pos+1,rule); + children{end+1}=child; %#ok +end +node=struct('type','and','id','','children',{children}); +end + + +function [node,pos]=parseFactor(toks,pos,rule) +% factor := gene | '(' expr ')' +if pos>numel(toks) + error('RAVEN:badGrRule',['grRule ends after an operator: ' rule]); +end +tok=toks{pos}; +if isOp(tok,'(') + [node,pos]=parseExpr(toks,pos+1,rule); + if pos>numel(toks) || ~isOp(toks{pos},')') + error('RAVEN:badGrRule',['Unbalanced parentheses in grRule: ' rule]); + end + pos=pos+1; + return +end +if isstruct(tok) + node=struct('type','gene','id',tok.gene,'children',{{}}); + pos=pos+1; + return +end +error('RAVEN:badGrRule',['Expected a gene but found "' tokenText(tok) '" in grRule: ' rule]); +end + + +function tf=isOp(tok,op) +tf=ischar(tok) && strcmp(tok,op); +end + + +function txt=tokenText(tok) +if isstruct(tok) + txt=tok.gene; +else + txt=tok; +end +end