fix: give the same result as geckopy where the parity scenarios found differences

src/change_model/addNewRxnsToEC.m

@@ -35,7 +35,7 @@
 % - rxns : cell array with unique strings that identifies each reaction.
 % - rxnNames : cell array with the names of each reaction.
 % - equations : cell array with equation strings. Decimal coefficients are expressed as "1.2". Reversibility is indicated by "<=>" or "=>".
-% - grRules : cell array with the gene-reaction relationship for each reaction. For example "(A and B) or (C)" means that the reaction could be catalyzed by a complex between A & B or by C on its own. All the genes have to be present in model.genes. Add genes with addGenesRaven before calling this function if needed (opt, default '').
+% - grRules : cell array with the gene-reaction relationship for each reaction. For example "(A and B) or (C)" means that the reaction could be catalyzed by a complex between A & B or by C on its own. A rule with brackets is split into its isozymes in disjunctive normal form, so "A and (B or C)" is the two isozymes "A and B" and "A and C". All the genes have to be present in model.genes. Add genes with addGenesRaven before calling this function if needed (opt, default '').
 %
 % The newEnzymes structure has the following fields:
 %
@@ -136,11 +136,17 @@
 
     % 3. Add new reactions
 
+    % The operators of a gene rule are case insensitive: write them in lower
+    % case, as the rest of the function and RAVEN's addRxns expect. Gene ids
+    % that contain the letters (Sand1, andA, ORF1) are left as they are.
+    newRxns.grRules = regexprep(convertCharArray(newRxns.grRules), ...
+        '(?<![^\s()])(and|or)(?![^\s()])', '${lower($0)}', 'ignorecase');
+
     % Before add the reactions validate that all the genes in grRules are
     % present in the model.
-    genesInRules = strjoin(convertCharArray(newRxns.grRules));
-    genesInRules = regexp(genesInRules,' |)|(|and|or','split'); % Remove all grRule punctuation
-    genesInRules = genesInRules(~cellfun(@isempty,genesInRules));  % Remove spaces and empty genes
+    genesInRules = strjoin(newRxns.grRules);
+    genesInRules = regexp(genesInRules,'[^\s()]+','match'); % Gene ids and the words 'and' and 'or', without spaces and brackets
+    genesInRules = genesInRules(~ismember(lower(genesInRules),{'and','or'}));  % Only the gene ids
     genesInRules = setdiff(unique(genesInRules),model.genes); % Only keep new genes
 
     % If there is any gene which is not present, as to be pass as
@@ -168,13 +174,14 @@
     end
 
     % Get reactions to separate isozymes
-    toSplit = find(cellfun(@(x) contains(x,'or'),newRxns.grRules));
+    isozymes = cellfun(@gprToDnf, newRxns.grRules, 'UniformOutput', false);
+    toSplit = find(cellfun(@numel,isozymes) > 1);
     if ~isempty(toSplit)
         splitRxns = {}; splitRxnNames = {}; splitEquations = {}; splitGrRules = {};
         for i=1:numel(toSplit)
             idx = toSplit(i);
-            newRules = split(newRxns.grRules(idx),'or');
-            newRules = strtrim(erase(newRules, {'(' ')'}));
+            % One reaction per AND-clause of the gene rule in disjunctive normal form
+            newRules = cellfun(@(x) strjoin(x,' and '), isozymes{idx}, 'UniformOutput', false);
 
             for j=1:numel(newRules)
                 splitRxns{end+1} = [newRxns.rxns{idx} '_EXP_' num2str(j)];
@@ -249,8 +256,8 @@
         model.ec.eccodes(end+1) = {''};
 
         % Get rules
-        genes = newRxns.grRules(rxnIdx);
-        genes = strtrim(erase(split(genes,'and'), {'(' ')'}));
+        genes = regexp(newRxns.grRules{rxnIdx},'[^\s()]+','match'); % Gene ids and the word 'and', without spaces and brackets
+        genes = genes(~ismember(lower(genes),{'and','or'}));
         [~,enzIdx] = ismember(genes, model.ec.genes);
 
         % Update the enzyme rxns matrix
@@ -260,3 +267,59 @@
     rxnsAdded = newRxns.rxns(:);
 end
 end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Function that expands a gene rule into disjunctive normal form: a cell array
