fix: parsing grRules; addMets and addRxns; importModel annotations; r…

…eplaceMets and changeGrRules options (#571)

* fix: addMets if metNames and not mets are given

* feat: changeGrRules multiple rxns same grRule

* fix: addRxns subSystems correct format

* fix: subSystems if cell array of cell arrays

* fix: importModel if identifiers.org/name:value

* refactor: checkModelStruct use getGenesFromGrRules

* feat: replaceMets can use identifiers instead

* refactor: use getGenesFromGrRules

* doc: checkRxn clarify meaning of output

core/addExchangeRxns.m

@@ -62,7 +62,11 @@
     model.eccodes=[model.eccodes;filler];
 end
 if isfield(model,'subSystems')
-    model.subSystems=[model.subSystems;filler];
+    fillerSub = filler;
+    if iscell(model.subSystems(1,1))
+        fillerSub = repmat({fillerSub},numel(J),1);
+    end
+    model.subSystems=[model.subSystems;fillerSub];
 end
 if isfield(model,'grRules')
     model.grRules=[model.grRules;filler];

core/addMets.m

@@ -70,6 +70,7 @@
 
 %Check some stuff regarding the required fields
 if ~isfield(metsToAdd,'mets')
+    metsToAdd.metNames=convertCharArray(metsToAdd.metNames);
     metsToAdd.mets=generateNewIds(newModel,'mets',prefix,numel(metsToAdd.metNames));
 else
     metsToAdd.mets=convertCharArray(metsToAdd.mets);

core/addRxns.m

@@ -363,7 +363,7 @@
     if ~isfield(newModel,'subSystems')
         newModel.subSystems=celllargefiller;
     end
-    newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems];
+    newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems(:)];
 else
     %Fill with standard if it doesn't exist
     if isfield(newModel,'subSystems')

core/changeGrRules.m

@@ -23,22 +23,22 @@
 rxns=convertCharArray(rxns);
 grRules=convertCharArray(grRules);
 
+if isscalar(grRules) && ~isscalar(rxns)
+    grRules = repmat(grRules,1,numel(rxns));
+end
 if ~(numel(grRules)==numel(rxns))
     error('Number of rxns and grRules should be identical')
 end
 
-for i=1:length(rxns)
-    % Add genes to model
-    geneList=transpose(cell(unique(regexp(grRules{i},'[)(]*|( and )*|( or )*','split')))); % Extract individual, unique genes from the geneAssoc provided
-    geneList=geneList(~cellfun(@isempty, geneList));
-    genesToAdd.genes=setdiff(geneList,model.genes); % Only keep the genes that are not yet part of the model.genes.
-    if ~isempty(genesToAdd.genes)
-        model=addGenesRaven(model,genesToAdd); % Add genes
-    end
-
-    % Find reaction and gene indices
-    idx=getIndexes(model,rxns,'rxns');
+% Add genes to model
+geneList = getGenesFromGrRules(grRules);
+genesToAdd.genes=setdiff(geneList,model.genes); % Only keep the genes that are not yet part of the model.genes.
+if ~isempty(genesToAdd.genes)
+    model=addGenesRaven(model,genesToAdd); % Add genes
 end
+
+% Find reaction and gene indices
+idx=getIndexes(model,rxns,'rxns');
 
 % Change gene associations
 if replace==true % Replace old gene associations

core/checkModelStruct.m

@@ -124,9 +124,7 @@ function checkModelStruct(model,throwErrors,trimWarnings)
         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 = getGenesFromGrRules(model.grRules);
         geneList = setdiff(unique(geneList),model.genes);
         if ~isempty(geneList)
             problemGrRules = model.rxns(contains(model.grRules,geneList));

core/checkRxn.m

@@ -14,11 +14,11 @@
 %
 %   report
 %       reactants   array with reactant indexes
-%       canMake     boolean array, true if the corresponding reactant can be
-%                   synthesized
+%       canMake     boolean array, true if the corresponding reactant can
+%                   be synthesized by the rest of the metabolic network
 %       products    array with product indexes
-%       canConsume  boolean array, true if the corresponding reactant can
-%                   be consumed
+%       canConsume  boolean array, true if the corresponding product can
+%                   be consumed by the rest of the metabolic network
 %
 % Usage: report=checkRxn(model,rxn,cutoff,revDir,printReport)
 

core/replaceMets.m

@@ -1,4 +1,4 @@
-function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose)
+function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose,identifiers)
 % replaceMets
 %   Replaces metabolite names and annotation with replacement metabolite
 %   that is already in the model. If this results in duplicate metabolites,
@@ -13,42 +13,59 @@
 %   verbose         logical whether to print the ids of reactions that
 %                   involve the replaced metabolite (optional, default
 %                   false)
