fix: getBlast, importModel, replaceMets, ravenCobraWrapper / feat: ge…

…tIndexes (#551)

* refactor: speedup readYAMLmodel

* fix: ravenCobraWrapper rxnReferences cobra=>raven

* feat: getIndexes ec-fields

* fix: replaceMets only contracts rxns with replMets

* fix: getBlast handles spaces in path

also more informative error message

* fix: importModel with empty confidenceScores

* refactor: checkInstallation duplicate warning

* fix: KEGG model empty subSystems

* doc: getIndexes error if model has duplicate IDs

* fix: getDiamond fails unittest

* fix: getIndexes move ec-field check

* fix: constructS avoid warning

"Colon operands must be real scalars"

* fix: randomSampling outputs full matrix

sparse matrix can give weird summary results when applying functions like min, max, mean etc.

* fix: getBlast if temppath has space char

core/constructS.m

@@ -120,7 +120,7 @@
                 coeff=1;
                 name=strtrim(metabolites{j});
             else
-                name=strtrim(metabolites{j}(space+1:end));
+                name=strtrim(metabolites{j}(space(1)+1:end));
             end
         end
 

core/contractModel.m

@@ -25,7 +25,7 @@
 %   NOTE: This code might not work for advanced grRules strings
 %         that involve nested expressions of 'and' and 'or'.
 %
-% Usage: [reducedModel, removedRxns, indexedDuplicateRxns]=contractModel(model,distReverse)
+% Usage: [reducedModel, removedRxns, indexedDuplicateRxns]=contractModel(model,distReverse,mets)
 
 if nargin<2
     distReverse=true;

core/getIndexes.m

@@ -2,23 +2,29 @@
 % getIndexes
 %   Retrieves the indexes for a list of reactions or metabolites
 %
-%   Input:
+% Input:
 %   model           a model structure
 %   objects         either a cell array of IDs, a logical vector with the
 %                   same number of elements as metabolites in the model,
 %                   of a vector of indexes
 %   type            'rxns', 'mets', or 'genes' depending on what to retrieve
 %                   'metnames' queries metabolite names, while 'metcomps'
 %                   allows to provide specific metabolites and their
-%                   compartments in the format metaboliteName[comp]
+%                   compartments in the format metaboliteName[comp]. If a
+%                   model.ec structure exists (GECKO 3), then also
+%                   'ecenzymes', 'ecrxns' and 'ecgenes' are allowed
 %   returnLogical   Sets whether to return a logical array or an array with
 %                   the indexes (optional, default false)
 %
-%   Output:
+% Output:
 %   indexes         can be a logical array or a double array depending on
 %                   the value of returnLogical
 %
-% 	Usage: indexes=getIndexes(model, objects, type, returnLogical)
+% Note: If 'ecenzymes', 'ecrxns' or 'ecgenes' are used with a GECKO 3
+% model, then the indexes are from the model.ec.enzymes, model.ec.rxns or
+% model.ec.genes fields, respectively.
+% 
+% Usage: indexes=getIndexes(model, objects, type, returnLogical)
 
 if nargin<4
     returnLogical=false;
@@ -62,12 +68,27 @@
             indexes=tempIndexes;
         end
         return %None of the remaining function needs to run if metcomps
+    case 'ecrxns'
+        if ~isfield(model,'ec')
+            error('Type %s cannot be used if no model.ec structure is present.',type)
+        end
+        searchIn=model.ec.rxns;
+    case 'ecenzymes'
+        if ~isfield(model,'ec')
+            error('Type %s cannot be used if no model.ec structure is present.',type)
+        end
+        searchIn=model.ec.enzymes;
+    case 'ecgenes'
+        if ~isfield(model,'ec')
+            error('Type %s cannot be used if no model.ec structure is present.',type)
+        end
+        searchIn=model.ec.genes;
     otherwise
-        if contains(objects,{'rxns','mets','metnames','metcomps','genes'})
-            error('The second and third parameter provided to run getIndexes are likely switched. Note that the second parameter is the object to find the index for, while the third parameter is the object type (''rxns'', ''mets'', ''metnames'', ''metcomps'' or ''genes'').')
+        if contains(objects,{'rxns','mets','metnames','metcomps','genes','ecgenes','ecenzymes','ecrxns'})
+            error('The second and third parameter provided to run getIndexes are likely switched. Note that the second parameter is the object to find the index for, while the third parameter is the object type (''rxns'', ''mets'', ''metnames'', ''metcomps'', ''genes'', ''ecgenes'', ''ecenzymes'' or ''ecrxns'').')
         else
-            error('Incorrect value of the ''type'' parameter. Allowed values are ''rxns'', ''mets'', ''metnames'', ''metcomps'' or ''genes''.');
-        end
+            error('Incorrect value of the ''type'' parameter. Allowed values are ''rxns'', ''mets'', ''metnames'', ''metcomps'', ''genes'', ''ecgenes'', ''ecenzymes'' or ''ecrxns''.');
+       end
 end
 
 if iscell(objects)
@@ -76,6 +97,9 @@
         if strcmpi(type,'metnames')
             indexes{i}=index;
         elseif ~isempty(index)
+            if length(index) > 1
+                error('There are multiple instances of object "%s" in the model, while "%s" type should be unique', objects{i}, type)
+            end
             indexes(i)=index;
         else
             error(['Could not find object ''' objects{i} ''' in the model']);

core/randomSampling.m

@@ -166,7 +166,6 @@
     [~, I]=ismember(model.rxns,originalRxns);
     solutions=zeros(numel(originalRxns),nSamples);
     solutions(I,:)=sols;
-    solutions=sparse(solutions);
 else
     solutions=[];
 end

core/replaceMets.m

@@ -1,21 +1,28 @@
-function model=replaceMets(model,metabolite,replacement,verbose)
+function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose)
 % replaceMets
 %   Replaces metabolite names and annotation with replacement metabolite
 %   that is already in the model. If this results in duplicate metabolites,
 %   the replacement metabolite will be kept, while the S matrix is updated
-%   to use the replacement metabolite instead.
+%   to use the replacement metabolite instead. At the end, contractModel is
+%   run to remove any duplicate reactions that might have occured.
 %
-%   model               a model structure
-%   metabolite          string with name of metabolite to be replace
-%   replacement         string with name of replacement metabolite
-%   verbose             logical whether to print the ids of reactions that
-%                       involve the replaced metabolite (optional, default
-%                       false)
+% Input:
+%   model           a model structure
+%   metabolite      string with name of metabolite to be replace
+%   replacement     string with name of replacement metabolite
+%   verbose         logical whether to print the ids of reactions that
+%                   involve the replaced metabolite (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
 %
-%   This function is useful when the model contains both 'oxygen' and 'o2'
-%   as metabolites.
+% Note: This function is useful when the model contains both 'oxygen' and
+% 'o2' as metabolites. 
 %
-% Usage: model=replaceMets(model,metabolite,replacement,verbose)
+% Usage: [model, removedRxns, idxDuplRxns] = replaceMets(model, metabolite, replacement, verbose)
 
 metabolite=char(metabolite);
 replacement=char(replacement);
@@ -116,5 +123,5 @@
 end
 
 % This could now have created duplicate reactions. Contract model.
-model=contractModel(model);
+model=contractModel(model,[],repIdx);
 end

doc/core/constructS.html

@@ -187,7 +187,7 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
 0120                 coeff=1;
 0121                 name=strtrim(metabolites{j});
 0122             <span class="keyword">else</span>
-0123                 name=strtrim(metabolites{j}(space+1:end));
+0123                 name=strtrim(metabolites{j}(space(1)+1:end));
 0124             <span class="keyword">end</span>
 0125         <span class="keyword">end</span>
 0126         

doc/core/contractModel.html

@@ -53,7 +53,7 @@ <h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="
    NOTE: This code might not work for advanced grRules strings
          that involve nested expressions of 'and' and 'or'.
 
- Usage: [reducedModel, removedRxns, indexedDuplicateRxns]=contractModel(model,distReverse)</pre></div>
+ Usage: [reducedModel, removedRxns, indexedDuplicateRxns]=contractModel(model,distReverse,mets)</pre></div>
 
 <!-- crossreference -->
 <h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
@@ -62,7 +62,7 @@ <h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^
 <li><a href="dispEM.html" class="code" title="function dispEM(string,throwErrors,toList,trimWarnings)">dispEM</a>	dispEM</li><li><a href="getIndexes.html" class="code" title="function indexes=getIndexes(model, objects, type, returnLogical)">getIndexes</a>	getIndexes</li><li><a href="removeReactions.html" class="code" title="function reducedModel=removeReactions(model,rxnsToRemove,removeUnusedMets,removeUnusedGenes,removeUnusedComps)">removeReactions</a>	removeReactions</li><li><a href="sortModel.html" class="code" title="function model=sortModel(model,sortReversible,sortMetName,sortReactionOrder)">sortModel</a>	sortModel</li><li><a href="standardizeGrRules.html" class="code" title="function [grRules,rxnGeneMat,indexes2check] = standardizeGrRules(model,embedded)">standardizeGrRules</a>	standardizeGrRules</li></ul>
 This function is called by:
 <ul style="list-style-image:url(../matlabicon.gif)">
-<li><a href="mergeCompartments.html" class="code" title="function [model, deletedRxns, duplicateRxns]=mergeCompartments(model,keepUnconstrained,deleteRxnsWithOneMet,distReverse)">mergeCompartments</a>	mergeCompartments</li><li><a href="replaceMets.html" class="code" title="function model=replaceMets(model,metabolite,replacement,verbose)">replaceMets</a>	replaceMets</li><li><a href="simplifyModel.html" class="code" title="function [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,deleteUnconstrained, deleteDuplicates, deleteZeroInterval, deleteInaccessible, deleteMinMax, groupLinear, constrainReversible, reservedRxns, suppressWarnings)">simplifyModel</a>	simplifyModel</li></ul>
+<li><a href="mergeCompartments.html" class="code" title="function [model, deletedRxns, duplicateRxns]=mergeCompartments(model,keepUnconstrained,deleteRxnsWithOneMet,distReverse)">mergeCompartments</a>	mergeCompartments</li><li><a href="replaceMets.html" class="code" title="function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose)">replaceMets</a>	replaceMets</li><li><a href="simplifyModel.html" class="code" title="function [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,deleteUnconstrained, deleteDuplicates, deleteZeroInterval, deleteInaccessible, deleteMinMax, groupLinear, constrainReversible, reservedRxns, suppressWarnings)">simplifyModel</a>	simplifyModel</li></ul>
 <!-- crossreference -->
 
 <h2><a name="_subfunctions"></a>SUBFUNCTIONS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
@@ -97,7 +97,7 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
 0025 <span class="comment">%   NOTE: This code might not work for advanced grRules strings</span>
 0026 <span class="comment">%         that involve nested expressions of 'and' and 'or'.</span>
 0027 <span class="comment">%</span>
-0028 <span class="comment">% Usage: [reducedModel, removedRxns, indexedDuplicateRxns]=contractModel(model,distReverse)</span>
+0028 <span class="comment">% Usage: [reducedModel, removedRxns, indexedDuplicateRxns]=contractModel(model,distReverse,mets)</span>
 0029 
 0030 <span class="keyword">if</span> nargin&lt;2
 0031     distReverse=true;