SysBioChalmers / GECKO

Toolbox for including enzyme constraints on a genome-scale model.
http://sysbiochalmers.github.io/GECKO/
MIT License
66 stars 51 forks source link

Flux comparison using GLPK & CPLEX #332

Closed Yu-sysbio closed 1 year ago

Yu-sysbio commented 1 year ago

I used GLPK and CPLEX (via COBRA) to run Step 29 but got quite different plots (see attached plots.zip). The deviation comes from the comparison between extremely low fluxes. While the plots look different, the general conclusion still holds: the simulated objective fluxes are the same, and there is no big difference between mapped fluxes simulated by light and full ecModels.

I have an issue with the license of GUROBI currently and thus do not have the corresponding result for comparison. But I suppose that the plot in the manuscript was generated using GUROBI, which looks much better than the other two solvers. Maybe we could mention in the manuscript that all the anticipated results are obtained with GUROBI (I also got a slightly different FVA plot using GLPK). And maybe only plot high absolute fluxes (e.g. > 1e-6) as the extremely low fluxes make no sense.

Below I also put what I got in the Command Window using these two solvers:

GLPK:

Comparison of duration light vs. full ecModel ecModel reconstruction: 91% (11 vs 12 seconds) WARNING: Could not solve the problem of minimizing the sum of fluxes. Uses output from original problem WARNING: Could not solve the problem of minimizing the sum of fluxes. Uses output from original problem FBA: 90% (1.34 vs 1.50 seconds) Mapping fluxes: 60% (0.037 vs 0.063 seconds) Growth rate that is reached: 0.0252 vs 0.0252

CPLEX:

Comparison of duration light vs. full ecModel ecModel reconstruction: 78% (11 vs 14 seconds) FBA: 68% (0.61 vs 0.90 seconds) Mapping fluxes: 39% (0.038 vs 0.099 seconds) Growth rate that is reached: 0.0252 vs 0.0252

I hereby confirm that I have:

edkerk commented 1 year ago

It might be good to double-check which reactions have those very low fluxes.

But it indeed seems to be an artifact of the solver and not the true difference of the flux distributions. I can modify plotlightVSfull to set these fluxes that are below the solver cut-off to 0.

In the Materials section is mentioned that Gurobi is strongly recommended.

edkerk commented 1 year ago

Please review #331 to see if the plots look more similar.

Yu-sysbio commented 1 year ago

It now looks similar using GLPK. But using CPLEX via COBRA got the error:

The logical indices in position 1 contain a true value outside of the array bounds. Error in mapRxnsToConv (line 29) fluxes(revRxns,:) = -fluxes(revRxns,:); Error in plotlightVSfull (line 51) tic; solF = mapRxnsToConv(ecFull,model,solFull.x);

I am not sure if this should be fixed as CPLEX is not mentioned in the Materials section. But this means that there must be some bug in the code.

edkerk commented 1 year ago

I imagine this would happen if CPLEX failed to find a solution (solFull.x is empty), not sure why this would happen though. I'll include a check in mapRxnsToConv that throws an error if provided with empty flux vector.