<jupyter_start><jupyter_text>NEW FUNCTIONS<jupyter_code>def gapfilling(model, template, integer_threshold, it = 1):
    """
    Calls GapFiller class from cobrapy.
    """
    # gurobi max integer threshold value -> 1e-9
    gapfiller = GapFiller(model, template, integer_threshold = integer_threshold,
                         demand_reactions = False, exchange_reactions=False)
    # not valid for gurobi:
    # gapfiller.model.solver.configuration.tolerances.feasability = 1e-10
    return gapfiller.fill(iterations = it)

def add_exchange_reactions(model, template):
    """
    Adds all exchange reactions from a template matching model metabolites. 
    """
    EX = [template.reactions.get_by_id(i.id) for i in template.reactions if i.id.startswith("EX_")]
    for reaction in EX:
        if reaction not in model.reactions and str(list(reaction.metabolites.keys())[0]) in [i.id for i in model.metabolites]:
            model.add_reaction(reaction.copy())
    return model
# add_exchange_reactions test
x = read_sbml_model("BiGG_files/e_coli_core.xml")
x<jupyter_output>--------------------------------------------
Warning: your license will expire in 7 days
--------------------------------------------

Using license file /home/fco/gurobi.lic
Academic license - for non-commercial use only - expires 2021-03-12<jupyter_text>95 reactions<jupyter_code>universal = cobra.Model("universal_reactions")
for i in [i.id for i in x.metabolites.glc__D_e.reactions]:
    reaction = x.reactions.get_by_id(i)
    universal.add_reaction(reaction.copy())
    x.remove_reactions([reaction])
universal.reactions
x<jupyter_output><empty_output><jupyter_text>93 reactions; GLCpts and EX_glc__D_e have been removed from the model.<jupyter_code>add_exchange_reactions(x, universal)<jupyter_output><empty_output><jupyter_text>Unlike GLCpts, EX_glc__D_e has been added to the model.<jupyter_code>x.reactions.EX_glc__D_e<jupyter_output><empty_output><jupyter_text>HOMOLOGY GAP FILLING FUNCTION<jupyter_code>def homology_gapfilling(model, templates, model_obj = None, template_obj = None, use_all_templates = False,
                       integer_threshold = 1e-6, force_exchange = False, force_transport = True):
    """
    Performs gap filling on a model using homology models as templates.
    """
    model = load_query_model(model, obj = model_obj)
    model.solver = 'gurobi'
    templates = load_template_models(templates, obj = template_obj)
    # this dict will store used models, genes and reactions
    added_reactions = {}
    # initial flux value
    value = model.optimize().objective_value
    if value == None:
        value = 0.0
    if use_all_templates == False:
        for template in templates:
            # adding exchange reactions
            if force_exchange == True:
                add_exchange_reactions(model, template)
            # adding transport reactions
            if force_transport == True:
                add_transport(model, template)
            template.solver = 'gurobi'
            try:
                # result variable will store the reactions ids
                result = gapfilling(model, template, integer_threshold = integer_threshold)
                # reactions "log"
                log = []
                for reaction in result[0]:
                    if reaction.id.startswith("EX_"):
                        log.append((reaction.id, "Exchange reaction"))
                    else:
                        log.append((reaction.id, [str(list(reaction.genes)[i]) 
                                              for i in range(len(list(reaction.genes)))]))
                added_reactions[str(template)] = log
                # Adding reactions to the model
                [model.add_reaction(reaction.copy()) for reaction in result[0]]
                # Flux will be evaluated here
                new_value = model.optimize().objective_value
                if new_value != None and new_value > value:
                    value = new_value
                elif new_value == None:
                    continue
                elif new_value != None and new_value == value:
                    break
            except RuntimeError:
                print("\n" + str(template) + ": failed to validate gapfilled model, try lowering the integer_threshold")
            except exceptions.Infeasible:
                print("\n" + str(template) + ": gapfilling optimization failed (infeasible)")
        return model, added_reactions
    else:
        for template in templates:
            # adding exchange reactions
            if force_exchange == True:
                add_exchange_reactions(model, template)
            # adding transport reactions
            if force_transport == True:
                add_transport(model, template)
            template.solver = 'gurobi'
            try:
                result = gapfilling(model, template, integer_threshold = integer_threshold)
                log = []
                for reaction in result[0]:
                    if reaction.id.startswith("EX_"):
                        log.append((reaction.id, "Exchange reaction"))
                    else:
                        log.append((reaction.id, [str(list(reaction.genes)[i]) 
                                              for i in range(len(list(reaction.genes)))]))
                added_reactions[str(template)] = log
                [model.add_reaction(reaction.copy()) for reaction in result[0]]
            except RuntimeError:
                print("\n" + str(template) + ": failed to validate gapfilled model, try lowering the integer_threshold")
            except exceptions.Infeasible:
                print("\n" + str(template) + ": gapfilling optimization failed (infeasible)") 
        return model, added_reactions<jupyter_output><empty_output><jupyter_text>LOADING MODELS<jupyter_code># Pseudomonas putida KT2440
