# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
# ---

# +
# -*- coding: utf-8 -*-
import sys
import numpy as np
import pandas as pd
import cobra

print('Python version:', sys.version)
print('numpy version:', np.__version__)
print('pandas version:', pd.__version__)
print('cobrapy version:', cobra.__version__)


# -

def AddRxn(model, newRxnFile):
    """Function of adding new reactions to the model."""
    n1 = len(model.reactions)
    AllAddRxn = pd.read_csv(newRxnFile, sep=',', index_col='RxnID', skipinitialspace=True)
    n2 = len(AllAddRxn)
    for i in range(n2):
        ID = AllAddRxn.index.values[i]
        addRxn = cobra.Reaction(ID)
        model.add_reactions([addRxn])
        addRxnInf = model.reactions[n1 + i]
        addRxnInf.name = AllAddRxn.loc[ID, 'RxnName']
        addRxnInf.reaction = AllAddRxn.loc[ID, 'RxnFormula']
        addRxnInf.subsystem = AllAddRxn.loc[ID, 'Subsystem']
        addRxnInf.lower_bound = AllAddRxn.loc[ID, 'LowerBound']
        addRxnInf.upper_bound = AllAddRxn.loc[ID, 'UpperBound']
    return model



def KORxn(model, rxns2KO):
    """Function for knocking out reactions."""
    for ID in rxns2KO:
        reaction = model.reactions.get_by_id(ID)
        reaction.knock_out()



def flux2file(model, product, psw, save=False):
    """Function of exporting flux data."""
    if save:
        n = len(model.reactions)
        modelMatrix = np.empty([n, 9], dtype = object)
        for i in range(len(model.reactions)):
            x = model.reactions[i]
            modelMatrix[i, 0] = i + 1
            modelMatrix[i, 1] = x.id
            modelMatrix[i, 2] = x.name
            modelMatrix[i, 3] = x.reaction
            modelMatrix[i, 4] = x.subsystem
            modelMatrix[i, 5] = x.lower_bound
            modelMatrix[i, 6] = x.upper_bound
            modelMatrix[i, 7] = x.flux
            modelMatrix[i, 8] = abs(x.flux)
            
        df = pd.DataFrame(data = modelMatrix, 
                          columns = ['N', 'RxnID', 'RxnName', 'Reaction', 'SubSystem', 
                          'LowerBound', 'UpperBound', 'Flux-core', 'abs(Flux)'])
        df.to_excel('{}_{}.xlsx'.format(product, psw), index=False)



def change_objective(model, product):
    """change objective function"""
    if product in ['accoa_c', 'succoa_c']:
        rxn = cobra.Reaction(id=f"EX_{product}", lower_bound=-1000, upper_bound=1000)
        model.add_reaction(rxn)
        rxn.add_metabolites({product: -1, 'coa_c': 1})
    elif product in ['g6p_c', 'f6p_c', 'e4p_c', 'r5p_c', 'g3p_c', '3pg_c', 'pep_c']:
        rxn = cobra.Reaction(id=f"EX_{product}", lower_bound=-1000, upper_bound=1000)
        model.add_reaction(rxn)
        rxn.add_metabolites({product: -1, 'pi_c': 1})
    elif f"EX_{product}" not in model.reactions:
        rxn = cobra.Reaction(id=f"EX_{product}", lower_bound=-1000, upper_bound=1000)
        model.add_reaction(rxn)
        rxn.add_metabolites({product: -1})
    model.objective = f"EX_{product}"



def prodFBA(model, psw, product, save=False):
    solution = model.optimize()
    # print("%s: " %product, round(solution.objective_value, 3))
    # export results
    flux2file(model, product, psw, save)
    return round(solution.objective_value, 3)
    


