Python cobra.Model类代码示例

 def test_inequality(self):
     # The space enclosed by the constraints is a 2D triangle with
     # vertexes as (3, 0), (1, 2), and (0, 1)
     solver = self.solver
     # c1 encodes y - x > 1 ==> y > x - 1
     # c2 encodes y + x < 3 ==> y < 3 - x
     c1 = Metabolite("c1")
     c2 = Metabolite("c2")
     x = Reaction("x")
     x.lower_bound = 0
     y = Reaction("y")
     y.lower_bound = 0
     x.add_metabolites({c1: -1, c2: 1})
     y.add_metabolites({c1: 1, c2: 1})
     c1._bound = 1
     c1._constraint_sense = "G"
     c2._bound = 3
     c2._constraint_sense = "L"
     m = Model()
     m.add_reactions([x, y])
     # test that optimal values are at the vertices
     m.objective = "x"
     self.assertAlmostEqual(solver.solve(m).f, 1.0)
     self.assertAlmostEqual(solver.solve(m).x_dict["y"], 2.0)
     m.objective = "y"
     self.assertAlmostEqual(solver.solve(m).f, 3.0)
     self.assertAlmostEqual(solver.solve(m, objective_sense="minimize").f,
                            1.0)

 def test_solve_mip(self):
     solver = self.solver
     if not hasattr(solver, "_SUPPORTS_MILP") or not solver._SUPPORTS_MILP:
         self.skipTest("no milp support")
     cobra_model = Model('MILP_implementation_test')
     constraint = Metabolite("constraint")
     constraint._bound = 2.5
     x = Reaction("x")
     x.lower_bound = 0.
     x.objective_coefficient = 1.
     x.add_metabolites({constraint: 2.5})
     y = Reaction("y")
     y.lower_bound = 0.
     y.objective_coefficient = 1.
     y.add_metabolites({constraint: 1.})
     cobra_model.add_reactions([x, y])
     float_sol = solver.solve(cobra_model)
     # add an integer constraint
     y.variable_kind = "integer"
     int_sol = solver.solve(cobra_model)
     self.assertAlmostEqual(float_sol.f, 2.5)
     self.assertAlmostEqual(float_sol.x_dict["y"], 2.5)
     self.assertEqual(int_sol.status, "optimal")
     self.assertAlmostEqual(int_sol.f, 2.2)
     self.assertAlmostEqual(int_sol.x_dict["y"], 2.0)

def test_model_history(tmp_path):
    """Testing reading and writing of ModelHistory."""
    model = Model("test")
    model._sbml = {
        "creators": [{
            "familyName": "Mustermann",
            "givenName": "Max",
            "organisation": "Muster University",
            "email": "[email protected]",
        }]
    }

    sbml_path = join(str(tmp_path), "test.xml")
    with open(sbml_path, "w") as f_out:
        write_sbml_model(model, f_out)

    with open(sbml_path, "r") as f_in:
        model2 = read_sbml_model(f_in)

    assert "creators" in model2._sbml
    assert len(model2._sbml["creators"]) is 1
    c = model2._sbml["creators"][0]
    assert c["familyName"] == "Mustermann"
    assert c["givenName"] == "Max"
    assert c["organisation"] == "Muster University"
    assert c["email"] == "[email protected]"

class TestCobraSolver(TestCase):
    def setUp(self):
        self.model = create_test_model()
        initialize_growth_medium(self.model, "MgM")
        self.old_solution = 0.320064
        self.infeasible_model = Model()
        metabolite_1 = Metabolite("met1")
        # metabolite_2 = Metabolite("met2")
        reaction_1 = Reaction("rxn1")
        reaction_2 = Reaction("rxn2")
        reaction_1.add_metabolites({metabolite_1: 1})
        reaction_2.add_metabolites({metabolite_1: 1})
        reaction_1.lower_bound = 1
        reaction_2.upper_bound = 2
        self.infeasible_model.add_reactions([reaction_1, reaction_2])

 def test_change_coefficient(self):
     solver = self.solver
     c = Metabolite("c")
     c._bound = 6
     x = Reaction("x")
     x.lower_bound = 1.
     y = Reaction("y") 
     y.lower_bound = 0.
     x.add_metabolites({c: 1})
     #y.add_metabolites({c: 1})
     z = Reaction("z")
     z.add_metabolites({c: 1})
     z.objective_coefficient = 1
     m = Model("test_model")
     m.add_reactions([x, y, z])
     # change an existing coefficient
     lp = solver.create_problem(m)
     solver.solve_problem(lp)
     sol1 = solver.format_solution(lp, m)
     solver.change_coefficient(lp, 0, 0, 2)
     solver.solve_problem(lp)
     sol2 = solver.format_solution(lp, m)
     self.assertAlmostEqual(sol1.f, 5.0)
     self.assertAlmostEqual(sol2.f, 4.0)
     # change a new coefficient
     z.objective_coefficient = 0.
     y.objective_coefficient = 1.
     lp = solver.create_problem(m)
     solver.change_coefficient(lp, 0, 1, 2)
     solver.solve_problem(lp)
     solution = solver.format_solution(lp, m)
     self.assertAlmostEqual(solution.x_dict["y"], 2.5)

