Metabolism in soy mature leaves.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modelling metabolism in soy mature leaves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating SOY core model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra import io\n",
    "from cobra.core import Metabolite\n",
    "from cobra import flux_analysis\n",
    "\n",
    "#import sbml file\n",
    "fname = \"Data/Soy_core_model_GA.sbml\"\n",
    "model = io.sbml.read_sbml_model(fname)\n",
    "\n",
    "#change Biomass_tx to Protein_biomass\n",
    "rxn = model.reactions.get_by_id(\"Biomass_tx\")\n",
    "rxn.id = \"Protein_biomass_tx\"\n",
    "for met in [\"K_b\",\"Ca_b\",\"Mg_b\"]:\n",
    "    met  = model.metabolites.get_by_id(met)\n",
    "    coeff = rxn.metabolites.get(met)\n",
    "    rxn.add_metabolites({met:-1*coeff})\n",
    "met = Metabolite(\"Protein_b\",name=\"Protein_biomass\")\n",
    "met.compartment=\"b\"\n",
    "formula_dict = rxn.check_mass_balance()\n",
    "met.formula = \"\".join([atom+str(formula_dict[atom]*-1) for atom in formula_dict.keys() if atom != \"charge\"])\n",
    "met.charge = formula_dict[\"charge\"]*-1\n",
    "rxn.add_metabolites({met:1})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2098"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(model.genes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing transcript data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/sanu/.local/lib/python3.6/site-packages/IPython/core/interactiveshell.py:3051 \u001b[1;31mDtypeWarning\u001b[0m: Columns (0,2,3,13,14) have mixed types. Specify dtype option on import or set low_memory=False.\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "df = pd.read_csv(\"/home/sanu/Downloads/gmax_leaf_CO2_normExprn_1.txt\",sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Glyma2.0', 'Query_Length', 'Uniport_ID', 'Hit_Alias', 'Hit_Length', 'E-Value', 'Alignment_Length', 'Query_Coverage', 'Percent_Identity', 'Identities_Covered_percent', 'Number_of_Matches', 'Bit_Score', 'HSPs', 'gene name/function', 'EC', 'ProbeID', 'Glyma 1.1', 'X474_E_Amb.CEL', 'X554_E_Amb.CEL', 'X560_E_Amb.CEL', 'X561_E_Amb.CEL', 'X580_E_Amb.CEL', 'X582_E_Amb.CEL', 'X583_E_Amb.CEL', 'X585_E_Amb.CEL']\n",
      "['Glyma2.0', 'Query_Length', 'Uniport_ID', 'Hit_Alias', 'Hit_Length', 'E-Value', 'Alignment_Length', 'Query_Coverage', 'Percent_Identity', 'Identities_Covered_percent', 'Number_of_Matches', 'Bit_Score', 'HSPs', 'gene name/function', 'EC', 'ProbeID', 'Glyma 1.1', 'X460_E_Ele.CEL', 'X467_E_Ele.CEL', 'X468_E_Ele.CEL', 'X469_E_Ele.CEL', 'X556_E_Ele.CEL', 'X557_E_Ele.CEL', 'X558_E_Ele.CEL', 'X559_E_Ele.CEL']\n"
     ]
    }
   ],
   "source": [
    "amb_list = list(df.columns[0:17])\n",
    "ele_list = list(df.columns[0:17])\n",
    "for x in df.columns:\n",
    "    if \"Amb\" in x:\n",
    "        amb_list.append(x)\n",
    "    elif \"Ele\" in x:\n",
    "        ele_list.append(x)\n",
    "    else:\n",
    "        continue\n",
    "\n",
    "print(amb_list)\n",
    "print(ele_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_amb = df[amb_list]\n",
    "df_ele = df[ele_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "modelGenes = dict()\n",
    "for gene in model.genes:\n",
    "    geneID = gene.id.split(\".\")[0]+\".\"+gene.id.split(\".\")[1]\n",
    "    modelGenes[geneID]=gene\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "modelGenesNotInTranscript = set(modelGenes.keys()) - set(df[\"Glyma2.0\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1094"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(modelGenesNotInTranscript)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "fout = open(\"/home/sanu/SoyGenesInModel.csv\",\"w\")\n",
    "for gene in model.genes:\n",
    "    geneID = gene.id.split(\".\")[0]+\".\"+gene.id.split(\".\")[1]\n",
    "    fout.write(geneID+\"\\n\")\n",
    "fout.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "fout = open(\"/home/sanu/SoyGenesInModelMinusTranscriptomics.csv\",\"w\")\n",
    "for gene in modelGenesNotInTranscript:\n",
    "    geneID = gene.split(\".\")[0]+\".\"+gene.split(\".\")[1]\n",
    "    fout.write(geneID+\"\\n\")\n",
    "fout.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "401\n"
     ]
    }
   ],
   "source": [
    "i=0\n",
    "for rxn in model.reactions:\n",
    "    if len(rxn.genes)!=0:\n",
    "        i=i+1\n",
    "print(i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use biomass composition to generate biomass equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'CPD_17291'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-13-41a3f2818a75>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0mrxn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mReaction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Biomass_leaf_tx\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mmet\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mdf_biomass\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"Unnamed: 0\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m     \u001b[0mrxn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_metabolites\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmetabolites\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_by_id\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmet\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdf_biomass\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdf_biomass\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"Unnamed: 0\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0mmet\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"leaf\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     11\u001b[0m \u001b[0mrxn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlower_bound\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0mrxn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupper_bound\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1000\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/cobra/core/dictlist.py\u001b[0m in \u001b[0;36mget_by_id\u001b[0;34m(self, id)\u001b[0m\n\u001b[1;32m     56\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget_by_id\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mid\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     57\u001b[0m         \u001b[0;34m\"\"\"return the element with a matching id\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 58\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__getitem__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_dict\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mid\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     59\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     60\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mlist_attr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mattribute\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: 'CPD_17291'"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df_biomass = pd.read_csv(\"Data/biomass_soy.csv\")\n",
    "\n",
    "\n",
    "from cobra.core import Reaction\n",
    "\n",
    "rxn = Reaction(\"Biomass_leaf_tx\")\n",
    "for met in df_biomass[\"Unnamed: 0\"]:\n",
    "    rxn.add_metabolites({model.metabolites.get_by_id(met):-1*float(df_biomass[df_biomass[\"Unnamed: 0\"]==met][\"leaf\"])})\n",
    "rxn.lower_bound = 0\n",
    "rxn.upper_bound = 1000\n",
    "model.add_reaction(rxn)\n",
    "\n",
    "rxn = Reaction(\"Biomass_stem_tx\")\n",
    "for met in df_biomass[\"Unnamed: 0\"]:\n",
    "    rxn.add_metabolites({model.metabolites.get_by_id(met):-1*float(df_biomass[df_biomass[\"Unnamed: 0\"]==met][\"stem\"])})\n",
    "rxn.lower_bound = 0\n",
    "rxn.upper_bound = 1000\n",
    "model.add_reaction(rxn)\n",
    "\n",
    "rxn = Reaction(\"Biomass_root_tx\")\n",
    "for met in df_biomass[\"Unnamed: 0\"]:\n",
    "    rxn.add_metabolites({model.metabolites.get_by_id(met):-1*float(df_biomass[df_biomass[\"Unnamed: 0\"]==met][\"root\"])})\n",
    "rxn.lower_bound = 0\n",
    "rxn.upper_bound = 1000\n",
    "model.add_reaction(rxn)\n",
    "\n",
    "rxn = Reaction(\"Biomass_seed_tx\")\n",
    "for met in df_biomass[\"Unnamed: 0\"]:\n",
    "    rxn.add_metabolites({model.metabolites.get_by_id(met):-1*float(df_biomass[df_biomass[\"Unnamed: 0\"]==met][\"seed\"])})\n",
    "rxn.lower_bound = 0\n",
    "rxn.upper_bound = 1000\n",
    "model.add_reaction(rxn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run stem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra import flux_analysis\n",
    "stem_model = model.copy()\n",
    "stem_model.reactions.GLC_tx.upper_bound = 0\n",
    "stem_model.reactions.GLC_tx.lower_bound = 0\n",
    "stem_model.reactions.Sucrose_tx.upper_bound = 1\n",
    "stem_model.reactions.Sucrose_tx.lower_bound = 1\n",
    "stem_model.reactions.Photon_tx.upper_bound = 0\n",
    "stem_model.reactions.Photon_tx.lower_bound = 0\n",
    "stem_model.reactions.Biomass_stem_tx.objective_coefficient=1\n",
    "flux_analysis.parsimonious.pfba(stem_model)\n",
    "print(\"Biomass flux =\"+str(stem_model.reactions.Biomass_stem_tx.flux))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra import flux_analysis\n",
    "root_model = model.copy()\n",
    "root_model.reactions.GLC_tx.upper_bound = 0\n",
    "root_model.reactions.GLC_tx.lower_bound = 0\n",
    "root_model.reactions.Sucrose_tx.upper_bound = 1\n",
    "root_model.reactions.Sucrose_tx.lower_bound = 1\n",
    "root_model.reactions.Photon_tx.upper_bound = 0\n",
    "root_model.reactions.Photon_tx.lower_bound = 0\n",
    "root_model.reactions.Biomass_root_tx.objective_coefficient=1\n",
    "flux_analysis.parsimonious.pfba(root_model)\n",
    "print(\"Biomass flux =\"+str(root_model.reactions.Biomass_root_tx.flux))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra import flux_analysis\n",
    "seed_model = model.copy()\n",
    "seed_model.reactions.GLC_tx.upper_bound = 0\n",
    "seed_model.reactions.GLC_tx.lower_bound = 0\n",
    "seed_model.reactions.Sucrose_tx.upper_bound = 1\n",
    "seed_model.reactions.Sucrose_tx.lower_bound = 1\n",
    "seed_model.reactions.Photon_tx.upper_bound = 0\n",
    "seed_model.reactions.Photon_tx.lower_bound = 0\n",
    "seed_model.reactions.Biomass_seed_tx.objective_coefficient=1\n",
    "sol=flux_analysis.parsimonious.pfba(seed_model)\n",
    "print(\"Biomass flux =\"+str(seed_model.reactions.Biomass_seed_tx.flux))\n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate diel leaf model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "########################################################\n",
    "#This function was used to set up a C3 leaf diel model #\n",
    "########################################################\n",
    "def setupC3DielModel(core_model,transferMets,starch_sucrose_ratio=90.0/100):\n",
    "    from cobra.core import Metabolite, Reaction\n",
    "    import re\n",
    "\n",
    "    #create two copies of model elements for day and night\n",
    "    cobra_model2 = core_model.copy()\n",
    "    for met in cobra_model2.metabolites:\n",
    "        met.id = met.id+\"1\"\n",
    "        met.compartment = met.compartment+\"1\"\n",
    "    for rxn in cobra_model2.reactions:\n",
    "        rxn.id = rxn.id+\"1\"\n",
    "\n",
    "    cobra_model3 = core_model.copy()\n",
    "    for met in cobra_model3.metabolites:\n",
    "        met.id = met.id+\"2\"\n",
    "        met.compartment = met.compartment+\"2\"\n",
    "    for rxn in cobra_model3.reactions:\n",
    "        rxn.id = rxn.id+\"2\"\n",
    "\n",
    "    #merge the day and night model\n",
    "    cobra_model = cobra_model2+cobra_model3\n",
    "    for met in cobra_model3.metabolites:\n",
    "        if not cobra_model.metabolites.__contains__(met.id):\n",
    "            cobra_model.add_metabolites(met.copy())\n",
    "\n",
    "    met1 = Metabolite(\"X_Phloem_contribution_t1\",name=\"Phloem output during the day\")\n",
    "    cobra_model.reactions.get_by_id(\"Phloem_output_tx1\").add_metabolites({met1:1})\n",
    "    met2 = Metabolite(\"X_Phloem_contribution_t2\",name=\"Phloem output during at night\")\n",
    "    cobra_model.reactions.get_by_id(\"Phloem_output_tx2\").add_metabolites({met2:1})\n",
    "\n",
    "    rxn = Reaction(\"diel_biomass\")\n",
    "    rxn.add_metabolites({met1:-3,met2:-1})\n",
    "    rxn.lower_bound = 0\n",
    "    rxn.upper_bound = 1000\n",
    "    cobra_model.add_reaction(rxn)\n",
    "\n",
    "    #Adding reactions to allow for day-night metabolite accumulations\n",
    "    tmfile = open(transferMets,\"r\")\n",
    "    tmset=set()\n",
    "    for line in tmfile:\n",
    "        tmset.add(line.replace(\"\\n\",\"\"))\n",
    "\n",
    "    for met in tmset:\n",
    "        if met == \"AMMONIUM_v\" or met==\"FRUCTAN_v\":\n",
    "            continue\n",
    "        tempRxn = Reaction(met+\"_dielTransfer\")\n",
    "        tempRxn.add_metabolites({cobra_model.metabolites.get_by_id(met+\"1\"):-1,cobra_model.metabolites.get_by_id(met+\"2\"):1})\n",
    "        tempRxn.lower_bound=-1000\n",
    "        if not ((met == \"STARCH_p\") or (met == \"SUCROSE_v\") or (met == \"MAL_v\") or (met == \"aMAL_v\") or (met == \"NITRATE_v\") or (met == \"CIT_v\") or (met == \"aCIT_v\") or (met == \"PROTON_v\")):\n",
    "            tempRxn.lower_bound=0\n",
    "        tempRxn.upper_bound=1000\n",
    "        cobra_model.add_reaction(tempRxn)\n",
    "\n",
    "    fractionMets=dict()\n",
    "    for rxn in cobra_model.reactions:\n",
    "        for met in rxn.metabolites.keys():\n",
    "            prefix=\"\"\n",
    "            a=re.search(\"^a{1,3}\",met.id)\n",
    "            anion=\"\"\n",
    "            if a:\n",
    "                anion=a.group(0)\n",
    "                prefix=anion\n",
    "            b=re.search(\"^b{1,3}\",met.id)\n",
    "            basic=\"\"\n",
    "            if b:\n",
    "                basic=b.group(0)\n",
    "                prefix=basic\n",
    "            if ((not prefix == \"\") and met.compartment == \"v1\"):\n",
    "                fractionMets[met]=prefix\n",
    "\n",
    "    temp=cobra_model.copy()\n",
    "    for met in fractionMets.keys():\n",
    "        for rxn in met.reactions:\n",
    "            if rxn.id.__contains__(\"_dielTransfer\"):\n",
    "                continue\n",
    "            else:\n",
    "                mainMet = met.id[len(fractionMets[met]):]\n",
    "                coeff1 = temp.reactions.get_by_id(rxn.id).metabolites.get(temp.metabolites.get_by_id(mainMet))\n",
    "                coeff2 = temp.reactions.get_by_id(rxn.id).metabolites.get(temp.metabolites.get_by_id(met.id))\n",
    "                if not coeff1:\n",
    "                    coeff1=0\n",
    "                if not coeff2:\n",
    "                    coeff2=0\n",
    "                total = coeff1 + coeff2\n",
    "                coeff1 = float(coeff1)/total\n",
    "                coeff2 = float(coeff2)/total\n",
    "                if cobra_model.reactions.has_id(met.id[0:len(met.id)-1]+\"_dielTransfer\"):\n",
    "                    ub = temp.reactions.get_by_id(met.id[0:len(met.id)-1]+\"_dielTransfer\").upper_bound\n",
    "                    lb = temp.reactions.get_by_id(met.id[0:len(met.id)-1]+\"_dielTransfer\").lower_bound\n",
    "                    temp.reactions.get_by_id(met.id[0:len(met.id)-1]+\"_dielTransfer\").remove_from_model()\n",
    "                    temp.reactions.get_by_id(mainMet[0:len(mainMet)-1]+\"_dielTransfer\").remove_from_model()\n",
    "                    Reac = Reaction(mainMet[0:len(mainMet)-1]+\"_dielTransfer\",name=mainMet+\"_dielTransfer\")\n",
    "                    Reac.add_metabolites({temp.metabolites.get_by_id(met.id[0:len(met.id)-1]+\"1\"):-coeff2,temp.metabolites.get_by_id(met.id[0:len(met.id)-1]+\"2\"):coeff2,temp.metabolites.get_by_id(mainMet[0:len(mainMet)-1]+\"1\"):-coeff1,temp.metabolites.get_by_id(mainMet[0:len(mainMet)-1]+\"2\"):coeff1})\n",
    "                    Reac.lower_bound=lb\n",
    "                    Reac.upper_bound=ub\n",
    "                    temp.add_reaction(Reac)\n",
    "                    print Reac.reaction\n",
    "                break\n",
    "    ####ADD CONSTRAINTS TO MODEL####\n",
    "    cobra_model = temp.copy()\n",
    "\n",
    "    #objective function\n",
    "    cobra_model.reactions.get_by_id(\"diel_biomass\").objective_coefficient=1\n",
    "    #Leaves - light\n",
    "    cobra_model.reactions.get_by_id(\"Sucrose_tx1\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"Sucrose_tx1\").upper_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"GLC_tx1\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"GLC_tx1\").upper_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"CO2_tx1\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"NH4_tx1\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"NH4_tx1\").upper_bound=0\n",
    "    #Leaves - dark\n",
    "    cobra_model.reactions.get_by_id(\"Sucrose_tx2\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"Sucrose_tx2\").upper_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"GLC_tx2\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"GLC_tx2\").upper_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"Photon_tx2\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"Photon_tx2\").upper_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"NH4_tx2\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"NH4_tx2\").upper_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"CO2_tx2\").upper_bound=0\n",
    "\n",
    "    #Set pG6P transporter to 0\n",
    "    cobra_model.reactions.get_by_id(\"G6P_Pi_pc1\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"G6P_Pi_pc1\").upper_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"G6P_Pi_pc2\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"G6P_Pi_pc2\").upper_bound=0\n",
    "\n",
    "    #Turn off PTOX\n",
    "    cobra_model.reactions.get_by_id(\"Plastoquinol_Oxidase_p1\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"Plastoquinol_Oxidase_p1\").upper_bound=0\n",
    "\n",
    "    #nitrate uptake constrain\n",
    "    Nitrate_balance = Metabolite(\"Nitrate_bal_c\", name = \"Weights to balance nitrate uptake\", compartment = \"c1\")\n",
    "    cobra_model.reactions.get_by_id(\"Nitrate_ec1\").add_metabolites({Nitrate_balance:-2})\n",
    "    cobra_model.reactions.get_by_id(\"Nitrate_ec2\").add_metabolites({Nitrate_balance:3})\n",
    "\n",
    "    #Rubisco balance\n",
    "    Rubisco_balance = Metabolite(\"rubisco_bal_p1\", name = \"Weights to balance RuBP carboxygenase oxygenase balance\", compartment = \"p1\")\n",
    "    cobra_model.reactions.get_by_id(\"RXN_961_p1\").add_metabolites({Rubisco_balance:3})\n",
    "    cobra_model.reactions.get_by_id(\"RIBULOSE_BISPHOSPHATE_CARBOXYLASE_RXN_p1\").add_metabolites({Rubisco_balance:-1})\n",
    "\n",
    "    #generic ATPase and NADPH oxidase\n",
    "    Maintenance_constraint = Metabolite(\"ATPase_NADPHoxidase_constraint_c1\",name =  \"ATPase_NADPHoxidase_constraint_c1\", compartment = \"c1\")\n",
    "    Maintenance_constraint2 = Metabolite(\"ATPase_NADPHoxidase_constraint_c2\",name =  \"ATPase_NADPHoxidase_constraint_c2\", compartment = \"c2\")\n",
    "    Maintenance_constraint3 = Metabolite(\"Light_dark_maintainence_constraint\",name =  \"Light_dark_maintainence_constraint\", compartment = \"c1\")\n",
    "    cobra_model.reactions.get_by_id(\"ATPase_tx1\").add_metabolites({Maintenance_constraint:1,Maintenance_constraint3:1})\n",
    "    cobra_model.reactions.get_by_id(\"ATPase_tx2\").add_metabolites({Maintenance_constraint2:1,Maintenance_constraint3:-1})\n",
    "    cobra_model.reactions.get_by_id(\"NADPHoxc_tx1\").add_metabolites({Maintenance_constraint:-3})\n",
    "    cobra_model.reactions.get_by_id(\"NADPHoxc_tx2\").add_metabolites({Maintenance_constraint2:-3})\n",
    "    cobra_model.reactions.get_by_id(\"NADPHoxm_tx1\").add_metabolites({Maintenance_constraint:-3})\n",
    "    cobra_model.reactions.get_by_id(\"NADPHoxm_tx2\").add_metabolites({Maintenance_constraint2:-3})\n",
    "    cobra_model.reactions.get_by_id(\"NADPHoxp_tx1\").add_metabolites({Maintenance_constraint:-3})\n",
    "    cobra_model.reactions.get_by_id(\"NADPHoxp_tx2\").add_metabolites({Maintenance_constraint2:-3})\n",
    "\n",
    "    ##constrain sucrose and starch storage\n",
    "    Sucorse_starch_balance = Metabolite(\"sucrose_starch_bal_c\", name = \"Weights to balance sucrose-starch uptake\", compartment = \"c1\")\n",
    "    cobra_model.reactions.get_by_id(\"SUCROSE_v_dielTransfer\").add_metabolites({Sucorse_starch_balance:-1*starch_sucrose_ratio})\n",
    "    cobra_model.reactions.get_by_id(\"STARCH_p_dielTransfer\").add_metabolites({Sucorse_starch_balance:1})\n",
    "\n",
    "    #Plastid enolase was not detected in Arabidopsis mesophyll tissue\n",
    "    cobra_model.reactions.get_by_id(\"2PGADEHYDRAT_RXN_p1\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"2PGADEHYDRAT_RXN_p1\").upper_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"2PGADEHYDRAT_RXN_p2\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"2PGADEHYDRAT_RXN_p2\").upper_bound=0\n",
    "\n",
    "    #Setting chloroplastic NADPH dehydrogenase to 0  ((Yamamoto et al., 2011)\n",
    "    cobra_model.reactions.get_by_id(\"NADPH_Dehydrogenase_p1\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"NADPH_Dehydrogenase_p1\").upper_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"NADPH_Dehydrogenase_p2\").lower_bound=0\n",
    "    cobra_model.reactions.get_by_id(\"NADPH_Dehydrogenase_p2\").upper_bound=0\n",
    "\n",
    "    #ATP_ADP_Pi constrained to 0 because while there is evidence for its existance, it does not carry flux during the day\n",
    "    cobra_model.reactions.get_by_id(\"ATP_ADP_Pi_pc1\").lower_bound = 0\n",
    "    cobra_model.reactions.get_by_id(\"ATP_ADP_Pi_pc1\").upper_bound = 0\n",
    "\n",
    "    return cobra_model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5 CIT_v1 + 0.5 aCIT_v1 <=> 0.5 CIT_v2 + 0.5 aCIT_v2\n",
      "0.7 MAL_v1 + 0.3 aMAL_v1 <=> 0.7 MAL_v2 + 0.3 aMAL_v2\n",
      "bHIS_v1 --> bHIS_v2\n",
      "106.674498326\n",
      "1500\n",
      "10.1351\n",
      "21.6455696203\n",
      "Phloem output = \n",
      "0.452844044741\n",
      "1500\n",
      "10.1351\n",
      "24.4303797468\n",
      "Phloem output = \n",
      "0.521458572285\n"
     ]
    }
   ],
   "source": [
    "leaf_model = model.copy()\n",
    "\n",
    "leaf_model = setupC3DielModel(leaf_model,\"MetabolitesToTransfer.txt\",starch_sucrose_ratio=5)\n",
    "\n",
    "for rxn in leaf_model.reactions:\n",
    "    if rxn.lower_bound==-1000:\n",
    "        rxn.lower_bound=-1500\n",
    "    if rxn.upper_bound==1000:\n",
    "        rxn.upper_bound=1500\n",
    "\n",
    "\n",
    "PPFD = 1500\n",
    "ATPase = (0.0049*PPFD) + 2.7851\n",
    "leaf_model.reactions.get_by_id(\"ATPase_tx1\").lower_bound = ATPase\n",
    "leaf_model.reactions.get_by_id(\"ATPase_tx1\").upper_bound = ATPase\n",
    "leaf_model.reactions.diel_biomass.objective_coefficient=1\n",
    "sol = flux_analysis.parsimonious.pfba(leaf_model)\n",
    "print(leaf_model.reactions.CO2_tx1.flux)\n",
    "\n",
    "from sweetlovegroup import analysis\n",
    "temp = leaf_model.copy()\n",
    "out_400 = analysis.estimateOutputFromNetCO2(temp,21.645569620253173)\n",
    "leaf_model_400 = leaf_model.copy()\n",
    "leaf_model_400.reactions.diel_biomass.lower_bound = out_400\n",
    "leaf_model_400.reactions.diel_biomass.upper_bound = out_400\n",
    "sol = flux_analysis.parsimonious.pfba(leaf_model_400)\n",
    "print(\"Phloem output = \")\n",
    "print(leaf_model_400.reactions.diel_biomass.flux)\n",
    "\n",
    "\n",
    "temp = leaf_model.copy()\n",
    "out_550 = analysis.estimateOutputFromNetCO2(temp,24.43037974683546)\n",
    "leaf_model_550 = leaf_model.copy()\n",
    "leaf_model_550.reactions.diel_biomass.lower_bound = out_550\n",
    "leaf_model_550.reactions.diel_biomass.upper_bound = out_550\n",
    "sol = flux_analysis.parsimonious.pfba(leaf_model_550)\n",
    "print(\"Phloem output = \")\n",
    "print(leaf_model_550.reactions.diel_biomass.flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "leaf_model.compartments={'b1': 'Biomass_1',\"b2\":\"Biomass_2\"}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b1\n",
      "Biomass_1\n",
      "b2\n",
      "Biomass_2\n",
      "mc2\n",
      "\n",
      "mc1\n",
      "\n",
      "m1\n",
      "\n",
      "m2\n",
      "\n",
      "i1\n",
      "\n",
      "i2\n",
      "\n",
      "v1\n",
      "\n",
      "v2\n",
      "\n",
      "x2\n",
      "\n",
      "c2\n",
      "\n",
      "c1\n",
      "\n",
      "x1\n",
      "\n",
      "e1\n",
      "\n",
      "e2\n",
      "\n",
      "p2\n",
      "\n",
      "p1\n",
      "\n",
      "r1\n",
      "\n",
      "r2\n",
      "\n",
      "l2\n",
      "\n",
      "l1\n",
      "\n",
      "mi1\n",
      "\n",
      "mi2\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for comp in leaf_model.compartments.keys():\n",
    "    print comp\n",
    "    print leaf_model.compartments[comp]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gene Glyma.06G074900for rxn ACETYLGLUTKIN_RXN_p1 present\n",
      "Gene Glyma.13G140700for rxn DIAMINOPIMDECARB_RXN_p1 present\n",
      "Gene Glyma.02G222400for rxn F16ALDOLASE_RXN_p1 present\n",
      "Gene Glyma.02G222400for rxn F16ALDOLASE_RXN_c1 present\n",
      "Gene Glyma.20G079400for rxn NADH_DEHYDROG_A_RXN_mi1 present\n",
      "Gene Glyma.13G097100for rxn RIBULOSE_BISPHOSPHATE_CARBOXYLASE_RXN_p1 present\n",
      "Gene Glyma.14G053300for rxn IGPSYN_RXN_p1 present\n",
      "Gene Glyma.10G016200for rxn FUMHYDR_RXN_m1 present\n",
      "Gene Glyma.13G204700for rxn GLYCOGENSYN_RXN_p1 present\n",
      "Gene Glyma.01G026700for rxn ALANINE_AMINOTRANSFERASE_RXN_c1 present\n",
      "Gene Glyma.03G173400for rxn DIAMINOPIMEPIM_RXN_p1 present\n",
      "Gene Glyma.10G064000for rxn THREDEHYD_RXN_p1 present\n",
      "Gene Glyma.20G148100for rxn ACSERLY_RXN_c1 present\n",
      "Gene Glyma.17G159700for rxn 2_ISOPROPYLMALATESYN_RXN_p1 present\n",
      "Gene Glyma.14G217300for rxn 3_ISOPROPYLMALISOM_RXN_p1 present\n",
      "Gene Glyma.12G042400for rxn GLUC1PADENYLTRANS_RXN_p1 present\n",
      "Gene Glyma.20G148100for rxn ACSERLY_RXN_p1 present\n",
      "Gene Glyma.12G193800for rxn ACONITATEHYDR_RXN_m1 present\n",
      "Gene Glyma.13G144900for rxn ISOCITRATE_DEHYDROGENASE_NAD_RXN_m1 present\n",
      "Gene Glyma.03G027500for rxn 1TRANSKETO_RXN_p1 present\n",
      "Gene Glyma.11G158400for rxn IMPCYCLOHYDROLASE_RXN_p1 present\n",
      "Gene Glyma.07G212800for rxn FERREDOXIN_NITRITE_REDUCTASE_RXN_p1 present\n",
      "Gene Glyma.16G049300for rxn ASPARTATEKIN_RXN_p1 present\n",
      "Gene Glyma.01G082200for rxn 3_DEHYDROQUINATE_DEHYDRATASE_RXN_p1 present\n",
      "Gene Glyma.06G015900for rxn 1_PERIOD_2_PERIOD_1_PERIOD_13_RXN_p1 present\n",
      "Gene Glyma.18G002900for rxn HOMOSERDEHYDROG_RXN_NADP_p1 present\n",
      "Gene Glyma.09G040100for rxn CYTOCHROME_C_OXIDASE_RXN_mi1 present\n",
      "Gene Glyma.17G101700for rxn RIBULP3EPIM_RXN_p1 present\n",
      "Gene Glyma.02G254000for rxn AICARSYN_RXN_p1 present\n",
      "Gene Glyma.17G261700for rxn CATAL_RXN_x1 present\n",
      "Gene Glyma.15G263100for rxn IMP_DEHYDROG_RXN_c1 present\n",
      "Gene Glyma.02G254000for rxn AMPSYN_RXN_p1 present\n",
      "Gene Glyma.10G197700for rxn MALATE_DEHYDROGENASE_NADP_RXN_p1 present\n",
      "Gene Glyma.13G336200for rxn TRIOSEPISOMERIZATION_RXN_c1 present\n",
      "Gene Glyma.12G193800for rxn ACONITATEDEHYDR_RXN_c1 present\n",
      "Gene Glyma.13G336200for rxn TRIOSEPISOMERIZATION_RXN_p1 present\n",
      "Gene Glyma.11G254700for rxn ADENOSYLHOMOCYSTEINASE_RXN_c1 present\n",
      "Gene Glyma.11G000500for rxn 2_PERIOD_4_PERIOD_1_PERIOD_111_RXN_c1 present\n",
      "Gene Glyma.16G087200for rxn 1_PERIOD_10_PERIOD_2_PERIOD_2_RXN_mi1 present\n",
      "Gene Glyma.11G241900for rxn HISTCYCLOHYD_RXN_p1 present\n",
      "Gene Glyma.07G182300for rxn F16BDEPHOS_RXN_c1 present\n",
      "Gene Glyma.07G182300for rxn F16BDEPHOS_RXN_p1 present\n",
      "Gene Glyma.05G240700for rxn ASPARTATE_SEMIALDEHYDE_DEHYDROGENASE_RXN_p1 present\n",
      "Gene Glyma.10G241000for rxn PHOSPHOGLUCMUT_RXN_p1 present\n",
      "Gene Glyma.10G201100for rxn PEPDEPHOS_RXN_p1 present\n",
      "Gene Glyma.03G240100for rxn BRANCHED_CHAINAMINOTRANSFERILEU_RXN_p1 present\n",
      "Gene Glyma.09G270500for rxn 3PGAREARR_RXN_c1 present\n",
      "Gene Glyma.10G241000for rxn PHOSPHOGLUCMUT_RXN_c1 present\n",
      "Gene Glyma.05G005300for rxn ISOCITDEH_RXN_c1 present\n",
      "Gene Glyma.05G046700for rxn ARGSUCCINSYN_RXN_p1 present\n",
      "Gene Glyma.05G181000for rxn PREPHENATE_TRANSAMINE_RXN_p1 present\n",
      "Gene Glyma.16G028900for rxn DIHYDROPICRED_RXN_NADP_p1 present\n",
      "Gene Glyma.10G242600for rxn SULFATE_ADENYLYLTRANS_RXN_p1 present\n",
      "Gene Glyma.14G165000for rxn NITRATE_REDUCTASE_NADH_RXN_c1 present\n",
      "Gene Glyma.11G078700for rxn 3_DEHYDROQUINATE_SYNTHASE_RXN_p1 present\n",
      "Gene Glyma.16G058800for rxn ANTHRANSYN_RXN_p1 present\n",
      "Gene Glyma.13G049500for rxn RXN0_2381_p1 present\n",
      "Gene Glyma.05G183700for rxn SHIKIMATE_KINASE_RXN_p1 present\n",
      "Gene Glyma.08G195000for rxn INORGPYROPHOSPHAT_RXN_x1 present\n",
      "Gene Glyma.19G131500for rxn PYRROLINECARBREDUCT_RXN_NAD_c1 present\n",
      "Gene Glyma.17G100600for rxn MALATE_DEH_RXN_c1 present\n",
      "Gene Glyma.18G025500for rxn GLUC1PURIDYLTRANS_RXN_c1 present\n",
      "Gene Glyma.20G009400for rxn N_ACETYLGLUTPREDUCT_RXN_p1 present\n",
      "Gene Glyma.05G181000for rxn ASPAMINOTRANS_RXN_p1 present\n",
      "Gene Glyma.17G100600for rxn MALATE_DEH_RXN_p1 present\n",
      "Gene Glyma.04G094600for rxn ARGSUCCINLYA_RXN_p1 present\n",
      "Gene Glyma.17G100600for rxn MALATE_DEH_RXN_m1 present\n",
      "Gene Glyma.08G195000for rxn INORGPYROPHOSPHAT_RXN_c1 present\n",
      "Gene Glyma.07G250600for rxn GLUTATHIONE_REDUCT_NADPH_RXN_p1 present\n",
      "Gene Glyma.05G181000for rxn ASPAMINOTRANS_RXN_c1 present\n",
      "Gene Glyma.08G195000for rxn INORGPYROPHOSPHAT_RXN_p1 present\n",
      "Gene Glyma.17G100600for rxn MALATE_DEH_RXN_x1 present\n",
      "Gene Glyma.04G092400for rxn PGLUCISOM_RXN_c1 present\n",
      "Gene Glyma.04G092400for rxn PGLUCISOM_RXN_p1 present\n",
      "Gene Glyma.11G084100for rxn SERINE_O_ACETTRAN_RXN_c1 present\n",
      "Gene Glyma.18G244900for rxn GPH_RXN_p1 present\n",
      "Gene Glyma.17G109700for rxn SUCROSE_PHOSPHATE_SYNTHASE_RXN_c1 present\n",
      "Gene Glyma.08G302600for rxn SERINE_GLYOXYLATE_AMINOTRANSFERASE_RXN_x1 present\n",
      "Gene Glyma.11G084100for rxn SERINE_O_ACETTRAN_RXN_p1 present\n",
      "Gene Glyma.12G210000for rxn ACETOLACTREDUCTOISOM_RXN_p1 present\n",
      "Gene Glyma.08G180600for rxn RXN_7737_p1 present\n",
      "Gene Glyma.06G172600for rxn GAPOXNPHOSPHN_RXN_c1 present\n",
      "Gene Glyma.11G079000for rxn 3_ISOPROPYLMALDEHYDROG_RXN_p1 present\n",
      "Gene Glyma.19G088000for rxn ADENYL_KIN_RXN_p1 present\n",
      "Gene Glyma.18G093200for rxn GLUTAMIDOTRANS_RXN_p1 present\n",
      "Gene Glyma.19G088000for rxn ADENYL_KIN_RXN_m1 present\n",
      "Gene Glyma.07G081900for rxn HYDROXYPYRUVATE_REDUCTASE_RXN_NAD_x1 present\n",
      "Gene Glyma.15G027000for rxn THRESYN_RXN_p1 present\n",
      "Gene Glyma.04G009800for rxn HISTIDPHOS_RXN_p1 present\n",
      "Gene Glyma.11G226900for rxn SEDOHEPTULOSE_BISPHOSPHATASE_RXN_p1 present\n",
      "Gene Glyma.01G026700for rxn GLYCINE_AMINOTRANSFERASE_RXN_x1 present\n",
      "Gene Glyma.14G184500for rxn GLYOHMETRANS_RXN_m1 present\n",
      "Gene Glyma.16G155800for rxn HISTAMINOTRANS_RXN_p1 present\n",
      "Gene Glyma.09G168900for rxn GLUTDECARBOX_RXN_c1 present\n",
      "Gene Glyma.12G199400for rxn PLASTOQUINOL_PLASTOCYANIN_REDUCTASE_RXN_p1 present\n",
      "Gene Glyma.11G241900for rxn HISTPRATPHYD_RXN_p1 present\n",
      "Gene Glyma.20G019700for rxn GCVMULTI_RXN_m1 present\n",
      "Gene Glyma.07G118300for rxn ACETYLORNTRANSAM_RXN_p1 present\n",
      "Gene Glyma.18G061100for rxn ASNSYNA_RXN_c1 present\n",
      "Gene Glyma.15G131700for rxn RXN_8001_p1 present\n",
      "Gene Glyma.19G130800for rxn GLUTAMATE_SYNTHASE_FERREDOXIN_RXN_p1 present\n",
      "Gene Glyma.17G262400for rxn PSII_RXN_p1 present\n",
      "Gene Glyma.14G065600for rxn GLUTSEMIALDEHYDROG_RXN_c1 present\n",
      "Gene Glyma.19G239800for rxn SHIKIMATE_5_DEHYDROGENASE_RXN_p1 present\n",
      "Gene Glyma.15G130100for rxn 1_PERIOD_18_PERIOD_1_PERIOD_2_RXN_p1 present\n",
      "Gene Glyma.08G309300for rxn PRPPSYN_RXN_p1 present\n",
      "Gene Glyma.02G244000for rxn GLUTAMINESYN_RXN_c1 present\n",
      "Gene Glyma.02G244000for rxn GLUTAMINESYN_RXN_m1 present\n",
      "Gene Glyma.04G256100for rxn 1_PERIOD_5_PERIOD_1_PERIOD_20_RXN_c1 present\n",
      "Gene Glyma.07G186000for rxn BRANCHED_CHAINAMINOTRANSFERLEU_RXN_p1 present\n",
      "Gene Glyma.13G207900for rxn DIHYDROXYISOVALDEHYDRAT_RXN_p1 present\n",
      "Gene Glyma.14G102700for rxn CHORISMATEMUT_RXN_p1 present\n",
      "Gene Glyma.15G102400for rxn GUANYL_KIN_RXN_c1 present\n",
      "Gene Glyma.18G026400for rxn PRAISOM_RXN_p1 present\n",
      "Gene Glyma.10G086600for rxn SUCROSE_PHOSPHATASE_RXN_c1 present\n",
      "Gene Glyma.18G221700for rxn DIHYDRODIPICSYN_RXN_p1 present\n",
      "Gene Glyma.20G145900for rxn IMIDPHOSDEHYD_RXN_p1 present\n",
      "Gene Glyma.05G045200for rxn GDPKIN_RXN_p1 present\n",
      "Gene Glyma.03G212700for rxn METHENYLTHFCYCLOHYDRO_RXN_c1 present\n",
      "Gene Glyma.03G240100for rxn BRANCHED_CHAINAMINOTRANSFERVAL_RXN_p1 present\n",
      "Gene Glyma.02G311600for rxn HOMOSERKIN_RXN_p1 present\n",
      "Gene Glyma.03G176400for rxn ATPPHOSPHORIBOSYLTRANS_RXN_p1 present\n",
      "Gene Glyma.06G169700for rxn ACETOLACTSYN_RXN_p1 present\n",
      "Gene Glyma.12G229400for rxn PEPCARBOX_RXN_c1 present\n",
      "Gene Glyma.01G010200for rxn PHOSPHORIBULOKINASE_RXN_p1 present\n",
      "Gene Glyma.05G045200for rxn GDPKIN_RXN_c1 present\n",
      "Gene Glyma.09G153900for rxn 2PGADEHYDRAT_RXN_c1 present\n",
      "Gene Glyma.04G184900for rxn PRTRANS_RXN_p1 present\n",
      "Gene Glyma.14G176600for rxn DAHPSYN_RXN_p1 present\n",
      "Gene Glyma.18G261600for rxn CYSPH_RXN_p1 present\n",
      "Gene Glyma.01G139600for rxn 2_PERIOD_5_PERIOD_1_PERIOD_19_RXN_p1 present\n",
      "Gene Glyma.06G030600for rxn ORNCARBAMTRANSFER_RXN_p1 present\n",
      "Gene Glyma.15G131700for rxn HISTOLDEHYD_RXN_p1 present\n",
      "Gene Glyma.13G131900for rxn GLUTAMATE_N_ACETYLTRANSFERASE_RXN_p1 present\n",
      "Gene Glyma.11G158400for rxn AICARTRANSFORM_RXN_p1 present\n",
      "Gene Glyma.16G165900for rxn CELLULOSE_SYNTHASE_GDP_FORMING_RXN_c1 present\n",
      "Gene Glyma.03G246300for rxn RIB5PISOM_RXN_p1 present\n",
      "Gene Glyma.19G142400for rxn ADENYLOSUCCINATE_SYNTHASE_RXN_p1 present\n",
      "Gene Glyma.05G045200for rxn UDPKIN_RXN_c1 present\n",
      "Gene Glyma.01G031500for rxn RXN66_3_m1 present\n",
      "Gene Glyma.09G178100for rxn RXN_10773_c1 present\n",
      "Gene Glyma.14G217300for rxn RXN_8991_p1 present\n",
      "Gene Glyma.08G113600for rxn ACETATE_COA_LIGASE_RXN_p1 present\n",
      "Gene Glyma.13G049500for rxn RXN0_2382_p1 present\n",
      "Gene Glyma.15G267000for rxn ATP_CITRATE_PRO_S_LYASE_RXN_c1 present\n",
      "Gene Glyma.03G055900for rxn AMP_DEAMINASE_RXN_c1 present\n",
      "Gene Glyma.13G319000for rxn CARBOXYCYCLOHEXADIENYL_DEHYDRATASE_RXN_p1 present\n",
      "Gene Glyma.08G015000for rxn GMP_SYN_GLUT_RXN_c1 present\n",
      "Gene Glyma.03G181900for rxn ADENOSINE_KINASE_RXN_c1 present\n",
      "Gene Glyma.11G093200for rxn SULFITE_REDUCTASE_FERREDOXIN_RXN_p1 present\n",
      "Gene Glyma.14G065600for rxn GLUTKIN_RXN_c1 present\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gene Glyma.09G167000for rxn SUCROSE_SYNTHASE_RXN_c1 present\n",
      "Gene Glyma.10G211400for rxn CHORISMATE_SYNTHASE_RXN_p1 present\n",
      "Gene Glyma.03G129700for rxn CYSTATHIONINE_BETA_LYASE_RXN_p1 present\n",
      "Gene Glyma.13G153700for rxn PRIBFAICARPISOM_RXN_p1 present\n",
      "Gene Glyma.03G212700for rxn METHYLENETHFDEHYDROG_NADP_RXN_c1 present\n",
      "Gene Glyma.02G222400for rxn F16ALDOLASE_RXN_p2 present\n",
      "Gene Glyma.02G222400for rxn F16ALDOLASE_RXN_c2 present\n",
      "Gene Glyma.20G079400for rxn NADH_DEHYDROG_A_RXN_mi2 present\n",
      "Gene Glyma.11G101400for rxn RXN_9663_p2 present\n",
      "Gene Glyma.02G015700for rxn FUMHYDR_RXN_m2 present\n",
      "Gene Glyma.02G015700for rxn FUMHYDR_RXN_c2 present\n",
      "Gene Glyma.11G248700for rxn UREASE_RXN_c2 present\n",
      "Gene Glyma.01G026700for rxn ALANINE_AMINOTRANSFERASE_RXN_c2 present\n",
      "Gene Glyma.17G159700for rxn 2_ISOPROPYLMALATESYN_RXN_p2 present\n",
      "Gene Glyma.14G217300for rxn 3_ISOPROPYLMALISOM_RXN_p2 present\n",
      "Gene Glyma.20G148100for rxn ACSERLY_RXN_p2 present\n",
      "Gene Glyma.11G101400for rxn RXN_9661_p2 present\n",
      "Gene Glyma.12G193800for rxn ACONITATEHYDR_RXN_m2 present\n",
      "Gene Glyma.08G072200for rxn RXN_6883_mc2 present\n",
      "Gene Glyma.03G027500for rxn 1TRANSKETO_RXN_p2 present\n",
      "Gene Glyma.11G158400for rxn IMPCYCLOHYDROLASE_RXN_p2 present\n",
      "Gene Glyma.03G121100for rxn MALTODEG_RXN_c2 present\n",
      "Gene Glyma.17G131300for rxn ARGINASE_RXN_m2 present\n",
      "Gene Glyma.19G239800for rxn 3_DEHYDROQUINATE_DEHYDRATASE_RXN_p2 present\n",
      "Gene Glyma.13G362900for rxn CYTOCHROME_C_OXIDASE_RXN_mi2 present\n",
      "Gene Glyma.17G101700for rxn RIBULP3EPIM_RXN_p2 present\n",
      "Gene Glyma.05G029200for rxn PYRROLINECARBDEHYDROG_RXN_NADP_m2 present\n",
      "Gene Glyma.17G261700for rxn CATAL_RXN_x2 present\n",
      "Gene Glyma.02G254000for rxn AMPSYN_RXN_p2 present\n",
      "Gene Glyma.13G336200for rxn TRIOSEPISOMERIZATION_RXN_c2 present\n",
      "Gene Glyma.12G193800for rxn ACONITATEDEHYDR_RXN_c2 present\n",
      "Gene Glyma.16G087200for rxn 1_PERIOD_10_PERIOD_2_PERIOD_2_RXN_mc2 present\n",
      "Gene Glyma.11G254700for rxn ADENOSYLHOMOCYSTEINASE_RXN_c2 present\n",
      "Gene Glyma.16G087200for rxn 1_PERIOD_10_PERIOD_2_PERIOD_2_RXN_mi2 present\n",
      "Gene Glyma.11G241900for rxn HISTCYCLOHYD_RXN_p2 present\n",
      "Gene Glyma.10G039000for rxn SUCCCOASYN_RXN_m2 present\n",
      "Gene Glyma.19G009000for rxn 1_PERIOD_2_PERIOD_1_PERIOD_2_RXN_m2 present\n",
      "Gene Glyma.18G265300for rxn ACETYL_COA_CARBOXYLTRANSFER_RXN_p2 present\n",
      "Gene Glyma.10G241000for rxn PHOSPHOGLUCMUT_RXN_p2 present\n",
      "Gene Glyma.03G002300for rxn RXN_9549_p2 present\n",
      "Gene Glyma.10G201100for rxn PEPDEPHOS_RXN_c2 present\n",
      "Gene Glyma.10G201100for rxn PEPDEPHOS_RXN_p2 present\n",
      "Gene Glyma.06G133500for rxn 6PGLUCONOLACT_RXN_p2 present\n",
      "Gene Glyma.09G270500for rxn 3PGAREARR_RXN_c2 present\n",
      "Gene Glyma.10G241000for rxn PHOSPHOGLUCMUT_RXN_c2 present\n",
      "Gene Glyma.05G005300for rxn ISOCITDEH_RXN_c2 present\n",
      "Gene Glyma.13G354900for rxn MALIC_NADP_RXN_p2 present\n",
      "Gene Glyma.11G078700for rxn 3_DEHYDROQUINATE_SYNTHASE_RXN_p2 present\n",
      "Gene Glyma.13G354900for rxn MALIC_NADP_RXN_c2 present\n",
      "Gene Glyma.06G185100for rxn SUCCINATE_DEHYDROGENASE_UBIQUINONE_RXN_mi2 present\n",
      "Gene Glyma.05G183700for rxn SHIKIMATE_KINASE_RXN_p2 present\n",
      "Gene Glyma.08G195000for rxn INORGPYROPHOSPHAT_RXN_x2 present\n",
      "Gene Glyma.17G100600for rxn MALATE_DEH_RXN_c2 present\n",
      "Gene Glyma.18G025500for rxn GLUC1PURIDYLTRANS_RXN_c2 present\n",
      "Gene Glyma.05G181000for rxn ASPAMINOTRANS_RXN_p2 present\n",
      "Gene Glyma.17G100600for rxn MALATE_DEH_RXN_p2 present\n",
      "Gene Glyma.17G100600for rxn MALATE_DEH_RXN_m2 present\n",
      "Gene Glyma.05G181000for rxn ASPAMINOTRANS_RXN_x2 present\n",
      "Gene Glyma.05G181000for rxn ASPAMINOTRANS_RXN_c2 present\n",
      "Gene Glyma.08G195000for rxn INORGPYROPHOSPHAT_RXN_p2 present\n",
      "Gene Glyma.17G100600for rxn MALATE_DEH_RXN_x2 present\n",
      "Gene Glyma.04G092400for rxn PGLUCISOM_RXN_c2 present\n",
      "Gene Glyma.04G092400for rxn PGLUCISOM_RXN_p2 present\n",
      "Gene Glyma.11G084100for rxn SERINE_O_ACETTRAN_RXN_c2 present\n",
      "Gene Glyma.17G109700for rxn SUCROSE_PHOSPHATE_SYNTHASE_RXN_c2 present\n",
      "Gene Glyma.11G084100for rxn SERINE_O_ACETTRAN_RXN_p2 present\n",
      "Gene Glyma.12G210000for rxn ACETOLACTREDUCTOISOM_RXN_p2 present\n",
      "Gene Glyma.02G249300for rxn RXN_1381_p2 present\n",
      "Gene Glyma.06G172600for rxn GAPOXNPHOSPHN_RXN_c2 present\n",
      "Gene Glyma.11G079000for rxn 3_ISOPROPYLMALDEHYDROG_RXN_p2 present\n",
      "Gene Glyma.19G088000for rxn ADENYL_KIN_RXN_p2 present\n",
      "Gene Glyma.18G093200for rxn GLUTAMIDOTRANS_RXN_p2 present\n",
      "Gene Glyma.06G172600for rxn GAPOXNPHOSPHN_RXN_p2 present\n",
      "Gene Glyma.05G056400for rxn TRANSALDOL_RXN_p2 present\n",
      "Gene Glyma.15G027000for rxn THRESYN_RXN_p2 present\n",
      "Gene Glyma.11G101400for rxn RXN_9662_p2 present\n",
      "Gene Glyma.04G009800for rxn HISTIDPHOS_RXN_p2 present\n",
      "Gene Glyma.20G183100for rxn RXN0_5260_m2 present\n",
      "Gene Glyma.11G226900for rxn SEDOHEPTULOSE_BISPHOSPHATASE_RXN_p2 present\n",
      "Gene Glyma.14G184500for rxn GLYOHMETRANS_RXN_m2 present\n",
      "Gene Glyma.08G334000for rxn RXN0_5184_c2 present\n",
      "Gene Glyma.14G184500for rxn GLYOHMETRANS_RXN_p2 present\n",
      "Gene Glyma.16G155800for rxn HISTAMINOTRANS_RXN_p2 present\n",
      "Gene Glyma.09G168900for rxn GLUTDECARBOX_RXN_c2 present\n",
      "Gene Glyma.11G241900for rxn HISTPRATPHYD_RXN_p2 present\n",
      "Gene Glyma.20G019700for rxn GCVMULTI_RXN_m2 present\n",
      "Gene Glyma.15G131700for rxn RXN_8001_p2 present\n",
      "Gene Glyma.01G082200for rxn SHIKIMATE_5_DEHYDROGENASE_RXN_p2 present\n",
      "Gene Glyma.08G309300for rxn PRPPSYN_RXN_p2 present\n",
      "Gene Glyma.02G244000for rxn GLUTAMINESYN_RXN_c2 present\n",
      "Gene Glyma.02G244000for rxn GLUTAMINESYN_RXN_m2 present\n",
      "Gene Glyma.04G256100for rxn 1_PERIOD_5_PERIOD_1_PERIOD_20_RXN_c2 present\n",
      "Gene Glyma.07G186000for rxn BRANCHED_CHAINAMINOTRANSFERLEU_RXN_p2 present\n",
      "Gene Glyma.13G207900for rxn DIHYDROXYISOVALDEHYDRAT_RXN_p2 present\n",
      "Gene Glyma.13G059400for rxn CHORISMATEMUT_RXN_p2 present\n",
      "Gene Glyma.02G072200for rxn GLUTAMATE_DEHYDROGENASE_RXN_m2 present\n",
      "Gene Glyma.10G086600for rxn SUCROSE_PHOSPHATASE_RXN_c2 present\n",
      "Gene Glyma.20G145900for rxn IMIDPHOSDEHYD_RXN_p2 present\n",
      "Gene Glyma.05G045200for rxn GDPKIN_RXN_p2 present\n",
      "Gene Glyma.03G212700for rxn METHENYLTHFCYCLOHYDRO_RXN_c2 present\n",
      "Gene Glyma.03G240100for rxn BRANCHED_CHAINAMINOTRANSFERVAL_RXN_p2 present\n",
      "Gene Glyma.19G177100for rxn ATPPHOSPHORIBOSYLTRANS_RXN_p2 present\n",
      "Gene Glyma.06G169700for rxn ACETOLACTSYN_RXN_p2 present\n",
      "Gene Glyma.09G217000for rxn RXN_1827_p2 present\n",
      "Gene Glyma.02G288500for rxn CITSYN_RXN_x2 present\n",
      "Gene Glyma.09G153900for rxn 2PGADEHYDRAT_RXN_c2 present\n",
      "Gene Glyma.14G176600for rxn DAHPSYN_RXN_p2 present\n",
      "Gene Glyma.18G261600for rxn CYSPH_RXN_p2 present\n",
      "Gene Glyma.01G139600for rxn 2_PERIOD_5_PERIOD_1_PERIOD_19_RXN_p2 present\n",
      "Gene Glyma.15G131700for rxn HISTOLDEHYD_RXN_p2 present\n",
      "Gene Glyma.02G288500for rxn CITSYN_RXN_m2 present\n",
      "Gene Glyma.11G158400for rxn AICARTRANSFORM_RXN_p2 present\n",
      "Gene Glyma.03G246300for rxn RIB5PISOM_RXN_p2 present\n",
      "Gene Glyma.10G128200for rxn 2_PERIOD_7_PERIOD_1_PERIOD_90_RXN_c2 present\n",
      "Gene Glyma.19G142400for rxn ADENYLOSUCCINATE_SYNTHASE_RXN_p2 present\n",
      "Gene Glyma.07G115900for rxn UDPKIN_RXN_c2 present\n",
      "Gene Glyma.01G031500for rxn RXN66_3_m2 present\n",
      "Gene Glyma.14G217300for rxn RXN_8991_p2 present\n",
      "Gene Glyma.08G113600for rxn ACETATE_COA_LIGASE_RXN_p2 present\n",
      "Gene Glyma.15G267000for rxn ATP_CITRATE_PRO_S_LYASE_RXN_c2 present\n",
      "Gene Glyma.16G063200for rxn GLU6PDEHYDROG_RXN_p2 present\n",
      "Gene Glyma.03G181900for rxn ADENOSINE_KINASE_RXN_c2 present\n",
      "Gene Glyma.20G121200for rxn PHOSPHATIDATE_PHOSPHATASE_RXN_p2 present\n",
      "Gene Glyma.09G167000for rxn SUCROSE_SYNTHASE_RXN_c2 present\n",
      "Gene Glyma.12G013400for rxn RXN_14351_pc2 present\n",
      "Gene Glyma.10G211400for rxn CHORISMATE_SYNTHASE_RXN_p2 present\n",
      "Gene Glyma.03G129700for rxn CYSTATHIONINE_BETA_LYASE_RXN_p2 present\n",
      "Gene Glyma.15G020900for rxn 6PFRUCTPHOS_RXN_p2 present\n",
      "Gene Glyma.13G153700for rxn PRIBFAICARPISOM_RXN_p2 present\n",
      "Gene Glyma.03G212700for rxn METHYLENETHFDEHYDROG_NADP_RXN_c2 present\n"
     ]
    }
   ],
   "source": [
    "ActiveRxnNotInTranscriptomicsData=set()\n",
    "for rxn in leaf_model.reactions:\n",
    "    if abs(rxn.flux) > 0 and len(rxn.genes)!=0:\n",
    "        flag=False\n",
    "        for gene in rxn.genes:\n",
    "            if gene.id in set(df[\"Glyma2.0\"]):\n",
    "                flag=True\n",
    "                print \"Gene \"+gene.id+\"for rxn \"+rxn.id+\" present\"\n",
    "                break\n",
    "        if not flag:\n",
    "            ActiveRxnNotInTranscriptomicsData.add(rxn.id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'1_ACYLGLYCEROL_3_P_ACYLTRANSFER_RXN_p2',\n",
       " 'GLUCOKIN_RXN_c2',\n",
       " 'GLY3KIN_RXN_p1',\n",
       " 'PHOSGLYPHOS_RXN_c1',\n",
       " 'PHOSGLYPHOS_RXN_c2',\n",
       " 'PHOSGLYPHOS_RXN_p1',\n",
       " 'PHOSGLYPHOS_RXN_p2',\n",
       " 'THREONINE_ALDOLASE_RXN_c1'}"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ActiveRxnNotInTranscriptomicsData"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.7630348276875"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_amb[df_amb[\"Glyma2.0\"]==\"Glyma.14G053300\"][df_amb.columns[17:]].mean().mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "totalTranscript_400 = dict()\n",
    "for rxn in leaf_model.reactions:\n",
    "    x=0\n",
    "    for gene in rxn.genes:\n",
    "        #print gene\n",
    "        y = df_amb[df_amb[\"Glyma2.0\"]==gene.id][df_amb.columns[17:]].mean().mean()\n",
    "        #print y\n",
    "        if pd.isna(y):\n",
    "            x = x + 0\n",
    "        else:\n",
    "            x = x + y\n",
    "    totalTranscript_400[rxn.id]=x\n",
    "    \n",
    "totalTranscript_550 = dict()\n",
    "for rxn in leaf_model.reactions:\n",
    "    x=0\n",
    "    for gene in rxn.genes:\n",
    "        #print gene\n",
    "        y = df_ele[df_ele[\"Glyma2.0\"]==gene.id][df_ele.columns[17:]].mean().mean()\n",
    "        #print y\n",
    "        if pd.isna(y):\n",
    "            x = x + 0\n",
    "        else:\n",
    "            x = x + y\n",
    "    totalTranscript_550[rxn.id]=x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INORGPYROPHOSPHAT_RXN_p1\n",
      "1.101657328803641\n",
      "1.08903172312\n",
      "1.0115934232333557\n",
      "NADH_DEHYDROG_A_RXN_mi2\n",
      "1.0046849867090821\n",
      "1.01675560787\n",
      "0.9881282964439743\n",
      "GAPOXNPHOSPHN_RXN_c2\n",
      "1.0193884551113048\n",
      "1.01191025142\n",
      "1.007390184735083\n",
      "PEPDEPHOS_RXN_c2\n",
      "1.0215797817779837\n",
      "1.01178551139\n",
      "1.0096801844678556\n",
      "FUMHYDR_RXN_m2\n",
      "0.9982047427477271\n",
      "1.00721935887\n",
      "0.9910499971563462\n",
      "TRIOSEPISOMERIZATION_RXN_p1\n",
      "1.0072391861039773\n",
      "1.01857028924\n",
      "0.9888754823771729\n",
      "RIB5PISOM_RXN_p2\n",
      "0.979800200721292\n",
      "0.994290325557\n",
      "0.9854266661725176\n",
      "PHOSPHOGLUCMUT_RXN_p2\n",
      "1.0003436901517928\n",
      "0.994290325557\n",
      "1.0060881258110699\n",
      "CITSYN_RXN_m2\n",
      "1.0367447835558694\n",
      "1.02202578301\n",
      "1.0144017898517175\n",
      "SUCCINATE_DEHYDROGENASE_UBIQUINONE_RXN_mi2\n",
      "1.0057316875838567\n",
      "1.00717613236\n",
      "0.9985658468948113\n",
      "TRANSALDOL_RXN_p2\n",
      "1.0012607681881907\n",
      "0.994223426223\n",
      "1.0070782298827117\n",
      "MALATE_DEH_RXN_c2\n",
      "1.0152293808837947\n",
      "1.01191025142\n",
      "1.0032800630897292\n",
      "MALATE_DEH_RXN_m2\n",
      "1.0152293808837947\n",
      "1.01737953002\n",
      "0.9978865810950349\n",
      "SUCCCOASYN_RXN_m2\n",
      "1.010030296025145\n",
      "1.00717613236\n",
      "1.0028338277406335\n",
      "CYTOCHROME_C_OXIDASE_RXN_mi2\n",
      "1.0130306629446855\n",
      "1.01600150255\n",
      "0.9970759495974426\n",
      "6PGLUCONOLACT_RXN_p2\n",
      "0.9881540574178065\n",
      "0.994245732325\n",
      "0.9938730690919103\n",
      "GLU6PDEHYDROG_RXN_p2\n",
      "1.0071355981789694\n",
      "0.994245732325\n",
      "1.012964466866148\n"
     ]
    }
   ],
   "source": [
    "rxnList = list()\n",
    "transcriptList = list()\n",
    "fluxList = list()\n",
    "import numpy as np\n",
    "\n",
    "for rxnID in totalTranscript_550.keys():\n",
    "    rxnList.append(rxnID)\n",
    "    a=0\n",
    "    b=0\n",
    "    if round(totalTranscript_400[rxnID],3)==0 and round(totalTranscript_550[rxnID],3)==0:\n",
    "        transcriptList.append(0)\n",
    "    elif round(totalTranscript_400[rxnID],3)==0:\n",
    "        transcriptList.append(np.nan)\n",
    "    else:\n",
    "        a = totalTranscript_550[rxnID]/totalTranscript_400[rxnID]\n",
    "        transcriptList.append(a)\n",
    "    if round(leaf_model_550.reactions.get_by_id(rxnID).flux,3)==0 and round(leaf_model_400.reactions.get_by_id(rxnID).flux,3)==0:\n",
    "        fluxList.append(0)\n",
    "    elif round(leaf_model_400.reactions.get_by_id(rxnID).flux,3)==0:\n",
    "        fluxList.append(np.nan)\n",
    "    else:\n",
    "        b = leaf_model_550.reactions.get_by_id(rxnID).flux/leaf_model_400.reactions.get_by_id(rxnID).flux\n",
    "        fluxList.append(b)\n",
    "    if not (a==0 and b==0):\n",
    "        if b == 0:\n",
    "            continue\n",
    "        if round(a/b,2)>0.98 and round(a/b,2)<1.02:\n",
    "            print rxnID\n",
    "            print a\n",
    "            print b\n",
    "            print a/b\n",
    "\n",
    "df_ratios = pd.DataFrame(data={\"rxn ID\":rxnList,\"Transcript ratio\":transcriptList,\"Flux ratio\":fluxList})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "111.53797876447916"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "totalTranscript_400[\"INORGPYROPHOSPHAT_RXN_p1\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "leaf_model.reactions.L_CITRULLINE_pc1.genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "leaf_model_sink = model.copy()\n",
    "\n",
    "leaf_model_sink = setupC3DielModel(leaf_model_sink,\"MetabolitesToTransfer.txt\",starch_sucrose_ratio=0.6)\n",
    "\n",
    "met = Metabolite(\"X_Biomass_contribution_t1\")\n",
    "leaf_model_sink.reactions.Biomass_leaf_tx1.add_metabolites({met:1})\n",
    "leaf_model_sink.metabolites.X_Phloem_contribution_t1.remove_from_model()\n",
    "leaf_model_sink.reactions.diel_biomass.add_metabolites({met:-3})\n",
    "\n",
    "met = Metabolite(\"X_Biomass_contribution_t2\")\n",
    "leaf_model_sink.reactions.Biomass_leaf_tx2.add_metabolites({met:1})\n",
    "leaf_model_sink.metabolites.X_Phloem_contribution_t2.remove_from_model()\n",
    "leaf_model_sink.reactions.diel_biomass.add_metabolites({met:-1})\n",
    "\n",
    "leaf_model_sink.reactions.diel_biomass.objective_coefficient=1\n",
    "flux_analysis.parsimonious.pfba(leaf_model_sink)\n",
    "print(\"Biomass accumulation = \")\n",
    "print(leaf_model_sink.reactions.diel_biomass.flux)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}