# ## Model background
#
#  * Using E. coli full model iML1515.
#  * Modified the transhydrogenase reaction (THD2pp) from 2 protons translocation to 1.
#  * Changed HSDy (homoserine DH) to be irrversible, towards homoserine.
#  * Changed MMM (Methylmalonyl-CoA mutase) to be reversible.
#  * Base deletions: PFL, OBTFL, POR5 (pyruvate synthase, Ferredoxin), FDH4pp, FDH5pp, GLYCK (glycerate kinase, to 3pg), ATPM.
#  * FolD is reversible, MTHFC and MTHFD.
#  * FTL reaction is existed, FTHFLi (irrversible)
#  * The model has part of the Ethylmalonyl-CoA pathway, ACACT1r, HACD1 and ECOAH1,
#  propanoyl-CoA degradation, MMM and MMCD already.

# +
model = cobra.io.load_json_model(r'..\0_ecoli_models\iML1515.json')

model.reactions.THD2pp.add_metabolites({"h_p": 1, "h_c": -1})
model.reactions.HSDy.bounds = (-1000, 0)
model.reactions.MMM.bounds = (-1000, 1000)
# -

# ## base deletions:
#
# * POR5, GLYCK, FDH4pp, FDH5pp;
# * ATPM;
# * aerobic condition: PFL, OBTFL;
# * GART makes confusion, since it has alternative way and the deletion of it does not 
#  change the FBA result.
# * in the case of AcCoA production (maybe also acetate) of RuMP, 
# most product is from r5p --> acald --> ; DRPA is in the pathway.
# And deletion of it does not change the biomass result of WT on glc.
# * The same for PAI2T.

# +
KORxn_base = ['POR5', 'GLYCK', 'FDH4pp', 'FDH5pp', 'ATPM', 
              'PFL', 'OBTFL', 'GART', 'DRPA', 'PAI2T']

KORxn(model, KORxn_base)

# +
# Set medium: MeOH as sole the carbon source.
for rxn in model.boundary:
    if "C" in rxn.check_mass_balance(): # knock out all other carbon-related transporters
        rxn.bounds = (0.0, 1000)
model.reactions.get_by_id('EX_meoh_e').bounds = (-10, 1000)
model.reactions.get_by_id('EX_co2_e').bounds = (-1000, 1000)

model.reactions.get_by_id('EX_h2s_e').bounds = (0.0, 0.0) # avoid using S as electron accepeter
model.reactions.get_by_id('EX_fe3_e').bounds = (0.0, 0.0) # avoid using Fe2 as electron donor

# -

psws = {'RuMP': 'NewRxns4Full_RuMP.csv',
        'SerCyc': 'NewRxns4Full_SerCyc.csv',
        'MserCyc': 'NewRxns4Full_MSerCyc.csv',
        'HomSer': 'NewRxns4Full_HomSer.csv',}
# Serine cycle has ethylmalonyl-coa pathway reactions (EtMaCoA)


# ### Try first with the 12 precursors (DOI: 10.1016/j.cell.2016.02.004) + citrate

precursors = ['g6p_c', 'f6p_c', 'e4p_c', 'r5p_c', 'g3p_c', 
              '3pg_c', 'pep_c', 'pyr_c', 'accoa_c', 'akg_c', 
              'succoa_c', 'oaa_c', 'cit_c']

# +
# saveif = True
saveif = False
precursors_results = pd.DataFrame()

for psw, pswrxns in psws.items():
    with model:
        model = AddRxn(model, pswrxns)
        if psw == 'RuMP':
            AddKO = ['FTHFLi', 'THRA', 'THRD',]
        elif psw == 'HomSer':
            AddKO = ['FTHFLi', ]
        elif psw == 'SerCyc':
            AddKO = ['THRA', 'THRD', 'EtMaCoA']  # no EtMaCoA to be fair
        elif psw == 'MserCyc':
            AddKO = ['THRA', 'THRD',]
        KORxn(model, AddKO)
        for precursor in precursors:
            with model:
                change_objective(model, precursor)
                sol = prodFBA(model, psw, precursor, saveif)
                precursors_results.loc[precursor, psw] = sol

precursors_results
# -

import matplotlib.pyplot as plt 
plt.rc('axes', axisbelow=True)
print('matplotlib version:', matplotlib.__version__)

