61_load_LTEE_fitness_data.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)\n",
    "\n",
    "We follow the procedures from Wiser et al. 2013 [here](https://doi.org/10.1126/science.1243357). From the Supplemental Material, we are given the following information. \n",
    "\n",
    "\n",
    "- Summarizing statistical procedures to fit the two models\n",
    "- Models were fit to fitness trajectories using the ‘nls’ package in r. \n",
    "- Model fits were compared using the BIC information criterion scores. These were then converted into an odds ratio. \n",
    "    - Table S1 shows the BIC scores and odds ratios for fits to subsets of the data: a) all 12 populations and all time points, b) excluding 3 populations with incomplete trajectories and c) excluding 6 populations that evolved hypermutability\n",
    "    - Table S2 summarizes BIC scores for fits to individual populations. This also indicates if the population was truncated or a hypermutator \n",
    "    - Table S4 lists the estimated parameters for the power law fit\n",
    "\n",
    "On the bigger picture, there is also the talk from 2013 by Wiser on [Youtube](https://www.youtube.com/watch?v=CmyBn5Cezy4) with 127 views as of September 2022. "
   ]
  },
  {
   "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\n",
    "\n",
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "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_fit/'\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('./data/Concatenated.LTEE.data.all.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": [
    "df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "### split into two subsets, according to the Ara-marker of evolving population\n",
    "is_Ara_positive = np.array(['+' in v for v in df['Population'].values])\n",
    "\n",
    "df_pos = df[is_Ara_positive]\n",
    "df_neg = df[~is_Ara_positive] ### only use Ara negative lineages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## treat Ara-posative population\n",
    "assert all(df_pos['Red.Pop'] == '606') # wild-type is always the red population\n",
    "\n",
    "### re-construct population sizes\n",
    "df_pos['Nwt.0'] = df_pos['Red.0']*df_pos['D.0']\n",
    "df_pos['Nmut.0'] = df_pos['White.0']*df_pos['D.0']\n",
    "df_pos['Nwt.1'] = df_pos['Red.1']*df_pos['D.1']\n",
    "df_pos['Nmut.1'] = df_pos['White.1']*df_pos['D.1']\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## treat Ara-negative population\n",
    "assert all(df_neg['White.Pop'] == '607') # wild-type is always the white population\n",
    "\n",
    "### re-construct population sizes\n",
    "df_neg['Nwt.0'] = df_neg['White.0']*df_neg['D.0']\n",
    "df_neg['Nmut.0'] = df_neg['Red.0']*df_neg['D.0']\n",
    "df_neg['Nwt.1'] = df_neg['White.1']*df_neg['D.1']\n",
    "df_neg['Nmut.1'] = df_neg['Red.1']*df_neg['D.1']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### join\n",
    "\n",
    "df = df_pos.append(df_neg)\n",
    "\n",
    "## reconstruct frequencies\n",
    "df['xmut.0'] = df['Nmut.0']/(df['Nwt.0'] + df['Nmut.0'])\n",
    "df['xmut.1'] = df['Nmut.1']/(df['Nwt.1']  + df['Nmut.1'])\n",
    "\n",
    "## reconstruct log fold-changes\n",
    "df['logfc_mut'] = np.log(df['Nmut.1']/df['Nmut.0'])\n",
    "df['logfc_wt']  =  np.log(df['Nwt.1']/df['Nwt.0'])\n",
    "\n",
    "## reconstuct selection coefficients based on logit\n",
    "df['logit_percycle_from_freq'] = np.log(df['xmut.1']/(1-df['xmut.1'])) - np.log(df['xmut.0']/(1-df['xmut.0'])) \n",
    "df['logit_percycle']  = df['logfc_mut'] - df['logfc_wt']\n",
    "df['logit_pergen']  = df['logit_percycle']/df['logfc_wt']\n",
    "df['W'] = df['logit_pergen'] + 1\n",
    "\n",
    "## reconstruct selection coefficients based on log\n",
    "df['log_percycle']  =  np.log(df['xmut.1']) - np.log(df['xmut.0']) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### check that  my number is consistent with existing value for 'Fitness' in the dataset\n",
    "np.allclose(df['W'], df['Fitness'],equal_nan=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## manual check\n",
    "df['gap'] = df['W'] - df['Fitness']\n",
    "print(df['gap'].max())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Investigate truncation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop2truncation = dict()\n",
    "list_pop = set(df['Population'])\n",
    "for pop in list_pop:\n",
    "    truncation = df.loc[df['Population']==pop,'Generation'].max()\n",
    "    pop2truncation[pop] = truncation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop2truncation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare the data for plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "### split into two subsets, according to the Ara-marker of evolving population\n",
    "is_Ara_positive = np.array(['+' in v for v in df['Population'].values])\n",
    "\n",
    "df_pos = df[is_Ara_positive]\n",
    "df_neg = df[~is_Ara_positive] ### only use Ara negative lineages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## treat Ara-posative population\n",
    "assert all(df_pos['Red.Pop'] == '606') # wild-type is always the red population\n",
    "\n",
    "### re-construct population sizes\n",
    "df_pos['Nwt.0'] = df_pos['Red.0']*df_pos['D.0']\n",
    "df_pos['Nmut.0'] = df_pos['White.0']*df_pos['D.0']\n",
    "df_pos['Nwt.1'] = df_pos['Red.1']*df_pos['D.1']\n",
    "df_pos['Nmut.1'] = df_pos['White.1']*df_pos['D.1']\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## treat Ara-negative population\n",
    "assert all(df_neg['White.Pop'] == '607') # wild-type is always the white population\n",
    "\n",
    "### re-construct population sizes\n",
    "df_neg['Nwt.0'] = df_neg['White.0']*df_neg['D.0']\n",
    "df_neg['Nmut.0'] = df_neg['Red.0']*df_neg['D.0']\n",
    "df_neg['Nwt.1'] = df_neg['White.1']*df_neg['D.1']\n",
    "df_neg['Nmut.1'] = df_neg['Red.1']*df_neg['D.1']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### join\n",
    "\n",
    "df = df_pos.append(df_neg)\n",
    "\n",
    "## reconstruct frequencies\n",
    "df['xmut.0'] = df['Nmut.0']/(df['Nwt.0'] + df['Nmut.0'])\n",
    "df['xmut.1'] = df['Nmut.1']/(df['Nwt.1']  + df['Nmut.1'])\n",
    "\n",
    "## reconstruct log fold-changes\n",
    "df['logfc_mut'] = np.log(df['Nmut.1']/df['Nmut.0'])\n",
    "df['logfc_wt']  =  np.log(df['Nwt.1']/df['Nwt.0'])\n",
    "\n",
    "## reconstuct selection coefficients based on logit\n",
    "df['logit_percycle_from_freq'] = np.log(df['xmut.1']/(1-df['xmut.1'])) - np.log(df['xmut.0']/(1-df['xmut.0'])) \n",
    "df['logit_percycle']  = df['logfc_mut'] - df['logfc_wt']\n",
    "df['logit_pergen']  = df['logit_percycle']/df['logfc_wt']\n",
    "df['W'] = df['logit_pergen'] + 1\n",
    "\n",
    "## reconstruct selection coefficients based on log\n",
    "df['log_percycle']  =  np.log(df['xmut.1']) - np.log(df['xmut.0']) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## test\n",
    "fig, axes = plt.subplots(1,2, figsize = (2*FIGHEIGHT_TRIPLET, FIGHEIGHT_TRIPLET))\n",
    "axes[0].scatter(df['logit_percycle_from_freq'], df['logit_percycle'])\n",
    "axes[1].scatter(df['logit_percycle_from_freq'], df['logit_percycle_from_freq'] -df['logit_percycle'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### check that  my number is consistent with existing value for 'Fitness' in the dataset\n",
    "np.allclose(df['W'], df['Fitness'],equal_nan=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## manual check\n",
    "df['gap'] = df['W'] - df['Fitness']\n",
    "print(df['gap'].max())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### ok, can drop these auxiliary columns\n",
    "df = df.drop(['Fitness','logit_percycle_from_freq', 'gap' ], axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## save\n",
    "df.to_csv('./output/LTEE_all_data.csv',index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Handle the replication"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have different numbers of replicates for the different population, we cannot pool all the replicates into one fit. Instead, we calculate the average for each populationa at each timepoint. This way, each population contributes exactly one datapoint at each timepoint (except for later timepoints, where some populations are truncated and do not contribute at all.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will calculate averages of the fitness values, frequencies and LFC values. For the later analysis, we will use the average at the level of fitness values. This ensures that the grand mean (average over all populations) is the same in our procedure as in the original analysis (Wiser et al 2013). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## drop metavariables to simplify the averaging\n",
    "columns_to_drop = ['Red.Pop', 'White.Pop', # we do not want to distinguish these case\n",
    "                   'Rep', # we do not need that marker variable, we do it manually\n",
    "                   'Complete','Mutator.Ever', # these are strings, cannot be averaged\n",
    "                   'Red.0', 'White.0', 'D.0', 'Red.1', 'White.1', 'D.1'] # superfluous information for our analysis\n",
    "\n",
    "## drop columns\n",
    "df_simple = df.drop(columns_to_drop, axis = 1)\n",
    "## reset index\n",
    "df_simple = df_simple.set_index(['Generation', 'Population'])\n",
    "## show shape\n",
    "df_simple.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## create dataframe for avg\n",
    "index = df_simple.index.drop_duplicates() #only on entry per population per timepoint\n",
    "df_avg = pd.DataFrame(index=index, columns = df_simple.columns)\n",
    "\n",
    "## \n",
    "for v in df_avg.index: \n",
    "    data = df_simple.loc[v] ## get all replicate entries\n",
    "    df_avg.loc[v] = data.mean(axis = 0, skipna = True)\n",
    "    df_avg.at[v,'#Rep'] = data.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_avg['#Rep'].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## show degree of replication\n",
    "ax = df_avg['#Rep'].hist(log = True, bins = np.arange(-0.5,12.5, step=1))\n",
    "ax.set_xlabel('no. replicates')\n",
    "ax.set_ylabel('no. evolved lineages\\n(specific timepoint & replicate)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## save\n",
    "df_avg.to_csv('./output/LTEE_averaged_data.csv',index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### alternative script for averaging\n",
    "\n",
    "\n",
    "### remove populations \n",
    "df_subset = df.copy(deep = True)\n",
    "\n",
    "### we create a new dataframe where each timepoint and population is only represented once\n",
    "df_subset = df_subset.sort_values(by = ['Population', 'Generation', 'Rep']) # first sort for nice look\n",
    "df_averaged = df_subset.drop_duplicates([ 'Population', 'Generation']).copy(deep=True)\n",
    "df_averaged = df_averaged.drop(['Rep'], axis = 1)\n",
    "df_averaged= df_averaged.reset_index()\n",
    "\n",
    "## we average across the number of replicates\n",
    "df_averaged['no_replicates'] = -1 # as a collateral statistic, we count the number of replicates\n",
    "\n",
    "for i in df_averaged.index:\n",
    "    row = df_averaged.loc[i]\n",
    "\n",
    "    pop = row['Population']\n",
    "    gen = row['Generation']\n",
    "\n",
    "    is_gen = np.array([v == gen for v in df['Generation'].values])\n",
    "    is_pop = np.array([v == pop for v in df['Population'].values])\n",
    "    df_replicates = df.loc[is_gen & is_pop]\n",
    "    df_averaged.at[i,'no_replicates'] = df_replicates.shape[0]\n",
    "\n",
    "    for v in ['xmut.0', 'xmut.1', 'logit_percycle', 'logit_pergen', 'log_percycle']:\n",
    "        df_averaged.at[i, v] = df_replicates[~df_replicates[v].isna()][v].mean()\n",
    "\n",
    "\n",
    "## shift data points for alternative statistics, based at fitness = 1\n",
    "df_averaged['logit_percycle+1'] = df_averaged['logit_percycle'] +1\n",
    "df_averaged['log_percycle+1'] = df_averaged['log_percycle'] +1\n",
    "df_averaged['logit_pergen+1'] = df_averaged['logit_pergen'] +1\n",
    "\n"
   ]
  }
 ],
 "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
}