+%   identifiers     true if 'metabolite' and 'replacement' refer to
+%                   metabolite identifiers instead of metabolite names
+%                   (optional, default false)
 % 
 % Output:
 %   model           model structure with selected metabolites replaced
 %   removedRxns     identifiers of duplicate reactions that were removed
 %   idxDuplRxns     index of removedRxns in original model
 %
 % Note: This function is useful when the model contains both 'oxygen' and
-% 'o2' as metabolites. 
+% 'o2' as metabolite names. If 'oxygen' and 'o2' are identifiers instead,
+% then the 'identifiers' flag should be set to true.
 %
 % Usage: [model, removedRxns, idxDuplRxns] = replaceMets(model, metabolite, replacement, verbose)
 
 metabolite=char(metabolite);
 replacement=char(replacement);
 
-if nargin<4
+if nargin<4 || isempty(verbose)
     verbose=false;
 end
+if nargin<5
+    identifiers = false;
+end
 
 % Find occurence of replacement metabolites. Annotation will be taken from
-% first metabolite found. Metabolite ID from replacement will be used where
-% possible.
-repIdx = find(strcmp(replacement,model.metNames));
+% first metabolite found.
+if identifiers
+    repIdx = find(strcmp(replacement,model.mets));
+else
+    repIdx = find(strcmp(replacement,model.metNames));
+end
 if isempty(repIdx)
-    error('The replacement metabolite name cannot be found in the model.');
+    error('The replacement metabolite cannot be found in the model.');
 end
 
 % Change name and information from metabolite to replacement metabolite
-metIdx = find(strcmp(metabolite,model.metNames));
+if identifiers
+    metIdx = find(strcmp(metabolite,model.mets));
+else
+    metIdx = find(strcmp(metabolite,model.metNames));
+end
 if isempty(metIdx)
-    error('The to-be-replaced metabolite name cannot be found in the model.');
+    error('The to-be-replaced metabolite cannot be found in the model.');
 end
+
+rxnsWithMet = find(model.S(metIdx,:));
 if verbose==true
-    fprintf('\n\nThe following reactions contain the replaced metabolite as reactant:\n')
-    fprintf(strjoin(model.rxns(find(model.S(metIdx,:))),'\n'))
+    fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n')
+    fprintf(strjoin(model.rxns(rxnsWithMet),'\n'))
     fprintf('\n')
 end
+
 model.metNames(metIdx) = model.metNames(repIdx(1));
 if isfield(model,'metFormulas')
     model.metFormulas(metIdx) = model.metFormulas(repIdx(1));
@@ -68,24 +85,32 @@
 if isfield(model,'metSmiles')
     model.metSmiles(metIdx) = model.metSmiles(repIdx(1));
 end
-% Run through replacement metabolites and their compartments. If any of the
-% to-be-replaced metabolites is already present (checked by
-% metaboliteName[compartment], then the replacement metabolite is kept and
-% the to-be-replace metabolite ID deleted.
-
-% Build list of metaboliteName[compartment]
-metCompsN =cellstr(num2str(model.metComps));
-map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps);
-metCompsN = map.values(metCompsN);
-metCompsN = strcat(lower(model.metNames),'[',metCompsN,']');
 
 idxDelete=[];
-for i = 1:length(repIdx)
-    metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN));
-    if length(metCompsNidx)>1
-        for j = 2:length(metCompsNidx)
-            model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:);
-            idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete
+if identifiers
+    originalStoch = model.S(metIdx,rxnsWithMet);
+    model.S(repIdx,rxnsWithMet) = originalStoch;
+    model.S(metIdx,rxnsWithMet) = 0;
+    idxDelete = metIdx;
+else
+    % Run through replacement metabolites and their compartments. If any of the
+    % to-be-replaced metabolites is already present (checked by
+    % metaboliteName[compartment], then the replacement metabolite is kept and
+    % the to-be-replace metabolite ID deleted.
+
+    % Build list of metaboliteName[compartment]
+    metCompsN =cellstr(num2str(model.metComps));
+    map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps);
+    metCompsN = map.values(metCompsN);
+    metCompsN = strcat(lower(model.metNames),'[',metCompsN,']');
+
+    for i = 1:length(repIdx)
+        metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN));
+        if length(metCompsNidx)>1
+            for j = 2:length(metCompsNidx)
+                model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:);
+                idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete
+            end
         end
     end
 end

doc/core/addExchangeRxns.html

@@ -126,37 +126,41 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
 0062     model.eccodes=[model.eccodes;filler];
 0063 <span class="keyword">end</span>
 0064 <span class="keyword">if</span> isfield(model,<span class="string">'subSystems'</span>)
