Applying Shlomi et al's biomarker prediction method to PKU on RECON2

**Authors**: Thierry D.G.A Mondeel, Stefania Astrologo, Ewelina Weglarz-Tomczak & Hans V. Westerhoff
University of Amsterdam
2017

The original publication looked at RECON2's predecessor RECON1. We will reproduce their analysis on RECON2 instead.



In [ ]:

    
import cobra
from utils import findBiomarkers
import pandas as pd

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

M = cobra.io.load_json_model('models/recon_2_2_simple_medium.json')
model = M.copy() # this way we can edit model but leave M unaltered

**Assignment (10 min):** Take the previous tutorial's last command cell as a template (we already copy-pasted it for you) to do the same analysis here on the full RECON2 model and the PKU disease state.

Tips

Don't reinvent the wheel, the idea is the same as the last tutorial.
Instead of giving 'R1' as the disease reaction you should now give the PKU gene HGNC:8582 or the two reactions it catalyzes as input.
Also think about the fact that there are two equivalent reactions that you have to account for not just one.
The number of exchange reactions here is much bigger than in the example. The FVA computation will take quite a bit longer as a result. It may take 2 minutes or so.



In [ ]:

    
exchanges = [ rxn for rxn in model.reactions if rxn.products == [] and 'EX_' in rxn.id ]

# Shlomi et al suggested using the following medium: everything in (-1) with a few exceptions 
# everything out (unlimited)
# This was actually first proposed in (Sahoo et al, 2012)
for rxn in exchanges:
    rxn.lower_bound = -1
    rxn.upper_bound = 999999
        
# specifics
M.reactions.EX_o2_e.lower_bound = -40

for rxn in ['EX_h2o_e','EX_h_e','EX_co2_e','EX_nh4_e','EX_pi_e','EX_hco3_e','EX_so4_e']:
    M.reactions.get_by_id(rxn).lower_bound = - 100

# to reduce computation time we check all amino acids + a couple neurotransmitters
biomarkers_to_check = ['EX_his_L_e','EX_ile_L_e','EX_leu_L_e','EX_lys_L_e','EX_met_L_e',
                       'EX_phe_L_e','EX_thr_L_e','EX_trp_L_e','EX_val_L_e','EX_cys_L_e',
                       'EX_glu_L_e','EX_tyr_L_e','EX_ala_L_e','EX_asp_L_e','EX_gly_e',
                       'EX_arg_L_e','EX_gln_L_e','EX_pro_L_e','EX_ser_L_e','EX_asn_L_e',
                        'EX_dopa_e','EX_adrnl_e','EX_srtn_e']

# UNCOMMENT & FIX THE LINE BELOW
findBiomarkers(model,fvaRxns=biomarkers_to_check,mods=['HGNC:8582'],synchronous=True)



In [ ]:

    
len(exchanges)



In [ ]:

    
findBiomarkers(model,fvaRxns=exchanges,mods=['HGNC:8582'],synchronous=True)

**Assignment (3 min):** If you see the prediction of the biofluids/tissue as the brain tissue: does the model correctly predict issues with neurotransmitters in the brain?

**Assignment (10 min):** What biomarkers are predicted when you focus on blocking the cofactor biopterin recycling reactions that also produce PKU? To get you started we included some code below.



In [ ]:

    
model.reactions.DHPR.gene_reaction_rule, model.reactions.DHPR.reaction

model.reactions.DHPR2.gene_reaction_rule,model.reactions.DHPR2.reaction

model.reactions.r0398.gene_reaction_rule,model.reactions.r0398.reaction



In [ ]:

    
# UNCOMMENT & FIX THE LINE BELOW
# findBiomarkers(model,fvaRxns=biomarkers_to_check,mods=[ADD THE RELEVANT GENES HERE],synchronous=True)

**Assignment (3 min):** Is this what you expected? If not, is the gene annotation perhaps an issue? Check what would happen if you gave 'findBiomarkers' the full list of reactions as input.



In [ ]:

    
# Copy what you did above but give reactions as input

**Assignment (3 min):** Are there any differences between the predictions for the two different ways to get PKU?



In [ ]: