Gapfilling Algorthim Failing

[jccvila]

So i'm trying to create strain-specific models of p.putida similar to nogales et al 2020. I'm starting with the p.putida iJN1463 mode and have a list of genes to delete to recreate each strain. I then want to gapfill using the original model as a reference to add back the reactions that would be needed for growth. However the gap-filling algorithm in cobrapy is unable to find a solution (even though one has to exist).

Code Sample

I've created a minimum reproducible example of this in which i delete a single essential gene that controls 1 reaction and show that the gap-filling still fails. Is there something obvious that i am doing wrong?

#Load Pputida model
P =c.io.read_sbml_model('../Data/Bigg_Model/iJN1463.xml')
#confirm it can grow
P.optimize()
P.summary()
#create copy
P2 = P.copy()
# delete genes from copy
genes_to_delete = ['PP_0380']
c.manipulation.delete_model_genes(P2,genes_to_delete)
#confirm it can't groww
P2.optimize()
P2.summary()
# gapfill mutant using original putida model
solution = c.flux_analysis.gapfill(P2, P,demand_reactions=False)

The error

N_FLUXES IN_FLUXES OUT_FLUXES OUT_FLUXES      OBJECTIVES     OBJECTIVES
    ID       FLUX       ID        FLUX             ID            FLUX   
     o2_e    11.4      h2o_e      27.8     BIOMASS_KT2440_WT3   0.586   
    nh4_e    6.24      co2_e      12.4                            nan   
 glc__D_e       6        h_e      5.74                            nan   
     pi_e   0.539                  nan                            nan   
Read LP format model from file /tmp/tmp0wxnzc5e.lp
Reading time = 0.01 seconds
: 2153 rows, 5854 columns, 23272 nonzeros
IN_FLUXES IN_FLUXES OUT_FLUXES OUT_FLUXES      OBJECTIVES     OBJECTIVES
    ID       FLUX       ID        FLUX             ID            FLUX   
 glc__D_e    0.92     glcn_e      0.92     BIOMASS_KT2440_WT3      0    
     o2_e    0.46        h_e      0.92                           nan    
Read LP format model from file /tmp/tmp49ip_w5z.lp
Reading time = 0.01 seconds
: 2153 rows, 5854 columns, 23272 nonzeros
Read LP format model from file /tmp/tmp4o0bofdu.lp
Reading time = 0.02 seconds
: 2153 rows, 5854 columns, 23272 nonzeros
---------------------------------------------------------------------------
Infeasible                                Traceback (most recent call last)
<ipython-input-147-d9a73a74eeb9> in <module>
     13 print(P2.summary())
     14 # gapfill mutant using original putida model
---> 15 solution = c.flux_analysis.gapfill(P2, P,demand_reactions=False)

~/anaconda3/lib/python3.7/site-packages/cobra/flux_analysis/gapfilling.py in gapfill(model, universal, lower_bound, penalties, demand_reactions, exchange_reactions, iterations)
    311                           demand_reactions=demand_reactions,
    312                           exchange_reactions=exchange_reactions)
--> 313     return gapfiller.fill(iterations=iterations)

~/anaconda3/lib/python3.7/site-packages/cobra/flux_analysis/gapfilling.py in fill(self, iterations)
    232         for i in range(iterations):
    233             self.model.slim_optimize(error_value=None,
--> 234                                      message='gapfilling optimization failed')
    235             solution = [self.model.reactions.get_by_id(ind.rxn_id)
    236                         for ind in self.indicators if

~/anaconda3/lib/python3.7/site-packages/cobra/core/model.py in slim_optimize(self, error_value, message)
   1053             return error_value
   1054         else:
-> 1055             assert_optimal(self, message)
   1056 
   1057     def optimize(self, objective_sense=None, raise_error=False):

~/anaconda3/lib/python3.7/site-packages/cobra/util/solver.py in assert_optimal(model, message)
    439         exception_cls = OPTLANG_TO_EXCEPTIONS_DICT.get(
    440             status, OptimizationError)
--> 441         raise exception_cls("{} ({})".format(message, status))
    442 
    443 

Infeasible: gapfilling optimization failed (infeasible)

[gregmedlock]

Hi Jean,

The delete_model_genes function just closes the upper and lower bounds for reactions that are no longer functional with the gene deleted; so the issue is likely that your model with the deletion still has the original reactions (with bounds of 0,0), but reactions with the same id cannot be added from the universal reaction set. Thus, a feasible solution can’t be found because the reactions of interest are already in the model but with bounds that make the model infeasible. Try again with either 1) a full removal of the target reactions or 2) renaming the reactions in the model copy that you’re using as the universal.