def minimized_shuffle(small_model):
    model = small_model.copy()
    chosen = sample(list(set(model.reactions) - set(model.exchanges)), 10)
    new = Model("minimized_shuffle")
    new.add_reactions(chosen)
    LOGGER.debug("'%s' has %d metabolites, %d reactions, and %d genes.",
                 new.id, new.metabolites, new.reactions, new.genes)
    return new

def gapfillfunc(model, database, runs):
    """ This function gapfills the model using the growMatch algorithm that is built into cobrapy

    Returns a dictionary which contains the pertinent information about the gapfilled model such as
    the reactions added, the major ins and outs of the system and the objective value of the gapfilled
    model.
    This function calls on two other functions sort_and_deduplicate to assure the uniqueness of the solutions
    and findInsAndOuts to find major ins and outs of the model when gapfilled when certain reactions
    Args:
        model: a model in SBML format
        database: an external database database of reactions to be used for gapfilling
        runs: integer number of iterations the gapfilling algorithm will run through
    """
    # Read model from SBML file and create Universal model to add reactions to
    func_model = create_cobra_model_from_sbml_file(model)
    Universal = Model("Universal Reactions")
    f = open(database, 'r')
    next(f)
    rxn_dict = {}
    # Creates a dictionary of the reactions from the tab delimited database, storing their ID and the reaction string
    for line in f:
        rxn_items = line.split('\t')
        rxn_dict[rxn_items[0]] = rxn_items[6], rxn_items[1]
    # Adds the reactions from the above dictionary to the Universal model
    for rxnName in rxn_dict.keys():
        rxn = Reaction(rxnName)
        Universal.add_reaction(rxn)
        rxn.reaction = rxn_dict[rxnName][0]
        rxn.name = rxn_dict[rxnName][1]

    # Runs the growMatch algorithm filling gaps from the Universal model
    result = flux_analysis.growMatch(func_model, Universal, iterations=runs)
    resultShortened = sort_and_deduplicate(uniq(result))
    rxns_added = {}
    rxn_met_list = []
    print resultShortened
    for x in range(len(resultShortened)):
        func_model_test = deepcopy(func_model)
        # print func_model_test.optimize().f
        for i in range(len(resultShortened[x])):
            addID = resultShortened[x][i].id
            rxn = Reaction(addID)
            func_model_test.add_reaction(rxn)
            rxn.reaction = resultShortened[x][i].reaction
            rxn.reaction = re.sub('\+ dummy\S+', '', rxn.reaction)
            rxn.name = resultShortened[x][i].name
            mets = re.findall('cpd\d{5}_c0|cpd\d{5}_e0', rxn.reaction)
            for met in mets:
                y = func_model_test.metabolites.get_by_id(met)
                rxn_met_list.append(y.name)
        obj_value = func_model_test.optimize().f
        fluxes = findInsAndOuts(func_model_test)
        sorted_outs = fluxes[0]
        sorted_ins = fluxes[1]
        rxns_added[x] = resultShortened[x], obj_value, sorted_ins, sorted_outs, rxn_met_list
        rxn_met_list = []
    return rxns_added

