fluxVariability

PURPOSE ^

fluxVariability Performs flux variablity analysis

SYNOPSIS ^

function [minFlux,maxFlux,Vmin,Vmax] = fluxVariability(model,optPercentage,osenseStr,rxnNameList,verbFlag, allowLoops)

DESCRIPTION ^

fluxVariability Performs flux variablity analysis

 [minFlux,maxFlux] = fluxVariability(model,optPercentage,osenseStr,rxnNameList,verbFlag, allowLoops)

INPUT
 model             COBRA model structure

OPTIONAL INPUTS
 optPercentage     Only consider solutions that give you at least a certain
                   percentage of the optimal solution (Default = 100
                   or optimal solutions only)
 osenseStr         Objective sense ('min' or 'max') (Default = 'max')
 rxnNameList       List of reactions for which FVA is performed
                   (Default = all reactions in the model)
 verbFlag          Verbose output (opt, default false)
 allowLoops        Whether loops are allowed in solution. (Default = true)
                   See optimizeCbModel for description

OUTPUT
 minFlux           Minimum flux for each reaction
 maxFlux           Maximum flux for each reaction

OPTIONAL OUTPUT
 Vmin          Matrix of column flux vectors, where each column is a 
               separate minimization.
 Vmax          Matrix of column flux vectors, where each column is a 
               separate maximization.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [minFlux,maxFlux,Vmin,Vmax] = fluxVariability(model,optPercentage,osenseStr,rxnNameList,verbFlag, allowLoops)
