Analyzing Features¶

import impact as impt
# import cobra
# import cobra.test
# import cobra.io
import numpy as np

# import matplotlib.pyplot as plt
% matplotlib inline

from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
from plotly.graph_objs import Bar, Layout, Figure, Scatter
init_notebook_mode()

# We include this to ensure the js is loaded when viewed online
from IPython.display import HTML
HTML('<script src="https://cdn.plot.ly/plotly-latest.min.js"></script>')

# Let's grab the iJO1366 E. coli model, cobra's test module has a copy
model = cobra.test.create_test_model("ecoli")

# Simulate anaerobic conditions by prevent oxygen uptake
model.reactions.get_by_id('EX_o2_e').knock_out()

# Optimize the model
cobra.flux_analysis.optimize_minimal_flux(model)
model.summary()

# Let's consider one substrate and five products
biomass_keys = ['Ec_biomass_iJO1366_core_53p95M']
substrate_keys = ['EX_glc_e']
product_keys = ['EX_for_e','EX_ac_e','EX_etoh_e','EX_succ_e']
analyte_keys = biomass_keys+substrate_keys+product_keys
# We'll use numpy to generate an arbitrary time vector

# The initial conditions (mM) [biomass, substrate,
#                              product1, product2, ..., product_n]
y0 = [0.05, 100, 0, 0, 0, 0]
t = np.linspace(0,16,8)

# Returns a dictionary of the profiles
from impact.helpers.synthetic_data import generate_data
dFBA_profiles = generate_data(y0, t, model,
                                      biomass_keys, substrate_keys,
                                      product_keys, plot = True)

# We can query the rates instead.
for keys in [biomass_keys,substrate_keys,product_keys]:
#     plt.figure()
    data = []
    for analyte in keys:
        trace = Scatter(x=t,
            y=single_trial.analyte_dict[analyte].gradient,
            name = analyte)
        data.append(trace)
    fig = Figure(data=data,
                 layout=Layout(title="Calculated rates",
                               yaxis = {'title':'Production rate [mmol /  h]'})
                 )
    iplot(fig)
#         plt.plot(t,single_trial.analyte_dict[analyte].gradient)
#     plt.legend(keys)

exchange_keys = biomass_keys + substrate_keys + product_keys

median_specific_productivity = {}
model_specific_productivity = {}
for analyte in exchange_keys:
    model_specific_productivity[analyte] = model.solution.x_dict[analyte]
    median_specific_productivity[analyte] = np.median(single_trial.analyte_dict[analyte].specific_productivity.data)

trace = Bar(x=[analyte for analyte in median_specific_productivity],
            y=[median_specific_productivity[analyte] for analyte in median_specific_productivity],
           name = 'Calculated')
trace2 = Bar(x=[analyte for analyte in model_specific_productivity],
            y=[model_specific_productivity[analyte] for analyte in model_specific_productivity],
            name = 'Model')
data = [trace, trace2]
fig = Figure(data=data,
             layout=Layout(title="Calculated vs Model Specific Productivity",
                           yaxis = {'title':'Specific Productivity [mmol / (gdw h)]'})
             )
iplot(fig)

# Let's take the model and add bounds for the known reactions
for analyte in analyte_keys:
    with_noise = np.median(single_trial.analyte_dict[analyte].specific_productivity.data)
    print(analyte,' ',with_noise)
    model.reactions.get_by_id(analyte).lower_bound = with_noise
    model.reactions.get_by_id(analyte).upper_bound = with_noise
print('\n\n\n')
solution = model.optimize()
model.summary()

# Convert mmol/gdw - hr to g / gdw - hr

# Create a dictionary for the mapping of flux reactions to metabolite names to get the molar mass
metabolite_export_reaction_dict = {
    'EX_for_e':'for_e',
    'EX_ac_e':'ac_e',
    'EX_etoh_e':'etoh_e',
    'EX_glc_e':'glc__D_e',
    'EX_succ_e':'succ_e'
}

mass_flux = {}
for analyte in product_keys+substrate_keys:
    mass_flux[analyte] = model.solution.x_dict[analyte] \
    * model.metabolites.get_by_id(metabolite_export_reaction_dict[analyte]).formula_weight