iJN746 = read_sbml_model("BiGG_files/iJN746.xml")
iJN746
# Extended model
iJN1463 = read_sbml_model("BiGG_files/iJN1463.xml")
iJN1463
# Chlamy<NAME>
iRC1080 = read_sbml_model("BiGG_files/iRC1080.xml")
iRC1080
"""
# adding an objective
iRC1080 = load_query_model(iRC1080, obj = "biomass")
"""
# <NAME>
# Core-reduced
e_coli_core = read_sbml_model("BiGG_files/e_coli_core.xml")
e_coli_core
# Strain C
iEC1344_C = read_sbml_model("BiGG_files/iEC1344_C.xml")
iEC1344_C
# Strain W
iEC1364_W = read_sbml_model("BiGG_files/iEC1364_W.xml")
iEC1364_W
# Bacillus subtilis
iYO844 = read_sbml_model("BiGG_files/iYO844.xml")
iYO844 
# Saccharomyces cerevisiae
iMM904 = read_sbml_model("BiGG_files/iMM904.xml")
iMM904
iND750 = read_sbml_model("BiGG_files/iND750.xml")
iND750
# Mus musculus
iMM1415 = read_sbml_model("BiGG_files/iMM1415.xml")
iMM1415<jupyter_output>No objective coefficients in model. Unclear what should be optimized<jupyter_text>TESTING<jupyter_code># Case 1a: gapfill a simple model with itself (biomass)
e_coli_core2 = read_sbml_model("BiGG_files/e_coli_core.xml")
e_coli_core2
# moving essential reactions to a new model
universal = cobra.Model("universal_reactions")
for i in [i.id for i in e_coli_core2.metabolites.glc__D_e.reactions]:
    reaction = e_coli_core2.reactions.get_by_id(i)
    universal.add_reaction(reaction.copy())
    e_coli_core2.remove_reactions([reaction])
e_coli_core2.optimize().objective_value
e_coli_core2_gp, reactions1 = homology_gapfilling(e_coli_core2, templates=[universal])
reactions1
e_coli_core2_gp.summary()
# Case 1b: gapfill a model with itself (reaction)
e_coli_core3 = read_sbml_model("BiGG_files/e_coli_core.xml")
e_coli_core3.reactions.CO2t
# We set an objective other than biomass
e_coli_core3 = load_query_model(e_coli_core3, obj='CO2t')
e_coli_core3.summary()
# we will eliminate reactions which don't contain co2_e or co2_c
universal = cobra.Model("universal_reactions")
for i in [i.id for i in e_coli_core3.metabolites.succ_e.reactions]:
    reaction = e_coli_core3.reactions.get_by_id(i)
    universal.add_reaction(reaction.copy())
    e_coli_core3.remove_reactions([reaction])
e_coli_core3.optimize().objective_value
universal.reactions
e_coli_core3.optimize()
# standard procedure
e_coli_core3_gp, reactions2 = homology_gapfilling(e_coli_core3, templates = [universal], template_obj='CO2t')
e_coli_core3_gp.optimize().objective_value
reactions2<jupyter_output><empty_output><jupyter_text>It doesn't add the succinate reactions although they would improve the flux<jupyter_code># forcing succinate addition
with e_coli_core3:
    e_coli_core3.objective = e_coli_core3.add_boundary(e_coli_core3.metabolites.succ_e, type='demand')
    solution = gapfill(e_coli_core3, universal)
    print([reaction for reaction in solution[0]])
    [e_coli_core3.add_reaction(reaction.copy()) for reaction in solution[0]]
e_coli_core3.optimize().objective_value<jupyter_output><empty_output><jupyter_text>No reactions have been added from universal; it is a viable model but it should have add them to maximize fluxes Another case is the following: removing certain reactions from a model, for example EX_CO2, and try to optimize that reaction. That gap filling would be infeasible.<jupyter_code># Case 1c: gapfill a complex model with itself (biomass)
iMM904.summary()
universal = cobra.Model("universal_reactions")
for i in [i.id for i in iMM904.metabolites.g6p_r.reactions]:
    reaction = iMM904.reactions.get_by_id(i)
    universal.add_reaction(reaction.copy())
    iMM904.remove_reactions([reaction])
