pFBA and fva : infinite ?

[Ludo Cottret]

Hi,

Two weeks ago, I've launched pFBA and fva on the iAF1260 E.coli model
and it's still running... I'm quite surprised because fva with OptFlux
takes only a couple of hours.
Do you think that GPLK is again the source of the problem here ? Does
it exist a free alternative solver that works better with the Cobra
toolbox or do I need to buy another solver to be able to run the
Cobra toolbox efficiently ?

Ludo

[Ronan Fleming]

Try Gurobi
http://www.gurobi.com/

--
Mr. Ronan MT Fleming B.V.M.S. Dip. Math. Ph.D.
-----------------------------------------------------------------
Research Scientist,
Science Institute & Center for Systems Biology,
University of Iceland.
http://systemsbiology.hi.is/
http://www.hi.is/~rfleming
Ph:  +354 618 6245
-----------------------------------------------------------------

[Ludo Cottret]

Is the Cobra toolbox directly compatible with Gurobi ?

2011/1/6 Ronan Fleming <ronan.mt...@gmail.com>:

[Ben]

Hi Ludo,

I've found that it is easy to use Gurobi with the Cobra toolbox, once
you have all the pieces together. I've also found that Gurobi is much
faster than GLPK for OptKnock - I suspect it would be for OptFlux,
too.

There are a few steps to get Gurobi working with matlab:
- instructions for Gurobi are here: http://www.gurobi.com/html/academic.html
- instructions for interfacing Gurobi and MATLAB are here:
http://www.convexoptimization.com/wikimization/index.php/Gurobi_Mex:_A_MATLAB_interface_for_Gurobi
- and a free compiler (microsoft visual C++) is available here:
http://www.microsoft.com/express/Downloads/)

Once you have all the pieces, you compile a gurobi mex file in matlab
with a command like:

mex -O -I"C:\Gurobi400\win32\include" "C:\Documents and Settings
\bdh32\Desktop\gurobi_mex_v1.45\gurobi_mex_v1.45\gurobi_mex.c" "C:
\Gurobi400\win32\lib\gurobi40.lib" "C:\Program Files\MATLAB\R2010a
\extern\lib\win32\microsoft\libut.lib"

(you'll have to adjust the paths to get it to work for you).

Gurobi includes some sample code so you can test it out before trying
it with the COBRA toolbox.

Then, you can use the COBRA toolbox command

changeCobraSolver('gurobi', 'MILP')

and life should be good.

I've found it useful to edit my initcobratoolbox script to move gurobi
up so that it is initialized first for the MILP solver.

I hope that helps!
-Ben

On Jan 7, 5:00 am, Ludo Cottret <l.cott...@gmail.com> wrote:
> Is the Cobra toolbox directly compatible with Gurobi ?
>

> 2011/1/6 Ronan Fleming <ronan.mt.flem...@gmail.com>:
>
> > Try Gurobi
> >http://www.gurobi.com/

[dr scientist]

Ludo,

It is highly possible that the culprit is the COBRA Toolboxes
interface to GLPK. I use GLPK extensively in my non-MATLAB research
and it can handle the calculations. One of the goals of The openCOBRA
project is to find and fix these type of bugs in the core COBRA
Toolbox code. As the community grows and we acquire dedicated
resources for the project, we'll be able to squash bugs more rapidly.

If you could submit this as a bug (http://sourceforge.net/projects/
opencobra/support) to the bug tracker along with the steps and files
to reproduce then we can see about fixing this.

Better yet, if you figure out where the problem is you can submit a
patch (http://sourceforge.net/projects/opencobra/support) and as long
as it doesn't break any of the exisiting features it'll make its way
into the code.


As an alternative, Gurobi which was mentioned by Ronan and Ben is an
excellent solver free to academics and its interface in the COBRA
Toolbox has received more love than GLPK so it might do the job.
Compiling the mex on GNU/Linux and Mac OS X is trivial and Ben's
instructions for MS Windows look just as simple.

On Jan 7, 5:55 am, Ben <bheav...@gmail.com> wrote:
> Hi Ludo,
>
> I've found that it is easy to use Gurobi with the Cobra toolbox, once
> you have all the pieces together. I've also found that Gurobi is much
> faster than GLPK for OptKnock - I suspect it would be for OptFlux,
> too.
>
> There are a few steps to get Gurobi working with matlab:
> - instructions for Gurobi are here:http://www.gurobi.com/html/academic.html

> - instructions for interfacing Gurobi and MATLAB are here:http://www.convexoptimization.com/wikimization/index.php/Gurobi_Mex:_...

[SaraS]

I tried gurobi but it still says inconsistent objective values and does not stop running. I would really appreciate your thoughts

[Ronan M.T. Fleming]

Hi Sara,

1) Please make sure you are using the most recent version

git clone https://github.com/opencobra/cobratoolbox.git

https://git-scm.com/book/en/v2/Getting-Started-Git-Basics

2) Provide more detail in your questions and people will more probably respond. When you take one line to write an email, potential responders might look at it an think, well if the question is not important enough to the person asking for them to write it out in full, then why should I take the time to respond?