def tiny_toy_model():
    tiny = Model("Toy Model")
    m1 = Metabolite("M1")
    d1 = Reaction("ex1")
    d1.add_metabolites({m1: -1})
    d1.upper_bound = 0
    d1.lower_bound = -1000
    tiny.add_reactions([d1])
    tiny.objective = 'ex1'
    return tiny

 def test_quadratic(self):
     solver = self.solver
     if not hasattr(solver, "set_quadratic_objective"):
         self.skipTest("no qp support")
     c = Metabolite("c")
     c._bound = 2
     x = Reaction("x")
     x.objective_coefficient = -0.5
     x.lower_bound = 0.
     y = Reaction("y")
     y.objective_coefficient = -0.5
     y.lower_bound = 0.
     x.add_metabolites({c: 1})
     y.add_metabolites({c: 1})
     m = Model()
     m.add_reactions([x, y])
     lp = self.solver.create_problem(m)
     quadratic_obj = scipy.sparse.eye(2) * 2
     solver.set_quadratic_objective(lp, quadratic_obj)
     solver.solve_problem(lp, objective_sense="minimize")
     solution = solver.format_solution(lp, m)
     self.assertEqual(solution.status, "optimal")
     # Respecting linear objectives also makes the objective value 1.
     self.assertAlmostEqual(solution.f, 1.)
     self.assertAlmostEqual(solution.x_dict["y"], 1.)
     self.assertAlmostEqual(solution.x_dict["y"], 1.)
     # When the linear objectives are removed the objective value is 2.
     solver.change_variable_objective(lp, 0, 0.)
     solver.change_variable_objective(lp, 1, 0.)
     solver.solve_problem(lp, objective_sense="minimize")
     solution = solver.format_solution(lp, m)
     self.assertEqual(solution.status, "optimal")
     self.assertAlmostEqual(solution.f, 2.)
     # test quadratic from solve function
     solution = solver.solve(m, quadratic_component=quadratic_obj,
                             objective_sense="minimize")
     self.assertEqual(solution.status, "optimal")
     self.assertAlmostEqual(solution.f, 1.)
     c._bound = 6
     z = Reaction("z")
     x.objective_coefficient = 0.
     y.objective_coefficient = 0.
     z.lower_bound = 0.
     z.add_metabolites({c: 1})
     m.add_reaction(z)
     solution = solver.solve(m, quadratic_component=scipy.sparse.eye(3),
                             objective_sense="minimize")
     # should be 12 not 24 because 1/2 (V^T Q V)
     self.assertEqual(solution.status, "optimal")
     self.assertAlmostEqual(solution.f, 6)
     self.assertAlmostEqual(solution.x_dict["x"], 2, places=6)
     self.assertAlmostEqual(solution.x_dict["y"], 2, places=6)
     self.assertAlmostEqual(solution.x_dict["z"], 2, places=6)

 def test_remove_breaks(self):
     model = Model("test model")
     A = Metabolite("A")
     r = Reaction("r")
     r.add_metabolites({A: -1})
     r.lower_bound = -1000
     r.upper_bound = 1000
     model.add_reaction(r)
     convert_to_irreversible(model)
     model.remove_reactions(["r"])
     with pytest.raises(KeyError):
         revert_to_reversible(model)