0002 %fluxVariability Performs flux variablity analysis
0003 %
0004 % [minFlux,maxFlux] = fluxVariability(model,optPercentage,osenseStr,rxnNameList,verbFlag, allowLoops)
0005 %
0006 %INPUT
0007 % model             COBRA model structure
0008 %
0009 %OPTIONAL INPUTS
0010 % optPercentage     Only consider solutions that give you at least a certain
0011 %                   percentage of the optimal solution (Default = 100
0012 %                   or optimal solutions only)
0013 % osenseStr         Objective sense ('min' or 'max') (Default = 'max')
0014 % rxnNameList       List of reactions for which FVA is performed
0015 %                   (Default = all reactions in the model)
0016 % verbFlag          Verbose output (opt, default false)
0017 % allowLoops        Whether loops are allowed in solution. (Default = true)
0018 %                   See optimizeCbModel for description
0019 %
0020 %OUTPUT
0021 % minFlux           Minimum flux for each reaction
0022 % maxFlux           Maximum flux for each reaction
0023 %
0024 %OPTIONAL OUTPUT
0025 % Vmin          Matrix of column flux vectors, where each column is a
0026 %               separate minimization.
0027 % Vmax          Matrix of column flux vectors, where each column is a
0028 %               separate maximization.
0029 %
0030 
0031 % Markus Herrgard  8/21/06 Original code.
0032 % Ronan Fleming   01/20/10 Take the extremal flux from the flux vector,
0033 %                          not from the objective since this is invariant
0034 %                          to the value and sign of the coefficient
0035 % Ronan Fleming   27/09/10 Vmin,Vmax
0036 
0037 if (nargin < 2)
0038     optPercentage = 100;
0039 end
0040 if (nargin < 3)
0041     osenseStr = 'max';
0042 end
0043 if (nargin < 4)
0044     rxnNameList = model.rxns;
0045 end
0046 if (nargin < 5)
0047     verbFlag = false;
0048 end
0049 if (nargin < 6)
0050     allowLoops = true;
0051 end
0052 if (isempty(optPercentage))
0053     optPercentage = 100;
0054 end
0055 if (isempty(osenseStr))
0056     osenseStr = 'max';
0057 end
0058 if (isempty(rxnNameList))
0059     rxnNameList = model.rxns;
0060 end
0061 
0062 % LP solution tolerance
0063 global CBT_LP_PARAMS
0064 if (exist('CBT_LP_PARAMS', 'var'))
0065     if isfield(CBT_LP_PARAMS, 'objTol')
0066         tol = CBT_LP_PARAMS.objTol;
0067     else
0068         tol = 1e-6;
0069     end
0070     if isfield(CBT_LP_PARAMS, 'minNorm')
0071         minNorm = CBT_LP_PARAMS.minNorm;
0072     else
0073         minNorm = 0;
0074     end
0075 else
0076     tol = 1e-6;
0077     minNorm = 0;
0078 end
0079 
0080 % Determine constraints for the correct space (0-100% of the full space)
0081 if (sum(model.c ~= 0) > 0)
0082     hasObjective = true;
0083     optSol = optimizeCbModel(model,osenseStr, 0, allowLoops);
0084     if (optSol.stat > 0)
0085         objRxn = model.rxns(model.c~=0);
0086         if (strcmp(osenseStr,'max'))
0087             objValue = floor(optSol.f/tol)*tol*optPercentage/100;
0088         else
0089             objValue = ceil(optSol.f/tol)*tol*optPercentage/100;
0090         end
0091     else
0092         error('Infeasible problem - no optimal solution!');
0093     end
0094 else
0095     hasObjective = false;
0096 end
0097 
0098 if (verbFlag == 1)  
0099     h = waitbar(0,'Flux variability analysis in progress ...');
0100 end
0101 if (verbFlag > 1)
0102     fprintf('%4s\t%4s\t%10s\t%9s\t%9s\n','No','Perc','Name','Min','Max');
0103 end
0104 
0105 if (~isfield(model,'b'))
0106     model.b = zeros(size(model.S,1),1);
0107 end
0108 % Set up the general problem
0109 [nMets,nRxns] = size(model.S);
0110 rxnListFull = model.rxns;
0111 LPproblem.c = model.c;
0112 LPproblem.lb = model.lb;
0113 LPproblem.ub = model.ub;
0114 LPproblem.csense(1:nMets) = 'E';
0115 LPproblem.csense = LPproblem.csense';
0116 if hasObjective
0117     LPproblem.A = [model.S;columnVector(model.c)'];
0118     LPproblem.b = [model.b;objValue];
0119     if (strcmp(osenseStr,'max'))
0120         LPproblem.csense(end+1) = 'G';
0121     else
0122         LPproblem.csense(end+1) = 'L';
0123     end
0124 else
0125     LPproblem.A = model.S;
0126     LPproblem.b = model.b;
0127 end
0128 
0129 
0130 % %solve to generate initial basis
0131 LPproblem.osense = -1;
0132 tempSolution = solveCobraLP(LPproblem);
0133 LPproblem.basis = tempSolution.basis;
0134 
0135 % Loop through reactions
0136 maxFlux = zeros(length(rxnNameList), 1);
0137 minFlux = zeros(length(rxnNameList), 1);
0138 
0139 if length(minNorm)> 1 || minNorm > 0
0140     Vmin=zeros(nRxns,nRxns);
0141     Vmax=zeros(nRxns,nRxns);
0142     %minimizing the Euclidean norm gets rid of the loops, so there
0143     %is no need for a second slower MILP approach
0144     allowLoops=1;
0145 else
0146     Vmin=[];
0147     Vmax=[];
0148 end
0149 
0150 solutionPool = zeros(length(model.lb), 0); 
0151 if ~exist('matlabpool') || (matlabpool('size') == 0) %aka nothing is active
0152     m = 0;
0153     for i = 1:length(rxnNameList)
0154         if mod(i,10) == 0, clear mex, end
0155         if (verbFlag == 1),fprintf('iteration %d.  skipped %d\n', i, round(m));end
0156         LPproblem.c = zeros(nRxns,1);
0157         rxnBool=strcmp(rxnListFull,rxnNameList{i});
0158         LPproblem.c(rxnBool) = 1; %no need to set this more than 1
0159         % do LP always
0160         LPproblem.osense = -1;
0161         LPsolution = solveCobraLP(LPproblem);
0162         %take the maximum flux from the flux vector, not from the obj -Ronan
0163         maxFlux(i) = LPsolution.full(LPproblem.c~=0);
0164         
0165         %minimise the Euclidean norm of the optimal flux vector to remove
0166         %loops -Ronan
0167         if length(minNorm)> 1 || minNorm > 0
0168             QPproblem=LPproblem;
0169             QPproblem.lb(LPproblem.c~=0)=maxFlux(i)-1e-12;
0170             QPproblem.ub(LPproblem.c~=0)=maxFlux(i)+1e12;
0171             QPproblem.c(:)=0;
0172             %Minimise Euclidean norm using quadratic programming
0173             if length(minNorm)==1
0174                 minNorm=ones(nRxns,1)*minNorm;
0175             end
0176             QPproblem.F = spdiags(minNorm,0,nRxns,nRxns);
0177             %quadratic optimization
0178             solution = solveCobraQP(QPproblem);
0179             if isempty(solution.full)
0180                 pause(eps)
0181             end
0182             Vmax(:,rxnBool)=solution.full(1:nRxns,1);
0183         end
0184         
0185         LPproblem.osense = 1;
0186         LPsolution = solveCobraLP(LPproblem);
0187         %take the maximum flux from the flux vector, not from the obj -Ronan
0188         minFlux(i) = LPsolution.full(LPproblem.c~=0);
0189         
0190         %minimise the Euclidean norm of the optimal flux vector to remove
0191         %loops
0192         %minimise the Euclidean norm of the optimal flux vector to remove
0193         %loops
0194         if length(minNorm)> 1 || minNorm > 0
0195             QPproblem=LPproblem;
0196             QPproblem.lb(LPproblem.c~=0)=maxFlux(i)-1e-12;
0197             QPproblem.ub(LPproblem.c~=0)=maxFlux(i)+1e12;
0198             QPproblem.c(:)=0;
0199             QPproblem.F = spdiags(minNorm,0,nRxns,nRxns);
0200             %Minimise Euclidean norm using quadratic programming
0201             if length(minNorm)==1
0202                 minNorm=ones(nRxns,1)*minNorm;
0203             end
0204             QPproblem.F = spdiags(minNorm,0,nRxns,nRxns);
0205             %quadratic optimization
0206             solution = solveCobraQP(QPproblem);
0207             Vmin(:,rxnBool)=solution.full(1:nRxns,1);
0208         end
0209 
0210         
0211         if ~allowLoops
0212             if any( abs(LPproblem.c'*solutionPool - maxFlux(i)) < tol) % if any previous solutions are good enough.
0213                 % no need to do anything.
0214                 m = m+.5;
0215             else
0216                 LPproblem.osense = -1;
0217                 LPsolution = solveCobraMILP(addLoopLawConstraints(LPproblem, model));
0218                 maxFlux(i) = LPsolution.obj/1000;
0219               end
0220             if any( abs(LPproblem.c'*solutionPool - minFlux(i)) < tol)
0221                 m = m+.5;
0222                 % no need to do anything.
0223             else
0224                 LPproblem.osense = 1;
0225                 LPsolution = solveCobraMILP(addLoopLawConstraints(LPproblem, model));
0226                 minFlux(i) = LPsolution.obj/1000;
0227             end
0228         end
0229         if (verbFlag == 1)
0230             waitbar(i/length(rxnNameList),h);
0231         end
0232         if (verbFlag > 1)
0233             fprintf('%4d\t%4.0f\t%10s\t%9.3f\t%9.3f\n',i,100*i/length(rxnNameList),rxnNameList{i},minFlux(i),maxFlux(i));
0234         end
0235     end
0236 else % parallel job.  pretty much does the same thing.
0237     parfor i = 1:length(rxnNameList)
0238         %if mod(i,10) == 0, clear mex, end
0239         %if (verbFlag == 1),fprintf('iteration %d.  skipped %d\n', i, round(m));end
0240         c = zeros(nRxns,1);
0241         c(strcmp(rxnListFull,rxnNameList{i})) = 1000;
0242         if allowLoops % do LP
0243             LPsolution = solveCobraLP(struct(...
0244                 'A', LPproblem.A,...
0245                 'b', LPproblem.b,...
0246                 'lb', LPproblem.lb,...
0247                 'ub', LPproblem.ub,...
0248                 'csense', LPproblem.csense,...
0249                 'c',c,...
0250                 'osense',-1, ...
0251                 'basis', LPproblem.basis ...
0252             ));
0253             %take the maximum flux from the flux vector, not from the obj -Ronan
0254             maxFlux(i) = LPsolution.full(c~=0);
0255             %LPproblemb.osense = 1;
0256             LPsolution = solveCobraLP(struct(...
0257                 'A', LPproblem.A,...
0258                 'b', LPproblem.b,...
0259                 'lb', LPproblem.lb,...
0260                 'ub', LPproblem.ub,...
0261                 'csense', LPproblem.csense,...
0262                 'c',c,...
0263                 'osense',1, ... %only part that's different.
0264                 'basis', LPproblem.basis ...
0265             ));
0266             minFlux(i) = LPsolution.full(c~=0);
0267         else
0268             LPsolution = solveCobraMILP(addLoopLawConstraints(struct(...
0269                 'A', LPproblem.A,...
0270                 'b', LPproblem.b,...
0271                 'lb', LPproblem.lb,...
0272                 'ub', LPproblem.ub,...
0273                 'csense', LPproblem.csense,...
0274                 'c',c,...
0275                 'osense',-1 ...
0276             ), model));
0277             maxFlux(i) = LPsolution.obj/1000;
0278             
0279             LPsolution = solveCobraMILP(addLoopLawConstraints(struct(...
0280                 'A', LPproblem.A,...
0281                 'b', LPproblem.b,...
0282                 'lb', LPproblem.lb,...
0283                 'ub', LPproblem.ub,...
0284                 'csense', LPproblem.csense,...
0285                 'c',c,...
0286                 'osense',1 ...
0287             ), model));%
0288             minFlux(i) = LPsolution.obj/1000;
0289         end
0290     end
0291 end
0292     
0293     
0294 if (verbFlag == 1)
0295     if ( regexp( version, 'R20') )
0296             close(h);
0297     end
0298 end
0299 
0300 maxFlux = columnVector(maxFlux);
0301 minFlux = columnVector(minFlux);

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003