fix: predictLocalization, setParam and copyToComps (#596)

* fix: predictLocalization with CELLO output

* fix: setParam bug

* WIP: mergeModels copyToComps

* fix: copyToComps keeps ...From fields

* chore: fitTasks correct flags to fillGaps

* fix: mergeModels and copyToComps ...From fields

* fix: addRxns make subSystems nested cell arrays

* refactor: exportModel speedup subSystem parsing

* feat: runRAVENtests for quick local tests

Author: edkerk

core/addRxns.m

@@ -350,15 +350,19 @@
 end
 
 if isfield(rxnsToAdd,'subSystems')
+    if ischar(rxnsToAdd.subSystems)
+        rxnsToAdd.subSystems = {{rxnsToAdd.subSystems}};
+    else
+        for i=1:numel(rxnsToAdd.subSystems)
+            if ischar(rxnsToAdd.subSystems{i})
+                rxnsToAdd.subSystems{i}=rxnsToAdd.subSystems(i);
+            end
+        end
+    end
     if numel(rxnsToAdd.subSystems)~=nRxns
         EM='rxnsToAdd.subSystems must have the same number of elements as rxnsToAdd.rxns';
         dispEM(EM);
     end
-    for i=1:numel(rxnsToAdd.subSystems)
-        if ischar(rxnsToAdd.subSystems{i})
-            rxnsToAdd.subSystems{i}=rxnsToAdd.subSystems(i);
-        end
-    end
     %Fill with standard if it doesn't exist
     if ~isfield(newModel,'subSystems')
         newModel.subSystems=celllargefiller;

core/copyToComps.m

@@ -25,24 +25,26 @@
 %
 % Usage: model=copyToComps(model,toComps,rxns,deleteOriginal,compNames,compOutside)
 
-if nargin<3
-    rxns=model.rxns;
-elseif ~islogical(rxns) && ~isnumeric(rxns)
-    rxns=convertCharArray(rxns);
+arguments
+    model (1,1) struct
+    toComps {emptyOrTextOrCellOfText}
+    rxns = model.rxns
+    deleteOriginal {emptyOrLogicalScalar} = false
+    compNames {emptyOrTextOrCellOfText} = toComps
+    compOutside {emptyOrTextOrCellOfText} = '';
 end
-if nargin<4
-    deleteOriginal=false;
+
+if nargin >= 3 && ~islogical(rxns) && ~isnumeric(rxns)
+    rxns = convertCharArray(rxns);
 end
-if nargin<5
-    compNames=toComps;
-else
+if nargin >= 5
     compNames=convertCharArray(compNames);
 end
-if nargin<6
-    compOutside=cell(numel(toComps),1);
-    compOutside(:)={''};
-else
+if nargin >= 6
     compOutside=convertCharArray(compOutside);
+    if length(compOutside) ~= length(compNames)
+        error('compOutside and compNames should be of equal size.');
+    end
 end
 
 originalID=model.id;
@@ -79,15 +81,20 @@
         modelToAdd.compMiriams=modelToAdd.compMiriams(J);
     end
     modelToAdd.metComps=ones(numel(modelToAdd.mets),1);
-    
+    if isfield(modelToAdd,'metFrom')
+        modelToAdd = rmfield(modelToAdd,'metFrom');
+    end
+    if isfield(modelToAdd,'rxnFrom')
+        modelToAdd = rmfield(modelToAdd,'rxnFrom');
+    end
+    if isfield(modelToAdd,'geneFrom')
+        modelToAdd = rmfield(modelToAdd,'geneFrom');
+    end
+  
     %Merge the models
-    model=mergeModels({model;modelToAdd},'metNames');
+    model=mergeModels({model;modelToAdd},'metNames',[],true);
 end
 
-model=rmfield(model,'rxnFrom');
-model=rmfield(model,'metFrom');
-model=rmfield(model,'geneFrom');
-
 if deleteOriginal==true
     model=removeReactions(model,rxns,true,true,true); %Also delete unused compartments
 end

core/fitTasks.m

@@ -266,7 +266,7 @@
             %Only do gap-filling if it cannot be solved
             failed=false;
             try
-                [~, ~, newRxns, newModel, exitFlag]=fillGaps(tModel,refModel,false,true,supressWarnings,rxnScores,params);
+                [~, ~, 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'];
                     dispEM(EM,false);

core/mergeModels.m

@@ -1,53 +1,83 @@
-function model=mergeModels(models,metParam,supressWarnings)
+function model=mergeModels(models,metParam,supressWarnings,copyToComps)
 % mergeModels
 %   Merges models into one model structure. Reactions are added without any
 %   checks, so duplicate reactions might appear. Metabolites are matched by
 %   their name and compartment (metaboliteName[comp]), while genes are
 %   matched by their name.
 %
+% Input:
 %   models          a cell array with model structures
-%   metParam        string specifying whether to refer to metabolite name
-%                   (metNames) or ID (mets) for matching (default, metNames)
-%   supressWarnings true if warnings should be supressed (optional, default
-%                   false)
+%   metParam        string metabolite name ('metNames') or ID ('mets') are
+%                   used for matching (optional, default 'metNames')
+%   supressWarnings logical whether warnings should be supressed (optional,
+%                   default false)
+%   copyToComps     logical whether mergeModels is run via copyToComps
+%                   (optional, default false)
 %
-%   model     a model structure with the merged model. Follows the structure
-%             of normal models but also has 'rxnFrom/metFrom/geneFrom' fields
-%             to indicate from which model each reaction/metabolite/gene was
-%             taken
+% Output:
+%   model           a model structure with the merged model. Follows the
+%                   structure of normal models but also has 'rxnFrom/
+%                   metFrom/geneFrom' fields to indicate from which model
+%                   each reaction/metabolite/gene was taken. If the model
+%                   already has 'rxnFrom/metFrom/geneFrom' fields, then
+%                   these fields are not modified.
 %
 % Usage: model=mergeModels(models)
 
+arguments
+    models;
+    metParam {emptyOrTextScalar} = "metNames"
+    supressWarnings {emptyOrLogicalScalar} = false
+    copyToComps {emptyOrLogicalScalar} = false
+end
+
+metParam = char(metParam);
+
 %Just return the model
 if numel(models)<=1
     model=models{1};
     return;
 end
 
-if nargin<2
-    metParam='metNames';
-else
-    metParam=char(metParam);
-end
+hasMetFrom  = cellfun(@(s) isfield(s,'metFrom'), models);
+hasGeneFrom = cellfun(@(s) isfield(s,'geneFrom'), models);
+hasRxnFrom  = cellfun(@(s) isfield(s,'rxnFrom'), models);
 
-if nargin<3
-    supressWarnings=false;
+for i = 1:numel(models)
+    if copyToComps
+        if hasMetFrom(1)
+            models{2}.metFrom = repmat({''},numel(models{i}.mets),1);
+        end
+    elseif ~any(hasMetFrom)
+        models{i}.metFrom = repmat({models{i}.id},numel(models{i}.mets),1);
+    elseif ~hasMetFrom(i)
+        models{i}.metFrom = repmat({''},numel(models{i}.mets),1);
+    end
+    if copyToComps
+        if hasRxnFrom(1)
+            models{2}.rxnFrom = repmat({''},numel(models{i}.rxns),1);
+        end
+    elseif ~any(hasRxnFrom)
+        models{i}.rxnFrom = repmat({models{i}.id},numel(models{i}.rxns),1);
+    elseif ~hasRxnFrom(i)
+        models{i}.rxnFrom = repmat({''},numel(models{i}.rxns),1);
+    end
+    if copyToComps
+        if hasGeneFrom(1)
+            models{2}.geneFrom = repmat({''},numel(models{i}.genes),1);
+        end
+    elseif ~any(hasGeneFrom) && any(cellfun(@(s) isfield(s,'genes'), models))
+        models{i}.geneFrom = repmat({models{i}.id},numel(models{i}.genes),1);
+    elseif ~hasGeneFrom(i)
+        models{i}.geneFrom = repmat({''},numel(models{i}.genes),1);
+    end
 end
 
 %Add new functionality in the order specified in models
 model=models{1};
 model.id='MERGED';
 model.name='';
 
-model.rxnFrom=cell(numel(models{1}.rxns),1);
-model.rxnFrom(:)={models{1}.id};
-model.metFrom=cell(numel(models{1}.mets),1);
-model.metFrom(:)={models{1}.id};
-if isfield(models{1},'genes')
-    model.geneFrom=cell(numel(models{1}.genes),1);
-    model.geneFrom(:)={models{1}.id};
-end
-
 if isfield(model,'equations')
     model=rmfield(model,'equations');
 end
@@ -78,15 +108,15 @@
     end
 
     %Add all static stuff
-    rxnFrom=cell(numel(models{i}.rxns),1);
-    rxnFrom(:)={models{i}.id};
-    model.rxnFrom=[model.rxnFrom;rxnFrom];
-    model.rxns=[model.rxns;models{i}.rxns];
-    model.rxnNames=[model.rxnNames;models{i}.rxnNames];
-    model.lb=[model.lb;models{i}.lb];
-    model.ub=[model.ub;models{i}.ub];
-    model.c=[model.c;models{i}.c];
-    model.rev=[model.rev;models{i}.rev];
+    if any(hasRxnFrom) || (~copyToComps && ~any(hasRxnFrom))
+        model.rxnFrom  = [model.rxnFrom;  models{i}.rxnFrom];
+    end
+    model.rxns     = [model.rxns;     models{i}.rxns];
+    model.rxnNames = [model.rxnNames; models{i}.rxnNames];
+    model.lb       = [model.lb;       models{i}.lb];
+    model.ub       = [model.ub;       models{i}.ub];
+    model.c        = [model.c;        models{i}.c];
+    model.rev      = [model.rev;      models{i}.rev];
 
     if isfield(models{i},'subSystems')
         if isfield(model,'subSystems')
@@ -287,12 +317,12 @@
     end
 
     %Add static info on the metabolites
-    metFrom=cell(numel(metsToAdd),1);
-    metFrom(:)={models{i}.id};
-    model.metFrom=[model.metFrom;metFrom];
-    model.mets=[model.mets;models{i}.mets(metsToAdd)];
-    model.metNames=[model.metNames;models{i}.metNames(metsToAdd)];
-    model.b=[model.b;zeros(numel(metsToAdd),size(model.b,2))];
+    if any(hasMetFrom)
+        model.metFrom  = [model.metFrom;  models{i}.metFrom(metsToAdd)];
+    end
+    model.mets     = [model.mets;     models{i}.mets(metsToAdd)];
+    model.metNames = [model.metNames; models{i}.metNames(metsToAdd)];
+    model.b        = [model.b;        zeros(numel(metsToAdd),size(model.b,2))];
 
     if isfield(model,'unconstrained')
         if isfield(models{i},'unconstrained')
@@ -481,13 +511,13 @@
     if isfield(models{i},'genes')
         if ~isfield(model,'genes')
             %If there was no gene info before
-            model.genes=models{i}.genes;
-            model.rxnGeneMat=[sparse(numel(model.rxns),numel(models{i}.genes));models{i}.rxnGeneMat];
-            emptyGene=cell(numel(model.rxns),1);
-            emptyGene(:)={''};
-            model.grRules=[emptyGene;models{i}.grRules];
-            model.geneFrom=cell(numel(models{i}.genes),1);
-            model.geneFrom(:)={models{i}.id};
+            model.genes      = models{i}.genes;
+            model.rxnGeneMat = [sparse(numel(model.rxns),numel(models{i}.genes));models{i}.rxnGeneMat];
+            emptyGene        = repmat({''},numel(model.rxns),1);
+            model.grRules    = [emptyGene;models{i}.grRules];
+            if any(hasGeneFrom)
+                model.geneFrom   = models{i}.geneFrom;
+            end
 
             if isfield(models{i},'geneShortNames')
                 model.geneShortNames=models{i}.geneShortNames;
@@ -513,11 +543,9 @@
             %Only add extra gene info on new genes. This might not be
             %correct and should be changed later...
             if ~isempty(genesToAdd)
-                model.genes=[model.genes;models{i}.genes(genesToAdd)];
-                emptyGene=cell(numel(genesToAdd),1);
-                emptyGene(:)={models{i}.id};
-                model.geneFrom=[model.geneFrom;emptyGene];
-                model.rxnGeneMat=[model.rxnGeneMat sparse(size(model.rxnGeneMat,1),numel(genesToAdd))];
+                model.genes      = [model.genes;        models{i}.genes(genesToAdd)];
+                model.geneFrom   = [model.geneFrom;     models{i}.geneFrom(genesToAdd)];
+                model.rxnGeneMat = [model.rxnGeneMat    sparse(size(model.rxnGeneMat,1),numel(genesToAdd))];
 
                 if isfield(models{i},'geneShortNames')
                     if isfield(model,'geneShortNames')
@@ -587,7 +615,7 @@
             %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
             %now
-            [a, b]=ismember(models{i}.genes,model.genes);
+            a = ismember(models{i}.genes,model.genes);
 
             %Just a check
             if ~all(a)

core/parseTaskList.m

@@ -2,10 +2,11 @@
 % parseTaskList
 %   Parses a task list file.
 %
-%   inputFile       a task list in Excel format. The file must contain a
-%                   sheet named TASKS, which in turn may contain the
-%                   following column headers (note, all rows starting with
-%                   a non-empty cell are removed. The first row after that
+%   inputFile       a task list in either Excel (*.xlsx, with a sheet named
+%                   TASKS with all relevant content) or tab-delimited
+%                   (*.txt) format. The file may contain the following
+%                   column headers (note, all rows starting with a
+%                   non-empty cell are removed. The first row after that
 %                   is considered the headers):
 %                   ID
 %                       the only required header. Each task must have a

core/predictLocalization.m

@@ -621,7 +621,7 @@
         outModel.compNames(2)=GSS.compartments(1);
     end
 end
-outModel.compNames=[outModel.compNames;GSS.compartments(2:end)'];
+outModel.compNames=[outModel.compNames;GSS.compartments(2:end)];
 
 %Ugly little loop
 for i=1:numel(GSS.compartments)-1

core/setParam.m

@@ -68,7 +68,7 @@
     params(~Lia)=[];
     indexes(~Lia)=[];
     paramType(~Lia)=[];
-    dispEM('Reactions not present in model, will be ignored:',false,rxnLise(~Lia));
+    dispEM('Reactions not present in model, will be ignored:',false,rxnList(~Lia));
 end
 
 %Change the parameters