opencobra / cobrapy

COBRApy is a package for constraint-based modeling of metabolic networks.
http://opencobra.github.io/cobrapy/
GNU General Public License v2.0
463 stars 217 forks source link

Reaction.summary() broken for some models #1118

Open cdiener opened 3 years ago

cdiener commented 3 years ago

Problem description

Reconstructed a model for a course with CARVEME and was surprised that Reaction.summary() errors with a key error.

Code Sample

Create a minimal, complete, verifiable example.

In [1]: from cobra.io import read_sbml_model

In [2]: model = read_sbml_model("SRR9275468.asm1.xml")

In [3]: model.optimize()
Out[3]: <Solution 95.558 at 0x7f51eb05da90>

In [4]: print(model.reactions[0].summary())
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-4-b09348ecbd9e> in <module>
----> 1 print(model.reactions[0].summary())

~/code/cobrapy/src/cobra/summary/summary.py in __str__(self)
     73     def __str__(self) -> str:
     74         """Return a string representation of the summary."""
---> 75         return self.to_string()
     76 
     77     def _repr_html_(self) -> str:

~/code/cobrapy/src/cobra/summary/reaction_summary.py in to_string(self, names, threshold, float_format, column_width)
    195             header = shorten(self._reaction.id, width=column_width, placeholder="...")
    196 
--> 197         flux = self._string_flux(threshold, float_format)
    198 
    199         return dedent(

~/code/cobrapy/src/cobra/summary/reaction_summary.py in _string_flux(self, threshold, float_format)
    157         else:
    158             frame = self._flux.loc[self._flux["flux"].abs() >= threshold, :].copy()
--> 159             return f"{frame.at[self._reaction.id, 'flux']:{float_format}}"
    160 
    161     def to_string(

~/.local/lib/python3.9/site-packages/pandas/core/indexing.py in __getitem__(self, key)
   2273             return self.obj.loc[key]
   2274 
-> 2275         return super().__getitem__(key)
   2276 
   2277     def __setitem__(self, key, value):

~/.local/lib/python3.9/site-packages/pandas/core/indexing.py in __getitem__(self, key)
   2220 
   2221         key = self._convert_key(key)
-> 2222         return self.obj._get_value(*key, takeable=self._takeable)
   2223 
   2224     def __setitem__(self, key, value):

~/.local/lib/python3.9/site-packages/pandas/core/frame.py in _get_value(self, index, col, takeable)
   3577 
   3578         try:
-> 3579             loc = engine.get_loc(index)
   3580             return series._values[loc]
   3581         except AttributeError:

~/.local/lib/python3.9/site-packages/pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

~/.local/lib/python3.9/site-packages/pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

KeyError: '12DGR120tipp'

Here is the model that causes the crash: SRR9275468.asm1.xml.gz

Context

``` System Information ================== OS Linux OS-release 5.14.0-1-amd64 Python 3.9.7 Package Versions ================ appdirs~=1.4 not installed black 21.9b0 bumpversion 0.6.0 cobra 0.22.0 depinfo 1.7.0 diskcache~=5.0 not installed future 0.18.2 httpx~=0.14 not installed importlib_resources 3.3.0 isort; extra == "development" not installed numpy~=1.13 not installed optlang~=1.5 not installed pandas~=1.0 not installed pip 21.2.4 pydantic~=1.6 not installed python-libsbml 5.19.0 rich 10.11.0 ruamel.yaml~=0.16 not installed scipy 1.7.1 setuptools 58.1.0 swiglpk 5.0.3 tox 3.23.0 wheel 0.37.0 ```
cdiener commented 3 years ago

Oh I see, the problem is the filter in line 518. If the flux<threshold this will cause the key error. I don't get the logic though, what was the output supposed to be if the flux is basically zero? Shouldn't you still report the zero flux?

Midnighter commented 3 years ago

The summaries have always filtered low flux. I think the default was 1e-6. But indeed, need to handle the reactions not being present any longer.

cdiener commented 3 years ago

I think that makes sense for the model and metabolite summary but for the reaction I would just show the actual value. Setting it to zero wouldn't be correct either since we don't know if it's actually zero, so just showing the estimated flux seems okay to me.

Midnighter commented 3 years ago

So you wouldn't even set it to zero if the absolute value was below tolerance? I agree that it's weird to not show anything when you request the summary of a single reaction.

cdiener commented 3 years ago

Yeah, exactly. Because you don't know whether the flux is zero. Just that solver can't differentiate if the value is different from zero. An example would be if I set the reactions lower bound to 1e-7. In that case the summary would now show an infeasible flux even if the solution was feasible.

cdiener commented 3 years ago

Alternative suggestion: output 0 along with the tolerance range. So 0 ± 1e-6 for instance.