fix: randomSampling, importModel, setParam

INIT/ftINIT.m

@@ -54,9 +54,7 @@
 %   useScoresForTasks   true if the calculated reaction scored should be 
 %                       used as weights when fitting to tasks (optional, default
 %                       true)
-%   paramsFT            parameter structure as used by getMILPParams. This
-%                       is for the fitTasks step. For the INIT algorithm,
-%                       see params (optional, default [])
+%   paramsFT            *obsolete option*
 %   verbose             if true, the MILP progression will be shown. 
 %                       (optional, default false)
 %

INIT/ftINITFillGaps.m

@@ -23,7 +23,7 @@
 %                       reference tModel. 
 %                       The solver will try to maximize the sum of the
 %                       scores for the included reactions
-%   params              parameter structure as used by getMILPParams
+%   params              *obsolete option*
 %   verbose             if true, the MILP progression will be shown. 
 %
 %   addedRxns           the rxns added

INIT/ftINITFillGapsForAllTasks.m

@@ -18,7 +18,7 @@
 %                   reactions (optional, default is -1 for all reactions)
 %   taskStructure   structure with the tasks, as from parseTaskList. If
 %                   this is supplied then inputFile is ignored
-%   params          parameter structure as used by getMILPParams
+%   params          *obsolete option*
 %   verbose         if true, the MILP progression will be shown. 
 %
 %

INIT/ftINITFillGapsMILP.m

@@ -12,7 +12,7 @@
 %                 with the same number of elements as reactions in the model,
 %                 of a vector of indexes for the reactions that should be
 %                 minimized (optional, default model.rxns)
-%   params        parameter structure as used by getMILPParams (optional)
+%   params        *obsolete option*
 %   scores        vector of weights for the reactions. Negative scores
 %                 should not have flux. Positive scores are not possible in this
 %                 implementation, and they are changed to max(scores(scores<0)).

INIT/getINITModel.m

@@ -1,10 +1,10 @@
-function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=getINITModel(refModel, tissue, celltype, hpaData, arrayData, metabolomicsData, taskFile, useScoresForTasks, printReport, taskStructure, params, paramsFT)
+function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=getINITModel(refModel, tissue, celltype, hpaData, arrayData, metabolomicsData, taskFile, useScoresForTasks, printReport, taskStructure)
 % getINITModel_legacy
 %   Generates a model using the INIT algorithm, based on proteomics and/or
 %   transcriptomics and/or metabolomics and/or metabolic tasks. This is the original 
 %   implementation of tINIT, which is replaced by ftINIT.
 %
-%   Input:
+% Input:
 %   refModel            a model structure. The model should be in the
 %                       closed form (no exchange reactions open). Import
 %                       using import(filename,false). If the model is not
@@ -53,15 +53,8 @@
 %                       as an alternative way to define tasks when Excel
 %                       sheets are not suitable. Overrides taskFile (optional,
 %                       default [])
-%   params              parameter structure as used by getMILPParams. This is
-%                       for the INIT algorithm. For the the MILP problems
-%                       solved to fit tasks, see paramsFT (optional, default [])
-%   paramsFT            parameter structure as used by getMILPParams. This is
-%                       for the fitTasks step. For the INIT algorithm, see
-%                       params (optional, default [])
 %
-%
-%   Output:
+% Output:
 %   model                   the resulting model structure
 %   metProduction           array that indicates which of the
 %                           metabolites in metabolomicsData that could be

INIT/runINIT.m

@@ -1,4 +1,4 @@
-function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops,params)
+function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops)
 % 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,8 +38,6 @@
 %                   problem significantly more computationally intensive to
 %                   solve (two more integer constraints per reversible reaction)
 %                   (optional, default false)
-%   params          parameter structure as used by getMILPParams (optional,
-%                   default [])
 %
 %   outModel        the resulting model structure
 %   deletedRxns     reactions which were deleted by the algorithm

core/consumeSomething.m

