creat_compm_media.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import cobra\n",
    "import json"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_dict(data, name):\n",
    "    with open(name, 'w' ) as file:\n",
    "        json.dump( data, file )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating COMPM media\n",
    "\n",
    "We compute the competition inducing media by using the minimal minimum flux of all common exchange reactions. First we load all required stuff."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Growth:  0.2798991966652624\n",
      "Growth:  2.5586946126133867\n"
     ]
    }
   ],
   "source": [
    "model_DP = cobra.io.read_sbml_model(\"models/consistent_DP_SNM.xml\")\n",
    "model_SA = cobra.io.read_sbml_model(\"models/consistent_iYS854_SNM.xml\")\n",
    "print(\"Growth: \", model_DP.slim_optimize())\n",
    "print(\"Growth: \", model_SA.slim_optimize())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Compound</th>\n",
       "      <th>BiGG</th>\n",
       "      <th>ModelSeed</th>\n",
       "      <th>KEGG</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Alanine</td>\n",
       "      <td>ala__L</td>\n",
       "      <td>cpd00035</td>\n",
       "      <td>C00041</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Arginine</td>\n",
       "      <td>arg__L</td>\n",
       "      <td>cpd00051</td>\n",
       "      <td>C00062</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Cysteine</td>\n",
       "      <td>cys__L</td>\n",
       "      <td>cpd00084</td>\n",
       "      <td>C00097</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Glutamic acid</td>\n",
       "      <td>glu__L</td>\n",
       "      <td>cpd00023</td>\n",
       "      <td>C00025</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Glycine</td>\n",
       "      <td>gly</td>\n",
       "      <td>cpd00033</td>\n",
       "      <td>C00037</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        Compound    BiGG ModelSeed    KEGG\n",
       "0        Alanine  ala__L  cpd00035  C00041\n",
       "1       Arginine  arg__L  cpd00051  C00062\n",
       "2       Cysteine  cys__L  cpd00084  C00097\n",
       "3  Glutamic acid  glu__L  cpd00023  C00025\n",
       "4        Glycine     gly  cpd00033  C00037"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "snm3 = pd.read_csv(\"SNM3.csv\", sep =\";\")\n",
    "snm3.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of united exchange reactions:  111\n"
     ]
    }
   ],
   "source": [
    "# Union of exchange reactions\n",
    "ex_DP_ids = set(map(lambda x:x.id, model_DP.exchanges))\n",
    "ex_SA_ids = set(map(lambda x:x.id, model_SA.exchanges))\n",
    "ex_union_ids = ex_DP_ids.union(ex_SA_ids)\n",
    "print(\"Number of united exchange reactions: \", len(ex_union_ids))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "fva_DP = cobra.flux_analysis.flux_variability_analysis(model_DP)\n",
    "fva_SA = cobra.flux_analysis.flux_variability_analysis(model_SA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Computing COMPM media\n",
    "We compute the compm media, as following:\n",
    "* If the exchange reaction is common between the species, we take the minimum flux.\n",
    "* If the exchange reaction is not common, we take the minFVA value of the specific species \n",
    "\n",
    "Note that we make **two speperate media**, for each species one, as a **medium can only contain exchange reactions that are present in the model**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def compute_COMPM(fva_DP, fva_SA):\n",
    "    media_COMPM_DP = dict()\n",
    "    media_COMPM_SA = dict()\n",
    "    for ex in ex_union_ids:\n",
    "        if ex in fva_DP.index and ex in fva_SA.index:\n",
    "            # Take minimum of both minima if common\n",
    "            flux = min(fva_DP.loc[ex].minimum, fva_SA.loc[ex].minimum)\n",
    "            if flux < 0:\n",
    "                    media_COMPM_DP[ex] = abs(flux)\n",
    "                    media_COMPM_SA[ex] = abs(flux)\n",
    "        elif ex in fva_DP.index:\n",
    "            # Take just the minimum if not common\n",
    "            flux = fva_DP.loc[ex].minimum\n",
    "            if flux < 0:\n",
    "                    media_COMPM_DP[ex] = abs(flux)\n",
    "        elif ex in fva_SA.index:\n",
    "            # Take just the minimum if not common\n",
    "            flux = fva_SA.loc[ex].minimum\n",
    "            if flux < 0:\n",
    "                    media_COMPM_SA[ex] = abs(flux)\n",
    "        else:\n",
    "            print(\"error\")\n",
    "    return media_COMPM_DP, media_COMPM_SA       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "media_COMPM_DP, media_COMPM_SA = compute_COMPM(fva_DP,fva_SA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length of media:  36\n",
      "{'EX_cobalt2_e': 0.009894648965030426, 'EX_thr__L_e': 0.14751502550566087, 'EX_o2_e': 20.0, 'EX_so4_e': 0.010557173971642781, 'EX_zn2_e': 0.008915031973978057, 'EX_gly_e': 10.0, 'EX_nac_e': 7.953461359822521, 'EX_fum_e': 10.0, 'EX_arg__L_e': 10.0, 'EX_trp__L_e': 0.03925018692530497, 'EX_phe__L_e': 10.0, 'EX_ca2_e': 0.014549711222763873, 'EX_cys__L_e': 10.0, 'EX_his__L_e': 0.1541538579498835, 'EX_h2o_e': 10.0, 'EX_pi_e': 5.519506278294172, 'EX_cl_e': 0.016447785134853243, 'EX_mn2_e': 0.010614197937497382, 'EX_ni2_e': 0.009935074921301279, 'EX_mg2_e': 0.023991908100634878, 'EX_k_e': 0.47512400263709553, 'EX_ala__L_e': 10.0, 'EX_cu2_e': 0.009176397041292562, 'EX_orn_e': 10.0, 'EX_glc__D_e': 10.0, 'EX_leu__L_e': 0.8034579716308596, 'EX_fe2_e': 0.037761215092602184, 'EX_glu__L_e': 10.0, 'EX_na1_e': 5.090522931794308, 'EX_4abz_e': 0.00018133893906634825, 'EX_lys__L_e': 0.9223436956285925, 'EX_co2_e': 0.13166473758769873, 'EX_pro__L_e': 0.05991910190673909, 'EX_thm_e': 0.0011411777972255722, 'EX_ribflv_e': 0.0007637848411345633, 'EX_val__L_e': 10.0}\n"
     ]
    }
   ],
   "source": [
    "print(\"Length of media: \", len(media_COMPM_DP))\n",
    "print(media_COMPM_DP)\n",
    "save_dict(media_COMPM_DP, \"compm_DP.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length of media:  30\n",
      "{'EX_o2_e': 20.0, 'EX_so4_e': 0.010557173971642781, 'EX_zn2_e': 0.008915031973978057, 'EX_nac_e': 7.953461359822521, 'EX_arg__L_e': 10.0, 'EX_trp__L_e': 0.03925018692530497, 'EX_phe__L_e': 10.0, 'EX_cys__L_e': 10.0, 'EX_his__L_e': 0.1541538579498835, 'EX_h2o_e': 10.0, 'EX_pi_e': 5.519506278294172, 'EX_cl_e': 0.016447785134853243, 'EX_mn2_e': 0.010614197937497382, 'EX_mg2_e': 0.023991908100634878, 'EX_pyr_e': 10.0, 'EX_k_e': 0.47512400263709553, 'EX_orn_e': 10.0, 'EX_glc__D_e': 10.0, 'EX_leu__L_e': 0.8034579716308596, 'EX_urea_e': 10.0, 'EX_fe2_e': 0.037761215092602184, 'EX_glu__L_e': 10.0, 'EX_na1_e': 5.090522931794308, 'EX_ser__L_e': 10.0, 'EX_mobd_e': 1.7910862288293614e-05, 'EX_lys__L_e': 0.9223436956285925, 'EX_co2_e': 0.13166473758769873, 'EX_thm_e': 0.0011411777972255722, 'EX_ribflv_e': 0.0007637848411345633, 'EX_val__L_e': 10.0}\n"
     ]
    }
   ],
   "source": [
    "print(\"Length of media: \", len(media_COMPM_SA))\n",
    "print(media_COMPM_SA)\n",
    "save_dict(media_COMPM_SA, \"compm_SA.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SA growth:  2.5586946126133925\n",
      "DP growth:  0.2798991966652633\n"
     ]
    }
   ],
   "source": [
    "# Check for growth\n",
    "with model_SA:\n",
    "    model_SA.medium=media_COMPM_SA\n",
    "    print(\"SA growth: \",model_SA.slim_optimize())\n",
    "with model_DP:\n",
    "    model_DP.medium=media_COMPM_DP\n",
    "    print(\"DP growth: \",model_DP.slim_optimize())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we have optimal growth for both species, as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check that it is valid snm3\n",
    "for id,val in media_COMPM_DP.items():\n",
    "    idx = str(id[3:-2])\n",
    "    assert idx in list(snm3[\"BiGG\"])\n",
    "for id,val in media_COMPM_SA.items():\n",
    "    idx = str(id[3:-2])\n",
    "    assert idx in list(snm3[\"BiGG\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length of media:  36\n",
      "Length of media:  30\n"
     ]
    }
   ],
   "source": [
    "print(\"Length of media: \", len(media_COMPM_DP))\n",
    "print(\"Length of media: \", len(media_COMPM_SA))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) FVA bounds with variable thresholds:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1) FVA = 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "fva_DP_10 = cobra.flux_analysis.flux_variability_analysis(model_DP, fraction_of_optimum=0.1)\n",
    "fva_SA_10 = cobra.flux_analysis.flux_variability_analysis(model_SA, fraction_of_optimum=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "fva_DP_10.to_excel(\"fva_DP_010.xlsx\")\n",
    "fva_SA_10.to_excel(\"fva_SA_010.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SA growth:  2.4322218392467847\n",
      "DP growth:  0.2798991966652633\n"
     ]
    }
   ],
   "source": [
    "media_COMPM_DP_10, media_COMPM_SA_10 = compute_COMPM(fva_DP_10, fva_SA_10)\n",
    "# Check for growth\n",
    "with model_SA:\n",
    "    model_SA.medium=media_COMPM_SA_10\n",
    "    print(\"SA growth: \",model_SA.slim_optimize())\n",
    "with model_DP:\n",
    "    model_DP.medium=media_COMPM_DP_10\n",
    "    print(\"DP growth: \",model_DP.slim_optimize())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2) FVA = 0.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "fva_DP_40 = cobra.flux_analysis.flux_variability_analysis(model_DP, fraction_of_optimum=0.4)\n",
    "fva_SA_40 = cobra.flux_analysis.flux_variability_analysis(model_SA, fraction_of_optimum=0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "fva_DP_40.to_excel(\"fva_DP_040.xlsx\")\n",
    "fva_SA_40.to_excel(\"fva_SA_040.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SA growth:  2.432143714435963\n",
      "DP growth:  0.2798991966652633\n"
     ]
    }
   ],
   "source": [
    "media_COMPM_DP_40, media_COMPM_SA_40 = compute_COMPM(fva_DP_40, fva_SA_40)\n",
    "# Check for growth\n",
    "with model_SA:\n",
    "    model_SA.medium=media_COMPM_SA_40\n",
    "    print(\"SA growth: \",model_SA.slim_optimize())\n",
    "with model_DP:\n",
    "    model_DP.medium=media_COMPM_DP_40\n",
    "    print(\"DP growth: \",model_DP.slim_optimize())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3) FVA = 0.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "fva_DP_80 = cobra.flux_analysis.flux_variability_analysis(model_DP, fraction_of_optimum=0.8)\n",
    "fva_SA_80 = cobra.flux_analysis.flux_variability_analysis(model_SA, fraction_of_optimum=0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "fva_DP_80.to_excel(\"fva_DP_080.xlsx\")\n",
    "fva_SA_80.to_excel(\"fva_SA_080.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SA growth:  2.432221839246782\n",
      "DP growth:  0.2798991966652633\n"
     ]
    }
   ],
   "source": [
    "media_COMPM_DP_80, media_COMPM_SA_80 = compute_COMPM(fva_DP_80, fva_SA_80)\n",
    "# Check for growth\n",
    "with model_SA:\n",
    "    model_SA.medium=media_COMPM_SA_80\n",
    "    print(\"SA growth: \",model_SA.slim_optimize())\n",
    "with model_DP:\n",
    "    model_DP.medium=media_COMPM_DP_80\n",
    "    print(\"DP growth: \",model_DP.slim_optimize())"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}