Hope this works for you!

On Mon, Feb 3, 2020 at 5:57 PM Jean Vila notifications@github.com wrote:

So i'm trying to create strain-specific models of p.putida similar to nogales et al 2020. I'm starting with the p.putida iJN1463 mode and have a list of genes to delete to recreate each strain. I then want to gapfill using the original model as a reference to add back the reactions that would be needed for growth. However the gap-filling algorithm in cobrapy is unable to find a solution (even though one has to exist).

Code Sample

I've created a minimum reproducible example of this in which i delete a single essential gene that controls 1 reaction and show that the gap-filling still fails. Is there something obvious that i am doing wrong?

Load Pputida model

P =c.io.read_sbml_model('../Data/Bigg_Model/iJN1463.xml')

confirm it can grow

P.optimize() P.summary()

create copy

P2 = P.copy()

delete genes from copy

genes_to_delete = ['PP_0380'] c.manipulation.delete_model_genes(P2,genes_to_delete)

confirm it can't groww

P2.optimize() P2.summary()

gapfill mutant using original putida model

solution = c.flux_analysis.gapfill(P2, P,demand_reactions=False)

The error

N_FLUXES IN_FLUXES OUT_FLUXES OUT_FLUXES OBJECTIVES OBJECTIVES ID FLUX ID FLUX ID FLUX o2_e 11.4 h2o_e 27.8 BIOMASS_KT2440_WT3 0.586 nh4_e 6.24 co2_e 12.4 nan glc__D_e 6 h_e 5.74 nan pi_e 0.539 nan nan Read LP format model from file /tmp/tmp0wxnzc5e.lp Reading time = 0.01 seconds : 2153 rows, 5854 columns, 23272 nonzeros IN_FLUXES IN_FLUXES OUT_FLUXES OUT_FLUXES OBJECTIVES OBJECTIVES ID FLUX ID FLUX ID FLUX glc__D_e 0.92 glcn_e 0.92 BIOMASS_KT2440_WT3 0 o2_e 0.46 h_e 0.92 nan Read LP format model from file /tmp/tmp49ip_w5z.lp Reading time = 0.01 seconds : 2153 rows, 5854 columns, 23272 nonzeros Read LP format model from file /tmp/tmp4o0bofdu.lp Reading time = 0.02 seconds : 2153 rows, 5854 columns, 23272 nonzeros

Infeasible Traceback (most recent call last)

in 13 print(P2.summary()) 14 # gapfill mutant using original putida model ---> 15 solution = c.flux_analysis.gapfill(P2, P,demand_reactions=False) ~/anaconda3/lib/python3.7/site-packages/cobra/flux_analysis/gapfilling.py in gapfill(model, universal, lower_bound, penalties, demand_reactions, exchange_reactions, iterations) 311 demand_reactions=demand_reactions, 312 exchange_reactions=exchange_reactions) --> 313 return gapfiller.fill(iterations=iterations) ~/anaconda3/lib/python3.7/site-packages/cobra/flux_analysis/gapfilling.py in fill(self, iterations) 232 for i in range(iterations): 233 self.model.slim_optimize(error_value=None, --> 234 message='gapfilling optimization failed') 235 solution = [self.model.reactions.get_by_id(ind.rxn_id) 236 for ind in self.indicators if ~/anaconda3/lib/python3.7/site-packages/cobra/core/model.py in slim_optimize(self, error_value, message) 1053 return error_value 1054 else: -> 1055 assert_optimal(self, message) 1056 1057 def optimize(self, objective_sense=None, raise_error=False): ~/anaconda3/lib/python3.7/site-packages/cobra/util/solver.py in assert_optimal(model, message) 439 exception_cls = OPTLANG_TO_EXCEPTIONS_DICT.get( 440 status, OptimizationError) --> 441 raise exception_cls("{} ({})".format(message, status)) 442 443 Infeasible: gapfilling optimization failed (infeasible) — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub , or unsubscribe . -- Greg Medlock, PhD Postdoctoral Research Associate Moore Lab, Child Health Research Center Department of Pediatrics, Division of Gastroenterology & Nutrition University of Virginia School of Medicine

