fix: findRAVENroot, subSystems, keggRxns, YAML without genes (#601)

* fix: store ravenPath as MATLAB preference

* fix: subSystems either nested or not

* fix: regenerate keggRxns, bug in subSystems

* fix: fully support subSystems as nested or not

* fix: YAML file allowed to have no genes

* chore: updateDocumentation

* fix: subSystems

* fix: ravenCobraWrapper subSystems, init-params

* fix: mergeModels giving useless output

* fix: ravenCobraWrapper

Author: edkerk

INIT/runINIT.m

@@ -1,4 +1,4 @@
-function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops)
+function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops,params)
 % runINIT
 %	Generates a model using the INIT algorithm, based on proteomics and/or
 %   transcriptomics and/or metabolomics and/or metabolic tasks. This is the 
@@ -38,6 +38,7 @@
 %                   problem significantly more computationally intensive to
 %                   solve (two more integer constraints per reversible reaction)
 %                   (optional, default false)
+%   params          parameter structure for use by optimizeProb
 %
 %   outModel        the resulting model structure
 %   deletedRxns     reactions which were deleted by the algorithm

core/FSEOF.m

@@ -119,7 +119,12 @@
         A2=char(model.rxns(num));                         %enzyme ID
         A3=char(model.rxnNames(num));                     %enzyme Name
         if isfield(model,'subSystems') && ~isempty(model.subSystems{num});
-            A4=char(strjoin(model.subSystems{num,1},';'));                   %Subsystems
+            if ~any(cellfun(@(x) iscell(x), model.subSystems));
+                subSys = cellfun(@(x) {x}, model.subSystems, 'uni', 0);
+            else
+                subSys = model.subSystems;
+            end
+            A4=char(strjoin(subSys{num},';'));                   %Subsystems
         else
             A4='';
         end

core/addExchangeRxns.m

@@ -63,8 +63,8 @@
 end
 if isfield(model,'subSystems')
     fillerSub = filler;
-    if iscell(model.subSystems(1,1))
-        fillerSub = repmat({fillerSub},numel(J),1);
+    if any(cellfun(@(x) iscell(x), model.subSystems))
+        fillerSub(:)={{''}};
     end
     model.subSystems=[model.subSystems;fillerSub];
 end

core/addRxns.m

@@ -227,10 +227,12 @@
 nOldRxns=numel(model.rxns);
 filler=cell(nRxns,1);
 filler(:)={''};
-cellfiller=cellfun(@(x) cell(0,0),filler,'UniformOutput',false);
+cellfiller=cell(nRxns,1);
+cellfiller(:)={{''}};
 largeFiller=cell(nOldRxns,1);
 largeFiller(:)={''};
-celllargefiller=cellfun(@(x) cell(0,0),largeFiller,'UniformOutput',false);
+largeCellFiller=cell(nOldRxns,1);
+largeCellFiller(:)={{''}};
 
 %***Add everything to the model except for the equations.
 if numel(rxnsToAdd.equations)~=nRxns
@@ -350,28 +352,45 @@
 end
 
 if isfield(rxnsToAdd,'subSystems')
+    % Has to be cell array
     if ischar(rxnsToAdd.subSystems)
-        rxnsToAdd.subSystems = {{rxnsToAdd.subSystems}};
+        rxnsToAdd.subSystems = {rxnsToAdd.subSystems};
+    end
+    % If all nested cells are 1x1, then unnest
+    if all(cellfun(@(x) iscell(x) && isscalar(x), rxnsToAdd.subSystems))
+        rxnsToAdd.subSystems = transpose([rxnsToAdd.subSystems{:}]);
+    end
+    % Cell array should now be as simple as possible. Check if it is nested
+    subSysRxnsNested = any(cellfun(@(x) iscell(x), rxnsToAdd.subSystems));
+    if isfield(newModel,'subSystems')
+        subSysModelNested = any(cellfun(@(x) iscell(x), newModel.subSystems));
     else
-        for i=1:numel(rxnsToAdd.subSystems)
-            if ischar(rxnsToAdd.subSystems{i})
-                rxnsToAdd.subSystems{i}=rxnsToAdd.subSystems(i);
-            end
+        subSysModelNested = subSysRxnsNested;
+        if subSysRxnsNested
+            newModel.subSystems=largeCellFiller;
+        else
+            newModel.subSystems=largeFiller;
         end
     end
