fluxvariability.m

使用基于约束的化学计量模型进行代谢通量可变性分析的例程

function [minFlux,maxFlux] = fluxVariability(model,optPercentage,osenseStr,rxnNameList,verbFlag)
%fluxVariability Performs flux variablity analysis
%
% [minFlux,maxFlux] = fluxVariability(model,optPercentage,osenseStr,rxnNameList,verbFlag)
%
% model             Model structure
% optPercentage     Only consider solutions that give you at least a certain
%                   percentage of the optimal solution (opt, default = 100
%                   or optimal solutions only)
% osenseStr         Objective sense ('min' or 'max') (opt, default 'max')
% rxnNameList       List of reactions for which FVA is performed (opt,
%                   default all reactions in the model)
% verbFlag          Verbose output (opt, default false)
%
% minFlux   Minimum flux for each reaction
% maxFlux   Maximum flux for each reaction
%
% Markus Herrgard 8/21/06 

if (nargin < 2)
    optPercentage = 100;
end
if (nargin < 3)
    osenseStr = 'max';
end
if (nargin < 4)
    rxnNameList = model.rxns;
end
if (nargin < 5)
    verbFlag = false;
end

if (isempty(optPercentage))
    optPercentage = 100;
end
if (isempty(osenseStr))
    osenseStr = 'max';
end
if (isempty(rxnNameList))
    rxnNameList = model.rxns;
end

% LP solution tolerance
if (exist('CBTLPTOL'))
    tol = CBTLPTOL;
else
    tol = 1e-6;
end

% Determine constraints for the correct space (0-100% of the full space)
if (sum(model.c ~= 0) > 0)
    hasObjective = true;
    optSol = optimizeCbModel(model,osenseStr);
    if (optSol.stat > 0)
        objRxn = model.rxns(model.c~=0);
        if (strcmp(osenseStr,'max'))
            objValue = floor(optSol.f/tol)*tol*optPercentage/100;
        else
            objValue = ceil(optSol.f/tol)*tol*optPercentage/100;
        end
    else
        error('Infeasible problem - no optimal solution!');
    end
else
    hasObjective = false;
end

if (verbFlag == 1)  
    h = waitbar(0,'Flux variability analysis in progress ...');
end
if (verbFlag > 1)
    fprintf('%4s\t%4s\t%10s\t%9s\t%9s\n','No','Perc','Name','Min','Max');
end

if (~isfield(model,'b'))
    model.b = zeros(size(model.S,1),1);
end
% Set up the general problem
[nMets,nRxns] = size(model.S);
LPproblem.c = model.c;
LPproblem.lb = model.lb;
LPproblem.ub = model.ub;
LPproblem.csense(1:nMets) = 'E';

if hasObjective
    LPproblem.A = [model.S;columnVector(model.c)'];
    LPproblem.b = [model.b;objValue];
    if (strcmp(osenseStr,'max'))
        LPproblem.csense(end+1) = 'G';
    else
        LPproblem.csense(end+1) = 'L';
    end
else
    LPproblem.A = model.S;
    LPproblem.b = model.b;
end

% Loop through reactions
for i = 1:length(rxnNameList)
    LPproblem.c = zeros(nRxns,1);
    LPproblem.c(strcmp(model.rxns,rxnNameList{i})) = 1;
    LPproblem.osense = -1;
    LPsolution = solveCobraLP(LPproblem,true);
    maxFlux(i) = LPsolution.obj;
    LPproblem.osense = 1;
    LPsolution = solveCobraLP(LPproblem,true);
    minFlux(i) = LPsolution.obj;
    if (verbFlag == 1)
        waitbar(i/length(rxnNameList),h);
    end
    if (verbFlag > 1)
        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));    
    end 
end
if (verbFlag == 1)
    close(h);
end

maxFlux = columnVector(maxFlux);
minFlux = columnVector(minFlux);