31_plot_M3_growthcurves_warringer2003.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd           \n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as patches\n",
    "\n",
    "from scipy import integrate\n",
    "from scipy import stats\n",
    "import random\n",
    "\n",
    "from scipy.stats import spearmanr, pearsonr\n",
    "from latex_format import float2latex\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from selection_coefficient import Problem_M3, get_ODE_solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Update dependent parameters according to input\n",
    "import os\n",
    "import os.path\n",
    "from os import path\n",
    "\n",
    "## create export directory if necessary\n",
    "## foldernames for output plots/lists produced in this notebook\n",
    "import os\n",
    "FIG_DIR = f'./figures/M3_growthcurves/'\n",
    "os.makedirs(FIG_DIR, exist_ok=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "### execute script to load modules here\n",
    "exec(open('setup_aesthetics.py').read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATASET_COLOR = 'darkorange'\n",
    "palette = {'wild-type':'orange', 'knockout': 'dimgrey', 'wild-type median':'navy'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load wild-type traits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "INDEX_COL = [0,1,2,3,4]\n",
    "list_na_representations = ['not_present', 'failed_to_compute']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PCWS_TRAITS_WARRINGER = './output/df_M3_traits.csv'\n",
    "df_warringer = pd.read_csv(PCWS_TRAITS_WARRINGER, header = 0, index_col= INDEX_COL,\\\n",
    "                                  float_precision=None, na_values=list_na_representations)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### define default wild_type\n",
    "df_wildtypes = df_warringer[df_warringer['is_wildtype']==True]\n",
    "\n",
    "WILDTYPE = df_wildtypes.median(axis = 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load mutant data (averaged)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "PCWS_TRAITS_WARRINGER_AVERAGED = './output/df_M3_traits_averaged.csv'\n",
    "df_averaged = pd.read_csv(PCWS_TRAITS_WARRINGER_AVERAGED, header = 0, float_precision=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### assign wild-type label\n",
    "def is_wildtype(row):\n",
    "    genotype = row['genotype']\n",
    "    \n",
    "    if genotype == 'BY4741':\n",
    "        return True\n",
    "    else:\n",
    "        return False\n",
    "    \n",
    "\n",
    "row = df_averaged.iloc[0]\n",
    "is_wildtype(row)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_averaged['is_wildtype'] = df_averaged.apply(is_wildtype, axis =1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### append mutant values (averaged) to set of individual wild-type strains\n",
    "df_knockouts = df_averaged[~df_averaged['is_wildtype']]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set units of time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_knockouts['gmax'] = df_knockouts['gmax']*60 # change units to growth rate per hour\n",
    "df_knockouts['lag']  = df_knockouts['lag']/60 # change units to hour\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "WILDTYPE['gmax'] = WILDTYPE['gmax']*60 # change units to growth rate per hour\n",
    "WILDTYPE['lag']  = WILDTYPE['lag']/60 # change units to hour"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define initial conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### set initial resource concentrations\n",
    "\n",
    "CONCENTRATION_GLUCOSE = 20/180 * 1e3 # concentrations are recored  in milliMolar, to match the units of yield\n",
    "print(CONCENTRATION_GLUCOSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### define default initial_OD\n",
    "OD_START = 0.05  #df_warringer['od_start'].median()\n",
    "\n",
    "### compare to initial OD in the monoculture cycles\n",
    "fig, ax = plt.subplots(figsize = (FIGWIDTH_TRIPLET, FIGHEIGHT_TRIPLET))\n",
    "\n",
    "ax = df_warringer['od_start'].hist(bins=41, color = DATASET_COLOR, alpha = 0.6, log = True, rasterized = True)\n",
    "\n",
    "\n",
    "ax.axvline(OD_START, color = 'tab:red', label = f'median value: $N_0={OD_START:.3f}$')\n",
    "ax.legend()\n",
    "ax.set_xlabel('initial OD')\n",
    "ax.set_ylabel('no. growth curves')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot growth curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from selection_coefficient import Problem_M3, get_ODE_solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver_params = {'t_final': 20,\n",
    " 'timestep': 0.1,\n",
    " 'adaptive_timewindow': False,\n",
    " 'rtol': 1e-08,\n",
    " 'atol': 1e-12,\n",
    " 'scoeff_atol': 1e-08,\n",
    " 'scoeff_rtol': 1e-06}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot pairwise competition strain growth-curves using simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### choose a nice strain\n",
    "is_longer_lag = df_knockouts['lag'] > 1.5*WILDTYPE['lag']\n",
    "is_faster_growth = df_knockouts['gmax'] > 1.1*WILDTYPE['gmax']\n",
    "\n",
    "is_candidate = is_longer_lag & is_faster_growth\n",
    "list_candidates = np.arange(0,len(df_knockouts))[is_candidate]\n",
    "print(list_candidates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = 35\n",
    "l2, g2, Y2 = df_knockouts.loc[index, 'lag'], df_knockouts.loc[index, 'gmax'], df_knockouts.loc[index, 'yield']\n",
    "l1, g1, Y1 = WILDTYPE['lag'], WILDTYPE['gmax'], WILDTYPE['yield']\n",
    "\n",
    "strain_params = {'lam':[l1,l2], 'g':[g1,g2], 'Y':[Y1,Y2]}\n",
    "initial_conditions = {'R_0': CONCENTRATION_GLUCOSE, 'N_0': OD_START}\n",
    "problem = Problem_M3(**strain_params, **initial_conditions, x =0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = get_ODE_solution(problem, **solver_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (0.7*FIGWIDTH_TRIPLET,0.7*FIGHEIGHT_TRIPLET))\n",
    "\n",
    "ax.plot(sol.t,sol.y[0], color = palette['wild-type median'], lw = 3, zorder =54)\n",
    "ax.plot(sol.t,sol.y[1], color = palette['knockout'], lw = 3)\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylabel('absolute abundance [OD]')\n",
    "ax.set_xlabel('time [hours]')\n",
    "ax.set_xlim(0,sol.t[-1])\n",
    "ax.set_ylim(ymax = 20)\n",
    "\n",
    "fig.savefig(FIG_DIR + f\"pairwise_competition_{index}.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## store axis boundaries\n",
    "ymin,ymax = ax.get_ylim()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot single strain growth-curves using simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = Problem_M3(**strain_params, **initial_conditions, x =1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = get_ODE_solution(problem, **solver_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (0.7*FIGWIDTH_TRIPLET,0.7*FIGHEIGHT_TRIPLET))\n",
    "\n",
    "ax.plot(sol.t,sol.y[1], color = palette['knockout'], lw = 3)\n",
    "ax.set_ylim(ymin,ymax)\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylabel('absolute abundance [OD]')\n",
    "ax.set_xlabel('time [hours]')\n",
    "ax.set_xlim(0,sol.t[-1])\n",
    "\n",
    "fig.savefig(FIG_DIR + f\"knockout_averaged.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load raw growth curve data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exec(open('setup_plateau_finder_warringer2003.py').read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene = df_knockouts.loc[index, 'genotype']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## pick  curves to plot\n",
    "\n",
    "list_to_plot = []\n",
    "for name in df_trimmed.index:\n",
    "    if name[0]==gene:\n",
    "        list_to_plot.append(name)\n",
    "\n",
    "# list_to_plot = [0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot growth curve\n",
    "fig, ax = plt.subplots(figsize = (0.7*FIGWIDTH_TRIPLET,0.7* FIGHEIGHT_TRIPLET))\n",
    "\n",
    "# plot original growth curves\n",
    "for name in list_to_plot:\n",
    "\n",
    "    t_trimmed, od_trimmed = get_excess_od_timeseries_before_trim(name)\n",
    "    ax.scatter(t_trimmed/60, od_trimmed, color = 'tab:green', marker = 'o', s = 6)\n",
    "    \n",
    "# plot mean growth curve\n",
    "ax.plot(sol.t,sol.y[1], color = palette['knockout'], lw = 3, zorder = 3)\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlabel('time [hours]')\n",
    "ax.set_ylabel('absolute abundance [OD]')\n",
    "\n",
    "ax.set_xlim(0,sol.t[-1])\n",
    "ax.set_ylim(ymin,ymax)\n",
    "\n",
    "fig.savefig(FIG_DIR + f\"knockout_averaged_with_raw_data.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot wildtype growth-curves using simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = Problem_M3(**strain_params, **initial_conditions, x =0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = get_ODE_solution(problem, **solver_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (0.7*FIGWIDTH_TRIPLET,0.7*FIGHEIGHT_TRIPLET))\n",
    "\n",
    "ax.plot(sol.t,sol.y[0], color = palette['wild-type median'], lw = 3)\n",
    "ax.set_ylim(ymin,ymax)\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylabel('absolute abundance [OD]')\n",
    "ax.set_xlabel('time [hours]')\n",
    "ax.set_xlim(0,sol.t[-1])\n",
    "\n",
    "fig.savefig(FIG_DIR + f\"wildtype_averaged.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load raw growth curve data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exec(open('setup_plateau_finder_warringer2003.py').read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## pick wild-type curves to plot\n",
    "\n",
    "list_to_plot = []\n",
    "for name in df_wildtypes.index:\n",
    "    if name[1]==1:\n",
    "        list_to_plot.append(name)\n",
    "\n",
    "list_to_plot = df_wildtypes.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot growth curve\n",
    "fig, ax = plt.subplots(figsize = (0.7*FIGWIDTH_TRIPLET, 0.7*FIGHEIGHT_TRIPLET))\n",
    "\n",
    "# plot original growth curves\n",
    "for name in list_to_plot:\n",
    "\n",
    "    t_trimmed, od_trimmed = get_excess_od_timeseries_before_trim(name)\n",
    "    ax.plot(t_trimmed/60, od_trimmed, palette['wild-type'])\n",
    "    \n",
    "# plot mean growth curve\n",
    "ax.plot(sol.t,sol.y[0], color = palette['wild-type median'], lw = 3, zorder = 3)\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlabel('time [hours]')\n",
    "ax.set_ylabel('absolute abundance [OD]')\n",
    "\n",
    "ax.set_xlim(0,sol.t[-1])\n",
    "ax.set_ylim(ymin,ymax)\n",
    "\n",
    "fig.savefig(FIG_DIR + f\"wildtype_averaged_with_raw_data.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_knockouts.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}