+    if subSysRxnsNested && ~subSysModelNested
+        % Make all existing subSystems nested
+        newModel.subSystems = cellfun(@(x) {x}, newModel.subSystems, 'uni', 0);
+    elseif ~subSysRxnsNested && subSysModelNested
+        rxnsToAdd.subSystems = cellfun(@(x) {x}, rxnsToAdd.subSystems, 'uni', 0);
+    end
     if numel(rxnsToAdd.subSystems)~=nRxns
         EM='rxnsToAdd.subSystems must have the same number of elements as rxnsToAdd.rxns';
         dispEM(EM);
     end
-    %Fill with standard if it doesn't exist
-    if ~isfield(newModel,'subSystems')
-        newModel.subSystems=celllargefiller;
-    end
     newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems(:)];
 else
-    %Fill with standard if it doesn't exist
+    %Fill with standard if it does not exist
     if isfield(newModel,'subSystems')
-        newModel.subSystems=[newModel.subSystems;cellfiller];
+        if any(cellfun(@(x) iscell(x), newModel.subSystems))
+            newModel.subSystems=[newModel.subSystems;cellfiller];
+        else
+            newModel.subSystems=[newModel.subSystems;filler];
+        end
     end
 end
 if isfield(rxnsToAdd,'rxnMiriams')

core/addTransport.m

@@ -134,11 +134,15 @@
         model.eccodes=[model.eccodes;filler];
     end
     if isfield(model,'subSystems')
-        ssFiller=filler;
+        isNested = any(cellfun(@(x) iscell(x), model.subSystems));
+        ssFiller = filler;
         if isRev==1
-            ssFiller(:)={{['Transport between ' fromComp ' and ' toComps{i}]}};
+            ssFiller(:) = {['Transport between ' fromComp ' and ' toComps{i}]};
         else
-            ssFiller(:)={{['Transport from ' fromComp ' to ' toComps{i}]}};
+            ssFiller(:) = {['Transport from ' fromComp ' to ' toComps{i}]};
+        end
+        if isNested
+            ssFiller = cellfun(@(x) {x}, ssFiller, 'uni', 0);
         end
         model.subSystems=[model.subSystems;ssFiller];
     end

core/checkModelStruct.m

@@ -175,11 +175,11 @@ function checkModelStruct(model,throwErrors,trimWarnings)
     end
 end
 if isfield(model,'subSystems')
-    for i=1:numel(model.subSystems)
-        if ~iscell(model.subSystems{i,1})
-            EM='The "subSystems" field must be a cell array';
+    isNested  = any(cellfun(@(x) iscell(x), model.subSystems));
+    isCellStr = any(cellfun(@(x) ischar(x), model.subSystems));
+    if ~xor(isNested,isCellStr)
+        EM='The "subSystems" field must be a cell array of chars, *or* a cell array of cell arrays of chars';
             dispEM(EM,throwErrors);
-        end
     end
 end
 if isfield(model,'eccodes')

core/compareMultipleModels.m

@@ -92,16 +92,6 @@
 fprintf('*** Done \n\n')
 
 
-%% Flatten models' subSystems field
-% Convert from cell array of cells to cell array of strings
-% NOTE: this function currently only recognizes one subSystem per reaction;
-%       additional subSystems will be ignored!
-for i = 1:numel(models)
-    cells = cellfun(@iscell,models{i}.subSystems);
-    models{i}.subSystems(cells) = cellfun(@(s) s{1}, models{i}.subSystems(cells), 'UniformOutput', false);
-end
-
-
 %% Compare models structure & function based on high-dimensional methods
 % Compare number of reactions in each subsystem in each model using a heatmap
 field = 'subSystems';
@@ -110,6 +100,16 @@
     fprintf('\nWARNING: At least one model does not contain the field "subSystems". \n')
     fprintf('         Skipping subsystem comparison. \n\n')
 else
+    %% Flatten model subSystems field
+    % Convert from cell array of cells to cell array of strings
+    % NOTE: this function currently only recognizes one subSystem per reaction;
+    %       additional subSystems will be ignored!
+    for i = 1:numel(models)
+        if any(cellfun(@(x) iscell(x), models{i}.subSystems));
+            cells = cellfun(@iscell,models{i}.subSystems);
+            models{i}.subSystems(cells) = cellfun(@(s) s{1}, models{i}.subSystems(cells), 'UniformOutput', false);
+        end
+    end
     [id,compMat] = compareModelField(models,field);
     compStruct.subsystems.ID = id;
     compStruct.subsystems.matrix = compMat;

