62_plot_LTEE_fitness_trajectories.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Long-term fitness trajectories from the LTEE\n",
    "\n",
    "Data available in the Supplemental Material of Good et al. Nature 2017. \n",
    "Download possible from Ben Good's github repository [here](https://github.com/benjaminhgood/LTEE-metagenomic/blob/master/additional_data/Concatenated.LTEE.data.all.csv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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/LTEE_trajectories/'\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": [
    "#df = pd.read_csv('./output/LTEE_all_data.csv')\n",
    "df = pd.read_csv('./output/LTEE_averaged_data.csv')"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "### explanation of column headers and labels\n",
    "\n",
    "    607 is the wild-type strain 'REL607'\n",
    "    D.0 is the dilution factor applied to count colonies in initial inoculum. \n",
    "    D.1 is the dilution factor applied to count colonies in the saturated population. We expect D.1 = 100*D.0. \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_logit_percycle =  'tab:grey' \n",
    "color_logit_pergen = 'firebrick'\n",
    "color_log_percycle = 'navy'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## plot for single lineage\n",
    "\n",
    "fig, axes = plt.subplots(3,1, sharex=True, figsize = (2*FIGHEIGHT_TRIPLET,2*FIGHEIGHT_TRIPLET))\n",
    "\n",
    "ax = axes[0]\n",
    "ax2 =  axes[1]\n",
    "ax3 = axes[2]\n",
    "\n",
    "for pop in list(set(df['Population'].values)):\n",
    "    \n",
    "    df_subset = df[df['Population'] == pop]\n",
    "    ## sort by timepoint\n",
    "    df_subset = df_subset.sort_values('Generation', ascending = True)\n",
    "    \n",
    "\n",
    "    x = df_subset['Generation']\n",
    "    y = df_subset['logit_percycle']\n",
    "    ax.plot(x,y, color = color_logit_percycle, label = pop)\n",
    "\n",
    "\n",
    "    x = df_subset['Generation']\n",
    "    y = df_subset['logit_pergen']\n",
    "    ax2.plot(x,y, color = color_logit_pergen)\n",
    "\n",
    "\n",
    "\n",
    "    \n",
    "    x = df_subset['Generation']\n",
    "    y = df_subset['log_percycle']\n",
    "    ax3.plot(x,y, color = color_log_percycle)\n",
    "\n",
    "## set labels\n",
    "\n",
    "\n",
    "ax.set_ylabel('evolved fitness' + r'  $s^{\\mathrm{logit}}_{\\mathrm{cycle}}$')\n",
    "ax2.set_ylabel('evolved fitness' + r'  $s^{\\mathrm{logit}}_{\\mathrm{gen}}$')\n",
    "ax3.set_ylabel('evolved fitness' + r'  $s^{\\mathrm{log}}_{\\mathrm{cycle}}$')\n",
    "ax3.set_xlabel('time [#generations in evolution experiment]')\n",
    "\n",
    "for ax in axes:\n",
    "    ax.set_xlim(0,df_subset['Generation'].max())\n",
    "    ax.set_ylim(0)\n",
    "\n",
    "fig.savefig(FIG_DIR + \"timeseries_fitness_alternative_statistics.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list_generations = list(set(df['Generation'].values))\n",
    "list_generations.sort()\n",
    "\n",
    "df_stats = pd.DataFrame(index = list_generations)\n",
    "#df_stats['rho_Spearman'] = -1 # dummy default value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## calculate correlation coefficient at all timepoints as a vector\n",
    "for t in df_stats.index:\n",
    "    \n",
    "    df_to_calc = df[df['Generation'] == t]\n",
    "\n",
    "    ## correlation between logit_percycle and logit_pergen\n",
    "    x = df_to_calc['logit_percycle'].values\n",
    "    y = df_to_calc['logit_pergen'].values\n",
    "    rho, pval = stats.spearmanr(x,y)\n",
    "    df_stats.at[t,'rho_percycle-pergen'] = rho\n",
    "    r, pval = stats.pearsonr(x,y)\n",
    "    df_stats.at[t,'r_percycle-pergen'] = r\n",
    "    \n",
    "    ### correlation between s and delta_log\n",
    "    x = df_to_calc['logit_percycle'].values\n",
    "    z = df_to_calc['log_percycle'].values\n",
    "    rho, pval = stats.spearmanr(x,z)\n",
    "    df_stats.at[t,'rho_logit-log'] = rho\n",
    "    r, pval = stats.pearsonr(x,z)\n",
    "    df_stats.at[t,'r_logit-log'] = r\n",
    "    \n",
    "    ### maximum disranking\n",
    "    x = df_to_calc['logit_percycle'].rank()\n",
    "    y = df_to_calc['logit_pergen'].rank()\n",
    "    deltaranks = y-x\n",
    "    df_stats.at[t,'deltarank_max'] = deltaranks.max()\n",
    "    df_stats.at[t,'deltarank_abs_max']  = np.abs(deltaranks).max()\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "\n",
    "## compare different measures of  correlation\n",
    "\n",
    "fig, axes = plt.subplots(2,1, sharex=True, sharey=True,\n",
    "                         figsize = (2*FIGHEIGHT_TRIPLET,4/3*FIGHEIGHT_TRIPLET))\n",
    "\n",
    "ax = axes[0]\n",
    "x = df_stats.index\n",
    "y = df_stats['rho_percycle-pergen']\n",
    "ax.plot(x,y, color = color_logit_pergen)\n",
    "ax.set_ylabel(r'$s^{\\mathrm{logit}}_{\\mathrm{gen}}$ and $s^{\\mathrm{logit}}_{\\mathrm{cycle}}$:'+'\\n'\n",
    "             + r\"rank correlation $\\rho$\")\n",
    "\n",
    "ax = axes[1]\n",
    "x = df_stats.index\n",
    "y = df_stats['rho_logit-log']\n",
    "ax.plot(x,y, color = color_log_percycle)\n",
    "ax.set_ylabel(r'$s^{\\mathrm{log}}_{\\mathrm{cycle}}$ and $s^{\\mathrm{logit}}_{\\mathrm{cycle}}$:'+'\\n'\n",
    "             + r\"rank correlation $\\rho$\")\n",
    "\n",
    "ax.set_xlabel('time [#generations in evolution experiment]')\n",
    "ax.set_xlim(0,df_subset['Generation'].max())\n",
    "\n",
    "fig.savefig(FIG_DIR + \"timeseries_correlation_both.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### plot individual correlationat a fixed timepoint\n",
    "\n",
    "fig,ax = plt.subplots(figsize = (0.3*FIGHEIGHT_TRIPLET,0.3*FIGHEIGHT_TRIPLET))\n",
    "t = 15000\n",
    "\n",
    "\n",
    "df_to_calc = df[df['Generation'] == t]\n",
    "x = df_to_calc['logit_percycle']\n",
    "y = df_to_calc['log_percycle']\n",
    "\n",
    "ax.scatter(x,y, label = f'Generation\\nt={t:d}', color = color_log_percycle)\n",
    "\n",
    "## show diagonal\n",
    "xmin,xmax = ax.get_xlim()\n",
    "ymin,ymax = ax.get_ylim()\n",
    "#ax.plot([xmin,xmax],[xmin,xmax], label = 'y=x', color = 'black', ls = '--')\n",
    "ax.set_xlim(xmin,xmax)\n",
    "ax.set_ylim(ymin,ymax)\n",
    "\n",
    "#ax.legend( loc= 'upper left')\n",
    "\n",
    "rho = df_stats.at[t,'rho_logit-log'] \n",
    "r = df_stats.at[t,'r_logit-log']\n",
    "\n",
    "#title = fr'Pearson $r={r:.2f}$, Spearman $\\rho={rho:.2f}$'\n",
    "title = fr'$\\rho={rho:.2f}$' # short title\n",
    "ax.set_title(title, loc = 'center')\n",
    "ax.tick_params(labelleft=False, labelbottom = False, left =False, bottom=False)\n",
    "\n",
    "\n",
    "ax.set_xlabel(r'$s^{\\mathrm{logit}}_{\\mathrm{cycle}}$')\n",
    "ax.set_ylabel( r'$s^{\\mathrm{log}}_{\\mathrm{cycle}}$')\n",
    "\n",
    "\n",
    "fig.savefig(FIG_DIR + f\"snapshot_logit-vs-log_at_t={t}.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Focus on comparison per-time point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list_timepoints = [4000, 15000, 30000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "## compare different measures of  correlation\n",
    "\n",
    "fig, ax = plt.subplots(1,1, figsize = (2*FIGHEIGHT_TRIPLET,2/3*FIGHEIGHT_TRIPLET))\n",
    "\n",
    "x = df_stats.index\n",
    "y = df_stats['rho_percycle-pergen']\n",
    "ax.plot(x,y, color = 'silver', lw = 3)\n",
    "\n",
    "for t in list_timepoints:\n",
    "    ax.axvline(t, ls = 'dotted', color = 'black')\n",
    "ax.set_ylabel(r'$s^{\\mathrm{logit}}_{\\mathrm{gen}}$ and $s^{\\mathrm{logit}}_{\\mathrm{cycle}}$:'+'\\n'\n",
    "             + r\"rank correlation $\\rho$\")\n",
    "\n",
    "ax.set_xlabel('time [# generations in evolution experiment]')\n",
    "ax.set_xlim(0,x.max())\n",
    "\n",
    "fig.savefig(FIG_DIR + \"timeseries_correlation.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "## compare different measures of  correlation\n",
    "\n",
    "fig, ax = plt.subplots(1,1, figsize = (2*FIGHEIGHT_TRIPLET,2/3*FIGHEIGHT_TRIPLET))\n",
    "\n",
    "x = df_stats.index\n",
    "y = df_stats['deltarank_abs_max']\n",
    "ax.plot(x,y, color = 'silver', lw = 3)\n",
    "\n",
    "for t in list_timepoints:\n",
    "    ax.axvline(t, ls = 'dotted', color = 'black')\n",
    "ax.set_ylabel('max. rank difference\\n' + 'betw. $s^{\\mathrm{logit}}_{\\mathrm{gen}}$ and $s^{\\mathrm{logit}}_{\\mathrm{cycle}}$')\n",
    "\n",
    "ax.set_xlabel('time [# generations in evolution experiment]')\n",
    "ax.set_xlim(0,x.max())\n",
    "ax.set_ylim(0)\n",
    "\n",
    "fig.savefig(FIG_DIR + \"timeseries_maximum_mismatch.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "    \n",
    "### redefine colormap\n",
    "cmap = sns.color_palette('crest', as_cmap=True)\n",
    "norm = plt.Normalize(df['Generation'].min(), df['Generation'].max())\n",
    "sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot individual correlationat a fixed timepoint\n",
    "\n",
    "\n",
    "for t in list_timepoints[:]:\n",
    "    df_to_calc = df[df['Generation'] == t]  # reduce data to current timepoint\n",
    "    color_t = cmap(norm(t)) ## use colormap\n",
    "    \n",
    "    ## sort by misranking\n",
    "    df_sorted = df_to_calc.copy(deep=True)\n",
    "    df_sorted['deltarank_t'] = df_sorted['logit_pergen'].rank() - df_sorted['logit_percycle'].rank()\n",
    "    df_sorted['deltarank_t_abs'] = np.abs(df_sorted['deltarank_t'])\n",
    "    df_sorted = df_sorted.sort_values('deltarank_t_abs', ascending = True)\n",
    "    select = df_sorted.index[[-1]]\n",
    "    \n",
    "    ### plot selection coefficients\n",
    "    fig,ax = plt.subplots(figsize = (0.3*FIGHEIGHT_TRIPLET,0.3*FIGHEIGHT_TRIPLET))\n",
    "\n",
    "    x = df_to_calc['logit_percycle']\n",
    "    y = df_to_calc['logit_pergen']\n",
    "\n",
    "    ax.scatter(x,y, color = color_t) #color_logit_pergen)\n",
    "    \n",
    "    ## calculate correlation\n",
    "    rho = df_stats.at[t,'rho_percycle-pergen'] \n",
    "    r = df_stats.at[t,'r_percycle-pergen']\n",
    "\n",
    "    #title = fr'Pearson $r={r:.2f}$, Spearman $\\rho={rho:.2f}$'\n",
    "    title = fr'$\\rho={rho:.2f}$' # short title\n",
    "    ax.set_title(title, loc = 'center')\n",
    "    ax.tick_params(labelleft=False, labelbottom = False, left =False, bottom=False)\n",
    "\n",
    "    ax.set_xlabel(r'$s^{\\mathrm{logit}}_{\\mathrm{cycle}}$')\n",
    "    ax.set_ylabel(r'$s^{\\mathrm{logit}}_{\\mathrm{gen}}$')\n",
    "\n",
    "    fig.savefig(FIG_DIR + f\"snapshot_percycle-vs-pergen_at_t={t}.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n",
    "    \n",
    "    ## plot LFCS\n",
    "    fig,ax = plt.subplots(figsize = (0.8*FIGHEIGHT_TRIPLET,0.8*FIGHEIGHT_TRIPLET))\n",
    "    \n",
    "    x = df_to_calc['logfc_wt']\n",
    "    y = df_to_calc['logfc_mut']\n",
    "    \n",
    "\n",
    "    ax.scatter(x,y, color = color_t, zorder = 2)\n",
    "    \n",
    "    ## make plot square\n",
    "    #ax.set_xlim(lfc_xmin,lfc_xmax)\n",
    "    #ax.set_ylim(lfc_ymin,lfc_ymax)\n",
    "    \n",
    "    ### find axis limits\n",
    "    xmin, xmax = ax.get_xlim()\n",
    "    fcwt_vec = np.linspace(xmin,xmax, num = 100) \n",
    "\n",
    "\n",
    "    ## plot red cone for a select point\n",
    "    #select = range(df_to_calc.shape[0])\n",
    "    for i in select:\n",
    "    #for i in df_sorted.index:\n",
    "        A, B = float(df_sorted.loc[i, 'logfc_wt']), float(df_sorted.loc[i, 'logfc_mut'])\n",
    "        ax.scatter(A,B,s=200,color ='tab:blue', zorder = -1, alpha = 0.25)\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",
    "    ax.set_xlim(xmin,xmax)\n",
    "    ax.set_ylabel(r'mutant log fold-change $\\mathrm{LFC}_{\\mathrm{mut}}$')\n",
    "    ax.set_xlabel(r'wild-type log fold-change $\\mathrm{LFC}_{\\mathrm{wt}}$')\n",
    "    #ax.tick_params(labelleft=False, labelbottom = False, left =False, bottom=False)\n",
    "\n",
    "\n",
    "    \n",
    "    fig.savefig(FIG_DIR + f\"snapshot_logfc-wt-vs-logfc-mut_at_t={t}.pdf\", DPI = DPI, bbox_inches = 'tight', pad_inches = PAD_INCHES)\n",
    "    \n",
    "    "
   ]
  },
  {
   "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
}