diff --git a/core/convertCharArray.m b/core/convertCharArray.m old mode 100644 new mode 100755 diff --git a/core/findRAVENroot.m b/core/findRAVENroot.m old mode 100644 new mode 100755 diff --git a/core/getAllowedBounds.m b/core/getAllowedBounds.m index 6510f3f6..fc915f4a 100755 --- a/core/getAllowedBounds.m +++ b/core/getAllowedBounds.m @@ -3,25 +3,30 @@ % Returns the minimal and maximal fluxes through each reaction. % % Input: -% model a model structure -% rxns either a cell array of reaction IDs, a logical vector with the -% same number of elements as reactions in the model, or a vector -% of reaction indexes (opt, default model.rxns) -% runParallel make use of MATLAB parallel pool to speed up calculations. Not -% beneficial if only a limited number of reactions are simulated. -% (opt, default true) +% model a model structure +% rxns either a cell array of reaction IDs, a logical vector +% with the same number of elements as reactions in the +% model, or a vector of reaction indexes (opt, default +% model.rxns) +% runParallel speed up calculations by parallel processing. This is +% not beneficial if allowed bounds are calculated for +% only a few reactions, as the overhead of parallel +% processing will take longer. It requires MATLAB +% Parallel Computing Toolbox. If this is not installed, +% the calculations will not be parallelized, regardless +% what is indicated as runParallel. (opt, default true) % % Output: -% minFluxes minimal allowed fluxes -% maxFluxes maximal allowed fluxes -% exitFlags exit flags for min/max for each of the reactions. True if it was -% possible to calculate a flux +% minFluxes minimal allowed fluxes +% maxFluxes maximal allowed fluxes +% exitFlags exit flags for min/max for each of the reactions. True +% if it was possible to calculate a flux % % NOTE: In cases where no solution can be calculated, NaN is returned. % % Usage: [minFluxes, maxFluxes, exitFlags] = getAllowedBounds(model, rxns, runParallel) -if nargin<2 +if nargin<2 || isempty(rxns) rxns = 1:numel(model.rxns); elseif ~islogical(rxns) && ~isnumeric(rxns) rxns = convertCharArray(rxns); @@ -30,78 +35,39 @@ if nargin<3 runParallel = true; end -if runParallel - addonList = matlab.addons.installedAddons; - if ~any(strcmpi(addonList.Name,'Parallel Computing Toolbox')) - disp('Cannot find MATLAB Parallel Computing Toolbox, process is not parallelized.') - runParallel = false; - end -end + +[ps, oldPoolAutoCreateSetting] = parallelPoolRAVEN(runParallel); minFluxes = zeros(numel(rxns),1); maxFluxes = zeros(numel(rxns),1); exitFlags = zeros(numel(rxns),2); c = zeros(numel(model.rxns),1); -p = 1; -progressbar('Calculating the minimal and maximal fluxes') -if runParallel - D = parallel.pool.DataQueue; - afterEach(D, @updateProgressParallel); - parfor i = 1:numel(rxns) - tmpModel = model; - tmpModel.c = c; - - % Get minimal flux - tmpModel.c(rxns(i)) = -1; - solMin = solveLP(tmpModel); - if ~isempty(solMin.f) - minFluxes(i) = solMin.x(rxns(i)); - else - minFluxes(i) = NaN; - end +PB = ProgressBar2(numel(rxns),'Running getAllowedBounds','cli'); +parfor i = 1:numel(rxns) + count(PB) + tmpModel = model; + tmpModel.c = c; - % Get maximal flux - tmpModel.c(rxns(i)) = 1; - solMax=solveLP(tmpModel); - exitFlags(i,:) = [solMin.stat solMax.stat]; - if ~isempty(solMax.f) - maxFluxes(i) = solMax.x(rxns(i)); - else - maxFluxes(i) = NaN; - end - send(D, i); + % Get minimal flux + tmpModel.c(rxns(i)) = -1; + solMin = solveLP(tmpModel); + if ~isempty(solMin.f) + minFluxes(i) = solMin.x(rxns(i)); + else + minFluxes(i) = NaN; end -else - for i = 1:numel(rxns) - tmpModel = model; - tmpModel.c = c; - % Get minimal flux - tmpModel.c(rxns(i)) = -1; - solMin=solveLP(tmpModel); - if ~isempty(solMin.f) - minFluxes(i) = solMin.x(rxns(i)); - else - minFluxes(i) = NaN; - end - - % Get maximal flux - tmpModel.c(rxns(i)) = 1; - solMax = solveLP(tmpModel); - exitFlags(i,:) = [solMin.stat solMax.stat]; - if ~isempty(solMax.f) - maxFluxes(i) = solMax.x(rxns(i)); - else - maxFluxes(i) = NaN; - end - progressbar(i/numel(rxns)) + % Get maximal flux + tmpModel.c(rxns(i)) = 1; + solMax=solveLP(tmpModel); + exitFlags(i,:) = [solMin.stat solMax.stat]; + if ~isempty(solMax.f) + maxFluxes(i) = solMax.x(rxns(i)); + else + maxFluxes(i) = NaN; end end -progressbar(1) % Make sure it closes - - function updateProgressParallel(~) - progressbar(p/numel(rxns)) - p = p + 1; - end +% Reset original Parallel setting +ps.Pool.AutoCreate = oldPoolAutoCreateSetting; end diff --git a/core/getGenesFromGrRules.m b/core/getGenesFromGrRules.m old mode 100644 new mode 100755 index 7f84fa5f..ac97c771 --- a/core/getGenesFromGrRules.m +++ b/core/getGenesFromGrRules.m @@ -31,7 +31,7 @@ end % check if the grRules use written or symbolic boolean operators -if any(contains(grRules,{'&','|'})) +if any(contains(grRules,{' & ',' | '})) % fix some potential missing spaces between parentheses and &/| grRules = regexprep(grRules,'\)&',') &'); % ")&" -> ") &" grRules = regexprep(grRules,'&\(','& ('); % "&(" -> "& (" @@ -50,11 +50,12 @@ end % extract list of genes from each reaction -rxnGenes = cellfun(@(r) unique(regexp(r,'[^&|\(\) ]+','match')),grRules,'UniformOutput',false); +rxnGenes = cellfun(@(r) regexprep(unique(strsplit(r,{' | ',' & '})),'[\(\) ]+',''),grRules,'UniformOutput',false); % construct new gene list nonEmpty = ~cellfun(@isempty,rxnGenes); genes = unique([rxnGenes{nonEmpty}]'); +genes(cellfun(@isempty,genes)) = []; if ~isempty(originalGenes) if ~isequal(sort(originalGenes), sort(genes)) diff --git a/core/getMinNrFluxes.m b/core/getMinNrFluxes.m index 25deba8c..5b24530f 100755 --- a/core/getMinNrFluxes.m +++ b/core/getMinNrFluxes.m @@ -135,7 +135,7 @@ prob=rmfield(prob,{'blx','bux','blc','buc'}); % Optimize the problem -res = optimizeProb(prob,params); +res = optimizeProb(prob,params,false); isFeasible=checkSolution(res); if ~isFeasible diff --git a/core/getModelFromHomology.m b/core/getModelFromHomology.m index 6f4adeb0..30aca2ac 100755 --- a/core/getModelFromHomology.m +++ b/core/getModelFromHomology.m @@ -115,6 +115,11 @@ if isfield(models{i},'geneMiriams') models{i}=rmfield(models{i},'geneMiriams'); end + %The geneFrom field also loses meaning if the genes are replaced by + %orthologs + if isfield(models{i},'geneFrom') + models{i}=rmfield(models{i},'geneFrom'); + end end %Check that genes do not begin with ( or end with ), as this makes problematic grRules diff --git a/core/mapCompartments.m b/core/mapCompartments.m index 0e93eff0..dfb8034f 100755 --- a/core/mapCompartments.m +++ b/core/mapCompartments.m @@ -69,7 +69,7 @@ %Check that there are no compartments in the rules that are not in the %geneScoreStructure. uComps=upper(geneScoreStructure.compartments); -J=[uComps;{'OTHER'}]; +J=[uComps,{'OTHER'}]; if ~isempty(setdiff([toKeep;toMerge],J)) EM='There are compartment in the rules that are not in geneScoreStructure.compartments'; diff --git a/core/parallelPoolRAVEN.m b/core/parallelPoolRAVEN.m new file mode 100755 index 00000000..36018ce2 --- /dev/null +++ b/core/parallelPoolRAVEN.m @@ -0,0 +1,53 @@ +function [ps, oldPoolAutoCreate] = parallelPoolRAVEN(runParallel) +% handleParallelRAVEN +% Called by RAVEN functions that support parallel processing, to confirm +% whether the MATLAB Parallel Computing Toolbox is installed. +% - The toolbox is installed, and runParallel == true, ==> a parallel +% pool is started. +% - The toolbox is installed, but runParallel == false, ==> the auto- +% creation of a parallel pool is disabled, to prevent "parfor" in +% the target function to start a pool anyway. +% - The toolbox is not installed, and runParallel == true, ==> a warning +% is displayed that parallel computer is not possible. +% - The toolbox is not installed, and runParallel == false, ==> the +% target runs as intended, as "parfor" will automatically run in serial +% mode instead. +% +% Input: +% runParallel logical, whether the target function (which calls +% parallelPoolRAVEN) should be run in parallel (opt, +% default true) +% +% Output: +% ps parallel settings structure that will be used by +% the target function +% oldPoolAutoCreate logical, to reset the original ps.Pool.AutoCreate +% setting once the target function has finished +% +% Use: [ps, oldPoolAutoCreate] = parallelPoolRAVEN(runParallel) + +if nargin<1 || isempty(runParallel) + runParallel = true; +end + +addonList = matlab.addons.installedAddons; +ps = []; oldPoolAutoCreate = []; +if ~any(strcmpi(addonList.Name,'Parallel Computing Toolbox')) + if runParallel % User wants parallel, but will not be possible + disp('Cannot find MATLAB Parallel Computing Toolbox, process is not parallelized.') + end +else + if ~runParallel % User has Parallel toolbox, but does not want pool to start. + % If pool is already running, then it will use the max workers + % anyway, but this is probably okay. + ps = parallel.Settings; + oldPoolAutoCreate = ps.Pool.AutoCreate; + ps.Pool.AutoCreate = false; + else + pool = gcp('nocreate'); + if isempty(pool) + parpool('IdleTimeout',120) + end + end +end +end \ No newline at end of file diff --git a/core/predictLocalization.m b/core/predictLocalization.m index ccca5587..27da539d 100755 --- a/core/predictLocalization.m +++ b/core/predictLocalization.m @@ -1005,3 +1005,122 @@ transportCost=sum(transportCost(I)); score=geneScore-transportCost; end + +% To avoid dependency on stats toolbox, use this alternative implementation +% of randsample, source: +% https://github.com/gpeyre/numerical-tours/blob/dacee30081c04ef5f67b26b387ead85f2b193af9/matlab/toolbox_signal/randsample.m +function y = randsample(n, k, replace, w) +%RANDSAMPLE Random sample, with or without replacement. +% Y = RANDSAMPLE(N,K) returns Y as a vector of K values sampled uniformly +% at random, without replacement, from the integers 1:N. +% +% Y = RANDSAMPLE(POPULATION,K) returns K values sampled uniformly at +% random, without replacement, from the values in the vector POPULATION. +% +% Y = RANDSAMPLE(...,REPLACE) returns a sample taken with replacement if +% REPLACE is true, or without replacement if REPLACE is false (the default). +% +% Y = RANDSAMPLE(...,true,W) returns a weighted sample, using positive +% weights W, taken with replacement. W is often a vector of probabilities. +% This function does not support weighted sampling without replacement. +% +% Example: Generate a random sequence of the characters ACGT, with +% replacement, according to specified probabilities. +% +% R = randsample('ACGT',48,true,[0.15 0.35 0.35 0.15]) +% +% See also RAND, RANDPERM. + +% Copyright 1993-2008 The MathWorks, Inc. +% $Revision: 1.1.4.3 $ $Date: 2008/12/01 08:09:34 $ + +if nargin < 2 + error('stats:randsample:TooFewInputs','Requires two input arguments.'); +elseif numel(n) == 1 + population = []; +else + population = n; + n = numel(population); + if length(population)~=n + error('stats:randsample:BadPopulation','POPULATION must be a vector.'); + end +end + +if nargin < 3 + replace = false; +end + +if nargin < 4 + w = []; +elseif ~isempty(w) + if length(w) ~= n + if isempty(population) + error('stats:randsample:InputSizeMismatch',... + 'W must have length equal to N.'); + else + error('stats:randsample:InputSizeMismatch',... + 'W must have the same length as the population.'); + end + else + p = w(:)' / sum(w); + end +end + +switch replace + +% Sample with replacement +case {true, 'true', 1} + if isempty(w) + y = ceil(n .* rand(k,1)); + else + [dum, y] = histc(rand(k,1),[0 cumsum(p)]); + end + +% Sample without replacement +case {false, 'false', 0} + if k > n + if isempty(population) + error('stats:randsample:SampleTooLarge',... + 'K must be less than or equal to N for sampling without replacement.'); + else + error('stats:randsample:SampleTooLarge',... + 'K must be less than or equal to the population size.'); + end + end + + if isempty(w) + % If the sample is a sizeable fraction of the population, + % just randomize the whole population (which involves a full + % sort of n random values), and take the first k. + if 4*k > n + rp = randperm(n); + y = rp(1:k); + + % If the sample is a small fraction of the population, a full sort + % is wasteful. Repeatedly sample with replacement until there are + % k unique values. + else + x = zeros(1,n); % flags + sumx = 0; + while sumx < k + x(ceil(n * rand(1,k-sumx))) = 1; % sample w/replacement + sumx = sum(x); % count how many unique elements so far + end + y = find(x > 0); + y = y(randperm(k)); + end + else + error('stats:randsample:NoWeighting',... + 'Weighted sampling without replacement is not supported.'); + end +otherwise + error('stats:randsample:BadReplaceValue',... + 'REPLACE must be either true or false.'); +end + +if ~isempty(population) + y = population(y); +else + y = y(:); +end +end diff --git a/core/printFluxes.m b/core/printFluxes.m index 8b484aa0..5fc6d9d4 100755 --- a/core/printFluxes.m +++ b/core/printFluxes.m @@ -63,7 +63,10 @@ function printFluxes(model, fluxes, onlyExchange, cutOffFlux, outputFile,outputS else metaboliteList=convertCharArray(metaboliteList); end -if size(fluxes,1)~=numel(model.rxns) +if isempty(fluxes) + EM='Empty vector of fluxes, solveLP possibly returned infeasible'; + dispEM(EM); +elseif size(fluxes,1)~=numel(model.rxns) EM='The number of fluxes and the number of reactions must be the same'; dispEM(EM); end diff --git a/core/printOrange.m b/core/printOrange.m old mode 100644 new mode 100755 diff --git a/core/randomSampling.m b/core/randomSampling.m index ec6b5410..1edb9714 100755 --- a/core/randomSampling.m +++ b/core/randomSampling.m @@ -1,7 +1,9 @@ -function [solutions, goodRxns]=randomSampling(model,nSamples,replaceBoundsWithInf,supressErrors,showProgress,goodRxns,minFlux) +function [solutions, goodRxns]=randomSampling(model,nSamples,replaceBoundsWithInf,supressErrors,runParallel,goodRxns,minFlux) % randomSampling -% Returns a number of random solutions +% Performs random sampling of the solution space, as described in Bordel +% et al. (2010) PLOS Compt Biol (doi:10.1371/journal.pcbi.1000859). % +% Input: % model a model structure % nSamples the number of solutions to return % (opt, default 1000) @@ -19,9 +21,11 @@ % as unlimited glucose uptake) or too strict % (such as too many and too narrow constraints) % (opt, default false) -% showProgress if true, it will display in the command window -% how many iterations have been done (opt, default -% false) +% runParallel speed up calculations by parallel processing. +% Requires MATLAB Parallel Computing Toolbox. If +% this is not installed, the calculations will +% not be parallelized, regardless what is +% indicated as runParallel. (opt, default true) % goodRxns double vector of indexes of those reactions % that are not involved in loops and can be used % as random objective functions, as generated by @@ -34,6 +38,8 @@ % large number of samples are taken, but this is % not always the case (opt, default false) % +% +% Output: % solutions matrix with the solutions % goodRxns double vector of indexes of those reactions % that are not involved in loops and can be used @@ -43,25 +49,27 @@ % random set of three reactions. For reversible reactions it randomly % chooses between maximizing and minimizing. % -% Usage: solutions=randomSampling(model,nSamples,replaceBoundsWithInf,supressErrors,showProgress,goodRxns,minFlux) +% Usage: solutions=randomSampling(model,nSamples,replaceBoundsWithInf,supressErrors,runParallel,goodRxns,minFlux) -if nargin<2 +if nargin<2 | isempty(nSamples) nSamples=1000; end -if nargin<3 +if nargin<3 | isempty(replaceBoundsWithInf) replaceBoundsWithInf=true; end -if nargin<4 +if nargin<4 | isempty(supressErrors) supressErrors=false; end -if nargin<5 - showProgress=false; +if nargin<5 | isempty(runParallel) + runParallel=true; end -if nargin<7 +if nargin<7 | isempty(minFlux) minFlux=false; end -nRxns=2; %Number of reactions in the objective function in each iteration +[ps, oldPoolAutoCreateSetting] = parallelPoolRAVEN(runParallel); + +nObjRxns=2; %Number of reactions in the objective function in each iteration %Simplify the model to speed stuff up a little. Keep original mapping originalRxns=model.rxns; @@ -86,87 +94,78 @@ 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),1); - if showProgress - fprintf('Prepare goodRxns list of reactions not involved in loops... 0%% complete'); - end - for i=1:numel(model.rxns) - if showProgress && rem(i,100) == 0 - progress=num2str(floor(100*i/numel(model.rxns))); - progress=pad(progress,3,'left'); - fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b%s%% complete',progress); - end - if goodRxns(i)==true - testModel=setParam(model,'eq',model.rxns(i),1000); - sol=solveLP(testModel,0,[],hsSol); - if ~isempty(sol.f) - goodRxns(abs(sol.x)>999)=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) - goodRxns(abs(sol.x)>999)=false; - end + goodRxns = true(numel(model.rxns),numel(model.rxns)); + goodRxns = num2cell(goodRxns,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); + 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 end + count(PB); end - if showProgress - fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b100%% complete\n'); - end - goodRxns=find(goodRxns); + goodRxns = cell2mat(goodRxns); + goodRxns = find(prod(goodRxns,2)); end %Reserve space for a solution matrix -sols=zeros(numel(model.rxns),nSamples); +sols = zeros(numel(model.rxns),nSamples); +sols = num2cell(sols,1); %Main loop -counter=1; -badSolutions=0; -if showProgress - itrStr=num2str(nSamples); - fprintf('Performing random sampling: ready with%s0/%s iterations',repmat(' ',1,length(itrStr)),itrStr); -end -while counter<=nSamples - if showProgress && rem(counter,100) == 0 - fprintf([repmat('\b',1,length(itrStr)*2) '\b\b\b\b\b\b\b\b\b\b\b\b%s/%s iterations'],pad(num2str(counter),length(itrStr),'left'),itrStr); - end - rxns=randsample(numel(goodRxns),nRxns); - model.c=zeros(numel(model.rxns),1); - multipliers=randsample([-1 1],nRxns,true); - multipliers(model.rev(goodRxns(rxns))==0)=1; - model.c(goodRxns(rxns))=rand(nRxns,1).*multipliers; - if true(minFlux) - sol=solveLP(model,1,[],hsSol); - else - sol=solveLP(model,0,[],hsSol); - end - if any(sol.x) && abs(sol.f)>10^-8 - sols(:,counter)=sol.x; - counter=counter+1; - badSolutions=0; - else - badSolutions=badSolutions+1; - %If it only finds bad solutions then throw an error. - if badSolutions==100 && supressErrors==false - error('The program is having problems finding non-zero solutions (ignoring reactions that might be involved in loops). Review the constraints on your model. Set supressErrors to true to ignore this error'); +PB = ProgressBar2(nSamples,'Performing random sampling','cli'); +parfor i=1:nSamples + badSolutions = 1; + tmpModel = model; + 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; + tmpModel.c(goodRxns(rxns)) = rand(nObjRxns,1).*multipliers; + if true(minFlux) + sol=solveLP(tmpModel,1,[],hsSol); + else + sol=solveLP(tmpModel,0,[],hsSol); + end + if any(sol.x) && abs(sol.f)>10^-8 + sols{i} = sol.x; + badSolutions = 0; + else + badSolutions = badSolutions+1; + %If it only finds bad solutions then throw an error. + if badSolutions == 100 && supressErrors==false + error('The program is having problems finding non-zero solutions (ignoring reactions that might be involved in loops). Review the constraints on your model. Set supressErrors to true to ignore this error'); + end end end -end -if showProgress - fprintf('\n') + count(PB); end %Map to original model +sols = cell2mat(sols); [~, I]=ismember(model.rxns,originalRxns); solutions=zeros(numel(originalRxns),nSamples); solutions(I,:)=sols; solutions=sparse(solutions); +% Reset original Parallel setting +ps.Pool.AutoCreate = oldPoolAutoCreateSetting; end %To use instead of the normal Matlab randsample function. This is in order diff --git a/core/runDynamicFBA.m b/core/runDynamicFBA.m old mode 100644 new mode 100755 index 64a15616..cb776d3d --- a/core/runDynamicFBA.m +++ b/core/runDynamicFBA.m @@ -76,7 +76,7 @@ biomassVec = biomass; timeVec(1) = 0; [~,hsSol]=solveLP(model,1); -fprintf('Running dynamic FBA analysis... 0%% complete'); +PB = ProgressBar2(nSteps,'Running dynamic FBA analysis','cli'); for stepNo = 1:nSteps % Run FBA [sol,hsSol] = solveLP(model,1,[],hsSol); @@ -107,10 +107,8 @@ model.lb(excInd) = -uptakeBound; timeVec(stepNo+1) = stepNo*timeStep; - progress=pad(num2str(floor(stepNo/nSteps*100)),3,'left'); - fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b%s%% complete',progress); + count(PB); end -fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\bCOMPLETE\n'); selNonZero = any(concentrationMatrix>0,2); concentrationMatrix = concentrationMatrix(selNonZero,:); diff --git a/core/runPhenotypePhasePlane.m b/core/runPhenotypePhasePlane.m old mode 100644 new mode 100755 index 9b4a17ea..2b1d9f73 --- a/core/runPhenotypePhasePlane.m +++ b/core/runPhenotypePhasePlane.m @@ -49,10 +49,8 @@ shadowPrices2 = zeros(nPts); [~,hsSol] = solveLP(model); % Calulate points -fprintf('Running PhPP analysis... 0%% complete'); +PB = ProgressBar2(nPts,'Running PhPP analysis','cli'); for i = 1:nPts %ind1 - progress=pad(num2str(floor(i/nPts*100)),3,'left'); - fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b%s%% complete',progress); for j = 1:nPts %ind2 model1 = setParam(model,'eq',controlRxn1,-1*ind1(i)); model1 = setParam(model1,'eq',controlRxn2,-1*ind2(j)); @@ -63,6 +61,7 @@ shadowPrices2(j,i) = fbasol.sPrice(metID2); end end + count(PB); end fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\bCOMPLETE\n'); diff --git a/core/runProductionEnvelope.m b/core/runProductionEnvelope.m old mode 100644 new mode 100755 index 163da119..46eea9d6 --- a/core/runProductionEnvelope.m +++ b/core/runProductionEnvelope.m @@ -33,12 +33,10 @@ targetUpperBound = nan(1,numel(biomassValues)); targetLowerBound = nan(1,numel(biomassValues)); -fprintf('Creating production envelope... 0%% complete'); +PB = ProgressBar2(length(biomassValues),'Creating production envelope','cli'); % Max/min for target production model = setParam(model,'obj',targetRxn,1); for i = 1:length(biomassValues) - progress=pad(num2str(floor(i/numel(biomassValues)*100)),3,'left'); - fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b%s%% complete',progress); model1 = setParam(model,'eq',biomassRxn,biomassValues(i)); sol = solveLP(model1,0,[],hsSol); if (sol.stat > 0) @@ -53,8 +51,8 @@ else targetLowerBound(i) = NaN; end + count(PB); end -fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\bCOMPLETE\n'); % Plot results plot([biomassValues fliplr(biomassValues)],[targetUpperBound fliplr(targetLowerBound)],'blue','LineWidth',2); diff --git a/core/runRobustnessAnalysis.m b/core/runRobustnessAnalysis.m old mode 100644 new mode 100755 index c494f9aa..7504f3a7 --- a/core/runRobustnessAnalysis.m +++ b/core/runRobustnessAnalysis.m @@ -50,16 +50,15 @@ redCost = zeros(nPoints,1); controlFlux = linspace(solMin,solMax,nPoints)'; -fprintf('Running robustness analysis... 0%% complete'); + +PB = ProgressBar2(length(controlFlux),'Running robustness analysis','cli'); for i=1:length(controlFlux) - progress=pad(num2str(floor(i/numel(controlFlux)*100)),3,'left'); - fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b%s%% complete',progress); modelControlled = setParam(baseModel,'eq',controlRxnIdx,controlFlux(i)); solControlled = solveLP(modelControlled); objFlux(i) = solControlled.x(logical(modelControlled.c)); redCost(i) = solControlled.rCost(controlRxnIdx); + count(PB); end -fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\bCOMPLETE\n'); if plotRedCost yyaxis right diff --git a/core/runSimpleOptKnock.m b/core/runSimpleOptKnock.m old mode 100644 new mode 100755 diff --git a/core/setParam.m b/core/setParam.m index b5f00c67..c1ae9c79 100755 --- a/core/setParam.m +++ b/core/setParam.m @@ -6,11 +6,15 @@ % paramType the type of parameter to set: % 'lb' Lower bound % 'ub' Upper bound -% 'eq' Both upper and lower bound (equality -% constraint) +% 'eq' Both upper and lower bound (equality constraint) % 'obj' Objective coefficient -% 'rev' Reversibility +% 'rev' Reversibility (only changes the model.rev fields, +% does not affect model.lb and model.ub) % 'var' Variance around measured bound +% 'unc' Unconstrained, set lower and upper bound to the +% default values (-1000 and 1000, or any other values +% that are defined in model.annotation.defaultLB and +% .defaultUB) % rxnList a cell array of reaction IDs or a vector with their % corresponding indexes % params a vector of the corresponding values @@ -24,7 +28,7 @@ % Usage: model=setParam(model, paramType, rxnList, params) paramType=convertCharArray(paramType); -if ~any(strcmpi(paramType,{'lb','ub','eq','obj','rev','var'})) +if ~any(strcmpi(paramType,{'lb','ub','eq','obj','rev','var','unc'})) EM=['Incorrect parameter type: "' paramType '"']; dispEM(EM); end @@ -73,7 +77,7 @@ if ~isempty(indexes) if contains(paramType,'obj') model.c=zeros(numel(model.c),1); % parameter is changed, not added - end + end for j=1:length(paramType) if strcmpi(paramType{j},'eq') model.lb(indexes(j))=params(j); @@ -101,6 +105,20 @@ model.ub(indexes(j)) = params(j) * (1+var/200); end end + if strcmpi(paramType{j},'unc') + if isfield(model.annotation,'defaultLB') + lb = model.annotation.defaultLB; + else + lb = -1000; + end + if isfield(model.annotation,'defaultUB') + ub = model.annotation.defaultUB; + else + ub = 1000; + end + model.lb(indexes(j)) = lb; + model.ub(indexes(j)) = ub; + end end end end \ No newline at end of file diff --git a/core/simplifyModel.m b/core/simplifyModel.m index bf970416..86ccffe4 100755 --- a/core/simplifyModel.m +++ b/core/simplifyModel.m @@ -126,7 +126,7 @@ in=any(reducedModel.b<0,2); out=any(reducedModel.b>0,2); I=speye(numel(reducedModel.mets)); - revS=[reducedModel.S,reducedModel.S(:,reducedModel.rev~=0)*-1 I(:,in) I(:,out)*-1]; + revS=[reducedModel.S,reducedModel.S(:,reducedModel.lb<0)*-1 I(:,in) I(:,out)*-1]; metUsage=sum(abs(revS')>0); onlyProducts=sum(revS'>0) == metUsage; diff --git a/doc/core/dispEM.html b/doc/core/dispEM.html index 9553e946..fc08cb9a 100644 --- a/doc/core/dispEM.html +++ b/doc/core/dispEM.html @@ -49,7 +49,7 @@
0001 function [minFluxes, maxFluxes, exitFlags]=getAllowedBounds(model,rxns,runParallel) @@ -69,108 +72,74 @@SOURCE CODE % Returns the minimal and maximal fluxes through each reaction. 0004 % 0005 % Input: -0006 % model a model structure -0007 % rxns either a cell array of reaction IDs, a logical vector with the -0008 % same number of elements as reactions in the model, or a vector -0009 % of reaction indexes (opt, default model.rxns) -0010 % runParallel make use of MATLAB parallel pool to speed up calculations. Not -0011 % beneficial if only a limited number of reactions are simulated. -0012 % (opt, default true) -0013 % -0014 % Output: -0015 % minFluxes minimal allowed fluxes -0016 % maxFluxes maximal allowed fluxes -0017 % exitFlags exit flags for min/max for each of the reactions. True if it was -0018 % possible to calculate a flux -0019 % -0020 % NOTE: In cases where no solution can be calculated, NaN is returned. -0021 % -0022 % Usage: [minFluxes, maxFluxes, exitFlags] = getAllowedBounds(model, rxns, runParallel) -0023 -0024 if nargin<2 -0025 rxns = 1:numel(model.rxns); -0026 elseif ~islogical(rxns) && ~isnumeric(rxns) -0027 rxns = convertCharArray(rxns); -0028 rxns = getIndexes(model,rxns, 'rxns'); -0029 end -0030 if nargin<3 -0031 runParallel = true; -0032 end -0033 if runParallel -0034 addonList = matlab.addons.installedAddons; -0035 if ~any(strcmpi(addonList.Name,'Parallel Computing Toolbox')) -0036 disp('Cannot find MATLAB Parallel Computing Toolbox, process is not parallelized.') -0037 runParallel = false; -0038 end -0039 end +0006 % model a model structure +0007 % rxns either a cell array of reaction IDs, a logical vector +0008 % with the same number of elements as reactions in the +0009 % model, or a vector of reaction indexes (opt, default +0010 % model.rxns) +0011 % runParallel speed up calculations by parallel processing. This is +0012 % not beneficial if allowed bounds are calculated for +0013 % only a few reactions, as the overhead of parallel +0014 % processing will take longer. It requires MATLAB +0015 % Parallel Computing Toolbox. If this is not installed, +0016 % the calculations will not be parallelized, regardless +0017 % what is indicated as runParallel. (opt, default true) +0018 % +0019 % Output: +0020 % minFluxes minimal allowed fluxes +0021 % maxFluxes maximal allowed fluxes +0022 % exitFlags exit flags for min/max for each of the reactions. True +0023 % if it was possible to calculate a flux +0024 % +0025 % NOTE: In cases where no solution can be calculated, NaN is returned. +0026 % +0027 % Usage: [minFluxes, maxFluxes, exitFlags] = getAllowedBounds(model, rxns, runParallel) +0028 +0029 if nargin<2 || isempty(rxns) +0030 rxns = 1:numel(model.rxns); +0031 elseif ~islogical(rxns) && ~isnumeric(rxns) +0032 rxns = convertCharArray(rxns); +0033 rxns = getIndexes(model,rxns, 'rxns'); +0034 end +0035 if nargin<3 +0036 runParallel = true; +0037 end +0038 +0039 [ps, oldPoolAutoCreateSetting] = parallelPoolRAVEN(runParallel); 0040 0041 minFluxes = zeros(numel(rxns),1); 0042 maxFluxes = zeros(numel(rxns),1); 0043 exitFlags = zeros(numel(rxns),2); 0044 c = zeros(numel(model.rxns),1); 0045 -0046 p = 1; -0047 progressbar('Calculating the minimal and maximal fluxes') -0048 if runParallel -0049 D = parallel.pool.DataQueue; -0050 afterEach(D, @updateProgressParallel); -0051 parfor i = 1:numel(rxns) -0052 tmpModel = model; -0053 tmpModel.c = c; -0054 -0055 % Get minimal flux -0056 tmpModel.c(rxns(i)) = -1; -0057 solMin = solveLP(tmpModel); -0058 if ~isempty(solMin.f) -0059 minFluxes(i) = solMin.x(rxns(i)); -0060 else -0061 minFluxes(i) = NaN; -0062 end -0063 -0064 % Get maximal flux -0065 tmpModel.c(rxns(i)) = 1; -0066 solMax=solveLP(tmpModel); -0067 exitFlags(i,:) = [solMin.stat solMax.stat]; -0068 if ~isempty(solMax.f) -0069 maxFluxes(i) = solMax.x(rxns(i)); -0070 else -0071 maxFluxes(i) = NaN; -0072 end -0073 send(D, i); -0074 end -0075 else -0076 for i = 1:numel(rxns) -0077 tmpModel = model; -0078 tmpModel.c = c; -0079 -0080 % Get minimal flux -0081 tmpModel.c(rxns(i)) = -1; -0082 solMin=solveLP(tmpModel); -0083 if ~isempty(solMin.f) -0084 minFluxes(i) = solMin.x(rxns(i)); -0085 else -0086 minFluxes(i) = NaN; -0087 end -0088 -0089 % Get maximal flux -0090 tmpModel.c(rxns(i)) = 1; -0091 solMax = solveLP(tmpModel); -0092 exitFlags(i,:) = [solMin.stat solMax.stat]; -0093 if ~isempty(solMax.f) -0094 maxFluxes(i) = solMax.x(rxns(i)); -0095 else -0096 maxFluxes(i) = NaN; -0097 end -0098 progressbar(i/numel(rxns)) -0099 end -0100 end -0101 progressbar(1) % Make sure it closes -0102 -0103 function updateProgressParallel(~) -0104 progressbar(p/numel(rxns)) -0105 p = p + 1; -0106 end -0107 end