gurry_exvivo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a586201d-912a-4293-ad64-fb46960dfb3b",
   "metadata": {},
   "source": [
    "## This notebook shows the full workflow for building models, simulating growth and obtaining SCFA predictions from data collected by the _ex vivo_ study conducted by Gurry et al. 2021"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "f160f9ad-82b0-425a-affe-af74a1e0ba20",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np \n",
    "import micom as mm\n",
    "from plotnine import *\n",
    "import os\n",
    "import sys\n",
    "\n",
    "pd.options.mode.chained_assignment = None  # default='warn'\n",
    "\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30c83a66-2531-40fc-b0ed-c9e12e539e30",
   "metadata": {},
   "source": [
    "## First, we pull in the taxonomy table, matching each feature ID in the qiime2 output to a microbial taxa at the species level. We will build our models at the genus level, so collapse to this rank"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2335aca8-c11d-431c-94b6-12326dc3df96",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir('/proj/gibbons/nbohmann/exvivo/gurry1/data/qiime2/taxonomy/data')\n",
    "taxa = pd.read_csv('taxonomy.tsv',sep='\\t') # read table\n",
    "taxa.set_index('Feature ID',inplace=True)\n",
    "taxa = taxa.Taxon.str.split(';',expand=True) # split ranks\n",
    "taxa = taxa.rename(columns = {0:'Kingdom',1:'Phylum',2:'Class',3:'Order',4:'Family',5:'Genus',6:'Species'})\n",
    "taxa = taxa.dropna(subset = ['Genus']) # drop undefined columns\n",
    "taxa = taxa.drop(columns = ['Species']) # drop species column\n",
    "taxa = taxa.apply(lambda column: column.str.split('_').str[2]) # remove prefixes\n",
    "taxa = (taxa.apply(lambda row: \";\".join(row.str.capitalize().fillna(\"\")), axis=1)\n",
    "        .to_frame().rename(columns = {0:'taxon'})) # join columns into taxon identifier\n",
    "taxa"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e617ae1-e972-4cec-8ca4-7b6b3a4dfc6b",
   "metadata": {},
   "source": [
    "## Next we'll pull in the abundance table, with read counts for all present taxa. We'll drop those that aren't identified in the taxon list, and sum together duplicates. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8a72aa5-10db-4fbf-80d6-6b825a9a791d",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir('/proj/gibbons/nbohmann/exvivo/gurry1/data/qiime2/')\n",
    "unrarefied_table = q2.Artifact.load('table.qza') # read table\n",
    "abundance = unrarefied_table.view(pd.DataFrame).reset_index().rename(columns = {'index':'sample_id'})\n",
    "to_drop = (abundance[abundance.columns[1:]].columns[~abundance[abundance.columns[1:]].\n",
    "                                                    columns.isin(taxa.index)].to_list()) # taxa to drop\n",
    "abundance = abundance.drop(columns = to_drop)\n",
    "abundance = abundance.rename(columns = taxa['taxon'].to_dict())\n",
    "abundance = abundance.groupby(axis=1, level=0).sum() \n",
    "abundance = pd.melt(abundance, id_vars = 'sample_id', value_vars = abundance.columns[:-1],\n",
    "                    var_name = 'id', value_name = 'abundance') # melt into long form df \n",
    "abundance['genus'] = abundance['id'].str.split(';').str[-1] # need a genus column in df \n",
    "abundance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bd09f44-0488-4013-a28a-73ceea983b1b",
   "metadata": {},
   "source": [
    "## We also need a model database to pull our reconstructions from "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdce81ae-c0c7-41d3-900f-bb251e6a901e",
   "metadata": {},
   "outputs": [],
   "source": [
    "agora = ('/proj/gibbons/refs/micom_dbs/agora103_genus.qza')\n",
    "agora"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f3aca59-13a9-4a49-8304-78e6a2756f8f",
   "metadata": {},
   "source": [
    "## Now, we'll build our models, with cutoff of 0.001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e00c841b-17d9-44f7-8140-2bd3dd010cb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "models = mm.workflows.build(abundance,out_folder = '/proj/gibbons/nbohmann/exvivo/gurry1/micom/16S/16S_models/',\n",
    "                      model_db = agora, cutoff = 0.001, threads = 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c26d844-e526-49fd-b8d4-86c6c2f83d75",
   "metadata": {},
   "source": [
    "## We can take a look at the resulting model manifest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "260605d1-cf53-485f-87bf-bbcbd5ab8616",
   "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>sample_id</th>\n",
       "      <th>file</th>\n",
       "      <th>found_taxa</th>\n",
       "      <th>total_taxa</th>\n",
       "      <th>found_fraction</th>\n",
       "      <th>found_abundance_fraction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>H008</td>\n",
       "      <td>H008.pickle</td>\n",
       "      <td>30.0</td>\n",
       "      <td>50.0</td>\n",
       "      <td>0.600000</td>\n",
       "      <td>0.791525</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>H009</td>\n",
       "      <td>H009.pickle</td>\n",
       "      <td>33.0</td>\n",
       "      <td>44.0</td>\n",
       "      <td>0.750000</td>\n",
       "      <td>0.908978</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>H010</td>\n",
       "      <td>H010.pickle</td>\n",
       "      <td>38.0</td>\n",
       "      <td>56.0</td>\n",
       "      <td>0.678571</td>\n",
       "      <td>0.799867</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>H012</td>\n",
       "      <td>H012.pickle</td>\n",
       "      <td>28.0</td>\n",
       "      <td>43.0</td>\n",
       "      <td>0.651163</td>\n",
       "      <td>0.868381</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>H019</td>\n",
       "      <td>H019.pickle</td>\n",
       "      <td>22.0</td>\n",
       "      <td>33.0</td>\n",
       "      <td>0.666667</td>\n",
       "      <td>0.812679</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>H020</td>\n",
       "      <td>H020.pickle</td>\n",
       "      <td>26.0</td>\n",
       "      <td>34.0</td>\n",
       "      <td>0.764706</td>\n",
       "      <td>0.839140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>H021</td>\n",
       "      <td>H021.pickle</td>\n",
       "      <td>28.0</td>\n",
       "      <td>38.0</td>\n",
       "      <td>0.736842</td>\n",
       "      <td>0.840889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>H025</td>\n",
       "      <td>H025.pickle</td>\n",
       "      <td>27.0</td>\n",
       "      <td>43.0</td>\n",
       "      <td>0.627907</td>\n",
       "      <td>0.894737</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>H028</td>\n",
       "      <td>H028.pickle</td>\n",
       "      <td>20.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>0.800000</td>\n",
       "      <td>0.945193</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>H029</td>\n",
       "      <td>H029.pickle</td>\n",
       "      <td>37.0</td>\n",
       "      <td>58.0</td>\n",
       "      <td>0.637931</td>\n",
       "      <td>0.851163</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  sample_id         file  found_taxa  total_taxa  found_fraction  \\\n",
       "0      H008  H008.pickle        30.0        50.0        0.600000   \n",
       "1      H009  H009.pickle        33.0        44.0        0.750000   \n",
       "2      H010  H010.pickle        38.0        56.0        0.678571   \n",
       "3      H012  H012.pickle        28.0        43.0        0.651163   \n",
       "4      H019  H019.pickle        22.0        33.0        0.666667   \n",
       "5      H020  H020.pickle        26.0        34.0        0.764706   \n",
       "6      H021  H021.pickle        28.0        38.0        0.736842   \n",
       "7      H025  H025.pickle        27.0        43.0        0.627907   \n",
       "8      H028  H028.pickle        20.0        25.0        0.800000   \n",
       "9      H029  H029.pickle        37.0        58.0        0.637931   \n",
       "\n",
       "   found_abundance_fraction  \n",
       "0                  0.791525  \n",
       "1                  0.908978  \n",
       "2                  0.799867  \n",
       "3                  0.868381  \n",
       "4                  0.812679  \n",
       "5                  0.839140  \n",
       "6                  0.840889  \n",
       "7                  0.894737  \n",
       "8                  0.945193  \n",
       "9                  0.851163  "
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.chdir('/proj/gibbons/nbohmann/exvivo/gurry1/micom/16S/')\n",
    "manifest = pd.read_csv('16S_models/manifest.csv')\n",
    "manifest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "a64d4151",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3416f36f053f4357afc1295787ec2a21",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">[12:52:24] </span><span style=\"color: #800000; text-decoration-color: #800000\">WARNING </span> solver encountered an error infeasible                    <a href=\"file:///users/nbohmann/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/solution.py\" target=\"_blank\"><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">solution.py</span></a><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">:</span><a href=\"file:///users/nbohmann/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/solution.py#220\" target=\"_blank\"><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">220</span></a>\n",
       "</pre>\n"
      ],
      "text/plain": [
       "\u001b[2;36m[12:52:24]\u001b[0m\u001b[2;36m \u001b[0m\u001b[31mWARNING \u001b[0m solver encountered an error infeasible                    \u001b]8;id=151788;file:///users/nbohmann/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/solution.py\u001b\\\u001b[2msolution.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=431360;file:///users/nbohmann/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/solution.py#220\u001b\\\u001b[2m220\u001b[0m\u001b]8;;\u001b\\\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">           </span><span style=\"color: #800000; text-decoration-color: #800000\">WARNING </span> minimization of medium was unsuccessful                      <a href=\"file:///users/nbohmann/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/media.py\" target=\"_blank\"><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">media.py</span></a><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">:</span><a href=\"file:///users/nbohmann/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/media.py#236\" target=\"_blank\"><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">236</span></a>\n",
       "</pre>\n"
      ],
      "text/plain": [
       "\u001b[2;36m          \u001b[0m\u001b[2;36m \u001b[0m\u001b[31mWARNING \u001b[0m minimization of medium was unsuccessful                      \u001b]8;id=495143;file:///users/nbohmann/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/media.py\u001b\\\u001b[2mmedia.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=119359;file:///users/nbohmann/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/media.py#236\u001b\\\u001b[2m236\u001b[0m\u001b]8;;\u001b\\\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"color: #7fbfbf; text-decoration-color: #7fbfbf\">           </span><span style=\"color: #800000; text-decoration-color: #800000; font-weight: bold\">ERROR   </span> Could not get a minimal medium for sample H029.               <a href=\"file:///users/nbohmann/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/workflows/media.py\" target=\"_blank\"><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">media.py</span></a><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">:</span><a href=\"file:///users/nbohmann/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/workflows/media.py#36\" target=\"_blank\"><span style=\"color: #7f7f7f; text-decoration-color: #7f7f7f\">36</span></a>\n",
       "</pre>\n"
      ],
      "text/plain": [
       "\u001b[2;36m          \u001b[0m\u001b[2;36m \u001b[0m\u001b[1;31mERROR   \u001b[0m Could not get a minimal medium for sample H029.               \u001b]8;id=512130;file:///users/nbohmann/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/workflows/media.py\u001b\\\u001b[2mmedia.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=459068;file:///users/nbohmann/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/workflows/media.py#36\u001b\\\u001b[2m36\u001b[0m\u001b]8;;\u001b\\\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">\n",
       "</pre>\n"
      ],
      "text/plain": [
       "\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "ename": "OptimizationError",
     "evalue": "Could not find a growth medium that allows the specified growth rate for all taxa in all samples :(",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mOptimizationError\u001b[0m                         Traceback (most recent call last)",
      "Input \u001b[0;32mIn [66]\u001b[0m, in \u001b[0;36m<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0m min_med \u001b[38;5;241m=\u001b[39m \u001b[43mmm\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mworkflows\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmedia\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mminimal_media\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmanifest\u001b[49m\u001b[43m,\u001b[49m\u001b[43mmodel_folder\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43m16S_models\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msummarize\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmin_growth\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m0.1\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/miniconda3/envs/micom2022/lib/python3.8/site-packages/micom/workflows/media.py:61\u001b[0m, in \u001b[0;36mminimal_media\u001b[0;34m(manifest, model_folder, summarize, min_growth, threads)\u001b[0m\n\u001b[1;32m     59\u001b[0m results \u001b[38;5;241m=\u001b[39m workflow(_medium, args, threads)\n\u001b[1;32m     60\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28many\u001b[39m(r \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;28;01mfor\u001b[39;00m r \u001b[38;5;129;01min\u001b[39;00m results):\n\u001b[0;32m---> 61\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m OptimizationError(\n\u001b[1;32m     62\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mCould not find a growth medium that allows the specified \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m     63\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mgrowth rate for all taxa in all samples :(\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m     64\u001b[0m     )\n\u001b[1;32m     65\u001b[0m results \u001b[38;5;241m=\u001b[39m pd\u001b[38;5;241m.\u001b[39mconcat(results, axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m)\n\u001b[1;32m     66\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m summarize:\n",
      "\u001b[0;31mOptimizationError\u001b[0m: Could not find a growth medium that allows the specified growth rate for all taxa in all samples :("
     ]
    }
   ],
   "source": [
    "min_med = mm.workflows.media.minimal_media(manifest,model_folder = '16S_models', summarize = True, min_growth = 0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bdf43f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_med"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "359421f5-0215-413c-b4d9-77f5267dfdcd",
   "metadata": {},
   "source": [
    "## Load in the minimal media from the previous ex vivos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97172f5e-e989-47e0-8e14-8d4d7db38303",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir('/proj/gibbons/nbohmann/exvivo/')\n",
    "medium = pd.read_csv('minimal_media.csv',index_col =0)\n",
    "os.chdir('/proj/gibbons/nbohmann/exvivo/gurry1/micom/16S/')\n",
    "manifest = pd.read_csv('16S_models/manifest.csv')\n",
    "medium = mm.workflows.fix_medium(manifest,'16S_models',medium,community_growth=0.1,min_growth = 0.001,\n",
    "                                     minimize_components=False, weights='C',summarize=True)\n",
    "pectin_medium = medium.append({'reaction':'EX_pect_m','flux':300.00},ignore_index=True)\n",
    "inulin_medium = (medium[~medium.reaction.str.contains('inulin')].\n",
    "                 append({'reaction':'EX_inulin_m','flux':10.00},ignore_index=True))\n",
    "fos_medium = medium.append({'reaction':'EX_kesto_m','flux':100.00},ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5683ae7b-7a7a-46f9-844e-b45665695394",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir('/proj/gibbons/nbohmann/exvivo/gurry1/micom/16S/')\n",
    "tradeoff = mm.workflows.tradeoff(manifest, model_folder='16S_models',medium = fos_medium, \n",
    "                                tradeoffs = [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0],\n",
    "                                 presolve = True, threads = 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c676f98-b9f2-4bd4-870a-b3102e16188b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from micom.viz import plot_tradeoff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffba6a2e-98d4-4c27-a719-625121612e66",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir('/users/nbohmann/ex_vivo/gurry1')\n",
    "plot_tradeoff(tradeoff, filename='tradeoff_fos.html')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9afb1a7-37cf-4e59-a853-50339b6336cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir('/proj/gibbons/nbohmann/exvivo/gurry1/micom/16S/')\n",
    "ctrl_growth = mm.workflows.grow(manifest, model_folder='16S_models',medium = medium, \n",
    "                                tradeoff = 0.6, strategy = 'none', threads = 5)\n",
    "pectin_growth = mm.workflows.grow(manifest, model_folder='16S_models',medium = pectin_medium, \n",
    "                                   tradeoff = 0.6, strategy = 'none', threads = 5)\n",
    "inulin_growth = mm.workflows.grow(manifest, model_folder='16S_models/',medium = inulin_medium, \n",
    "                                   tradeoff = 0.6, strategy = 'none', threads = 5)\n",
    "fos_growth = mm.workflows.grow(manifest, model_folder='16S_models',medium = fos_medium, \n",
    "                                   tradeoff = 0.6, strategy = 'none', threads = 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eaa20a13-8292-4ed8-bf85-c324e56588bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_fluxes(growth, cond):\n",
    "    growth = growth.exchanges\n",
    "    growth = growth[growth.direction == \"export\"].groupby([\"sample_id\", \"metabolite\", \"reaction\"]).apply(lambda df: sum(df.flux * df.abundance)).reset_index()\n",
    "    growth = growth[(growth.reaction.str.startswith('EX_but(e)'))|(growth.reaction.str.startswith('EX_ppa(e)'))]\n",
    "    growth['index'] = growth['sample_id']+'_'+cond\n",
    "    return growth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a11a67c6-e3c1-4589-939c-dcf9b487ab83",
   "metadata": {},
   "outputs": [],
   "source": [
    "predicted = pd.DataFrame()\n",
    "predicted = pd.concat(predicted, get_fluxes(ctrl_growth, 'CTRL'))\n",
    "predicted = pd.concat(predicted, get_fluxes(inulin_growth, 'INUL'))\n",
    "predicted = pd.concat(predicted, get_fluxes(pectin_growth, 'PECT'))\n",
    "predicted = pd.concat(predicted, get_fluxes(fos_growth, 'FOS'))\n",
    "predicted = pd.pivot_table(predicted, index = 'index',columns = 'reaction', values = 0)\n",
    "but = predicted['EX_but(e)'].to_dict()\n",
    "ppa = predicted['EX_ppa(e)'].to_dict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a2310e2-7777-497e-af6f-a392e9f5f8bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "ctrl = ctrl_growth.exchanges\n",
    "ctrl = ctrl[ctrl.direction == \"export\"].groupby([\"sample_id\", \"metabolite\", \"reaction\"]).apply(lambda df: sum(df.flux * df.abundance)).reset_index()\n",
    "ctrl = ctrl[(ctrl.reaction.str.startswith('EX_but(e)'))|(ctrl.reaction.str.startswith('EX_ppa(e)'))]\n",
    "ctrl['index'] = ctrl['sample_id']+'_CTRL'\n",
    "pectin = pectin_growth.exchanges\n",
    "pectin = pectin[pectin.direction == \"export\"].groupby([\"sample_id\", \"metabolite\", \"reaction\"]).apply(lambda df: sum(df.flux * df.abundance)).reset_index()\n",
    "pectin = pectin[(pectin.reaction.str.startswith('EX_but(e)'))|(pectin.reaction.str.startswith('EX_ppa(e)'))]\n",
    "pectin['index'] = pectin['sample_id']+'_PECT'\n",
    "inulin = inulin_growth.exchanges\n",
    "inulin = inulin[inulin.direction == \"export\"].groupby([\"sample_id\", \"metabolite\", \"reaction\"]).apply(lambda df: sum(df.flux * df.abundance)).reset_index()\n",
    "inulin = inulin[(inulin.reaction.str.startswith('EX_but(e)'))|(inulin.reaction.str.startswith('EX_ppa(e)'))]\n",
    "inulin['index'] = inulin['sample_id']+'_INUL'\n",
    "fos = fos_growth.exchanges\n",
    "fos = fos[fos.direction == \"export\"].groupby([\"sample_id\", \"metabolite\", \"reaction\"]).apply(lambda df: sum(df.flux * df.abundance)).reset_index()\n",
    "fos = fos[(fos.reaction.str.startswith('EX_but(e)'))|(fos.reaction.str.startswith('EX_ppa(e)'))]\n",
    "fos['index'] = fos['sample_id']+'_FOS'\n",
    "predicted = ctrl.append(pectin,ignore_index = True)\n",
    "predicted = predicted.append(inulin,ignore_index = True)\n",
    "predicted = predicted.append(fos,ignore_index = True)\n",
    "predicted = pd.pivot_table(predicted, index = 'index',columns = 'reaction', values = 0)\n",
    "but = predicted['EX_but(e)'].to_dict()\n",
    "ppa = predicted['EX_ppa(e)'].to_dict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8b56cf5-92a1-44db-abc4-8e68c8b19592",
   "metadata": {},
   "outputs": [],
   "source": [
    "def flux_calculate(arg):\n",
    "    os.chdir('/proj/gibbons/nbohmann/exvivo/gurry1/data/gc_data/')\n",
    "    file = pd.read_csv(arg,index_col = 0)\n",
    "    file = file[['but','ppa']].dropna()\n",
    "    file = file[~file.index.str.contains(\"QC\")]\n",
    "    file['sample'] = file.index.str.split('-').str[0]\n",
    "    file['treatment'] = file.index.str.split('-').str[1]\n",
    "    file['timepoint'] = file.index.str.split('-').str[2]\n",
    "    file['replicate'] = file.index.str.split('-').str[4]\n",
    "    file = file.dropna()\n",
    "    baseline = file[file.timepoint.str.contains('0')]\n",
    "    baseline['treatment'] = 'INUL'\n",
    "    file = file.append(baseline)\n",
    "    baseline['treatment'] = 'PECT'\n",
    "    file = file.append(baseline)\n",
    "    baseline['treatment'] = 'FOS'\n",
    "    file = file.append(baseline)\n",
    "    file = file[(file.index.str.contains('CTRL'))|(file.index.str.contains('PECT'))\n",
    "                |(file.index.str.contains('INUL'))|(file.index.str.contains('FOS'))]\n",
    "    file = file[~file.timepoint.str.contains('2')]\n",
    "    file = file.sort_values(by=['sample','treatment','replicate','timepoint'])\n",
    "    file.set_index(['sample','treatment','replicate','timepoint'],inplace = True)\n",
    "    file = file.groupby(['sample','treatment','replicate']).diff().dropna().reset_index()\n",
    "    file = file.groupby(['sample','treatment']).mean().reset_index()\n",
    "    return file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb1c7d08-34b6-4cac-b46a-9163731ab885",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_list = ['H008-a.csv','H009-a.csv','H010-a.csv','H012-a.csv','H019-a.csv',\n",
    "               'H020-a.csv','H021-a.csv','H025-a.csv','H028-a.csv','H029-a.csv']\n",
    "os.chdir('/proj/gibbons/nbohmann/exvivo/gurry1/data/gc_data/')\n",
    "flux = pd.DataFrame([])\n",
    "for x in sample_list:\n",
    "    flux = flux.append(flux_calculate(x))\n",
    "flux.reset_index(inplace = True,drop = True)\n",
    "flux['index'] = flux['sample']+'_'+flux['treatment']\n",
    "flux.set_index('index',inplace = True)\n",
    "flux['predicted_but'] = flux.index.map(but)\n",
    "flux['predicted_ppa'] = flux.index.map(ppa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0c0ab53-d2b9-4548-8996-2e45156f99c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt1=(\n",
    "    ggplot(\n",
    "        flux,aes(x='ppa',y='predicted_ppa'))\n",
    "        +geom_point(aes(color='treatment'),size=5)\n",
    "        +geom_smooth(aes(group = 'treatment',color = 'treatment'),method='lm',linetype='--')\n",
    "#         +geom_text(aes(label = 'sample'),nudge_y = 1)\n",
    "        +labs(title='Propionate',x='Measured Propionate Production ($\\dfrac{mmol}{L*h}$)',\n",
    "              y = 'Predicted Propionate ($\\dfrac{mmol}{gDCW*h}$)')\n",
    "        +theme(text = element_text(size=15),panel_background=element_rect(fill = \"white\",\n",
    "                                colour = \"white\",size = 0.5, linetype = \"solid\"),\n",
    "                                panel_grid=element_line(size = .2, linetype = \"solid\",colour = \"gray\"),\n",
    "                                axis_line = element_line(size = 2, linetype = \"solid\",colour = \"black\"),\n",
    "                                legend_title=element_blank(),\n",
    "                                legend_position='right'))\n",
    "plt1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ab8c2c2-70dd-48a7-aaa6-f6ab386b5e85",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt1=(\n",
    "    ggplot(\n",
    "        flux,aes(x='but',y='predicted_but'))\n",
    "        +geom_point(aes(color='treatment'),size=5)\n",
    "        +geom_smooth(aes(group = 'treatment',color = 'treatment'),method='lm',linetype='--')\n",
    "#         +geom_text(aes(label = 'sample'),nudge_y = 1)\n",
    "        +ylim(0,80)\n",
    "        +labs(title='Butyrate',x='Measured Butyrate Flux ($\\dfrac{mmol}{L*h}$)',\n",
    "             y = 'Predicted Butyrate ($\\dfrac{mmol}{gDCW*h}$)')\n",
    "        +theme(text = element_text(size=15),panel_background=element_rect(fill = \"white\",\n",
    "                                colour = \"white\",size = 0.5, linetype = \"solid\"),\n",
    "                                panel_grid=element_line(size = .2, linetype = \"solid\",colour = \"gray\"),\n",
    "                                axis_line = element_line(size = 2, linetype = \"solid\",colour = \"black\"),\n",
    "                                legend_title=element_blank(),\n",
    "                                legend_position='right'))\n",
    "plt1"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}