fix: various functions

core/checkModelStruct.m

@@ -120,6 +120,21 @@ function checkModelStruct(model,throwErrors,trimWarnings)
         EM='The "grRules" field must be a cell array of strings';
         dispEM(EM,throwErrors);
     end
+    if ~isfield(model,'genes')
+        EM='If "grRules" field exists, the model should also contain a "genes" field';
+        dispEM(EM,throwErrors);
+    else
+        geneList = strjoin(model.grRules);
+        geneList = regexp(geneList,' |)|(|and|or','split'); % Remove all grRule punctuation
+        geneList = geneList(~cellfun(@isempty,geneList));  % Remove spaces and empty genes
+        geneList = setdiff(unique(geneList),model.genes);
+        if ~isempty(geneList)
+            problemGrRules = model.rxns(contains(model.grRules,geneList));
+            problemGrRules = strjoin(problemGrRules(:),'; ');
+            EM=['The reaction(s) "' problemGrRules '" contain the following genes in its "grRules" field, but these are not in the "genes" field:'];
+            dispEM(EM,throwErrors,geneList);
+        end
+    end
 end
 if isfield(model,'rxnComps')
     if ~isnumeric(model.rxnComps)
@@ -262,7 +277,10 @@ function checkModelStruct(model,throwErrors,trimWarnings)
 EM='The following reactions have contradicting bounds:';
 dispEM(EM,throwErrors,model.rxns(model.lb>model.ub),trimWarnings);
 EM='The following reactions have bounds contradicting their reversibility:';
-dispEM(EM,throwErrors,model.rxns(model.lb<0 & model.rev==0),trimWarnings);
+contradictBound = (model.lb < 0 & model.ub > 0 & model.rev==0) | ... % Reversible bounds, irreversible label
+                  (model.lb < 0 & model.ub <= 0 & model.rev==1) | ... % Negative bounds, reversible label
+                  (model.lb >= 0 & model.ub > 0 & model.rev==1); % Positive bounds, reversible label
+dispEM(EM,throwErrors,model.rxns(contradictBound),trimWarnings);
 
 %Multiple or no objective functions not allowed in SBML L3V1 FBCv2
 if numel(find(model.c))>1

core/constructS.m

@@ -152,7 +152,8 @@
         strjoin(unique(metsToS(~metsPresent)),', ')],'')
     else
         missingMet = find(~metsPresent);
-        missingMet = char(strcat(metsToS(missingMet),' (reaction:',rxns(rxnsToS(missingMet)),')\n'));
+        missingMet = strcat(metsToS(missingMet),' (reaction:',rxns(rxnsToS(missingMet)),')\n');
+        missingMet = strjoin(missingMet,'');
         error(['Could not find the following metabolites (reaction indicated) in the metabolite list: \n' ...
             missingMet '%s'],'');
     end

core/getExchangeRxns.m

@@ -1,47 +1,82 @@
 function [exchangeRxns, exchangeRxnsIndexes]=getExchangeRxns(model,reactionType)
 % getExchangeRxns
-%   Retrieves the exchange reactions from a model
+%   Retrieves the exchange reactions from a model. Exchange reactions are
+%   identified by having either no substrates or products.
 %
+% Input:
 %   model               a model structure
-%   reactionType        retrieve all reactions ('both'), only production
-%                       ('out'), or only consumption ('in') (optional, default
-%                       'both')
+%   reactionType        which exchange reactions should be returned
+%                       'all'     all reactions, irrespective of reaction
+%                                 bounds
+%                       'uptake'  reactions with bounds that imply that
+%                                 only uptake are allowed. Reaction
+%                                 direction, upper and lower bounds are
+%                                 all considered
+%                       'excrete' reactions with bounds that imply that
+%                                 only excretion are allowed. Reaction
+%                                 direction, upper and lower bounds are
+%                                 all considered
+%                       'reverse' reactions with non-zero upper and lower
+%                                 bounds that imply that both uptake and
+%                                 excretion are allowed
+%                       'blocked' reactions that have zero upper and lower
+%                                 bounds, not allowing any flux
+%                       'in'      reactions where the boundary metabolite
+%                                 is the substrate of the reaction, a
+%                                 positive flux value would imply uptake,
+%                                 but reaction bounds are not considered
+%                       'out'     reactions where the boundary metabolite
+%                                 is the substrate of the reaction, a
+%                                 positive flux value would imply uptake,
+%                                 but reaction bounds are not considered.
 %
