35_plot_disranking_handcrafted_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluate relative fitness from trait distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters of this notebook "
   ]
  },
  {
   "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"
   ]
  },
  {
   "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/handcrafted_correlation/'\n",
    "os.makedirs(FIG_DIR, exist_ok=True)\n",
    "print(\"All  plots will be stored in: \\n\" + FIG_DIR)"
   ]
  },
  {
   "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'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SUFFIX_DATASET = 'pleiotropic_handcrafted'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define initial condition for bulk growth cycle"
   ]
  },
  {
   "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": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### define default initial_OD\n",
    "OD_FACTOR = 1\n",
    "OD_START = 0.01 *OD_FACTOR  #want OD_START to match the warringer data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "INITIAL_FREQ = 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create trait distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### handcrafted example scaled to larger number of points\n",
    "\n",
    "mean_wt = (1,2,1./CONCENTRATION_GLUCOSE) #gmax, lag ,yield\n",
    "NO_CURVES = 100\n",
    "def lag2gmax(lag):\n",
    "    return 0.1*(lag+mean_wt[1]) - 0.2 + mean_wt[0]\n",
    "\n",
    "\n",
    "effect_lag = np.linspace(-1.2, -0.8, num = NO_CURVES)\n",
    "mutant_lag = mean_wt[1] + effect_lag\n",
    "mutant_gmax = lag2gmax(mutant_lag)\n",
    "\n",
    "def gmax2yield(gmax):\n",
    "    ## define two constants\n",
    "    c = np.log(3)/0.01\n",
    "    A = np.exp(-c*1.1)*OD_FACTOR/CONCENTRATION_GLUCOSE # use same units as empirical data\n",
    "    return A*np.exp(c*(gmax))\n",
    "\n",
    "mutant_yield = gmax2yield(mutant_gmax)\n",
    "\n",
    "effect_sample = np.vstack([mutant_gmax, mutant_lag, mutant_yield]).T\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FIG_DIR_DATASET = FIG_DIR\n",
    "os.makedirs(FIG_DIR_DATASET, exist_ok=True)\n",
    "\n",
    "os.makedirs('./output/synthetic/', exist_ok=True)\n",
    "EXPORT_DATASET = './output/synthetic/' + SUFFIX_DATASET + '.csv'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### set up as dataframe\n",
    "\n",
    "## create dataframe with mutant effects\n",
    "df_effects = pd.DataFrame(data=effect_sample, columns = ['gmax', 'lag', 'yield'])\n",
    "df_effects['is_wildtype'] = False\n",
    "\n",
    "mean_wildtype = pd.Series(data =mean_wt, index =['gmax', 'lag', 'yield'] )\n",
    "mean_wildtype['is_wildtype'] = True\n",
    "\n",
    "df_effects=df_effects.append(mean_wildtype, ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def row2label(row):\n",
    "    if row['is_wildtype'] == True:\n",
    "        return 'wild-type'\n",
    "    else:\n",
    "        return 'knockout'\n",
    "\n",
    "df_effects['label']  = df_effects.apply(row2label,axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "df_traits = df_effects\n",
    "df_traits.to_csv(EXPORT_DATASET, index = False, float_format= '%.6e')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot trait distributionn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import spearmanr, pearsonr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from latex_format import float2latex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "n_datapoints = df_traits.shape[0]\n",
    "is_wildtype = df_traits['is_wildtype']==True\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "palette = {'wild-type':'orange', 'knockout': 'dimgrey'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### prepare data\n",
    "\n",
    "is_wildtype =df_traits['is_wildtype'] == True \n",
    "df_knockouts = df_traits.loc[~is_wildtype]\n",
    "data = df_knockouts.mask(is_wildtype, other = np.nan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xvar = 'gmax'\n",
    "yvar = 'lag'\n",
    "\n",
    "ratio = 5\n",
    "grid = sns.jointplot(data=data, x=xvar, y=yvar, \n",
    "                     hue = 'label', palette = palette, legend = False,\n",
    "                     marginal_kws = {'multiple':'layer', 'fill':False},\n",
    "                     marginal_ticks= False, space = 0, ratio = ratio,\n",
    "                     height = (1+1/ratio)*FIGHEIGHT_TRIPLET, )\n",
    "\n",
    "ax = grid.ax_joint\n",
    "ax = sns.scatterplot(data=data, x=xvar, y=yvar, color = 'dimgrey', edgecolor = 'none', ax = ax, legend=False)\n",
    "\n",
    "## plot wild-type\n",
    "x = mean_wildtype[xvar]\n",
    "y = mean_wildtype[yvar]\n",
    "ax.scatter(x=x,y=y, color = palette['wild-type'], label = 'wild-type')\n",
    "#grid.ax_marg_y.axhline(y, color = palette['wild-type'])\n",
    "#grid.ax_marg_x.axvline(x, color = palette['wild-type'])\n",
    "\n",
    "\n",
    "### replot the marginal distributions in different colors\n",
    "ax = grid.ax_marg_x\n",
    "palette['knockout'] = 'dimgrey' ### fix the color for marginals in growth rate\n",
    "sns.kdeplot(data = df_traits, x=xvar, ax=ax, fill = True, hue = 'label', palette=palette, multiple = 'layer',\n",
    "       legend = False, cut = 0)\n",
    "\n",
    "### replot the marginal distributions in different colors\n",
    "ax = grid.ax_marg_y\n",
    "palette['knockout'] = 'dimgrey' ### fix the color for marginals in growth rate\n",
    "sns.kdeplot(data = df_traits, y=yvar, ax=ax, fill = True, hue = 'label', palette=palette, multiple = 'layer',\n",
    "       legend = False, cut = 0)\n",
    "\n",
    "## set label\n",
    "ax = grid.ax_joint\n",
    "ax.set_ylabel('lag time [hours]')\n",
    "ax.set_xlabel('growth rate [per hour]')\n",
    "## set legend\n",
    "#ax.legend(loc = 'lefleft', bbox_to_anchor = (-0.65,0.99), frameon=False,\n",
    "#          title = 'lag time [hours]', title_fontsize = MEDIUM_SIZE)\n",
    "#ax.legend(frameon = False)\n",
    "## set title\n",
    "\n",
    "title = f\"n={df_knockouts.shape[0]} mutants\"\n",
    "ax = grid.ax_marg_x\n",
    "#ax.annotate(title, (0.7,0.05), xycoords ='axes fraction') # right\n",
    "ax.annotate(title, (0.02,0.05), xycoords ='axes fraction') # left\n",
    "grid.fig.savefig(FIG_DIR_DATASET + f\"correlation_{xvar}-vs-{yvar}.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xvar = 'gmax'\n",
    "yvar = 'yield'\n",
    "\n",
    "ratio = 5\n",
    "grid = sns.jointplot(data=data, x=xvar, y=yvar, \n",
    "                     hue = 'label', palette = palette, legend = False,\n",
    "                     marginal_kws = {'multiple':'layer', 'fill':False},\n",
    "                     marginal_ticks= False, space = 0, ratio = ratio,\n",
    "                    height = (1+1/ratio)*FIGHEIGHT_TRIPLET, )\n",
    "\n",
    "\n",
    "## replot with coloring by lag time\n",
    "ax = grid.ax_joint\n",
    "ax = sns.scatterplot(data=data, x=xvar, y=yvar, color = 'dimgrey', edgecolor = 'none', ax = ax, legend=False)\n",
    "\n",
    "## plot wild-type\n",
    "x = mean_wildtype[xvar]\n",
    "y = mean_wildtype[yvar]\n",
    "ax.scatter(x=x,y=y, color = palette['wild-type'], label = 'wild-type')\n",
    "#grid.ax_marg_y.axhline(y, color = palette['wild-type'])\n",
    "#grid.ax_marg_x.axvline(x, color = palette['wild-type'])\n",
    "\n",
    "\n",
    "### replot the marginal distributions in different colors\n",
    "ax = grid.ax_marg_x\n",
    "palette['knockout'] = 'dimgrey' ### fix the color for marginals in growth rate\n",
    "sns.kdeplot(data = df_traits, x=xvar, ax=ax, fill = True, hue = 'label', palette=palette, multiple = 'layer',\n",
    "       legend = False, cut = 0)\n",
    "\n",
    "### replot the marginal distributions in different colors\n",
    "ax = grid.ax_marg_y\n",
    "palette['knockout'] = 'dimgrey' ### fix the color for marginals in growth rate\n",
    "sns.kdeplot(data = df_traits, y=yvar, ax=ax, fill = True, hue = 'label', palette=palette, multiple = 'layer',\n",
    "       legend = False, cut= 0)\n",
    "\n",
    "## fix y-axis limit, yield can only be positive\n",
    "ax.set_ylim(-0.0004)\n",
    "\n",
    "## set label\n",
    "ax = grid.ax_joint\n",
    "ax.set_ylabel('biomass yield [OD/mM glucose]')\n",
    "ax.set_xlabel('growth rate [per hour]')\n",
    "## set legend\n",
    "#ax.legend(loc = 'lefleft', bbox_to_anchor = (-0.65,0.99), frameon=False,\n",
    "#          title = 'lag time [hours]', title_fontsize = MEDIUM_SIZE)\n",
    "#ax.legend(frameon = False)\n",
    "## set title\n",
    "\n",
    "title = f\"n={df_knockouts.shape[0]} mutants\"\n",
    "ax = grid.ax_marg_x\n",
    "#ax.annotate(title, (0.7,0.05), xycoords ='axes fraction') # right\n",
    "ax.annotate(title, (0.02,0.05), xycoords ='axes fraction') # left\n",
    "grid.fig.savefig(FIG_DIR_DATASET + f\"correlation_{xvar}-vs-{yvar}.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load trait data into the standard form required by Michaels code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert notation, so we can reuse old code\n",
    "df_input = df_knockouts\n",
    "n_input = df_input.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create wildtype\n",
    "WILDTYPE = df_traits[df_traits['label']=='wild-type']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### growth rates\n",
    "gs = np.zeros(n_input+1)\n",
    "gs[0] = WILDTYPE['gmax']\n",
    "gs[1:] = df_input['gmax'].values\n",
    "\n",
    "### lag times\n",
    "ls = np.zeros(n_input+1)\n",
    "ls[0] = WILDTYPE['lag']\n",
    "ls[1:] = df_input['lag'].values\n",
    "\n",
    "### yield\n",
    "Ys = np.zeros(n_input+1)\n",
    "Ys[0] = WILDTYPE['yield']\n",
    "Ys[1:] = df_input['yield'].values\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate effective yield"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bulk_simulation_code import CalcRelativeYield"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### calculcate effective yields\n",
    "nus = CalcRelativeYield(Ys, R0 = CONCENTRATION_GLUCOSE, N0 = OD_START)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate pairwise competition growth cycles (scenario A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bulk_simulation_code import toPerGeneration, run_pairwise_experiment\n",
    "from bulk_simulation_code import CalcTotalSelectionCoefficientLogit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "xs_pair, xs_pair_final, tsats,fcs_both, fcs_wt, fcs_mut = run_pairwise_experiment(\n",
    "                                                                gs=gs,   ls=ls,   nus = nus, \n",
    "                                                                g1=gs[0],l1=ls[0],nu1=nus[0],\n",
    "                                                                x0 = INITIAL_FREQ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_percycle = CalcTotalSelectionCoefficientLogit(xs_pair,xs_pair_final)\n",
    "s_pergen = np.divide(s_percycle, np.log(fcs_wt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### store results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_output = df_input.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_output['logfc_wt'] = np.log(fcs_wt[1:])\n",
    "df_output['logfc_mut'] = np.log(fcs_mut[1:])\n",
    "\n",
    "df_output['logit_percycle'] = s_percycle[1:]\n",
    "df_output['logit_pergen'] = s_pergen[1:]\n",
    "\n",
    "df_output['logit_percycle_rank'] =df_output['logit_percycle'].rank()\n",
    "df_output['logit_pergen_rank'] = df_output['logit_pergen'].rank()\n",
    "\n",
    "df_output['deltarank'] = df_output['logit_pergen_rank'] - df_output['logit_percycle_rank'] \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare data for plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### sort by label prepare for plotting\n",
    "\n",
    "def row2label(row):\n",
    "    if row['is_wildtype'] == True:\n",
    "        return 'wild-type'\n",
    "    else:\n",
    "        return 'knockout'\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_output['label']  = df_output.apply(row2label,axis=1)\n",
    "\n",
    "df_output = df_output.sort_values('label')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## sort by misranking\n",
    "\n",
    "df_output['deltarank_abs'] = np.abs(df_output['deltarank'])\n",
    "df_sorted = df_output.sort_values('deltarank', ascending = True)\n",
    "select = df_sorted.index[[-1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_sorted.loc[select]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot misranking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot residuals\n",
    "\n",
    "palette = {'wild-type':'orange', 'knockout': 'dimgrey'}\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (FIGHEIGHT_TRIPLET, FIGHEIGHT_TRIPLET))\n",
    "\n",
    "x_var = 'logit_percycle'\n",
    "y_var = 'deltarank'\n",
    "data = df_output\n",
    "sns.scatterplot(data = data, x = x_var, y = y_var, rasterized = False, ax = ax, \n",
    "                color = 'dimgrey', edgecolor = 'none', legend =False)\n",
    "\n",
    "\n",
    "### plot select points\n",
    "is_labeled = True\n",
    "for i in select:\n",
    "    A, B = float(data.loc[i, x_var]), float(data.loc[i, y_var])\n",
    "    #ax.scatter(A-0.15,B,s=150,color ='tab:red', zorder = -1, marker = 5 )\n",
    "    if is_labeled == False: \n",
    "        label = 'max. disagreement' \n",
    "        is_labeled = True\n",
    "    else: label = None\n",
    "    #ax.scatter(A,B,s=200,color ='tab:blue', zorder = -1,label = label, alpha = 0.25)\n",
    "\n",
    "\n",
    "### plot horizontal line for orientation\n",
    "ax.axhline(0,ls = '--', color = 'black')\n",
    "\n",
    "### annotate\n",
    "ax.set_xlabel('relative fitness per-cycle:' + r'  $s^{\\mathrm{logit}}_{\\mathrm{cycle}}$')\n",
    "#ax.set_ylabel('rank difference to\\nrelative fitness per-generation [rank]')\n",
    "ax.set_ylabel('rank difference between fitness\\nper-generation and fitness per-cycle')\n",
    "ax.legend(loc = 'upper left', bbox_to_anchor = (-0.05,1.0), frameon=False) #inside\n",
    "\n",
    "title = f\"n = {sum(data['is_wildtype']==False)} mutants\"\n",
    "ax.set_title(title, loc = 'left')\n",
    "\n",
    "fig.savefig(FIG_DIR + f\"residuals_{x_var}_vs_{y_var}_x0={INITIAL_FREQ:.2f}.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot on foldchange phase diagram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eval_isocline_percycle(logfc_wt, level):\n",
    "    return logfc_wt + level\n",
    "\n",
    "def eval_isocline_pergen(logfc_wt, level):\n",
    "    return np.multiply((level+1),logfc_wt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "palette = {'wild-type':'orange', 'knockout': 'dimgrey'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot cloud of points\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (FIGHEIGHT_TRIPLET, FIGHEIGHT_TRIPLET))\n",
    "\n",
    "x_var = 'logfc_wt'\n",
    "y_var = 'logfc_mut'\n",
    "data = df_output\n",
    "sns.scatterplot(data = data, x = x_var, y = y_var, rasterized = False, ax = ax, \n",
    "                color = 'dimgrey', zorder = 2, edgecolor = 'none', legend = False)\n",
    "\n",
    "\n",
    "\n",
    "## find value limits\n",
    "\n",
    "fcmax = np.max([data[x_var].max(),data[y_var].max()])\n",
    "fcmin = np.min([data[x_var].min(),data[y_var].min()])\n",
    "assert fcmin > 0\n",
    "\n",
    "### set axis limits\n",
    "#xmin = 1\n",
    "#assert xmin < fcmin, 'We are cutting points from the dataplot!'\n",
    "#xmax = 5.5\n",
    "#assert xmax > fcmax, 'We are cutting points from the dataplot!'\n",
    "#ax.set_xlim(xmin,xmax)\n",
    "#ax.set_ylim(xmin,xmax)\n",
    "\n",
    "\n",
    "## take off axis spines\n",
    "#sns.despine(left=False, bottom = False, ax = ax)\n",
    "\n",
    "### find axis limits\n",
    "xmin, xmax = ax.get_xlim()\n",
    "fcwt_vec = np.linspace(xmin,xmax, num = 100) \n",
    "fcwt_vec = np.concatenate((-fcwt_vec,fcwt_vec))\n",
    "color_percycle = 'tab:grey'\n",
    "color_pergen = 'navy'\n",
    "\n",
    "### plot per cycle isoclines\n",
    "levels = np.outer([-1,1],np.linspace(0.01,8,num = 6)).flatten()\n",
    "levels.sort()\n",
    "\n",
    "for level in levels: \n",
    "\n",
    "    y = eval_isocline_percycle(fcwt_vec, level = level)\n",
    "    #ax.plot(fcwt_vec, y, color = color_percycle)\n",
    "    \n",
    "    \n",
    "## plot per generationa isoclines\n",
    "angles = np.linspace(0,np.pi/2 - 0.001, num = 6)\n",
    "levels = np.outer([-1,1],np.tan(angles)).flatten()\n",
    "\n",
    "for level in levels: \n",
    "    y = eval_isocline_pergen(fcwt_vec, level =level) \n",
    "    #ax.plot(fcwt_vec, y, color = color_pergen)\n",
    "    \n",
    "## plot diagonal \n",
    "xmin, xmax = ax.get_xlim()\n",
    "ymin, ymax = ax.get_ylim()\n",
    "ax.plot([-xmin,xmax],[-xmin,xmax], color = 'black', ls = '--')\n",
    "ax.set_xlim(xmin,xmax)\n",
    "ax.set_ylim(ymin,ymax)\n",
    "## add  legend items\n",
    "#ax.plot([],[], color = color_percycle, label = 'per-cycle $s$ isocline')\n",
    "#ax.plot([],[], color = color_pergen, label = 'per-generation $Q$ isocline')\n",
    "\n",
    "## plot red cone for a select point\n",
    "select = [0]\n",
    "for i in select:\n",
    "    A, B = float(data.loc[i, x_var]), float(data.loc[i, y_var])\n",
    "    #ax.scatter(A,B,s=70,color ='tab:red', zorder = 3)\n",
    "    #ax.scatter(A-0.05,B,s=150,color ='tab:red', zorder = -1, marker = 5 )\n",
    "    #ax.scatter(A,B,s=70,color ='tab:red', zorder = 3, marker = 'v')\n",
    "    #ax.scatter(A,B,s=200,color ='tab:blue', zorder = -1,label = label, alpha = 0.25)\n",
    "\n",
    "\n",
    "    x_fill = np.linspace(fcwt_vec[0],fcwt_vec[-1])\n",
    "    y_fill = B/A*x_fill\n",
    "\n",
    "    ax.fill_between(x_fill, (x_fill - A) + B, y_fill, color='tab:red', alpha=0.25)\n",
    "\n",
    "\n",
    "### annotate\n",
    "#ax.legend(loc = 'upper left', bbox_to_anchor = (1.3,1)) # outside\n",
    "#ax.legend(loc = 'upper left', bbox_to_anchor = (-0.05,0.25), frameon=False) #inside\n",
    "title = f\"n = {sum(data['is_wildtype']==False)} mutants\"\n",
    "#ax.set_title(title, loc = 'left')\n",
    "\n",
    "ax.set_xlabel(r\"wild-type log fold-change: $\\mathrm{LFC}_{\\mathrm{wt}}$\")\n",
    "ax.set_ylabel(r\"mutant log fold-change: $\\mathrm{LFC}_{\\mathrm{mut}}$\")\n",
    "\n",
    "title = f\"n = {sum(data['is_wildtype']==False)} mutants\"\n",
    "ax.set_title(title, loc = 'left')\n",
    "\n",
    "\n",
    "fig.savefig(FIG_DIR+ f'scatterplot_logfc_wt_vs_logfc_mut_x0={INITIAL_FREQ:.2f}.pdf', DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n",
    "              "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot DFEs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "palette = {'wild-type':'orange', 'knockout': 'dimgrey', 'wild-type median':'navy'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (FIGHEIGHT_TRIPLET, FIGHEIGHT_TRIPLET))\n",
    "\n",
    "data = df_output\n",
    "x_var = 'logit_percycle'\n",
    "#sns.histplot(data, x = 'logit_percycle', rasterized = False, ax = ax,\n",
    "#                hue = 'label', palette = palette)\n",
    "\n",
    "## need a context wrapper, else pandas throws a Future Warning\n",
    "## see https://stackoverflow.com/questions/15777951/how-to-suppress-pandas-future-warning\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "    # Warning-causing lines of code here\n",
    "    sns.kdeplot(data=data, x=x_var, hue=\"label\",common_norm = True,\n",
    "            palette = palette, multiple=\"layer\", ax = ax, fill = True, legend = False)\n",
    "    \n",
    "### plot selection coefficient zero for orientation\n",
    "#ax.axvline(0, ls = '--', color = 'black')\n",
    "#ax.legend_.set_title('')\n",
    "#ax.legend([],[])\n",
    "\n",
    "title = f\"n = {sum(data['is_wildtype']==False)} mutants\"\n",
    "ax.set_title(title, loc = 'left')\n",
    "\n",
    "ax.set_xlabel('relative fitness per-cycle:' + r'  $s^{\\mathrm{logit}}_{\\mathrm{cycle}}$')\n",
    "ax.set_ylabel('mutant density')\n",
    "\n",
    "\n",
    "fig.savefig(FIG_DIR+ f'dfeplot_logit_percycle_x0={INITIAL_FREQ:.2f}.pdf', DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (FIGHEIGHT_TRIPLET, FIGHEIGHT_TRIPLET))\n",
    "\n",
    "data = df_output\n",
    "x_var = 'logit_pergen'\n",
    "#sns.histplot(data, x = 'logit_percycle', rasterized = False, ax = ax,\n",
    "#                hue = 'label', palette = palette)\n",
    "\n",
    "## need a context wrapper, else pandas throws a Future Warning\n",
    "## see https://stackoverflow.com/questions/15777951/how-to-suppress-pandas-future-warning\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "    # Warning-causing lines of code here\n",
    "    sns.kdeplot(data=data, x=x_var, hue=\"label\",common_norm = True,\n",
    "            palette = palette, multiple=\"layer\", ax = ax, fill = True, legend = False)\n",
    "    \n",
    "### plot selection coefficient zero for orientation\n",
    "#ax.axvline(0, ls = '--', color = 'black')\n",
    "#ax.legend_.set_title('')\n",
    "#ax.legend([],[])\n",
    "\n",
    "title = f\"n = {sum(data['is_wildtype']==False)} mutants\"\n",
    "ax.set_title(title, loc = 'left')\n",
    "\n",
    "ax.set_xlabel('relative fitness per-generation:' + r'  $s^{\\mathrm{logit}}_{\\mathrm{gen}}$')\n",
    "ax.set_ylabel('mutant density')\n",
    "\n",
    "\n",
    "fig.savefig(FIG_DIR+ f'dfeplot_logit_pergen_x0={INITIAL_FREQ:.2f}.pdf', DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot global correlation in ranks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot correlation\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (FIGHEIGHT_TRIPLET, FIGHEIGHT_TRIPLET))\n",
    "\n",
    "x_var = 'logit_percycle_rank'\n",
    "y_var = 'logit_pergen_rank'\n",
    "data = df_output\n",
    "sns.scatterplot(data = data, x = x_var, y = y_var, rasterized = False, ax = ax, \n",
    "                color = 'dimgrey', edgecolor = 'none', legend = False)\n",
    "\n",
    "### plot select points\n",
    "is_labeled = True\n",
    "for i in select:\n",
    "    A, B = float(data.loc[i, x_var]), float(data.loc[i, y_var])\n",
    "    #ax.scatter(A-0.15,B,s=150,color ='tab:red', zorder = -1, marker = 5 )\n",
    "    if is_labeled == False: \n",
    "        label = 'max. disagreement' \n",
    "        is_labeled = True\n",
    "    else: label = None\n",
    "    ax.scatter(A,B,s=200,color ='tab:blue', zorder = -1,label = label, alpha = 0.25)\n",
    "\n",
    "\n",
    "### annotate\n",
    "ax.legend(loc = 'upper left', bbox_to_anchor = (-0.05,1.0), frameon=False) #inside\n",
    "\n",
    "ax.set_ylabel('relative fitness per-generation [rank]')\n",
    "ax.set_xlabel('relative fitness per-cycle [rank]')\n",
    "\n",
    "title = f\"n = {sum(data['is_wildtype']==False)} mutants\"\n",
    "ax.set_title(title, loc = 'left')\n",
    "\n",
    "fig.savefig(FIG_DIR + f\"scatterplot_{x_var}_vs_{y_var}_x0={INITIAL_FREQ:.2f}.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot global correlation in values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot correlation\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (FIGHEIGHT_TRIPLET, FIGHEIGHT_TRIPLET))\n",
    "\n",
    "x_var = 'logit_percycle'\n",
    "y_var = 'logit_pergen'\n",
    "data = df_output\n",
    "sns.scatterplot(data = data, x = x_var, y = y_var, rasterized = False, ax = ax,\n",
    "                color = 'dimgrey', edgecolor = 'none',  legend = False)\n",
    "\n",
    "### plot select points\n",
    "is_labeled = True\n",
    "for i in select:\n",
    "    A, B = float(data.loc[i, x_var]), float(data.loc[i, y_var])\n",
    "    #ax.scatter(A-0.15,B,s=150,color ='tab:red', zorder = -1, marker = 5 )\n",
    "    if is_labeled == False: \n",
    "        label = 'max. disagreement' \n",
    "        is_labeled = True\n",
    "    else: label = None\n",
    "    #ax.scatter(A,B,s=200,color ='tab:blue', zorder = -1,label = label, alpha = 0.25)\n",
    "\n",
    "\n",
    "### annotate\n",
    "ax.legend(loc = 'upper left', bbox_to_anchor = (-0.05,1.0), frameon=False) #inside\n",
    "\n",
    "ax.set_xlabel('relative fitness per-cycle:' + r'  $s^{\\mathrm{logit}}_{\\mathrm{cycle}}$')\n",
    "ax.set_ylabel('relative fitness per-generation:' + r'  $s^{\\mathrm{logit}}_{\\mathrm{gen}}$')\n",
    "\n",
    "title = f\"n = {sum(data['is_wildtype']==False)} mutants\"\n",
    "ax.set_title(title, loc = 'left')\n",
    "\n",
    "fig.savefig(FIG_DIR + f\"scatterplot_{x_var}_vs_{y_var}_x0={INITIAL_FREQ:.2f}.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check correlation with saturation time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot correlation\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (FIGHEIGHT_TRIPLET, FIGHEIGHT_TRIPLET))\n",
    "\n",
    "x_var = 'logfc_wt'\n",
    "y_var = 'logit_percycle'\n",
    "data = df_output\n",
    "sns.scatterplot(data = data, x = x_var, y = y_var, rasterized = False, ax = ax,\n",
    "                color = 'dimgrey', edgecolor = 'none', legend = False)\n",
    "\n",
    "### plot select points\n",
    "is_labeled = True\n",
    "for i in select:\n",
    "    A, B = float(data.loc[i, x_var]), float(data.loc[i, y_var])\n",
    "    #ax.scatter(A-0.15,B,s=150,color ='tab:red', zorder = -1, marker = 5 )\n",
    "    if is_labeled == False: \n",
    "        label = 'max. disagreement' \n",
    "        is_labeled = True\n",
    "    else: label = None\n",
    "    #ax.scatter(A,B,s=200,color ='tab:blue', zorder = -1,label = label, alpha = 0.25)\n",
    "\n",
    "\n",
    "### annotate\n",
    "#ax.legend(loc = 'upper left', bbox_to_anchor = (-0.05,1.0), frameon=False) #inside\n",
    "\n",
    "ax.set_ylabel('relative fitness per-cycle:' + r'  $s^{\\mathrm{logit}}_{\\mathrm{cycle}}$')\n",
    "ax.set_xlabel('wild-type fold-change:' + r'  $\\mathrm{LFC}_{\\mathrm{wt}}$')\n",
    "\n",
    "title = f\"n = {sum(data['is_wildtype']==False)} mutants\"\n",
    "ax.set_title(title, loc = 'left')\n",
    "\n",
    "fig.savefig(FIG_DIR + f\"scatterplot_{x_var}_vs_{y_var}_x0={INITIAL_FREQ:.2f}.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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
}