# +
x = np.arange(len(precursors))
witdh = 0.2

fig, ax = plt.subplots(figsize=(16, 6))
ax.set_xticks(x)
ax.set_xlim(-0.5, 12.5)
ax.set_xticklabels(precursors)
ax.set_yticklabels(map(str, np.arange(0, 7, 1)/10))
ax.set_ylabel('Maximal yeild (mol/mol-MeOH)')
plt.grid(which='both', axis='y')
ax.bar(x-1.5*witdh, precursors_results['RuMP'], witdh, label='RuMP')
ax.bar(x-0.5*witdh, precursors_results['SerCyc'], witdh, label = 'SerCyc')
ax.bar(x+0.5*witdh, precursors_results['MserCyc'], witdh, label = 'MserCyc')
ax.bar(x+1.5*witdh, precursors_results['HomSer'], witdh, label = 'HomSer')
plt.legend()
plt.savefig('13 precursors production figure.eps', dpi=300, bbox_inches='tight')
plt.show()
# -

# ### Short summary
# * The modified serine cycle has always the lowest yield.
# * RuMP, serine cycle and homoserine cycle have their own advantage over others for certain products.
# * When producing AcCoA and citrate, the homoserine cycle has the highest yield than the other pathways. 

# ### To produce some real chemicals. 
#

products = ['ac_e', 'etoh_e', 'acetone_e', 'but_e', 'buoh_e',
            'pyr_e', 'lac__D_e', 'ibtol_e', '23btdl_e',
            'cit_e', 'itacon_e', 'glu__L_e', 'asp__L_e', 'succ_e',
            'adpac_e', 'levu_e']


# ## Add new reactions for those pathways do not exist natively
# * acetone: DOI: 10.1016/j.ymben.2012.08.001 
# * butyrate: but_e; DOI: 10.1002/bit.24925
# * 1-butanol: buoh_e; DOI: 10.1186/1475-2859-7-36
# * isobutanol: ibtol_e; DOI: 10.1007/s00253-009-2085-6
# * 2,3-butanediol: 23btdl_e; DOI: 10.1021/acssynbio.5b00057
# * itaconate: itacon_e; DOI: 10.1007/s10529-017-2500-5
# * adipic acid: adpac_e; DOI: 10.1002/bit.25293
# * levulinic acid: levu_e; DOI: 10.1038/nbt.3505
#

model = AddRxn(model, 'NewRxn4Production.csv')

# +
# saveif = True
saveif = False
prod_results = pd.DataFrame()

for psw, pswrxns in psws.items():
    with model:
        model = AddRxn(model, pswrxns)
        # print('Pathway %s' % psw)
        if psw == 'RuMP':
            AddKO = ['FTHFLi', 'THRA', 'THRD',]
        elif psw == 'HomSer':
            AddKO = ['FTHFLi', ]
        elif psw == 'SerCyc':
            AddKO = ['THRA', 'THRD', 'EtMaCoA']  # no EtMaCoA to be fair
        elif psw == 'MserCyc':
            AddKO = ['THRA', 'THRD',]
        KORxn(model, AddKO)
        for product in products:
            with model:
                # print('Product: %s' % product)
                change_objective(model, product)
                sol = prodFBA(model, psw, product, saveif)
                prod_results.loc[product, psw] = sol

# -

prod_results

# +
x = np.arange(len(products))
witdh = 0.2

fig, ax = plt.subplots(figsize=(16, 6))
ax.set_xticks(x)
ax.set_xlim(-0.5, 15.5)
ax.set_xticklabels(products)
ax.set_yticklabels(map(str, np.arange(0, 7, 1)/10))
ax.set_ylabel('Maximal yeild (mol/mol-MeOH)')
plt.grid(which='both', axis='y')
ax.bar(x-1.5*witdh, prod_results['RuMP'], witdh, label='RuMP')
ax.bar(x-0.5*witdh, prod_results['SerCyc'], witdh, label = 'SerCyc')
ax.bar(x+0.5*witdh, prod_results['MserCyc'], witdh, label = 'MserCyc')
ax.bar(x+1.5*witdh, prod_results['HomSer'], witdh, label = 'HomSer')
plt.legend()
plt.savefig('chemical production from formladehyde figure.eps', dpi=300, bbox_inches='tight')
plt.show()