+% Output:
 %   exchangeRxns        cell array with the IDs of the exchange reactions
 %   exchangeRxnsIndexes vector with the indexes of the exchange reactions
 %
-%   Exchange reactions are defined as reactions which involve only products
-%   or only reactants. If the unconstrained field is present, then that is
-%   used instead.
+% Note:
+%   The union of 'in' and 'out' equals 'all'. Also, the union of 'uptake',
+%   'excrete', 'reverse' and 'blocked' equals all.
 %
 % Usage: [exchangeRxns,exchangeRxnsIndexes]=getExchangeRxns(model,reactionType)
 
 if nargin<2
-    reactionType='both';
+    reactionType='all';
 else
     reactionType=char(reactionType);
 end
 
-hasNoProducts=sparse(numel(model.rxns),1);
-hasNoReactants=sparse(numel(model.rxns),1);
-
-if isfield(model,'unconstrained')
-    if strcmpi(reactionType,'both') || strcmpi(reactionType,'out')
-        [~, I]=find(model.S(model.unconstrained~=0,:)>0);
-        hasNoProducts(I)=true;
-    end
-    if strcmpi(reactionType,'both') || strcmpi(reactionType,'in')
-        [~, I]=find(model.S(model.unconstrained~=0,:)<0);
-        hasNoReactants(I)=true;
-    end
+% Find exchange reactions
+if isfield(model, 'unconstrained')
+    [~, I]=find(model.S(model.unconstrained~=0,:)>0);
+    hasNoProd(I)=true;
+    [~, I]=find(model.S(model.unconstrained~=0,:)<0);
+    hasNoSubs(I)=true;
 else
-    if strcmpi(reactionType,'both') || strcmpi(reactionType,'out')
-        hasNoProducts=sum((model.S>0))==0;
-    end
-    if strcmpi(reactionType,'both') || strcmpi(reactionType,'in')
-        hasNoReactants=sum((model.S<0))==0;
-    end
+    hasNoProd = transpose(find(sum(model.S>0)==0));
+    hasNoSubs = transpose(find(sum(model.S<0)==0));
+end
+allExch   = [hasNoProd; hasNoSubs];
+
+switch reactionType
+    case {'both','all'} % For legacy reasons, 'both' is also allowed
+        exchangeRxnsIndexes = allExch;
+    case 'in'
+        exchangeRxnsIndexes = hasNoSubs;
+    case 'out'
+        exchangeRxnsIndexes = hasNoProd;
+    case 'blocked'
+        exchangeRxnsIndexes = allExch(model.lb(allExch) == 0 & model.ub(allExch) == 0);
+    case 'reverse'
+        exchangeRxnsIndexes = allExch(model.lb(allExch) < 0 & model.ub(allExch) > 0);
+    case 'uptake'
+        exchangeRxnsIndexes = allExch([(model.lb(hasNoSubs) >= 0 & model.ub(hasNoSubs) > 0); ...
+                              (model.lb(hasNoProd) < 0 & model.ub(hasNoProd) <= 0)]);
+    case 'excrete'
+        exchangeRxnsIndexes = allExch([(model.lb(hasNoSubs) < 0 & model.ub(hasNoSubs) <= 0); ...
+                              (model.lb(hasNoProd) >= 0 & model.ub(hasNoProd) > 0)]);
+    otherwise
+        error('Invalid reactionType specified')
 end
-exchangeRxnsIndexes=find(hasNoProducts(:) | hasNoReactants(:));
-exchangeRxns=model.rxns(exchangeRxnsIndexes);
+exchangeRxns = model.rxns(exchangeRxnsIndexes);
 end