def convert_to_cobra_model(the_network):
    """ Take a generic NAMpy model and convert to
    a COBRA model.  The model is assumed to be monopartite.
    You need a functional COBRApy for this.

    Arguments:
     the_network

    kwargs:
     flux_bounds


    """
    continue_flag = True
    try:
        from cobra import Model, Reaction, Metabolite
    except:
        print 'This function requires a functional COBRApy, exiting ...'

    if continue_flag:
        __default_objective_coefficient = 0
        
        if 'flux_bounds' in kwargs:
            flux_bounds = kwargs['flux_bounds'] 
        else:
            flux_bounds = len(the_nodes)

        metabolite_dict = {}
        for the_node in the_network.nodetypes[0].nodes: 
            the_metabolite = Metabolite(the_node.id)
            metabolite_dict.update({the_node.id: the_metabolite})

        cobra_reaction_list = []
        for the_edge in the_network.edges:
            the_reaction = Reaction(the_edge.id)
            cobra_reaction_list.append(the_reaction)
            the_reaction.upper_bound = flux_bounds
            the_reaction.lower_bound = -1 * flux_bounds
            cobra_metabolites = {}
            the_metabolite_id_1 = the_edge._nodes[0].id
            the_metabolite_id_2 = the_edge._nodes[1].id
            cobra_metabolites[metabolite_dict[the_metabolite_id_1]] = 1.
            cobra_metabolites[metabolite_dict[the_metabolite_id_2]] = -1.
            reaction.add_metabolites(cobra_metabolites)
            reaction.objective_coefficient = __default_objective_coefficient

        cobra_model = Model(model_id)
        cobra_model.add_reactions(cobra_reaction_list)

        return cobra_model
    else:
        return None

 def test_loopless(self):
     try:
         solver = get_solver_name(mip=True)
     except:
         self.skipTest("no MILP solver found")
     test_model = Model()
     test_model.add_metabolites(Metabolite("A"))
     test_model.add_metabolites(Metabolite("B"))
     test_model.add_metabolites(Metabolite("C"))
     EX_A = Reaction("EX_A")
     EX_A.add_metabolites({test_model.metabolites.A: 1})
     DM_C = Reaction("DM_C")
     DM_C.add_metabolites({test_model.metabolites.C: -1})
     v1 = Reaction("v1")
     v1.add_metabolites({test_model.metabolites.A: -1, test_model.metabolites.B: 1})
     v2 = Reaction("v2")
     v2.add_metabolites({test_model.metabolites.B: -1, test_model.metabolites.C: 1})
     v3 = Reaction("v3")
     v3.add_metabolites({test_model.metabolites.C: -1, test_model.metabolites.A: 1})
     DM_C.objective_coefficient = 1
     test_model.add_reactions([EX_A, DM_C, v1, v2, v3])
     feasible_sol = construct_loopless_model(test_model).optimize()
     v3.lower_bound = 1
     infeasible_sol = construct_loopless_model(test_model).optimize()
     self.assertEqual(feasible_sol.status, "optimal")
     self.assertEqual(infeasible_sol.status, "infeasible")

 def test_gene_knockout_computation(self):
     cobra_model = self.model
     delete_model_genes = delete.delete_model_genes
     get_removed = lambda m: {x.id for x in m._trimmed_reactions}
     gene_list = ['STM1067', 'STM0227']
     dependent_reactions = {'3HAD121', '3HAD160', '3HAD80', '3HAD140',
                            '3HAD180', '3HAD100', '3HAD181','3HAD120',
                            '3HAD60', '3HAD141', '3HAD161', 'T2DECAI',
                            '3HAD40'}
     delete_model_genes(cobra_model, gene_list)
     self.assertEqual(get_removed(cobra_model), dependent_reactions)
     # cumulative
     delete_model_genes(cobra_model, ["STM4221"],
                        cumulative_deletions=True)
     dependent_reactions.add('PGI')
     self.assertEqual(get_removed(cobra_model), dependent_reactions)
     # non-cumulative
     delete_model_genes(cobra_model, ["STM4221"],
                        cumulative_deletions=False)
     self.assertEqual(get_removed(cobra_model), {'PGI'})
     # make sure on reset that the bounds are correct
     reset_bound = cobra_model.reactions.get_by_id("T2DECAI").upper_bound
     self.assertEqual(reset_bound, 1000.)
     # test computation when gene name is a subset of another
     test_model = Model()
     test_reaction_1 = Reaction("test1")
     test_reaction_1.gene_reaction_rule = "eggs or (spam and eggspam)"
     test_model.add_reaction(test_reaction_1)
     delete.delete_model_genes(test_model, ["eggs"])
     self.assertEqual(get_removed(test_model), set())
     delete_model_genes(test_model, ["spam"], cumulative_deletions=True)
     self.assertEqual(get_removed(test_model), {'test1'})
     # test computation with nested boolean expression
     delete.undelete_model_genes(test_model)
     test_reaction_1.gene_reaction_rule = \
         "g1 and g2 and (g3 or g4 or (g5 and g6))"
     delete_model_genes(test_model, ["g3"], cumulative_deletions=False)
     self.assertEqual(get_removed(test_model), set())
     delete_model_genes(test_model, ["g1"], cumulative_deletions=False)
     self.assertEqual(get_removed(test_model), {'test1'})
     delete_model_genes(test_model, ["g5"], cumulative_deletions=False)
     self.assertEqual(get_removed(test_model), set())
     delete_model_genes(test_model, ["g3", "g4", "g5"],
                        cumulative_deletions=False)
     self.assertEqual(get_removed(test_model), {'test1'})

def model():
    A = Metabolite("A")
    B = Metabolite("B")
    C = Metabolite("C")
    r1 = Reaction("r1")
    r1.add_metabolites({A: -1, C: 1})
    r2 = Reaction("r2")
    r2.add_metabolites({B: -1, C: 1})
    r3 = Reaction("EX_A")
    r3.add_metabolites({A: 1})
    r4 = Reaction("EX_B")
    r4.add_metabolites({B: 1})
    r5 = Reaction("EX_C")
    r5.add_metabolites({C: -1})
    mod = Model("test model")
    mod.add_reactions([r1, r2, r3, r4, r5])
    conf = {"r1": 1, "r2": -1, "EX_A": 1, "EX_B": 1, "EX_C": 1}

    return (mod, conf)

 def test_solve_mip(self):
     solver = self.solver
     cobra_model = Model('MILP_implementation_test')
     constraint = Metabolite("constraint")
     constraint._bound = 2.5
     x = Reaction("x")
     x.lower_bound = 0.
     x.objective_coefficient = 1.
     x.add_metabolites({constraint: 2.5})
     y = Reaction("y")
     y.lower_bound = 0.
     y.objective_coefficient = 1.
     y.add_metabolites({constraint: 1.})
     cobra_model.add_reactions([x, y])
     float_sol = solver.solve(cobra_model)
     # add an integer constraint
     y.variable_kind = "integer"
     int_sol = solver.solve(cobra_model)
     self.assertAlmostEqual(float_sol.f, 2.5)
     self.assertAlmostEqual(float_sol.x_dict["y"], 2.5)
     self.assertAlmostEqual(int_sol.f, 2.2)
     self.assertAlmostEqual(int_sol.x_dict["y"], 2.0)