from cobra.flux_analysis.double_deletion import double_gene_deletion_fba, double_reaction_deletion_fba
from cobra.io import load_matlab_model, save_matlab_model
import pandas as pd

import os
os.chdir('/media/users/jmonk/helper/models/')

# Load the model
model = load_matlab_model('iJO1366.mat')

# Model is pre-defined with M9 media bounds - check its predicted growth rate:
model.optimize()

LBMEDIA = {}
rxn_exchange_dict = {}

LBMEDIA['EX_glc_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_ala_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_arg_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_asp_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_cys_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_glu_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_gly_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_his_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_ile_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_leu_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_lys_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_met_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_pro_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_thr_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_tyr_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_phe_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_ser_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_trp_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_val_DASH_L_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_pnto_DASH_R_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_nac_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_na1_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_cl_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_so4_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_k_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_pi_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_ca2_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_mg2_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_zn2_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_aso3_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_cd2_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_hg2_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_h_LPAREN_e_RPAREN_'] = -100
LBMEDIA['EX_h2o_LPAREN_e_RPAREN_'] = -100
LBMEDIA['EX_o2_LPAREN_e_RPAREN_'] = -18.5
LBMEDIA['EX_ins_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_hxan_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_dcyt_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_thymd_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_ura_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_uri_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_dad_DASH_2_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_adn_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_co2_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_cobalt2_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_cu2_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_fe2_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_fe3_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_mn2_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_mobd_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_tungs_LPAREN_e_RPAREN_'] = -1000
LBMEDIA['EX_cbl1_LPAREN_e_RPAREN_'] = -0.01
LBMEDIA['EX_fru_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_gal_LPAREN_e_RPAREN_'] = -5
LBMEDIA['EX_ni2_LPAREN_e_RPAREN_'] = -1000.0
LBMEDIA['EX_sel_LPAREN_e_RPAREN_'] = -1000.0
LBMEDIA['EX_slnt_LPAREN_e_RPAREN_'] = -1000.0

for r in model.reactions:
    if r.startswith('EX_'):
        r.lower_bound = 0
    if r.id in LBMEDIA.keys():
        r.lower_bound = LBMEDIA[r.id]

# Check model predicted growth rate on LB media
model.optimize()

dko_genes = double_gene_deletion_fba(model, return_frame=True)

dko_genes.head()

# Deleting RXNs (instead of genes)

dko_reactions = double_reaction_deletion_fba(model, return_frame=True)

dko_reactions.head()

def run_all(dkos, org):
    # get singly essential and synthetic lethal genes:
    essential = []
    syn_lethal = []
    for i in dkos.index:
        if sum(dkos.ix[i]) < 0.01:
            essential.append(i)
        else:
            for g in dkos.ix[i].index:
                if dkos[g].sum() > 0.01:
                    if dkos.ix[i][g] < 0.01:
                        syn_lethal.append(str(i)+','+str(g))
    print len(essential), 'essential genes'
    print len(syn_lethal), 'syn lethal genes'
    
    return {'syn':syn_lethal,'essential':essential}

eco = run_all(dko_genes, 'eco')

#print essential genes
eco['essential']

# print synthetic lethal genes
eco['syn']