# -

# add NADH oxidation reaction to check overflow reducing power
rxn = cobra.Reaction(id="NADHOx", name="NADH oxidation", lower_bound=0, upper_bound=1000)
model.add_reaction(rxn)
rxn.add_metabolites({"nadh_c": -1, "o2_c": -0.5, "h_c": -1, "nad_c": 1, "h2o_c": 1})

# +
prod_results2 = pd.DataFrame()
dissipation = pd.DataFrame()

for psw, pswrxns in psws.items():
    with model:
        model = AddRxn(model, pswrxns)
        # print('Pathway %s' % psw)
        if psw == 'RuMP':
            AddKO = ['FTHFLi', 'THRA', 'THRD',]
        elif psw == 'HomSer':
            AddKO = ['FTHFLi', ]
        elif psw == 'SerCyc':
            AddKO = ['THRA', 'THRD', 'EtMaCoA']  # no EtMaCoA to be fair
        elif psw == 'MserCyc':
            AddKO = ['THRA', 'THRD',]
        KORxn(model, AddKO)
        for product in products:
            with model:
                # print('Product: %s' % product)
                change_objective(model, product)
                # sol = prodFBA(model, psw, product, False)
                sol = model.optimize().objective_value
                prod_results2.loc[product, psw] = round(sol,3)
                flux2file(model, f"01_{product}_dis", psw, False)

                model.reactions.get_by_id(f"EX_{product}").bounds = (sol, sol)
                model.objective = 'NADHOx'
                sol2 = model.optimize()
                dissipation.loc[product, psw] = round(sol2.objective_value, 3)
                flux2file(model, f"02_{product}_dis", psw, False)
# -

prod_results2

dissipation

# * RuMP has overflow of electrons. 
# * The homoserine cycle uses PPC fixing CO2 to be a sink for electrons (reducing power).

yield2plot = pd.concat([precursors_results.loc[['g6p_c', 'r5p_c','g3p_c','3pg_c'],['RuMP','SerCyc','HomSer']],
                        prod_results.loc[['pyr_e','lac__D_e','ibtol_e','23btdl_e'],['RuMP','SerCyc','HomSer']],
                        precursors_results.loc[['accoa_c'],['RuMP','SerCyc','HomSer']],
                        prod_results.loc[['etoh_e','acetone_e','but_e','buoh_e','cit_e','itacon_e'],['RuMP','SerCyc','HomSer']],
                        precursors_results.loc[['akg_c'],['RuMP','SerCyc','HomSer']],
                        prod_results.loc[['succ_e','levu_e'],['RuMP','SerCyc','HomSer']]
])
yield2plot

# yield of RuMP as 100%
yield_norm = yield2plot.div(yield2plot.RuMP, axis=0) * 100
round(yield_norm,1)

# +
x = np.arange(len(yield_norm))
witdh = 0.25

fig, ax = plt.subplots(figsize=(18, 6))
ax.set_xticks(x)
ax.set_xlim(-0.6, 17.6)
ax.set_xticklabels(yield_norm.index.to_numpy())
# ax.set_ylim(0, 185)
ax.set_ylabel('')
plt.grid(which='both', axis='y')
ax.bar(x-1*witdh, yield_norm['RuMP'], witdh, label='RuMP')
ax.bar(x, yield_norm['SerCyc'], witdh, label = 'SerCyc')
ax.bar(x+1*witdh, yield_norm['HomSer'], witdh, label = 'HomSer')
plt.legend()
plt.savefig('chemical production comparison figure.eps', dpi=300, bbox_inches='tight')
plt.show()
# -



StarcoderdataJupyterScriptsDedupFiltered