Inconsistency in finding blocked reactions

[hvdinh16]

For the function cobra.flux_analysis.find_blocked_reactions:

If the function is implemented by the scripts below, the number of blocked reactions will be miscalculated, only specific for GLPK solver and not for Gurobi or CPLEX: import cobra model = cobra.io.read_sbml_model('./iRhtoC.xml') model.solver = 'glpk' result = cobra.flux_analysis.find_blocked_reactions(model, reaction_list=model.reactions, open_exchanges=True, zero_cutoff=1e-6)

If the function is implemented by the second method below, the number of blocked reactions is correct for GLPK solver: from cobra.flux_analysis import find_blocked_reactions from cobra.io import read_sbml_model config = cobra.Configuration() config.solver = 'glpk' model = read_sbml_model('iRhtoC.xml') result = find_blocked_reactions(model, open_exchanges=True, zero_cutoff=model.tolerance)

Scripts and model are attached below: blocked_reactions.ipynb.tar.gz iRhtoC.xml.zip

This issue propagates to memote, reported in https://github.com/opencobra/memote/issues/639

[cdiener]

Could you follow the issue template and provide more info? What exactly does "wrong" mean? Do you get more blocked reactions in the first case? You are setting a zero_cutoff that is an order of magnitude larger than the solver tolerance. So you will call reactions blocked even if they have evidence for non-zero flux.

However, there is a potential caveat here. Basically if your zero tolerance is stricter than your solver tolerance you get nonsense. We should throw a warning in that case.

[Midnighter]

However, there is a potential caveat here. Basically if your zero tolerance is stricter than your solver tolerance you get nonsense. We should throw a warning in that case.

That should actually already be handled adequately by the normalize_cutoff function. https://github.com/opencobra/cobrapy/blob/4808655418bc88f98c63eaaf76302cb13a93b818/cobra/flux_analysis/variability.py#L242

[hvdinh16]

Hi,

I have run the analysis again and the comparisons in the attached ipython notebook ( blocked_reactions-2019-06-27.zip) and shown below.

When switch to glpk, a warning of "solver status is 'infeasible'" is raised. find_blocked_reactions's results depend on the two steps of model.slim_optimize() and flux_variability_analysis. Because FVA results returned by glpk and cplex are comparable, we can infer there might be something wrong with model.slim_optimize() step.

For the actual results using the model attached in the first comment, glpk detected 48 more blocked reactions. These 48 reactions can carry flux in FVA (using flux_variability_analysis) even with glpk as the solver.

import cobra
from cobra.flux_analysis import find_blocked_reactions
from cobra.io import read_sbml_model

config = cobra.Configuration()
config.tolerance

1e-07

import cplex, optlang
print 'cobra:', cobra.__version__
print 'cplex:', cplex.__version__
print 'optlang:', optlang.__version__

cobra: 0.15.1
cplex: 12.8.0.0
optlang: 1.4.4

GLPK

model = read_sbml_model('./iRhtoC.xml')
model.solver = 'glpk'
model.tolerance

1e-07

brxns_glpk = find_blocked_reactions(model, open_exchanges=True, reaction_list=model.reactions,
                                    zero_cutoff=1e-7)
print len(brxns_glpk)

cobra/util/solver.py:416 UserWarning: solver status is 'infeasible'

725

CPLEX

model = read_sbml_model('./iRhtoC.xml')
model.solver = 'cplex'
model.tolerance

1e-07

brxns_cplex = find_blocked_reactions(model, open_exchanges=True, reaction_list=model.reactions,
                                zero_cutoff=1e-7)
print len(brxns_cplex)

677

Manual check using FVA with CPLEX

Without precheck using FBA

model = read_sbml_model('./iRhtoC.xml')
model.solver = 'cplex'
for rxn in model.exchanges:
    rxn.bounds = (min(rxn.lower_bound, -1000), max(rxn.upper_bound, 1000))
fva_cplex = cobra.flux_analysis.flux_variability_analysis(model,
            reaction_list=model.reactions, fraction_of_optimum=0)

