diff --git a/src/emll/test_antimony_utils.py b/src/emll/test_antimony_utils.py new file mode 100644 index 0000000..dd8fa5c --- /dev/null +++ b/src/emll/test_antimony_utils.py @@ -0,0 +1,41 @@ +"""Test file for Antimony to Cobra conversion""" + +import logging +from pathlib import Path + +import cobra +import pytest +import tellurium as te + +from emll.util import ant_to_cobra + +logging.getLogger("cobra").setLevel(logging.ERROR) + +HERE = Path(__file__).parent.resolve() +MODEL = HERE.joinpath("test_models/unlabeled18.ant") + + +@pytest.mark.parametrize( + "antimony_path", + [ + str(MODEL), + ], +) +def test_antimony_to_cobra_conversion(antimony_path): + """ + Compares the stoichiometric matrices of original antimony file + and cobra-compatible sbml file to see if they are the same. + + input: path to original antimony file (non-cobra-compatible) + """ + ant_to_cobra(antimony_path) + + r = te.loada(antimony_path) + N_te = r.getFullStoichiometryMatrix() + + file_name = antimony_path.split(".")[0] + model = cobra.io.read_sbml_model(file_name + "_cobra.xml") + N_cobra = cobra.util.create_stoichiometric_matrix(model) + + assert N_te.shape == N_cobra.shape + assert (N_te == N_cobra).all() diff --git a/src/emll/test_models/unlabeled18.ant b/src/emll/test_models/unlabeled18.ant new file mode 100644 index 0000000..b9a5345 --- /dev/null +++ b/src/emll/test_models/unlabeled18.ant @@ -0,0 +1,148 @@ + +v2: $A -> B; e_v2*(Vm2 / Km2 *(A-B/Keq2)) / (1 + A/Km2 + B/Km3) +v3: B -> R; e_v3*(Vm3 / Km4 *(B-R/Keq3)) / (1 + B/Km4 + R/Km5) +v13: R -> ; e_v13*(Vm13 / Km22 *(R/Keq13)) / (1 + R/Km22) + +v4: B -> C; e_v4*(Vm4 / Km6 *(B-C/Keq4)) / (1 + B/Km6 + C/Km7) +v5: C -> D; e_v5*(Vm5 / Km7 *(C-D/Keq5)) / (1 + C/Km7 + D/Km8) +v6: D -> E; e_v6*(Vm6 / Km9 *(D-E/Keq6)) / (1 + D/Km9 + E/Km10) +v7: E ->; e_v7*(Vm7 / Km11 *(E/Keq7)) / (1 + E/Km11) +v8: D -> F; e_v8*(Vm8 / Km12 *(D-F/Keq8)) / (1 + D/Km12 + F/Km13) +v9: D -> G; e_v9*(Vm9 / Km14 *(D-G/Keq9)) / (1 + D/Km14 + G/Km15) +v10: G -> H; e_v10*(Vm10 / Km16 *(G-H/Keq10)) / (1 + G/Km16 + H/Km17) +v11: H -> I; e_v11*(Vm11 / Km18 *(H-I/Keq11)) / (1 + H/Km18 + I/Km19) +v12: I -> $J; e_v12*(Vm12 / Km20 *(I-J/Keq12)) / (1 + I/Km20 + J/Km21) + +v14: F -> ; e_v14*(Vm14 / Km23 *(F/Keq14)) / (1 + F/Km23) +v15: F -> L; e_v15*(Vm15 / Km24 *(F-L/Keq15)) / (1 + F/Km24 + L/Km25) +v16: L -> M; e_v16*(Vm16 / Km26 *(L-M/Keq16)) / (1 + L/Km26 + M/Km27) +v17: M -> N; e_v17*(Vm17 / Km28 *(M-N/Keq17)) / (1 + M/Km28 + N/Km29) +v18: N-> O; e_v18*(Vm18 / Km30 *(N-O/Keq18)) / (1 + N/Km30 + O/Km31) +v19: O-> $P; e_v19*(Vm19 / Km31 *(O-P/Keq19)) / (1 + O/Km31 + P/Km32) + +v21: O -> $Q; e_v21*(Vm21 / Km34 *(O-Q/Keq21)) / (1 + O/Km34 + Q/Km35) + +e_v1 = 1; +e_v2 = 1; +e_v3 = 1; +e_v4 = 1; +e_v5 = 1; +e_v6 = 1; +e_v7 = 1; +e_v8 = 1; +e_v9 = 1; +e_v10 = 1; +e_v11 = 1; +e_v12 = 1; +e_v13 = 1; +e_v14 = 1; +e_v15 = 1; +e_v16 = 1; +e_v17 = 1; +e_v18 = 1; +e_v19 = 1; +e_v20 = 1; +e_v21 = 1; +e_v22 = 1; + +Vm1 = 1.448; +Vm2 = 4.881; +Vm3 = 3.556; +Vm4 = 9.404; +Vm5 = 7.653; +Vm6 = 7.487; +Vm7 = 9.037; +Vm8 = 0.834; +Vm9 = 5.522; +Vm10 = 5.845; +Vm11 = 9.619; +Vm12 = 2.921; +Vm13 = 2.408; +Vm14 = 1.003; +Vm15 = 0.164; +Vm16 = 9.295; +Vm17 = 6.699; +Vm18 = 7.852; +Vm19 = 2.817; +Vm20 = 5.864; +Vm21 = 0.64; +Vm22 = 4.856; + +Keq1 = 0.6875; +Keq2 = 2.155; +Keq3 = 9.474; +Keq4 = 7.309; +Keq5 = 2.539; +Keq6 = 2.133; +Keq7 = 5.182; +Keq8 = 0.257; +Keq9 = 2.075; +Keq10 = 4.247; +Keq11 = 3.742; +Keq12 = 4.636; +Keq13 = 2.776; +Keq14 = 5.868; +Keq15 = 8.639; +Keq16 = 1.175; +Keq17 = 5.174; +Keq18 = 1.321; +Keq19 = 7.169; +Keq20 = 3.961; +Keq21 = 5.654; +Keq22 = 1.833; + +Km1 = 0.398; +Km2 = 0.21; +Km3 = 0.186; +Km4 = 0.944; +Km5 = 0.74; +Km6 = 0.49; +Km7 = 0.227; +Km8 = 0.254; +Km9 = 0.058; +Km10 = 0.434; +Km11 = 0.312; +Km12 = 0.696; +Km13 = 0.378; +Km14 = 0.18; +Km15 = 0.025; +Km16 = 0.067; +Km17 = 0.679; +Km18 = 0.454; +Km19 = 0.537; +Km20 = 0.897; +Km21 = 0.99; +Km22 = 0.217; +Km23 = 0.663; +Km24 = 0.263; +Km25 = 0.021; +Km26 = 0.758; +Km27 = 0.32; +Km28 = 0.383; +Km29 = 0.588; +Km30 = 0.831; +Km31 = 0.629; +Km32 = 0.873; +Km33 = 0.274; +Km34 = 0.798; +Km35 = 0.186; +Km36 = 0.953; + +A = 10; +B = 0.19; +C = 0.511; +D = 0.224; +E = 0.098; +F = 0.862; +G = 0.973; +H = 0.961; +I = 0.907; +J = 3.32; +# K = 0.333; +L = 0.081; +M = 0.407; +N = 0.232; +O = 0.132; +P = 0.53; +Q = 2.26; +R = 0.011; diff --git a/src/emll/util.py b/src/emll/util.py index b729de9..145a5af 100644 --- a/src/emll/util.py +++ b/src/emll/util.py @@ -2,6 +2,12 @@ import scipy as sp import pytensor.tensor as at import pymc as pm +from pathlib import Path + +import tellurium as te +import libsbml + +import os def create_elasticity_matrix(model): @@ -25,11 +31,15 @@ def create_elasticity_matrix(model): array[r_ind(reaction), m_ind(metabolite)] = -np.sign(stoich) # Irrevesible in forward direction, only assign if met is reactant - elif (not reaction.reversibility) & (reaction.upper_bound > 0) & (stoich < 0): + elif ( + (not reaction.reversibility) & (reaction.upper_bound > 0) & (stoich < 0) + ): array[r_ind(reaction), m_ind(metabolite)] = -np.sign(stoich) # Irreversible in reverse direction, only assign in met is product - elif (not reaction.reversibility) & (reaction.lower_bound < 0) & (stoich > 0): + elif ( + (not reaction.reversibility) & (reaction.lower_bound < 0) & (stoich > 0) + ): array[r_ind(reaction), m_ind(metabolite)] = -np.sign(stoich) return array @@ -43,7 +53,8 @@ def create_Ey_matrix(model): boundary_indexes = [model.reactions.index(r) for r in model.medium.keys()] boundary_directions = [ - 1 if r.products else -1 for r in model.reactions.query(lambda x: x.boundary, None) + 1 if r.products else -1 + for r in model.reactions.query(lambda x: x.boundary, None) ] ny = len(boundary_indexes) Ey = np.zeros((len(model.reactions), ny)) @@ -151,7 +162,9 @@ def initialize_elasticity( name = "ex" if m_compartments is not None: - assert r_compartments is not None, "reaction and metabolite compartments must both be given" + assert ( + r_compartments is not None + ), "reaction and metabolite compartments must both be given" regulation_array = np.array( [[a in b for a in m_compartments] for b in r_compartments] @@ -220,3 +233,111 @@ def initialize_elasticity( E = flat_e_entries[flat_indexer].reshape(N.T.shape) return E + + +def ant_to_cobra(antimony_path): + """ + This method takes in an Antimony file and converts it to a Cobra- + friendly format by removing all boundary species and replacing all + rate laws with a constant. Both an Antimony file and an SBML are + produced, and the file name is returned. + """ + output_name = antimony_path.split("/")[-1] + output_name = output_name.split(".")[0] + output_path = Path(antimony_path).parent + + # load normal Antimony file + with open(antimony_path, "r") as file: + # reprint the antimony model into a temp file to ensure proper + # headers + # load antimony file + r = te.loada(antimony_path) + # print antimony file to temp + with open("tempA7K8L2P4W9.txt", "w") as f: + f.write(r.getCurrentAntimony()) + old_ant = "tempA7K8L2P4W9.txt" + + with open(old_ant, "r") as file: + lines = file.readlines() + section_indices = [] + c_index = -1 + for i, line in enumerate(lines): + if "//" in line: + section_indices.append(i) + if "// Compartment" in line and c_index == -1: + c_index = i + next_section = section_indices.index(c_index) + 1 + + with open(old_ant, "r") as file: + lines = file.readlines()[c_index : section_indices[next_section]] + + with open(f"{output_path}/{output_name}_cobra.ant", "w") as f: + f.write("") + + for line in lines: + line = line.strip() + if "$" not in line: + with open(f"{output_path}/{output_name}_cobra.ant", "a") as f: + f.write(line + "\n") + else: + no_bd_sp = [i for i in line.split(",") if "$" not in i] + no_bd_sp = [i for i in no_bd_sp if i != ""] + + line = ",".join(no_bd_sp) + if line != "" and line[-1] != ";": + line += ";" + if "species" not in line: + line = "species" + line + with open(f"{output_path}/{output_name}_cobra.ant", "a") as f: + f.write(line + "\n") + + r = te.loada(antimony_path) + doc = libsbml.readSBMLFromString(r.getSBML()) + model = doc.getModel() + + bd_sp = r.getBoundarySpeciesIds() + + reactants_list = [] + products_list = [] + + for n in range(len(r.getReactionIds())): + rxn = model.getReaction(n) + reactants = [] + products = [] + + for reactant in range(rxn.getNumReactants()): + stoich = rxn.getReactant(reactant).getStoichiometry() + if stoich == 1: + reactants.append(rxn.getReactant(reactant).species) + else: + reactants.append(str(stoich) + " " + rxn.getReactant(reactant).species) + reactants_list.append([i for i in reactants if i not in bd_sp]) + + for product in range(rxn.getNumProducts()): + stoich = rxn.getProduct(product).getStoichiometry() + if stoich == 1: + products.append(rxn.getProduct(product).species) + else: + products.append(str(stoich) + " " + rxn.getProduct(product).species) + products_list.append([i for i in products if i not in bd_sp]) + + for i in range(len(reactants_list)): + r1 = " + ".join(reactants_list[i]) + p1 = " + ".join(products_list[i]) + with open(f"{output_path}/{output_name}_cobra.ant", "a") as f: + f.write(r.getReactionIds()[i] + ": " + r1 + " -> " + p1 + "; 1;\n") + + with open(f"{output_path}/{output_name}_cobra.ant", "a") as f: + f.write("\n") + + for sp in r.getFloatingSpeciesIds(): + with open(f"{output_path}/{output_name}_cobra.ant", "a") as f: + f.write(sp + " = 1;\n") + + with open(f"{output_path}/{output_name}_cobra.xml", "w") as f: + f.write(te.loada(f"{output_path}/{output_name}_cobra.ant").getCurrentSBML()) + + if "tempA7K8L2P4W9.txt" in os.listdir(): + os.remove("tempA7K8L2P4W9.txt") + + return f"{output_path}/{output_name}_cobra"