+%with one cell per isozyme, each holding the gene ids that are ANDed in it.
+%Gene ids are whole words, so an id that contains 'or' or 'and' is kept
+%intact. An empty rule gives an empty cell array.
+function dnf = gprToDnf(rule)
+tokens = regexp(rule,'\(|\)|[^\s()]+','match');
+isOperator = strcmpi(tokens,'and') | strcmpi(tokens,'or');
+tokens(isOperator) = lower(tokens(isOperator)); % Operators are case insensitive
+if isempty(tokens)
+    dnf = {};
+    return
+end
+[dnf,pos] = parseOr(tokens,1,rule);
+if pos <= numel(tokens)
+    error('Cannot parse the gene rule "%s"',rule)
+end
+end
+
+function [dnf,pos] = parseOr(tokens,pos,rule)
+[dnf,pos] = parseAnd(tokens,pos,rule);
+while pos <= numel(tokens) && strcmp(tokens{pos},'or')
+    [rhs,pos] = parseAnd(tokens,pos+1,rule);
+    dnf = [dnf,rhs];
+end
+end
+
+function [dnf,pos] = parseAnd(tokens,pos,rule)
+[dnf,pos] = parseAtom(tokens,pos,rule);
+while pos <= numel(tokens) && strcmp(tokens{pos},'and')
+    [rhs,pos] = parseAtom(tokens,pos+1,rule);
+    product = cell(1,0);
+    for a = 1:numel(dnf)
+        for b = 1:numel(rhs)
+            product{end+1} = [dnf{a},rhs{b}]; %#ok<AGROW>
+        end
+    end
+    dnf = product;
+end
+end
+
+function [dnf,pos] = parseAtom(tokens,pos,rule)
+if pos > numel(tokens) || any(strcmp(tokens{pos},{')','and','or'}))
+    error('Cannot parse the gene rule "%s"',rule)
+elseif strcmp(tokens{pos},'(')
+    [dnf,pos] = parseOr(tokens,pos+1,rule);
+    if pos > numel(tokens) || ~strcmp(tokens{pos},')')
+        error('Cannot parse the gene rule "%s"',rule)
+    end
+    pos = pos + 1;
+else
+    dnf = {tokens(pos)};
+    pos = pos + 1;
+end
+end

src/change_model/applyComplexData.m

@@ -37,6 +37,9 @@
 %
 % Notes
 % -----