brxns_manual = []
for rxn in model.reactions:
    if abs(fva_cplex.minimum[rxn.id]) < 1e-7 and abs(fva_cplex.maximum[rxn.id]) < 1e-7:
        brxns_manual.append(rxn.id)

len(brxns_manual)

677

print 'Difference vs. find_blocked_reactions, GLPK:', len(set(brxns_glpk) - set(brxns_manual)),
print len(set(brxns_manual) - set(brxns_glpk))
print 'Difference vs. find_blocked_reactions, CPLEX:', len(set(brxns_cplex) - set(brxns_manual)),
print len(set(brxns_manual) - set(brxns_cplex))

Difference vs. find_blocked_reactions, GLPK: 48 0
Difference vs. find_blocked_reactions, CPLEX: 0 0

Compare: GLPK and CPLEX

print 'No. of rxns only obtained using GLPK:', len(set(brxns_glpk) - set(brxns_cplex))
print 'No. of rxns only obtained using CPLEX:', len(set(brxns_cplex) - set(brxns_glpk))

No. of rxns only obtained using GLPK: 48
No. of rxns only obtained using CPLEX: 0

model = read_sbml_model('./iRhtoC.xml')
model.solver = 'cplex'
for rxn in model.exchanges:
    rxn.bounds = (min(rxn.lower_bound, -1000), max(rxn.upper_bound, 1000))
fva_cplex = cobra.flux_analysis.flux_variability_analysis(model,
            reaction_list=list(set(brxns_glpk) - set(brxns_cplex)), fraction_of_optimum=0)

model = read_sbml_model('./iRhtoC.xml')
model.solver = 'glpk'
for rxn in model.exchanges:
    rxn.bounds = (min(rxn.lower_bound, -1000), max(rxn.upper_bound, 1000))
fva_glpk = cobra.flux_analysis.flux_variability_analysis(model,
            reaction_list=list(set(brxns_glpk) - set(brxns_cplex)), fraction_of_optimum=0)

for rxnid in set(brxns_glpk) - set(brxns_cplex):
    print rxnid, (fva_glpk.minimum[rxnid], fva_glpk.maximum[rxnid]), \
        (fva_cplex.minimum[rxnid], fva_cplex.maximum[rxnid])