Regards,

Ronan

--

---
You received this message because you are subscribed to the Google Groups "COBRA Toolbox" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cobra-toolbo...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

--
Mr. Ronan MT Fleming B.V.M.S. Dip. Math. Ph.D.

----------------------------------------------------------------------------
Senior research associate (EN) == Chercheur (FR),
Principal investigator,
Systems Biochemistry Group,
wwwen.uni.lu/lcsb/research/systems_biochemistry
Luxembourg Centre for Systems Biomedicine,
University of Luxembourg,
Campus Belval,
7, avenue des Hauts-Fourneaux
Esch-sur-Alzette,
Luxembourg,
L-4362.
&
National Centre of Excellence in Research on Parkinson’s disease
www.parkinson.lu

&

Adjunct Assistant Professor,

Division of Analytical Biosciences,

Leiden Academic Centre for Drug Research,

Faculty of Science,

University of Leiden.

http://analyticalbiosciences.leidenuniv.nl
----------------------------------------------------------------------------
Mobile:  +352 621 175 112
Office: +352 466 644 5528
Skype: ronan.fleming
----------------------------------------------------------------------------
(This message is confidential and may contain privileged information. It is intended for the named recipient only. If you receive it in error please notify me and permanently delete the original message and any copies.)

[Thierry Mondeel]

Hi Ludo and Sara,

This might be a cop-out answer but I usually do not want to perform FBA on all possible reactions in a large model. The fluxVariability function allows you to only calculate the FVA ranges on a select number of reactions through the rxnNameList argument. If this number is low enough it massively reduces computational time. For instance, just run FVA on a subsystem of interest or on all exchange reactions.

Thierry

[SaraS]

Hello,

I apologize for my short email. I had a similar problem as Ludo and was not sure if I should reiterate the issues. I tried using pFBA with three different solvers (ibm_cplex, gurobi, and glpk). I keep getting the error:

reduceModel.m: Inconsistent objective values 135.617 0 0 0

When I terminate the code, I get the following error:

Operation terminated by user during ismember>ismemberBuiltinTypes (line 314)

In ismember>ismemberR2012a (line 142)

            lia = ismemberBuiltinTypes(a,b);

In ismember (line 56)

    [varargout{1:max(1,nargout)}] = ismemberR2012a(A,B);

In verifyCobraProblem (line 106)

invalidCsense = find(~ismember(XPproblem.csense,['E','L','G']));

In optimizeCbModel (line 150)

if ~(verifyCobraProblem(LPproblem, [], [], false) == 1)

In reduceModel>checkConsistency (line 257)

    solRedMin = optimizeCbModel(modelRed,'min');

In reduceModel>expandBounds (line 240)

    modelOK = checkConsistency(model,tempModel,tol);

In reduceModel (line 211)

            modelRed = expandBounds(model,modelRed,tol);

In pFBA (line 89)

tempmodel = reduceModel(tempmodel,tol); % reduce the model to find blocked

reactions

-----------

I printed out the values of diffMax and diffMin in reduceModel.m

For diffMax, I get 135 which is much higher than the tolerance (1e-6).

Would performing pFBA on a select number of reactions possibly reduce the reliability of the answers I obtain?

Many thanks!

Sara

[Nathan Lewis]

Many of these errors come from some problem with reduceModel.m, which has been posing a problem in the recent past. I haven't taken the time to figure out its problems, so I don't recommend using that function until it's fixed. Thus, I have a different version of the pFBA function I typically use that doesn't use reduceModel.m

Here's the modified implementation of pFBA:

