Gapfilling failed in context manager (infeasible)
Opened this issue · 3 comments
Is there an existing issue for this?
- I have searched the existing issues
Problem description
What I tried to achieve:
- I have been trying to apply gapfill to context-specific human GEMs, but have been encountering an infeasibility error.
How I went about it:
- To showcase this error, I tried removing an essential reaction from Human1, and then applying gapfill using the original Human1 GEM as the universal model. I expect this to simply return the removed reaction.
- I also tried lowering the
model.solver.configuration.tolerances.feasibilitythreshold to1e-9, as previous error threads showed this helped. But this yields the same error.
Why the current behavior is a problem:
- Gapfill fails, and lowering the feasibility threshold does not resolve the issue.
Code sample
Code run:
# Load Human1 GEM
human1 = read_sbml_model('Human-GEM.xml')
# Try gapfilling an essential reaction (MAR06669)
with human1 as test_model:
rxn = test_model.reactions.get_by_id('MAR06669')
test_model.remove_reactions([rxn])
test_model.solver.configuration.tolerances.feasibility = 1e-9
solution = gapfill(test_model, human1, demand_reactions=False)
for reaction in solution[0]:
print(reaction.id)Traceback:
---------------------------------------------------------------------------
Infeasible Traceback (most recent call last)
[~\AppData\Local\Temp\ipykernel_26524\1945771835.py] in <module>
5 print(test_model.slim_optimize())
6 test_model.solver.configuration.tolerances.feasibility = 1e-9
----> 7 solution = gapfill(test_model, human1, demand_reactions=False)
8 for reaction in solution[0]:
9 print(reaction.id)
[~\anaconda3\lib\site-packages\cobra\flux_analysis\gapfilling.py] in gapfill(model, universal, lower_bound, penalties, demand_reactions, exchange_reactions, iterations)
405 exchange_reactions=exchange_reactions,
406 )
--> 407 return gapfiller.fill(iterations=iterations)
[~\anaconda3\lib\site-packages\cobra\flux_analysis\gapfilling.py] in fill(self, iterations)
294 used_reactions = []
295 for _ in range(iterations):
--> 296 self.model.slim_optimize(
297 error_value=None, message="gap filling optimization failed"
298 )
[~\anaconda3\lib\site-packages\cobra\core\model.py] in slim_optimize(self, error_value, message)
1198 return error_value
1199 else:
-> 1200 assert_optimal(self, message)
1201
1202 def optimize(
[~\anaconda3\lib\site-packages\cobra\util\solver.py] in assert_optimal(model, message)
588 if status != OPTIMAL:
589 exception_cls = OPTLANG_TO_EXCEPTIONS_DICT.get(status, OptimizationError)
--> 590 raise exception_cls(f"{message} ({status}).")
591
592
Infeasible: gap filling optimization failed (infeasible).
Environment
| Package | Version |
|---|---|
| cobra | 0.29.0 |
Dependency Information
| Package | Version |
|---|---|
| appdirs | 1.4.4 |
| black | 22.6.0 |
| bumpversion | missing |
| depinfo | 2.2.0 |
| diskcache | 5.4.0 |
| future | 0.18.2 |
| httpx | 0.25.2 |
| importlib-resources | 5.12.0 |
| isort | 5.9.3 |
| numpy | 1.21.6 |
| optlang | 1.8.1 |
| pandas | 1.3.5 |
| pydantic | 1.10.7 |
| python-libsbml | 5.19.7 |
| rich | 13.3.2 |
| ruamel.yaml | 0.17.21 |
| scipy | 1.7.3 |
| swiglpk | 5.0.10 |
| tox | missing |
Build Tools Information
| Package | Version |
|---|---|
| conda | 23.1.0 |
| pip | 24.0 |
| setuptools | 69.0.2 |
| wheel | 0.37.1 |
Platform Information
| Windows | 10-AMD64 |
| CPython | 3.9.13 |
Anything else?
No response
Hi, thanks for providing a reproducible example. In this case not an issue with the gapfiller (at least for GLPK) but rather an edge case due to the use of then context manager.
In this line
# ...
with human1 as test_model:
# ...human1 and test_model point to the same object/model. So when you remove the reaction from test_model you also remove it from human1. The with block only reverts all changes after the exiting but it does not take a copy of the model. Substituting this with a copy fixes it, at least for me:
In [70]: # Load Human1 GEM
...: human1 = read_sbml_model('Human-GEM.xml')
...: human1.solver = "glpk"
...:
...: # Try gapfilling an essential reaction (MAR06669)
...: with human1.copy() as test_model: # <-- note the copy here
...: rxn = test_model.reactions.get_by_id('MAR06669')
...: test_model.remove_reactions([rxn])
...: solution = gapfill(test_model, human1, demand_reactions=False)
...: for reaction in solution[0]:
...: print(reaction.id)
...:
MAR06669
I did notice an issue with CPLEX due to a low integer primal and the lack of an option to set this via gapfill. I will send a PR for this part.
Ah thank you for catching that! I am still having issues with trying to use gapfill however. My goal is to identify alternative reactions/pathways to an essential reaction in a context-specific model. My reasoning was this: if a reaction is essential in the context-specific model but not in the universal model, then there should be alternative reactions/pathways. To identify these, I removed the essential reaction from both the context-specific and universal models, then tried to gap-fill. But I still get the infeasibility error. Is there an error in this reasoning?
Code:
# Load Human1 GEM and context-specific model
human1 = read_sbml_model('Human-GEM.xml')
context_specific_model= read_sbml_model('context_specific_model.xml')
# find alternative reactions to compensate for MAR01820
with context_specific_model as pruned_model:
with human1 as univ_model:
pruned_model.remove_reactions([pruned_model.reactions.MAR01820])
univ_model.remove_reactions([univ_model.reactions.MAR01820])
test_model.solver.configuration.tolerances.feasibility = 1e-9
solution = gapfill(pruned_model, univ_model, demand_reactions=False)
for reaction in solution[0]:
print(reaction.id)
Context-specific model:
context_specific_model.zip
Also, the original Human1 check still breaks when I try to gapfill some other essential reactions (e.g. MAR06657, MAR00578, MAR06702).