[jccvila]

Hi greg.

Thanks for the response. So deleting one reaction rather than just using delete_model_genes still doest work though it gives me a slightly different error. I also get the same error when i delete the gene using (delete_model_genes) and then remove the reaction using (.delete())

New Code

#Load Pputida model
P = CAFBA_Model(c.io.read_sbml_model('../Data/Bigg_Model/iJN1463.xml'))
#confirm it can grow
P.optimize()
P.summary()
#create copy
P2 = P.copy()
# delete reaction
P2.reactions.PQQAC.delete()
P2.prune_model()
#confirm it can't groww
P2.optimize()
P2.summary()
# gapfill mutant using original putida model
solution = c.flux_analysis.gapfill(P2, P,demand_reactions=False)

Error

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-12-968017868380> in <module>
     15 print(P2.summary())
     16 # gapfill mutant using original putida model
---> 17 solution = c.flux_analysis.gapfill(P2, P,demand_reactions=False)

~/anaconda3/lib/python3.7/site-packages/cobra/flux_analysis/gapfilling.py in gapfill(model, universal, lower_bound, penalties, demand_reactions, exchange_reactions, iterations)
    311                           demand_reactions=demand_reactions,
    312                           exchange_reactions=exchange_reactions)
--> 313     return gapfiller.fill(iterations=iterations)

~/anaconda3/lib/python3.7/site-packages/cobra/flux_analysis/gapfilling.py in fill(self, iterations)
    237                         ind._get_primal() > self.integer_threshold]
    238             if not self.validate(solution):
--> 239                 raise RuntimeError('failed to validate gapfilled model, '
    240                                    'try lowering the integer_threshold')
    241             used_reactions.append(solution)

RuntimeError: failed to validate gapfilled model, try lowering the integer_threshold

[gregmedlock]

Hmm, the next thing I would check is the value of the objective returned in P.optimize(). gapfill() has a default lower bound through the former objective of 0.05, so if your model initially has a lower flux than that, it could raise the error you're getting.

[jccvila]

The value of the objective is 0.586 so i don't think that is the issue.

[gregmedlock]

I see--I just did a bit of digging around the gap filling module using the BIGG reconstruction you're using, and it looks like the primal value is not properly getting updated on the reaction that should be added, thus the list of reactions generated in GapFiller.fill() is empty, even though there is flux through PQQAC in the gap fill solution. I will investigate a bit further and make a PR to fix the issue.

[gregmedlock]

Hi Jean--finished with some more tinkering and found that it is likely a numerical issue. If you set the integer_threshold to 1e-12 or lower in the gapfill function, it should work. The integer_threshold parameter is not exposed in the gapfill() function, but it can be accessed through the GapFiller object used by gapfill(), e.g.,

gf = cobra.flux_analysis.gapfilling.GapFiller(P2,P, integer_threshold=1e-20)
gapfill_solution = gf.fill()

In some cases, lowering the integer_threshold could lead to false positive reaction additions (i.e., returned as part of the solution but not actually required for growth), but that does not seem to be the case with this specific model unless you go with something ludicrous like 1e-20.

[jccvila]

What solver/version are you using? With Gurobi i get

GurobiError: Unable to set parameter IntFeasTol to value 1e-12 (minimum is 1e-09)

Update: I tried switching the solver to cplex or glpk and for both your solution works (though i have to set demand_reactions=False). Thanks for the help.

[gregmedlock]

Great--glad it worked!