function [GeneClasses RxnClasses modelIrrevFM] = ModifiedpFBA(model, varargin)
% Parsimoneous enzyme usage Flux Balance Analysis - method that optimizes
% the user's objective function and then minimizes the flux through the
% model and subsequently classifies each gene by how it contributes to the
% optimal solution. See Lewis, et al. Mol Syst Bio doi:10.1038/msb.2010.47
%
% Inputs:
% model                 COBRA model
%
% varargin:
% 'geneoption'          1 = only minimize the sum of the flux through
%                       gene-associated fluxes (default), 0 = minimize the
%                       sum of all fluxes in the network
% 'tol'                 tolerance (default = 1e-6)
% 'map'                 map structure from readCbMap.m (no map written if empty
% 'mapoutname'          File Name for map
% 'skipclass'           1 = Don't classify genes and reactions. Only return
%                       model with the minimium flux set as upper bound.
%                       0 = classify genes and reactions (default).
%
% Output:
% GeneClasses           Structure with fields for each gene class
% RxnsClasses           Structure with fields for each reaction class
% modelIrrevFM          Irreversible model used for minimizing flux with
%                       the minimum flux set as a flux upper bound
%
%
%
% ** note on maps: Red (6) = Essential reactions, Orange (5) = pFBA optima
%       reaction, Yellow (4) = ELE reactions, Green (3) = MLE reactions,
%       blue (2) = zero flux reactions, purple (1) = blocked reactions,
%       black (0) = not classified
%
% Example:
% [GeneClasses RxnClasses modelIrrevFM] = pFBA(model, 'geneoption',0, 'tol',1e-7)
%
% by Nathan Lewis Aug 25, 2010

%

if nargin < 2
    tol = 1e-6;
    GeneOption = 1;
    map = []; % no map
    mapOutName = 'pFBA_map.svg';
    skipclass = 0;
end

if mod(length(varargin),2)==0
    for i=1:2:length(varargin)-1
        switch lower(varargin{i})
            case 'geneoption', GeneOption = varargin{i+1};
            case 'tol', tol = varargin{i+1};
            case 'map', map = varargin{i+1};
            case 'mapoutname', mapOutName = varargin{i+1};
            case 'skipclass', skipclass = varargin{i+1};
            otherwise, options.(varargin{i}) = varargin{i+1};
        end
    end
else
    error('Invalid number of parameters/values');
end

if ~exist('GeneOption','var'), GeneOption = 1;              end
if ~exist('tol','var'), tol = 1e-6;                         end
if ~exist('map','var'), map = [];                           end
if ~exist('mapOutName','var'), mapOutName = 'pFBA_map.svg'; end
if ~exist('skipclass','var'), skipclass = 0;                end

if skipclass % skip the model reduction and gene/rxn classification
    % minimize the network flux
FBAsoln = optimizeCbModel(model);
model.lb(model.c==1) = FBAsoln.f;
[ MinimizedFlux modelIrrevFM]= minimizeModelFlux_local(model,GeneOption);
modelIrrevFM = changeRxnBounds(modelIrrevFM,'netFlux',MinimizedFlux.f,'b');
GeneClasses = [];
RxnClasses = [];
else

% save a copy of the inputted model
model_sav = model;

% Remove all blocked reactions
[selExc,selUpt] = findExcRxns(model,0,0); % find and open up all exchanges
tempmodel = changeRxnBounds(model,model.rxns(selExc),-1000,'l');
tempmodel = changeRxnBounds(tempmodel,model.rxns(selExc),1000,'u');
[maxF minF] = fluxVariability(tempmodel,.001);
Blocked_Rxns = model.rxns(and(abs(maxF)<tol,abs(minF)<tol));

% tempmodel = reduceModel(tempmodel,tol); % reduce the model to find blocked reactions
% Blocked_Rxns_old = setdiff(model.rxns,regexprep(tempmodel.rxns,'_r$',''));
% model_old = removeRxns(model,setdiff(model.rxns,regexprep(tempmodel.rxns,'_r$',''))); % remove blocked reactions
% Ind2Remove_old  = find(~and(sum(full(model.rxnGeneMat),1),1));
% Blocked_genes_old  = model.genes(Ind2Remove_old );
% % Ind2Remove = find(~and(sum(full(model_sav.rxnGeneMat),1),1));
% % Blocked_genes = model_sav.genes(Ind2Remove);

model = removeRxns(model,Blocked_Rxns); % remove blocked reactions
Ind2Remove = find(~and(sum(full(model.rxnGeneMat),1),1));
Blocked_genes = model.genes(Ind2Remove);
model.genes(Ind2Remove)={'dead_end'}; % make sure genes that are unique to blocked reactions are tagged for removal

% find essential genes
grRatio = singleGeneDeletion(model);
grRatio(isnan(grRatio))=0;
pFBAEssential = model.genes(grRatio<tol);

if nargout > 1
    % find essential reactions
RxnRatio = singleRxnDeletion(model);
RxnRatio(isnan(RxnRatio))=0;
pFBAEssentialRxns = model.rxns(RxnRatio<tol);
end
   

% remove zero flux rxns
[maxF,minF] = fluxVariability(model,.001);
% tempmodel = reduceModel(model,tol);
% ZeroFluxRxns = setdiff(model.rxns,regexprep(tempmodel.rxns,'_r$',''));
ZeroFluxRxns = model.rxns(and(abs(maxF)<tol,abs(minF)<tol));
model = removeRxns(model,ZeroFluxRxns);

% find MLE reactions
FBAsoln = optimizeCbModel(model);
model.lb(model.c==1) = FBAsoln.f;
[minFlux,maxFlux] = fluxVariability(model,100);
for i = 1:length(minFlux)
    tmp(i,1) = max([abs(minFlux(i)) abs(maxFlux(i))])<tol;
end
MLE_Rxns = setdiff(model.rxns(tmp),ZeroFluxRxns);

% minimize the network flux
[ MinimizedFlux,modelIrrevFM]= minimizeModelFlux_local(model,GeneOption);

% separate pFBA optima rxns from ELE rxns
modelIrrevFM = changeRxnBounds(modelIrrevFM,'netFlux',MinimizedFlux.f,'b');
[minFlux,maxFlux] = fluxVariability(modelIrrevFM,100);
pFBAopt_Rxns = modelIrrevFM.rxns((abs(minFlux)+abs(maxFlux))>=tol);
ELE_Rxns = modelIrrevFM.rxns((abs(minFlux)+abs(maxFlux))<=tol);
pFBAopt_Rxns = unique(regexprep(pFBAopt_Rxns,'_[f|b]$',''));
pFBAopt_Rxns(ismember(pFBAopt_Rxns,MLE_Rxns))=[]; %%% removes non-gene associated reversible reactions that are only in the list because there is no constraint on them looping with the reverse reaction
ELE_Rxns = unique(regexprep(ELE_Rxns,'_[f|b]$',''));
ELE_Rxns = ELE_Rxns(~ismember(ELE_Rxns,pFBAopt_Rxns));
ELE_Rxns = ELE_Rxns(~ismember(ELE_Rxns,MLE_Rxns));

% determine pFBA optima genes
pFBAopt_Rxns(ismember(pFBAopt_Rxns,'netFlux'))=[];
[geneList]=findGenesFromRxns(model,pFBAopt_Rxns);
geneList2 = {};
for i = 1:length(geneList),
    geneList2(end+1:end+length(geneList{i}),1) = columnVector( geneList{i});
end
pFBAOptima = unique(geneList2);

% determine Zero Flux genes
Ind2Remove = find(~and(sum(full(model.rxnGeneMat),1),1));
ZeroFluxGenes = unique(model.genes(Ind2Remove));

% determine ELE genes
[geneList]=findGenesFromRxns(model,ELE_Rxns);
geneList2 = {};
for i = 1:length(geneList)
    geneList2(end+1:end+length(geneList{i}),1) = columnVector( geneList{i});
end
ELEGenes = unique(geneList2);
ELEGenes = setdiff(ELEGenes,[pFBAOptima;ZeroFluxGenes]);

% determine Met ineff genes
MLEGenes = setdiff(model.genes, [pFBAOptima;ZeroFluxGenes;ELEGenes]);

% clean up lists by removing non-genes
pFBAOptima(~cellfun('isempty',regexp(pFBAOptima,'dead_end')))=[];
ELEGenes(~cellfun('isempty',regexp(ELEGenes ,'dead_end')))=[];
MLEGenes(~cellfun('isempty',regexp(MLEGenes,'dead_end')))=[];
ZeroFluxGenes(~cellfun('isempty',regexp(ZeroFluxGenes,'dead_end')))=[];

pFBAOptima(cellfun('isempty',pFBAOptima))=[];
ELEGenes(cellfun('isempty',ELEGenes ))=[];
MLEGenes(cellfun('isempty',MLEGenes))=[];
ZeroFluxGenes(cellfun('isempty',ZeroFluxGenes))=[];

% filter out essential genes from pFBA optima
pFBAOptima(ismember(pFBAOptima,pFBAEssential))=[];

if nargout > 1
% filter out essential Rxns from pFBA optima
pFBAopt_Rxns(ismember(pFBAopt_Rxns,pFBAEssentialRxns))=[];
end

% prepare output variables
GeneClasses.pFBAEssential =pFBAEssential;
GeneClasses.pFBAoptima = pFBAOptima;
GeneClasses.ELEGenes = ELEGenes;
GeneClasses.MLEGenes = MLEGenes;
GeneClasses.ZeroFluxGenes = ZeroFluxGenes;
GeneClasses.Blockedgenes = Blocked_genes;

RxnClasses.Essential_Rxns = pFBAEssentialRxns;
RxnClasses.pFBAOpt_Rxns = pFBAopt_Rxns;
RxnClasses.ELE_Rxns = ELE_Rxns;
RxnClasses.MLE_Rxns = MLE_Rxns;
RxnClasses.ZeroFlux_Rxns = ZeroFluxRxns;
RxnClasses.Blocked_Rxns = Blocked_Rxns;

if ~isempty(map)
    MapVector = zeros(length(model_sav.rxns),1);
    MapVector(ismember(model_sav.rxns,Blocked_Rxns))= 1;
    MapVector(ismember(model_sav.rxns,ZeroFluxRxns))= 2;
    MapVector(ismember(model_sav.rxns,MLE_Rxns))= 3;
    MapVector(ismember(model_sav.rxns,ELE_Rxns))= 4;
    MapVector(ismember(model_sav.rxns,pFBAopt_Rxns))= 5;
    MapVector(ismember(model_sav.rxns,pFBAEssentialRxns))= 6;
   
    options.lb = 0;
    options.ub = 6;
    tmpCmap = hsv(18);
    tmpCmap = [tmpCmap([1,3,4,6,11,14],:); 0 0 0;];
    options.fileName = mapOutName;
    options.colorScale = flipud(round(tmpCmap*255));
   
    global CB_MAP_OUTPUT
    if ~exist('CB_MAP_OUTPUT', 'var') || isempty(CB_MAP_OUTPUT)
        changeCbMapOutput('svg');
    end
    drawFlux(map, model_sav, MapVector, options);
end
end
end

function [ MinimizedFlux modelIrrev]= minimizeModelFlux_local(model,GeneOption)
% This function finds the minimum flux through the network and returns the
% minimized flux and an irreversible model

% convert model to irrev
modelIrrev = convertToIrreversible(model);

% add pseudo-metabolite to measure flux through network
if nargin==1,GeneOption=0;
end
if GeneOption==0, % signal that you want to minimize the sum of all gene and non-gene associated fluxes
    modelIrrev.S(end+1,:) = ones(size(modelIrrev.S(1,:)));
elseif GeneOption==1, % signal that you want to minimize the sum of only gene-associated fluxes
    %find all reactions which are gene associated
    Ind=find(sum(modelIrrev.rxnGeneMat,2)>0);
    modelIrrev.S(end+1,:) = zeros(size(modelIrrev.S(1,:)));
    modelIrrev.S(end,Ind) = 1;
end
modelIrrev.b(end+1) = 0;
modelIrrev.mets{end+1} = 'fluxMeasure';

% add a pseudo reaction that measures the flux through the network
modelIrrev = addReaction(modelIrrev,'netFlux',{'fluxMeasure'},[-1],false,0,inf,0,'','');

% set the flux measuring demand as the objective
modelIrrev.c = zeros(length(modelIrrev.rxns),1);
modelIrrev = changeObjective(modelIrrev, 'netFlux');

% minimize the flux measuring demand (netFlux)
MinimizedFlux = optimizeCbModel(modelIrrev,'min');
end

[André Schultz]

I was getting a similar error to Sara's (reduceModel.m: Inconsistent objective values 135.617 0 0 0), and it did stem from reduceModel.m.
For me, it was because model.lb did not match model.rev. While some of the reactions had a lower bound lower than zero and were able to carry a negative flux, they were marked as irreversible and their minimum flux (and thus lower bound) was set to zero. This then led to the above error when checking model consistency.
I don't know if this is what was happening to you but I just thought I would mention.
André

[Ronan M.T. Fleming]

Hi Andre,

we are still looking for a volunteer to go clone the cobra toolbox, search for all the code that uses model.rev and replace it with code that reads directionality directly from model.lb and model.ub.

Regards,

Ronan