@@ -19,7 +19,7 @@
 %                   fluxes instead of the sum. Slower, but can be
 %                   used if the sum gives too many fluxes (optional, default
 %                   false)
-%   params          parameter structure as used by getMILPParams (optional)
+%   params          *obsolete option*
 %   ignoreIntBounds	true if internal bounds (including reversibility)
 %                   should be ignored. Exchange reactions are not affected.
 %                   This can be used to find unbalanced solutions which are

core/fillGaps.m

@@ -1,4 +1,4 @@
-function [newConnected, cannotConnect, addedRxns, newModel, exitFlag]=fillGaps(model,models,allowNetProduction,useModelConstraints,supressWarnings,rxnScores,params)
+function [newConnected, cannotConnect, addedRxns, newModel, exitFlag]=fillGaps(model,models,allowNetProduction,useModelConstraints,supressWarnings,rxnScores)
 % fillGaps
 %   Uses template model(s) to fill gaps in a model
 %
@@ -24,7 +24,6 @@
 %                       The solver will try to maximize the sum of the
 %                       scores for the included reactions (optional, default
 %                       is -1 for all reactions)
-%   params              parameter structure as used by getMILPParams (optional)
 %
 %   newConnected        cell array with the reactions that could be
 %                       connected. This is not calulated if

core/fitTasks.m

@@ -1,4 +1,4 @@
-function [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,rxnScores,taskStructure,params)
+function [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,rxnScores,taskStructure)
 % fitTasks
 %   Fills gaps in a model by including reactions from a reference model,
 %   so that the resulting model can perform all the tasks in a task list
@@ -17,7 +17,6 @@
 %                   reactions (optional, default is -1 for all reactions)
 %   taskStructure   structure with the tasks, as from parseTaskList. If
 %                   this is supplied then inputFile is ignored (optional)
-%   params          parameter structure as used by getMILPParams (optional)
 %
 %
 %   Output:

core/getMinNrFluxes.m

@@ -8,7 +8,7 @@
 %                 with the same number of elements as reactions in the model,
 %                 of a vector of indexes for the reactions that should be
 %                 minimized (optional, default model.rxns)
-%   params        parameter structure as used by getMILPParams (optional)
+%   params        *obsolete option*
 %   scores        vector of weights for the reactions. Negative scores
 %                 should not have flux. Positive scores are not possible in this
 %                 implementation, and they are changed to max(scores(scores<0)).

core/makeSomething.m

@@ -21,7 +21,7 @@
 %                   false)
 %   allowExcretion  allow for excretion of all other metabolites (optional,
 %                   default true)
-%   params          parameter structure as used by getMILPParams (optional)
+%   params          *obsolete option*
 %   ignoreIntBounds	true if internal bounds (including reversibility)
 %                   should be ignored. Exchange reactions are not affected.
 %                   This can be used to find unbalanced solutions which are

core/randomSampling.m

@@ -42,13 +42,18 @@
 % Output:
 %   solutions               matrix with the solutions
 %   goodRxns                double vector of indexes of those reactions
-%                           that are not involved in loops and can be used
-%                           as random objective functions
+%                           that are not involved in loops or always carry
+%                           zero flux and can be used as random objective
+%                           functions
 %
-% The solutions are generated by maximizing (with random weights) for a
-% random set of three reactions. For reversible reactions it randomly
+% Note: The solutions are generated by maximizing (with random weights) for
+% a random set of three reactions. For reversible reactions it randomly
 % chooses between maximizing and minimizing.
 %
+% If the model is a GECKO v3+ ecModel, then usage_prot reactions are not
+% selected for sampling, instead focusing on sampling from the metabolic
+% aspects that form the solution space.
+%
 % Usage: solutions = randomSampling(model, nSamples, replaceBoundsWithInf,...
 %                       supressErrors, runParallel, goodRxns, minFlux)
 
@@ -94,35 +99,39 @@
     warning('The model objective function cannot reach a non-zero value. This might be intended, so randomSampling will continue, but this could indicate problems with your model')
 end
 