core/findRAVENroot.m

@@ -1,23 +1,27 @@
 function [ravenPath, prevDir] = findRAVENroot()
 % findRAVENroot
-%   Finds the root of the RAVEN directory, by searching for the path to
+%   Finds the root of the RAVEN directory, first by  by searching for the path to
 %   RAVEN2.png. Can also record the current directory, in case a function will
 %   use the ravenPath to navigate to a precise folder, and it should return to
 %   the previous directory afterwards. See e.g. optimizeProb calling glpk.
 
 ST=dbstack('-completenames');
 prevDir = pwd();
-ravenPath = ST(strcmp({ST.name},'findRAVENroot')).file;
-rootFound = 0;
-while rootFound == 0
-    isRoot = isfile(fullfile(ravenPath,'RAVEN2.png'));
-    if isRoot
-        rootFound = 1;
-    else
-        ravenPathOld = ravenPath;
-        ravenPath = fileparts(ravenPath);
-        if strcmp(ravenPathOld,ravenPath)
-            error('Cannot find the RAVEN root directory. Make sure you have not removed the RAVEN2.png file from your RAVEN installation.')
+if ispref('RAVEN','ravenPath')
+    ravenPath = getpref('RAVEN','ravenPath');
+else
+    ravenPath = ST(strcmp({ST.name},'findRAVENroot')).file;
+    rootFound = 0;
+    while rootFound == 0
+        isRoot = isfile(fullfile(ravenPath,'RAVEN2.png'));
+        if isRoot
+            rootFound = 1;
+        else
+            ravenPathOld = ravenPath;
+            ravenPath = fileparts(ravenPath);
+            if strcmp(ravenPathOld,ravenPath)
+                error('Cannot find the RAVEN root directory. Make sure you have not removed the RAVEN2.png file from your RAVEN installation.')
+            end
         end
     end
 end

core/fitTasks.m

@@ -34,7 +34,7 @@
 %   the result.
 %
 % Usage: [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,...
-%           rxnScores,taskStructure,params)
+%           rxnScores,taskStructure)
 
 if nargin<4
     printOutput=true;
@@ -48,9 +48,6 @@
 if nargin<6
     taskStructure=[];
 end
-if nargin<7
-    params=[];
-end
 
 if isempty(taskStructure) && ~isfile(inputFile)
     error('Task file %s cannot be found',string(inputFile));
@@ -268,7 +265,7 @@
             try
                 [~, ~, newRxns, newModel, exitFlag]=fillGaps(tModel,refModel,false,true,supressWarnings,rxnScores);
                 if exitFlag==-2
-                    EM=['"[' taskStructure(i).id '] ' taskStructure(i).description '" was aborted before reaching optimality. Consider increasing params.maxTime\n'];
+                    EM=['"[' taskStructure(i).id '] ' taskStructure(i).description '" was aborted before reaching optimality.\n'];
                     dispEM(EM,false);
                 end
             catch

core/mergeModels.m

@@ -82,6 +82,11 @@
     model=rmfield(model,'equations');
 end
 
+for i=1:numel(models)
+    if isfield(models{i},'subSystems')
+        models{i}.subSystems = cellfun(@(x) {x}, models{i}.subSystems, 'uni', 0);
+    end
+end
 for i=2:numel(models)
     %Add the model id to the rxn id id it already exists in the model (id
     %have to be unique) This is because it makes a '[]' string if no new
@@ -283,7 +288,6 @@
     if strcmpi(metParam,'mets')
     %Get the new metabolites from matching the models. Metabolites are matched by metabolite ID (model.mets).
 
-        oldMetComps=model.comps(model.metComps);
         oldMets=model.mets;
 
         if ~isempty(models{i}.mets)
@@ -613,7 +617,7 @@
             end
 
             %Remap the genes from the new model. The same thing as with
-            %mets; this is a wasteful way to do it but I don't care right
+            %mets; this is a wasteful way to do it but I do not care right
             %now
             a = ismember(models{i}.genes,model.genes);
 
@@ -637,4 +641,10 @@
 [grRules,rxnGeneMat] = standardizeGrRules(model,true);
 model.grRules = grRules;
 model.rxnGeneMat = rxnGeneMat;
+%Flatten subSystems if possible
+if isfield(model,'subSystems')
+    if all(cellfun(@(x) iscell(x) && isscalar(x), model.subSystems))
+        model.subSystems = transpose([model.subSystems{:}]);
+    end
+end
 end