+% The subunit stoichiometry is assigned per protein ID: the proteins of a
+% complex need not be listed in the order of model.ec.enzymes.
+%
 % foundComplex and proposedComplex share the following columns: rxn
 % (reaction identifier), complexID and name (Complex Portal identifiers),
 % genes (gene names of the complex subunits), protID_model (matched protein
@@ -154,8 +157,11 @@
         foundComplex{foundComplexCount,7} = complexData(complexMatch).stochiometry;
 
         % Apply in rxnEnzMat
+        % The stoichiometry is taken per protein ID, as the proteins of the
+        % complex need not be listed in the order of model.ec.enzymes
         [~, modelProtsIdx] = ismember(modelProts, model.ec.enzymes);
-        model.ec.rxnEnzMat(i,modelProtsIdx) = complexData(complexMatch).stochiometry;
+        [~, complexPos] = ismember(modelProts, complexData(complexMatch).protID);
+        model.ec.rxnEnzMat(i,modelProtsIdx) = complexData(complexMatch).stochiometry(complexPos);
 
     elseif modComplexUnits > 1 % Only suggest potential complex if grRule > 1 protein
         moreUnitsInData = find(percMatch == 1 & totalMatch > 1);

src/change_model/applyCustomKcats.m

@@ -43,9 +43,13 @@
 % - gene_name : short gene name (optional, not used in matching).
 % - kcat : new kcat value (one per entry).
 % - rxns : reaction identifiers, multiple for the same kcat are separated by ',' (see further explanation below).
-% - notes : will be appended to model.ec.notes (optional).
+% - notes : will be appended to model.ec.notes, separated by ', ' from an existing note; an empty note changes nothing (optional).
 % - stoicho : complex stoichiometry, separated by ' + ' (examples: '1' or '3 + 1'), matching the order in proteins field.
 %
+% An entry with an empty kcat is skipped. Protein and reaction identifiers
+% are matched exactly, including their letter case. Light ecModels are
+% supported: the rows of model.ec are matched by their ec reaction ids.
+%
 % Matching order:
 %
 % 1. reactions are identified by .proteins and .rxns.
@@ -121,125 +125,141 @@
 enzInModel = cell(length(model.ec.rxns),1);
 ecRxnNoSuffix = regexprep(model.ec.rxns,'_EXP_\d+$','');
 
-% Implementation for full GECKO formulation
-if ~model.ec.geckoLight
-    for i = 1:numel(customKcats.proteins)
-        if isempty(customKcats.proteins{i})
-            %If only reaction ID(s) is/are specified (and no proteins),
-            %then apply the kcat to all isozymic reactions
-            rxns    = strtrim(strsplit(customKcats.rxns{i}, ','));
-            rxnIdxs = ismember(ecRxnNoSuffix,rxns);
-            rxnToUpdate(rxnIdxs) = 1;
-            model.ec.kcat(rxnIdxs) = customKcats.kcat(i);
-
-            % Add note mentioning manual kcat change
-            model.ec.source(rxnIdxs) = {'custom'};
-            if isfield(customKcats,'notes')
-                matchedIdxs = find(rxnIdxs);
-                for j = 1:numel(matchedIdxs)
-                    if isempty(model.ec.notes{matchedIdxs(j), 1}) && ~isempty(customKcats.notes{i})
-                        model.ec.notes{matchedIdxs(j), 1} = customKcats.notes{i};
-                    else
-                        model.ec.notes{matchedIdxs(j), 1} = [model.ec.notes{matchedIdxs(j), 1} ', ' customKcats.notes{i}];
-                    end
-                end
+for i = 1:numel(customKcats.proteins)
+    if isnan(customKcats.kcat(i))
+        printOrange(['WARNING: Entry ' num2str(i) ' has no kcat; it is skipped.']);
+        continue
+    end
+    if isempty(customKcats.proteins{i})
+        %If only reaction ID(s) is/are specified (and no proteins),
+        %then apply the kcat to all isozymic reactions
+        rxns    = strtrim(strsplit(customKcats.rxns{i}, ','));
+        rxnIdxs = ismember(ecRxnNoSuffix,rxns);
+        rxnToUpdate(rxnIdxs) = 1;
+        model.ec.kcat(rxnIdxs) = customKcats.kcat(i);
+
+        % Add note mentioning manual kcat change
+        model.ec.source(rxnIdxs) = {'custom'};
+        if isfield(customKcats,'notes')
+            matchedIdxs = find(rxnIdxs);
+            for j = 1:numel(matchedIdxs)
+                model.ec.notes{matchedIdxs(j), 1} = appendNote(model.ec.notes{matchedIdxs(j), 1}, customKcats.notes{i});
             end
+        end
+    else
+        prots = strtrim(strsplit(customKcats.proteins{i}, '+'));
+
+        % Find the index for the enzymes which kcat will be changed. In
+        % case the protein is not in the model generate a warning.
+        try
+            enzIdx = cellfun(@(x) find(strcmp(model.ec.enzymes, x)), prots);
+        catch
+            enzIdx = [];
+            printOrange(['WARNING: Protein(s) ' customKcats.proteins{i} ' were not found in the model.']);
+        end
+
+        % if not specific reactions are defined, find all the reaction
+        % index where the enzyme is used
+        if isempty(customKcats.rxns{i})
+            temp_rxnIdxs = arrayfun(@(x) find(model.ec.rxnEnzMat(:, x)), enzIdx, 'UniformOutput', false);
+            % otherwhise, If a set of reactions if defined, only get the index for those
+            % but ignore any _EXP_ suffixes
         else
-            prots = strtrim(strsplit(customKcats.proteins{i}, '+'));
-
-            % Find the index for the enzymes which kcat will be changed. In
-            % case the protein is not in the model generate a warning.
-            try
-                enzIdx = cellfun(@(x) find(strcmpi(model.ec.enzymes, x)), prots);
-            catch
-                enzIdx = [];
-                printOrange(['WARNING: Protein(s) ' customKcats.proteins{i} ' were not found in the model.']);
-            end
+            rxns = strtrim(strsplit(customKcats.rxns{i}, ','));
+            temp_rxnIdxs = cellfun(@(x) find(strcmp(ecRxnNoSuffix, x)), rxns, 'UniformOutput', false);
+        end
 
-            % if not specific reactions are defined, find all the reaction
-            % index where the enzyme is used
-            if isempty(customKcats.rxns{i})
-                temp_rxnIdxs = arrayfun(@(x) find(model.ec.rxnEnzMat(:, x)), enzIdx, 'UniformOutput', false);
-                % otherwhise, If a set of reactions if defined, only get the index for those
-                % but ignore any _EXP_ suffixes
-            else
-                rxns = strtrim(strsplit(customKcats.rxns{i}, ','));
-                temp_rxnIdxs = arrayfun(@(x) find(strcmpi(ecRxnNoSuffix, x)), rxns, 'UniformOutput', false);
+        if ~isempty(temp_rxnIdxs)
+            rxnIdxs = [];
+            for j = 1:numel(temp_rxnIdxs)
+                rxnIdxs = [rxnIdxs; temp_rxnIdxs{j}];
             end
 
-            if ~isempty(temp_rxnIdxs)
-                rxnIdxs = [];
-                for j = 1:numel(temp_rxnIdxs)
-                    rxnIdxs = [rxnIdxs; temp_rxnIdxs{j}];
-                end
+            % Check when multiple proteins are involved, since it can return same rxn n times
+            rxnIdxs = unique(rxnIdxs); %unique(rxnIdxs{1, :});
 
-                % Check when multiple proteins are involved, since it can return same rxn n times
-                rxnIdxs = unique(rxnIdxs); %unique(rxnIdxs{1, :});
+            for j = 1:numel(rxnIdxs)
+                % Get all the enzymes involved in the reaction
+                allEnzInRxn = find(model.ec.rxnEnzMat(rxnIdxs(j),:));
 
-                for j = 1:numel(rxnIdxs)
-                    % Get all the enzymes involved in the reaction
-                    allEnzInRxn = find(model.ec.rxnEnzMat(rxnIdxs(j),:));
+                C = intersect(enzIdx, allEnzInRxn);
 
-                    C = intersect(enzIdx, allEnzInRxn);
-
-                    % Determine the match percentage bewteen the rules
-                    if numel(C) == numel(enzIdx) && numel(C) == numel(allEnzInRxn)
-                        match = 1;
+                % Determine the match percentage bewteen the rules
+                if numel(C) == numel(enzIdx) && numel(C) == numel(allEnzInRxn)
+                    match = 1;
+                else
+                    if numel(enzIdx) < numel(allEnzInRxn)
+                        match = numel(C) / numel(allEnzInRxn);
                     else
-                        if numel(enzIdx) < numel(allEnzInRxn)
-                            match = numel(C) / numel(allEnzInRxn);
-                        else
-                            match = numel(C) / numel(enzIdx);
-                        end
+                        match = numel(C) / numel(enzIdx);
                     end
+                end
 
-                    % Update the kcat value stored in the model, if match 100%,
-                    % otherwhise if >= 0.5 inform to the user to be validated
-                    if match == 1
-                        rxnToUpdate(rxnIdxs(j)) = 1;
-                        model.ec.kcat(rxnIdxs(j)) = customKcats.kcat(i);
-
-                        % Add note mentioning manual kcat change
-                        model.ec.source{rxnIdxs(j),1} = 'custom';
-                        if isfield(customKcats,'notes')
-                            if isempty(model.ec.notes{rxnIdxs(j), 1}) && ~isempty(customKcats.notes{i})
-                                model.ec.notes{rxnIdxs(j), 1} = customKcats.notes{i};
-                            else
-                                model.ec.notes{rxnIdxs(j), 1} = [model.ec.notes{rxnIdxs(j), 1} ', ' customKcats.notes{i}];
-                            end
-                        end
-                    elseif match >= 0.5 && match < 1
-                        rxnNotMatch(rxnIdxs(j)) = 1;
-                        evaluatedRule{rxnIdxs(j), 1} = customKcats.proteins{i};
-                        enzInModel{rxnIdxs(j), 1} = strjoin(model.ec.enzymes(allEnzInRxn), ' + ');
+                % Update the kcat value stored in the model, if match 100%,
+                % otherwhise if >= 0.5 inform to the user to be validated
+                if match == 1
+                    rxnToUpdate(rxnIdxs(j)) = 1;
+                    model.ec.kcat(rxnIdxs(j)) = customKcats.kcat(i);
+
+                    % Add note mentioning manual kcat change
+                    model.ec.source{rxnIdxs(j),1} = 'custom';
+                    if isfield(customKcats,'notes')
+                        model.ec.notes{rxnIdxs(j), 1} = appendNote(model.ec.notes{rxnIdxs(j), 1}, customKcats.notes{i});
                     end
+                elseif match >= 0.5 && match < 1
+                    rxnNotMatch(rxnIdxs(j)) = 1;
+                    evaluatedRule{rxnIdxs(j), 1} = customKcats.proteins{i};
+                    enzInModel{rxnIdxs(j), 1} = strjoin(model.ec.enzymes(allEnzInRxn), ' + ');
                 end
             end
         end
     end
+end
 
-    % Apply the new kcat values to the model
-    if ~isempty(find(rxnToUpdate, 1))
-        model = applyKcatConstraints(model, rxnToUpdate);
+% Apply the new kcat values to the model
+if ~isempty(find(rxnToUpdate, 1))
+    if model.ec.geckoLight
+        %Light ecModels update by reaction of the model, not by ec row
+        updateRxns = ismember(model.rxns, extractAfter(model.ec.rxns(rxnToUpdate), 4));
     else
-        printOrange('WARNING: No matches found. Consider checking the IDs or proteins in customKcats.');
+        updateRxns = rxnToUpdate;
     end
+    model = applyKcatConstraints(model, updateRxns);
+else
+    printOrange('WARNING: No matches found. Consider checking the IDs or proteins in customKcats.');
+end
+
+rxnUpdated = model.ec.rxns(find(rxnToUpdate));
 
-    rxnUpdated = model.ec.rxns(find(rxnToUpdate));
-
-    % Remove from notMatch those reactions already updated
-    remove = and(rxnToUpdate, rxnNotMatch);
-    rxnNotMatch(remove) = 0;
-    evaluatedRule(remove) = '';
-    enzInModel(remove) = '';
-
-    idRxns = model.ec.rxns(find(rxnNotMatch));
-    fullIdx = cellfun(@(x) find(strcmpi(model.rxns, x)), idRxns);
-    rxnsNames = model.rxnNames(fullIdx);
-    evaluatedRule = evaluatedRule(~cellfun('isempty', evaluatedRule));
-    enzInModel = enzInModel(~cellfun('isempty', enzInModel));
-    rules = model.grRules(fullIdx);
-    notMatch = table(idRxns, rxnsNames, evaluatedRule, enzInModel, rules, ...
-        'VariableNames',{'rxns', 'name', 'custom enzymes', 'enzymes in model', 'rules'});
+% Remove from notMatch those reactions already updated
+remove = and(rxnToUpdate, rxnNotMatch);
+rxnNotMatch(remove) = 0;
+evaluatedRule(remove) = '';
+enzInModel(remove) = '';
+
+idRxns = model.ec.rxns(find(rxnNotMatch));
+if model.ec.geckoLight
+    modelRxnIds = extractAfter(idRxns, 4);
+else
+    modelRxnIds = idRxns;
+end
+fullIdx = cellfun(@(x) find(strcmpi(model.rxns, x), 1), modelRxnIds);
+rxnsNames = model.rxnNames(fullIdx);
+evaluatedRule = evaluatedRule(~cellfun('isempty', evaluatedRule));
+enzInModel = enzInModel(~cellfun('isempty', enzInModel));
+rules = model.grRules(fullIdx);
+notMatch = table(idRxns, rxnsNames, evaluatedRule, enzInModel, rules, ...
+    'VariableNames',{'rxns', 'name', 'custom enzymes', 'enzymes in model', 'rules'});
+end
+
+function notes = appendNote(notes, note)
+% An empty note changes nothing; a note is separated by ', ' from an existing one.
+note = strtrim(note);
+if isempty(note)
+    return
+elseif isempty(notes)
+    notes = note;
+else
+    notes = [notes ', ' note];
 end
 end