-[~,hsSol]=solveLP(model);
 nRxns = numel(model.rxns);
 %Reactions which can be involved in loops should not be optimized for.
 %Check which reactions reach an arbitary high upper bound
 if nargin<6 || isempty(goodRxns)
-    goodRxns = true(numel(model.rxns),numel(model.rxns));
-    goodRxns = num2cell(goodRxns,1);
+    revRxns = transpose(find(model.lb < 0));
+    fwdRxns = transpose(find(model.ub > 0));
+    rxnList = [fwdRxns,revRxns];
+    objList = [ones(numel(fwdRxns),1);-ones(numel(revRxns),1)];
+    testSol = zeros(numel(objList),1);
     PB = ProgressBar2(nRxns,'Prepare goodRxns not involved in loops','cli');
-    parfor (i=1:nRxns)
-        testModel=setParam(model,'eq',model.rxns(i),1000);
-        sol=solveLP(testModel,0,[],hsSol);
+
+    parfor i = 1:numel(objList)
+        testModel=setParam(model,'obj',rxnList(i),objList(i));
+        sol=solveLP(testModel,0);
         if ~isempty(sol.f)
-            notGood = abs(sol.x)>999;
-            goodRxns{i}(notGood)=false;
-        else
-            %If the reaction is reversible, also check in that direction
-            if model.rev(i)
-                testModel=setParam(model,'eq',model.rxns(i),-1000);
-                sol=solveLP(testModel,0,[],hsSol);
-                if ~isempty(sol.f)
-                    notGood = abs(sol.x)>999;
-                    goodRxns{i}(notGood)=false;
-                end
-            end
+            testSol(i) = sol.x(rxnList(i));
         end
         count(PB);
     end
-    goodRxns = cell2mat(goodRxns);
-    goodRxns = find(prod(goodRxns,2));
+    testSol  = testSol.*objList;
+    goodRxns = true(nRxns,1);
+    % Filter out reactions that can reach (-)1000 (= involved in loop)
+    goodRxns(fwdRxns(testSol(1:numel(fwdRxns)) > 999)) = false;
+    goodRxns(revRxns(testSol(numel(fwdRxns)+1:end) > 999)) = false;
+    testSol(revRxns) = testSol(revRxns) + testSol(numel(fwdRxns)+1:end);
+    % Filter out reactions that cannot carry flux
+    goodRxns(testSol(1:nRxns) == 0) = false;
+    % In ecModels do not sample from usage_prot reactions
+    if isfield(model,'ec')
+        usageRxns = startsWith(model.rxns,'usage_prot_');
+        goodRxns(usageRxns) = false;
+    end
+    goodRxns = find(goodRxns);
 end
 
 %Reserve space for a solution matrix
@@ -139,13 +148,16 @@
         while lt(0,badSolutions)
             rxns = randsample(numel(goodRxns),nObjRxns);
             tmpModel.c = zeros(numel(tmpModel.rxns),1);
+
             multipliers = randsample([-1 1],nObjRxns,true);
-            multipliers(tmpModel.rev(goodRxns(rxns))==0) = 1;
+            multipliers(tmpModel.ub(goodRxns(rxns))<0) = 1;
+
             tmpModel.c(goodRxns(rxns)) = rand(nObjRxns,1).*multipliers;
+
             if true(minFlux)
-                sol=solveLP(tmpModel,1,[],hsSol);
+                sol=solveLP(tmpModel,1);
             else
-                sol=solveLP(tmpModel,0,[],hsSol);
+                sol=solveLP(tmpModel,0);
             end
             if any(sol.x) && abs(sol.f)>10^-8
                 sols{i} = sol.x;

core/setParam.m

@@ -123,4 +123,9 @@
         end
     end
 end
+if any(ismember(paramType,{'lb','ub','unc'}))
+    invalidBound = model.lb(indexes) > model.ub(indexes);
+    if any(invalidBound)
+        error(['Invalid set of bounds for reaction(s): ', strjoin(model.rxns(indexes(invalidBound)),', '), '.'])
+    end
 end