SPMS_c (0.0, 1000.0) (0.0, 1000.0)
FACOAL183_x (0.0, 1000.0000000000002) (0.0, 1000.0)
ACOADS181_rm (0.0, 512.3656263407952) (0.0, 512.365626340795)
ACOAO141b_m (0.0, 333.3333333333333) (0.0, 333.3333333333333)
FAOX182lump_x (0.0, 111.1111111111113) (0.0, 111.1111111111111)
ACOAO100_m (0.0, 400.00000000000006) (0.0, 400.0)
GLUN_e (0.0, 1000.0) (0.0, 1000.0)
ACOAO160_m (0.0, 250.00000000000014) (0.0, 250.0)
CRN181t_m_x (0.0, 250.0000000000001) (0.0, 250.0)
ME2_c (0.0, 1000.0) (0.0, 1000.0)
ECOAI100a_m (0.0, 333.3333333333333) (0.0, 333.33333333333337)
FACOAE220_c (0.0, 499.9999999999998) (0.0, 500.0)
PAIL4Pt_en_rm (0.0, 1000.0000000000001) (0.0, 1000.0)
FACOA182tabc_x (0.0, 111.11111111111114) (0.0, 111.1111111111111)
FACOAL182_x (0.0, 1000.0) (0.0, 1000.0)
ACOAO161b_m (0.0, 333.3333333333333) (0.0, 333.33333333333337)
DM_5mta_c (0.0, 999.9999999999889) (0.0, 1000.0)
PEt_gm_rm (0.0, 333.38502439716325) (0.0, 333.385024397163)
CRNAT161_m (0.0, 333.33333333333314) (0.0, 333.3333333333333)
EX_caro_c (0.0, 114.54114999999999) (0.0, 114.54114999999999)
ASNN_e (0.0, 1000.0) (0.0, 1000.0)
FACOAE240_c (0.0, 1000.0) (0.0, 1000.0)
CRN161t_m_x (0.0, 333.3333333333333) (0.0, 333.33333333333326)
CRN160t_m_x (0.0, 249.99999999790646) (0.0, 250.0)
ACOAO161a_m (0.0, 250.00000000000006) (0.0, 250.0)
FACOAE240_x (0.0, 1000.0000000000005) (0.0, 1000.0)
TAGt_c_m (0.0, 333.3333333333341) (0.0, 333.33333333333337)
GLYCt_c_m (0.0, 333.3333333333326) (0.0, 333.3333333333333)
CRNAT181_m (0.0, 250.00000000000006) (0.0, 250.0)
ECOAI120a_m (0.0, 250.00000000000006) (0.0, 250.0)
SPRMS_c (0.0, 1000.0) (0.0, 1000.0)
FACOAE182_c (0.0, 1000.0) (0.0, 1000.0)
PAILt_rm_vm (0.0, 999.9999999999922) (0.0, 1000.0)
CRN140t_m_x (0.0, 333.3333333333333) (0.0, 333.3333333333333)
ACOAO120_m (0.0, 399.9999999999982) (0.0, 400.0)
FACOAE183_x (0.0, 1000.000000000004) (0.0, 1000.0)
ACOAO180_m (0.0, 200.00000000000003) (0.0, 200.0)
ACOAO181a_m (0.0, 249.99999999999997) (0.0, 250.0)
FACOAE182_x (0.0, 1000.0) (0.0, 1000.0)
MAGL_m (0.0, 333.3333333333333) (0.0, 333.3333333333333)
FACOAE183_c (0.0, 1000.0) (0.0, 1000.0)
ACOAO120a_m (0.0, 333.33333333333314) (0.0, 333.3333333333333)
CRN180t_m_x (0.0, 200.00000000000003) (0.0, 200.0)
ACOAO140_m (0.0, 333.3333333333331) (0.0, 333.3333333333333)
ACOAO60_m (0.0, 1000.0) (0.0, 1000.0)
ACOAO80_m (0.0, 1000.0) (0.0, 1000.0)
ACOAO141a_m (0.0, 249.99999999999994) (0.0, 250.0)
CRN120t_m_x (0.0, 399.99999999999983) (0.0, 400.0)

[cdiener]

Hi, thanks for the additional info. Unfortunately I can't reproduce your error even if I use the model you uploaded:

In [1]: import cobra

In [2]: mod = cobra.io.read_sbml_model("/Users/cdiener/Downloads/iRhtoC.xml")

In [3]: mod.tolerance
Out[3]: 1e-07

In [4]: mod.solver
Out[4]: <optlang.cplex_interface.Model at 0x12275e6d8>

In [5]: blocked = cobra.flux_analysis.find_blocked_reactions(
   ...:     mod, open_exchanges=True, reaction_list=mod.reactions,
   ...:     zero_cutoff=1e-7)

In [6]: print(len(blocked))
677

In [8]: mod.solver = "glpk"

In [9]: mod.tolerance
Out[9]: 1e-07

In [10]: blocked = cobra.flux_analysis.find_blocked_reactions(
    ...:     mod, open_exchanges=True, reaction_list=mod.reactions,
    ...:     zero_cutoff=1e-7)

In [11]: print(len(blocked))
677

The problem is definitely in the infeasible warning you get (and I don't). This should not happen and will screw up find_blocked_reactions so the differences don't surprise me. We should add an error message here.

Could you try to update cobrapy and all of its dependencies, especially swiglpk and libsbml? How did you install cobrapy and on what OS are you running it?

[hvdinh16]

Hi. The updates are for following packages: cobra (0.15.1 -> 0.15.3), numpy (1.16.2 -> 1.16.4), python-libsbml-experimental (5.17.2 -> 5.18.0), and ruamel.yaml (0.15.89 -> 0.15.97). The issue persists.

cobra was downloaded by cloning devel branch from github and installed by "python setup.py develop". I ran on ubuntu 18.04.2 (few months ago, it had been 16.04). All python packages are in a virtualenv.