iMM904.optimize().objective_value
universal.reactions
iMM904_gp, reactions3 = homology_gapfilling(iMM904, [universal])<jupyter_output>Read LP format model from file /tmp/tmpcf52haez.lp
Reading time = 0.01 seconds
: 1226 rows, 3150 columns, 13020 nonzeros
Read LP format model from file /tmp/tmpxoskyawq.lp
Reading time = 0.01 seconds
: 1226 rows, 3150 columns, 13020 nonzeros
Read LP format model from file /tmp/tmpihijg1gb.lp
Reading time = 0.00 seconds
: 6 rows, 4 columns, 14 nonzeros
Read LP format model from file /tmp/tmp7o_w29c9.lp
Reading time = 0.00 seconds
: 6 rows, 4 columns, 14 nonzeros
universal_reactions: failed to validate gapfilled model, try lowering the integer_threshold<jupyter_text>This time the error is related to the integer threshold of the Gapfiller class.<jupyter_code># lowering integer threshold
iMM904_gp, reactions3 = homology_gapfilling(iMM904, [universal], integer_threshold=1e-09)<jupyter_output>Read LP format model from file /tmp/tmp29z6qt5p.lp
Reading time = 0.01 seconds
: 1226 rows, 3150 columns, 13020 nonzeros
Read LP format model from file /tmp/tmp2r1nqly3.lp
Reading time = 0.01 seconds
: 1226 rows, 3150 columns, 13020 nonzeros
Read LP format model from file /tmp/tmpw3lv9bvs.lp
Reading time = 0.00 seconds
: 6 rows, 4 columns, 14 nonzeros
Read LP format model from file /tmp/tmp5d1r7p62.lp
Reading time = 0.00 seconds
: 6 rows, 4 columns, 14 nonzeros
universal_reactions: failed to validate gapfilled model, try lowering the integer_threshold<jupyter_text>Integer threshold must be lowered but 1e-9 is the minimum for gurobi<jupyter_code># adding exchange reactions
iMM904_gp, reactions3 = homology_gapfilling(iMM904, [universal], force_exchange=True, integer_threshold = 1e-09)
# Case 2a: Using other models as templates
e_coli_core4 = read_sbml_model("BiGG_files/e_coli_core.xml")
e_coli_core4.optimize().objective_value
t = load_template_models([iEC1344_C, iEC1364_W])
e_coli_core4_gp, reactions4 = homology_gapfilling(e_coli_core4, templates = t, use_all_templates=True)
e_coli_core4_gp.optimize().objective_value
reactions4<jupyter_output><empty_output><jupyter_text>No reactions added. That could mean that the core model has all the needed reactions for an optimal biomass flux.<jupyter_code># Case 2b: using more models as templates
e_coli_core5 = read_sbml_model("BiGG_files/e_coli_core.xml")
t = load_template_models([iEC1344_C, iEC1364_W, iJN1463, iJN746, iMM904, iND750])
e_coli_core5_gp, reactions5 = homology_gapfilling(e_coli_core5, t, use_all_templates=True, integer_threshold = 1e-9)
reactions5
# trying with simpler models
t = load_template_models([iEC1344_C, iEC1364_W, iJN1463, iJN746, iYO844])
e_coli_core5_gp, reactions5 = homology_gapfilling(e_coli_core5, t, use_all_templates=True, integer_threshold = 1e-9)
# Case 2c: using other models as templates (reaction)
e_coli_core6 = load_query_model(e_coli_core, obj='CO2t')
e_coli_core6.optimize().objective_value
e_coli_core6_gp, reactions6 = homology_gapfilling(e_coli_core6, t, use_all_templates=True, integer_threshold=1e-9)<jupyter_output>Read LP format model from file /tmp/tmprm4tymvp.lp
Reading time = 0.00 seconds
: 72 rows, 190 columns, 720 nonzeros
Read LP format model from file /tmp/tmpn0xt7cdh.lp
Reading time = 0.00 seconds
: 72 rows, 190 columns, 720 nonzeros
Read LP format model from file /tmp/tmp1kwpjhck.lp
Reading time = 0.01 seconds
: 1934 rows, 5452 columns, 21416 nonzeros
Read LP format model from file /tmp/tmp2ez4oeq8.lp
Reading time = 0.01 seconds
: 1934 rows, 5452 columns, 21416 nonzeros
Read LP format model from file /tmp/tmpl4jfuzus.lp
Reading time = 0.00 seconds
: 72 rows, 190 columns, 720 nonzeros
Read LP format model from file /tmp/tmpt1g6qgo1.lp
Reading time = 0.00 seconds
: 72 rows, 190 columns, 720 nonzeros
Read LP format model from file /tmp/tmp3oodpnhp.lp
Reading time = 0.01 seconds
: 1927 rows, 5528 columns, 21688 nonzeros
Read LP format model from file /tmp/tmpayi2hi6y.lp
Reading time = 0.01 seconds
: 1927 rows, 5528 columns, 21688 nonzeros
Read LP format model from file /tmp/tmp23hcxqjd.lp
R[...]
{'fraction_non_alphanumeric': 0.09215599051008304, 'fraction_numerical': 0.039961447212336895, 'mean_word_length': 3.732982456140351, 'pattern_counts': {'":': 4, '<': 46, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 15}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '371116', 'license': 'MIT', 'n_tokens_mistral': 4836, 'n_tokens_neox': 4281, 'n_words': 1101, 'path': 'gapfilling/testing/Homology_gapfilling_test.ipynb', 'repo_name': 'fcomnozz/GALLANT_gapfilling', 'size': '13511'}