Zero flux for modelSeed generated model

[marti...@isomerase.co.uk]

Hello everyone!

I'm trying to get to grips with some basics of producing my own genome scale model based on our own NGS assemblies. Being new to this at the moment I'm just following some prior examples of people doing something similar to what I need to do, and it seems the easiest way for a complete newbie such as myself is to upload my fasta sequence to modelSEED and produce a model there, download that and convert to COBRA compatible format with mackinac and then proceed ahead in COBRApy. I've found an example of this with a 2017 iGEM team project who followed these steps to generate a model of Paracoccus denitrificans strain DSM 413 and I've been using their script on their wiki (http://2017.igem.org/Team:Virginia/Model).

So my steps are to login to modelSEED and find the reference 266.10 assembly as they used, create a model from that, pull down the 266.10 model to my computer with the script below and run part of the script up to outputting an unmodified objective growth value.

I get an output as shown below:

1021 reactions initially
1022 metabolites initially
883 genes initially
Unmodified objective value: 0.000000000000000

Now the team's wiki says they got an output of 1550 reactions and 1556 metabolites so clearly something is different there so maybe I missed a vital step, but they claim to get an output of 224.324 and clearly I am getting just zero. So... what am I doing wrong at this very basic stage? I feel like I understand most of the basics of COBRA and adding/deleting reactions etc with a reference mocel such as iJO1366.json, but ultimately as I say I need to be able to model our own organisms and not just references.

My code is below and should be identical to what the iGEM team have. I've only modified to remove my modelSEED login details. Can anyone see anything wrong? Do you get the same number of reactions and a zero output or do you get something sensible?

Thanks so much for any help you can provide!!! :) :)

from __future__ import print_function
import mackinac
import cobra
from cobra import Model, Reaction, Metabolite
import os
from os.path import join

def boundary_list(model):
    """Creates and returns the list of boundary reactions (exchange, demand and sink-type reactions)
    """
    exList = []
    for rxn in model.reactions:
    if (rxn.boundary == True):
        exList.append(rxn)
    return exList

"""Defining parameters"""

modelseed_id = '266.10'
first_run = False
validate = False
perform_fva = True
export_to_escher = True
block_nitrite_EX = True
precision_fva = 0.999
rxns_KO = ['EX_cpd00644_e', 'EX_cpd10516_e', 'EX_cpd00060_e', 'rxn_cpd14010_c', 'rxn_cpd00109_c', 'rxn_cpd00570_c', 'rxn_cpd14202_c', 'rxn_cpd06111_c', 'rxn_cpd05625_c', 'EX_cpd00528_e']
simult_KO = ['rxn0625_c', 'EX_cpd00528_e']

data_dir = os.getcwd()
print("Script directory: ", data_dir)
 
token_username = ''
token_password = ''

if(first_run):
    mackinac.get_token(token_username, token_password, token_type = 'rast')
    mackinac.reconstruct_modelseed_model(modelseed_id)
    mackinac.gapfill_modelseed_model(modelseed_id)
    mackinac.optimize_modelseed_model(modelseed_id)

    MODEL = mackinac.create_cobra_model_from_modelseed_model(modelseed_id, validate = True)
    cobra.io.save_json_model(MODEL, "initMODEL.json")

if(not(first_run)):
    MODEL = cobra.io.load_json_model(join(data_dir, "initMODEL.json"))
    print("Model loaded from an existing file: ", join(data_dir, "initMODEL.json\n"))

if(validate):
    print("---- MODEL VALIDATION BEGINS ----")
    print(cobra.manipulation.validate.check_mass_balance(MODEL))
    print(cobra.manipulation.validate.check_metabolite_compartment_formula(MODEL))
    print(cobra.manipulation.validate.check_reaction_bounds(MODEL))
    print("---- END OF VALIDATION ----")

#SO NOW THE MODEL IS IMPORTED AS IS FROM MODELSEED AND CONVERTED TO COBRA MODEL

#ADD NEW HETEROLOGOUS REACTIONS
#nitric-oxide: ferricytomchrome-c oxireductase (EC 1.7.2.1)
#NO2 + H+ + cyt c2 <=> NO + H2O + cyt c3+
#cpd00075_c + cpd00067_c + cpd00110_c <=> cpd00418_c + cpd00001_c + cpd00109_c
cpd00418_c = MODEL.metabolites.get_by_id('cpd00418_c')
cpd00110_c = MODEL.metabolites.get_by_id('cpd00110_c')
cpd00109_c = MODEL.metabolites.get_by_id('cpd00109_c')
cpd00075_c = MODEL.metabolites.get_by_id('cpd00075_c')
cpd00067_c = MODEL.metabolites.get_by_id('cpd00067_c')
cpd00001_c = MODEL.metabolites.get_by_id('cpd00001_c')
cpd00659_c = MODEL.metabolites.get_by_id('cpd00659_c')
cpd00075_e = MODEL.metabolites.get_by_id('cpd00075_e')