-0065     model.subSystems=[model.subSystems;filler];
-0066 <span class="keyword">end</span>
-0067 <span class="keyword">if</span> isfield(model,<span class="string">'grRules'</span>)
-0068     model.grRules=[model.grRules;filler];
-0069 <span class="keyword">end</span>
-0070 <span class="keyword">if</span> isfield(model,<span class="string">'rxnFrom'</span>)
-0071     model.rxnFrom=[model.rxnFrom;filler];
-0072 <span class="keyword">end</span>
-0073 <span class="keyword">if</span> isfield(model,<span class="string">'rxnMiriams'</span>)
-0074     model.rxnMiriams=[model.rxnMiriams;filler];
-0075 <span class="keyword">end</span>
-0076 <span class="keyword">if</span> isfield(model,<span class="string">'rxnGeneMat'</span>)
-0077     model.rxnGeneMat=[model.rxnGeneMat;sparse(numel(J),numel(model.genes))];
-0078 <span class="keyword">end</span>
-0079 <span class="keyword">if</span> isfield(model,<span class="string">'rxnComps'</span>)
-0080     model.rxnComps=[model.rxnComps;ones(numel(J),1)];
-0081     fprintf(<span class="string">'NOTE: The exchange reactions are assigned to the first compartment\n'</span>);
+0065     fillerSub = filler;
+0066     <span class="keyword">if</span> iscell(model.subSystems(1,1))
+0067         fillerSub = repmat({fillerSub},numel(J),1);
+0068     <span class="keyword">end</span>
+0069     model.subSystems=[model.subSystems;fillerSub];
+0070 <span class="keyword">end</span>
+0071 <span class="keyword">if</span> isfield(model,<span class="string">'grRules'</span>)
+0072     model.grRules=[model.grRules;filler];
+0073 <span class="keyword">end</span>
+0074 <span class="keyword">if</span> isfield(model,<span class="string">'rxnFrom'</span>)
+0075     model.rxnFrom=[model.rxnFrom;filler];
+0076 <span class="keyword">end</span>
+0077 <span class="keyword">if</span> isfield(model,<span class="string">'rxnMiriams'</span>)
+0078     model.rxnMiriams=[model.rxnMiriams;filler];
+0079 <span class="keyword">end</span>
+0080 <span class="keyword">if</span> isfield(model,<span class="string">'rxnGeneMat'</span>)
+0081     model.rxnGeneMat=[model.rxnGeneMat;sparse(numel(J),numel(model.genes))];
 0082 <span class="keyword">end</span>
-0083 <span class="keyword">if</span> isfield(model,<span class="string">'rxnNotes'</span>)
-0084     model.rxnNotes=[model.rxnNotes;filler];
-0085 <span class="keyword">end</span>
-0086 <span class="keyword">if</span> isfield(model,<span class="string">'rxnReferences'</span>)
-0087     model.rxnReferences=[model.rxnReferences;filler];
-0088 <span class="keyword">end</span>
-0089 <span class="keyword">if</span> isfield(model,<span class="string">'rxnConfidenceScores'</span>)
-0090     model.rxnConfidenceScores=[model.rxnConfidenceScores;NaN(numel(J),1)];
-0091 <span class="keyword">end</span>
-0092 <span class="keyword">if</span> isfield(model,<span class="string">'rxnDeltaG'</span>)
-0093     model.rxnDeltaG=[model.rxnDeltaG;NaN(numel(J),1)];
-0094 <span class="keyword">end</span>
-0095 <span class="keyword">end</span></pre></div>
+0083 <span class="keyword">if</span> isfield(model,<span class="string">'rxnComps'</span>)
+0084     model.rxnComps=[model.rxnComps;ones(numel(J),1)];
+0085     fprintf(<span class="string">'NOTE: The exchange reactions are assigned to the first compartment\n'</span>);
+0086 <span class="keyword">end</span>
+0087 <span class="keyword">if</span> isfield(model,<span class="string">'rxnNotes'</span>)
+0088     model.rxnNotes=[model.rxnNotes;filler];
+0089 <span class="keyword">end</span>
+0090 <span class="keyword">if</span> isfield(model,<span class="string">'rxnReferences'</span>)
+0091     model.rxnReferences=[model.rxnReferences;filler];
+0092 <span class="keyword">end</span>
+0093 <span class="keyword">if</span> isfield(model,<span class="string">'rxnConfidenceScores'</span>)
+0094     model.rxnConfidenceScores=[model.rxnConfidenceScores;NaN(numel(J),1)];
+0095 <span class="keyword">end</span>
+0096 <span class="keyword">if</span> isfield(model,<span class="string">'rxnDeltaG'</span>)
+0097     model.rxnDeltaG=[model.rxnDeltaG;NaN(numel(J),1)];
+0098 <span class="keyword">end</span>
+0099 <span class="keyword">end</span></pre></div>
 <hr><address>Generated by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
 </body>
 </html>