mass_flux[biomass_keys[0]] = model.solution.x_dict[biomass_keys[0]]*1000
print(mass_flux)

# The balance is the sum of all metabolites, the uptake is already negative.
balance = sum(mass_flux[metabolite] for metabolite in mass_flux)
print(balance)
# The closure is the total substrate accounted for
percent_closure = balance/mass_flux['EX_glc_e']

print('\nThe mass balance is %f%% closed' % ((1-percent_closure)*100))
print(mass_flux)

# Two situations
# 1: balance is > 100% (must be some unaccounted substrate)
labels = [analyte.split('_')[1] for analyte in mass_flux]
print(percent_closure)
if percent_closure < 0:
    fig = {
    'data': [{'labels': labels + ['Unaccounted substrate'],
              'values': [mass_flux[metabolite] for metabolite in mass_flux]+[balance],
              'type': 'pie'}],
    'layout': {'title': 'Mass balance'}
     }
# 2: balance is < 100% (must be some unaccounted product)
else:
    fig = {
    'data': [{'labels': labels + ['Missing product'],
              'values': [mass_flux[metabolite] for metabolite in mass_flux]+[abs(balance)],
              'type': 'pie'}],
    'layout': {'title': 'Mass balance'}
         }
iplot(fig)

analyte_keys = substrate_keys + biomass_keys + product_keys

# Calculate the waterfall for the products and biomass
substrate_consuming_keys = biomass_keys+product_keys
running_total = balance
base = []
for i, key in enumerate(analyte_keys):
    if i == 0:
        running_total += -mass_flux[key]
        base.append(0)
    else:
        base.append(running_total)
        running_total += -mass_flux[key]


base_trace = Bar(x=analyte_keys,
                 y=base,
                 marker=dict(color='rgba(1,1,1, 0.0)',))

substrate_trace = Bar(x=[substrate_keys[0]],
                      y=[-mass_flux[substrate_keys[0]]],
                        marker=dict(color='rgba(55, 128, 191, 0.7)',
                                    line=dict(color='rgba(55, 128, 191, 1.0)',
                                              width=2,)))


biomass_trace = Bar(x=[biomass_keys[0]],
                    y=[-mass_flux[biomass_keys[0]]],
                       marker=dict(color='rgba(50, 171, 96, 0.7)',
                                   line=dict(color='rgba(50, 171, 96, 1.0)',
                                             width=2,)))

products_trace = Bar(x=product_keys,
                     y=[-mass_flux[metabolite] for metabolite in product_keys],
                      marker=dict(color='rgba(219, 64, 82, 0.7)',
                                  line=dict(color='rgba(219, 64, 82, 1.0)',
                                            width=2,)),
                    )

balance_base_trace = Bar(x=['Balance'],
                 y=[base[1]-balance],
                 marker=dict(color='rgba(1,1,1, 0.0)',))

balance_trace = Bar(x=['Balance'],
                     y=[balance],
                      marker=dict(color='rgba(0, 0, 0, 0.7)',
                                  line=dict(color='rgba(0, 0, 0, 1.0)',
                                            width=2,)),
                    )

data = [substrate_trace, balance_base_trace, base_trace, balance_trace, biomass_trace,products_trace]
layout = Layout(barmode='stack',
#                 paper_bgcolor='rgba(245, 246, 249, 1)',
#                 plot_bgcolor='rgba(245, 246, 249, 1)',
                yaxis={'title': 'Exchange rate (mmol/(gdw h))'},
                showlegend=False
            )
fig = Figure(data=data,layout=layout)
iplot(fig)

b = escher.Builder(map_name='iJO1366.Central metabolism',
                   reaction_data=solution.x_dict,
                   # color and size according to the absolute value
                   reaction_styles=['color', 'size', 'abs', 'text'],
                   # change the default colors
                   reaction_scale=[{'type': 'min', 'color': '#cccccc', 'size': 4},
                                   {'type': 'mean', 'color': '#0000dd', 'size': 20},
                                   {'type': 'max', 'color': '#ff0000', 'size': 40}],
                   # only show the primary metabolites
                   hide_secondary_metabolites=True)
b.display_in_notebook()