rxn05889_c = Reaction('rxn05889_c')
rxn05889_c.name = 'nitrite reductase (NO-forming)'
rxn05889_c.subsystem = 'Denitrification'
rxn05889_c.lower_bound = -1000.
rxn05889_c.upper_bound = 1000.
rxn05889_c.gene_reaction_rule = "( 266.10.peg.2397 )"

rxn05889_c.add_metabolites({
    cpd00075_c: -1.0,
    cpd00067_c: -1.0,
    cpd00110_c: -1.0,
    cpd00418_c: 1.0,
    cpd00001_c: 1.0,
    cpd00109_c: 1.0
})

cpd00528_c = MODEL.metabolites.get_by_id('cpd00528_c')

cpd00528_e = Metabolite(
'cpd00528_e',
formula='N2',
name='N2_e',
compartment='e')

rxn10577_c = Reaction('rxn10577_c')
rxn10577_c.name = 'Nitrogen exchange, diffusion'
rxn10577_c.subsystem = 'Transport'
rxn10577_c.lower_bound = -1000.
rxn10577_c.upper_bound = 1000.
rxn10577_c.add_metabolites({
    cpd00528_c: -1.0,
    cpd00528_e: 1.0
})

#add the exchange reactions got N2_e
MODEL.add_boundary(cpd00528_e, type="exchange", reaction_id="EX_cpd00528_e", ub=1000.)
MODEL.add_reactions([rxn05889_c, rxn10577_c])

#Excess NH3 inhibits NH-3 forming Nitrite reductase
MODEL.reactions.get_by_id('rxn00568_c').knock_out()
MODEL.reactions.get_by_id('rxn00569_c').knock_out()

#Block NO2- exchange, block one of the loop
#if(block_nitrite_EX):
 #   MODEL.reactions.get_by_id('cpd00075_e').knock_out()

#nitrate and nitrite reductatse only catalyze the forward reation in excess of ammonia
MODEL.reactions.get_by_id('rxn10121_c').bounds = (-10.,1000.)

print('')
print(len(MODEL.reactions), "reactions initially")
print(len(MODEL.metabolites), "metabolites initially")
print(len(MODEL.genes), "genes initially")

#Locally execute FBA on the model.
sol_1 = MODEL.optimize()
print('Unmodified objective value: %3.15f\n' % sol_1.objective_value)

#Performs FVA on the model. The result is a pandas.data_frame that contains the minimum and maximum flux range for exchnage reactions within 0.1% of the optimal solution
rxnList=boundary_list(MODEL)

print(sol_1)

[Moritz Beber]

Hi Martin,

Two general questions for you:

1. Since you have notes from another team performing this task, is there anyone on-site who could guide you through the process using their experience? (It's always easier to do these things in person.)
2. Have you explored KBase or carveme for the reconstruction workflow? I'm not saying their better but there are alternatives to modelSEED. (KBase maybe being a more visual and reproducible version of modelSEED.)

Best,

Moritz

--
You received this message because you are subscribed to the Google Groups "cobra pie" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cobra-pie+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[marti...@isomerase.co.uk]

Hi Moritz,

Thanks so much for your reply!

I'm working on this type of work alone at the moment. I work for a company and with the modelling/FBA I'm trying to train myself up to offer them new skills since we are beginning to carry out more NGS of organisms. Nobody else here is doing metabolic modelling so I'm rather on my own unfortunately... I haven't reached out to the Virginia iGEM team who produced the script and model since I am just using their example as it had the steps I was interested in trying to eventually replicate for my own purposes. From my own previous experience with iGEM in 2016 I know that once the jamboree ends the team disbands and the students return to their studies or go on to other things so as this project was from 2017 I wouldn't be hopeful of getting much anything back from them now. However, I can always try! I was more just wondering if anyone else here would be able to just replicate the steps and see if they also got a zero output to confirm to me it wasn't a problem at my end of things :)

As for using KBase well I'm up for anything! I've used modelSEED purely because I read a paper that said the type of bacteria we work on are annotated well by it, and as I want to keep everything as simple as possible - and I have no clue right at the moment of how to annotate and build the model from scratch from just the assembled genome - I thought it best to use an automated service to do this for me.

[Moritz Beber]

Hello,

I see KBase as somewhat of a continuation of ModelSEED since both concepts come from Chris Henry and there is large overlap in the teams. It takes a moment to create an account since it's US government infrastructure but take a look at this public narrative. Maybe you find it easier to generate a model this way.

https://narrative.kbase.us/#appcatalog/app/fba_tools/build_metabolic_model

--