{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic sample processing and analysis\n",
    "In this notebook we will analyse the outputs from the profiles after being run by taxpasta. \n",
    "\n",
    "First we perform some preprocessing to the data, such as simplifying the naming of the taxon hierarchy, or normalising the number of counts. \n",
    "We normalise the counts by dividing the number of reads in the FASTQ by the mean number of reads across all FASTQs. \n",
    "This correction helps us remove the bias of FASTQs with more reads. \n",
    "We use this, and not the number of mapped reads because we want a profiler-independent normalisation value; \n",
    "and we mapped to the reads that have NOT been filtered (before nf-core/rnaseq and Bowtie2 mapping) because the relevant number of reads is the one provided directly by the sequencer.\n",
    "\n",
    "\n",
    "Then, we will perform some QC cuts in the datasets. The cuts are in the following sense\n",
    "* Species that have a high CV across different profilers are flagged.\n",
    "* Species that have a low median number / sum of counts across profiles are flagged.\n",
    "* Species that have only counts in one or two profilers are removed. If they only have a record in one profiler, they are directly removed.\n",
    "\n",
    "With that, species that are flagged for 2 or more flags are removed, and the resulting table is created.\n",
    "\n",
    "Lastly, correlation across profilers are computed. For that, the Pearson correlation using normalised counts is computed, and then the mean value across all samples is computed."
   ]
  },
  {
   "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",
    "import seaborn as sns\n",
    "import os\n",
    "from datetime import datetime\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.options.display.float_format = '{:.3f}'.format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# VARS\n",
    "dir_diversity_output = '../results_diversity'\n",
    "dir_reads_fastq = '../data/EM_EVPools/control_sample/'\n",
    "fastq_basename = 'POOL.fastq.gz'\n",
    "dir_results_profiling = '../results_profiling/EM_EVPools'\n",
    "pools_file=\"../data/EM_EVPools/samples_profiling.txt\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "today = datetime.today().strftime('%Y-%m-%d')\n",
    "os.makedirs(f'{dir_diversity_output}/{today}', exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GENERAL VARIABLES\n",
    "POOL_list = !cat {pools_file}\n",
    "POOL_list_control = ['ACIDOLA', 'BLACTIS']\n",
    "profilers = ['kaiju', 'kraken_2', 'krakenuniq', 'centrifuge']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MINOR PROCESSING FUNCTIONS\n",
    "def process_df(df):\n",
    "    # Reverses the philogenetic order of the taxa, and removes \"root\" and \"cellular organisms\" labels\n",
    "    lineage_vals = df.lineage.values\n",
    "    new_lineage_vals = []\n",
    "    for lineage in lineage_vals:\n",
    "        val = ';'.join(lineage.split(';')[::-1]).replace('root;', '').replace('cellular organisms;', '')\n",
    "        new_lineage_vals.append(val)\n",
    "    df.lineage = new_lineage_vals\n",
    "    return df\n",
    "\n",
    "def get_FASTQ_len(POOL_list):\n",
    "    # Obtains a list of the size of the FASTQs before profiling.\n",
    "    dict_FASTQ_len = {}\n",
    "    \n",
    "    for POOL in POOL_list:\n",
    "        n_counts_fastq = !gzip -dc {dir_reads_fastq}/{fastq_basename.replace('POOL', POOL)} | wc -l\n",
    "        dict_FASTQ_len[POOL] = int(n_counts_fastq[0])\n",
    "    \n",
    "    return dict_FASTQ_len\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MAJOR PROCESSING FUNCTIONS\n",
    "\n",
    "def create_POOL_table(POOL, profilers, dict_FASTQ_len, cutoff_NA, cutoff_CV, cutoff_min_reads):\n",
    "    index_df_POOL = [] # index of taxonomy members to create a merged table with all methods\n",
    "    list_df_pools = []\n",
    "\n",
    "    # Loading the tables. We (1) reverse the phylogenetic lineage to start by root and and by species, (2) sort index by taxonomy and (3) rename count column to the method\n",
    "    for tax_method in profilers:\n",
    "        df_POOL_method = process_df(pd.read_csv(f'{dir_results_profiling}/{tax_method}/{POOL}.report.standardised', sep='\\t', index_col='taxonomy_id')).rename(columns={'count': tax_method})\n",
    "        list_df_pools.append(df_POOL_method)\n",
    "        index_df_POOL += df_POOL_method.index.tolist()\n",
    "    \n",
    "    index_df_POOL = list(set(index_df_POOL))\n",
    "    \n",
    "    # Creating the table and merging count columns\n",
    "    df_POOL = pd.DataFrame(index=index_df_POOL, columns= ['name', 'lineage'] + profilers)\n",
    "\n",
    "    for df_POOL_x, name_x in zip(list_df_pools, profilers):\n",
    "        df_POOL.loc[df_POOL_x.index, ['name', 'lineage', name_x]] = df_POOL_x.loc[df_POOL_x.index, ['name', 'lineage', name_x]]\n",
    "\n",
    "    # Apply normalisation based on other pools - this step takes the number of reads of the FASTQs and corrects the counts by dividing it by the mean number of counts, so that a a FASTQ with more total reads has fewer normalised counts\n",
    "    cols_POOL_norm = [f'{i}_norm' for i in profilers]\n",
    "    mean_FASTQ_len = np.median(np.array(list(dict_FASTQ_len.values())))\n",
    "    correction_factor = mean_FASTQ_len / dict_FASTQ_len[POOL]\n",
    "\n",
    "    for col_POOL in profilers:\n",
    "        df_POOL[f'{col_POOL}_norm'] = df_POOL[col_POOL] * correction_factor\n",
    "\n",
    "    # Calculate simple stats\n",
    "    df_POOL.loc[:, 'raw_mean'] = np.mean(df_POOL.loc[:, profilers], axis=1).astype(float)\n",
    "\n",
    "    df_POOL.loc[:, 'mean'] = np.mean(df_POOL.loc[:, cols_POOL_norm], axis=1).astype(float)\n",
    "    df_POOL.loc[:, 'std'] = np.std(df_POOL.loc[:, cols_POOL_norm], axis=1)\n",
    "\n",
    "\n",
    "    df_POOL = df_POOL[df_POOL['mean'] > 0]  # Filter step because in some cases it was 0 and CV y¡would yield NAN\n",
    "\n",
    "    df_POOL.loc[:, 'CV'] = df_POOL.loc[:, 'std'] / df_POOL.loc[:, 'mean']\n",
    "\n",
    "\n",
    "    # Calculate the number of reads in the fastq to obtain the relative abundance of the counts\n",
    "    df_POOL.loc[:, 'mean (%)'] = 100 * df_POOL.loc[:, 'mean'] / mean_FASTQ_len\n",
    "\n",
    "\n",
    "    df_POOL = df_POOL.sort_values(by='mean', ascending=False)\n",
    "\n",
    "\n",
    "    \n",
    "    # We use some quality metrics to flag and remove \"bad quality\" samples:\n",
    "    #   cutoff_CV to flag species that have very variable counts across profilers\n",
    "    #   cutoff_min_reads and cutoff_min_sum_reads to flag species that have a low count in one and in all profiling counts.\n",
    "    #   cutoff_NA to remove species that are only present in 2 or fewer samples\n",
    "    cutoff_median_reads = int(cutoff_min_reads * 2.5)\n",
    "    df_POOL[['quality_CV', 'quality_min_reads', 'quality_sum_reads', 'quality_NA']] = 0\n",
    "\n",
    "    df_POOL.loc[df_POOL.loc[:, 'CV'] > cutoff_CV, 'quality_CV'] += 1\n",
    "    df_POOL.loc[df_POOL.loc[:, cols_POOL_norm].min(1, skipna=True) < cutoff_min_reads, 'quality_min_reads'] +=  1\n",
    "    df_POOL.loc[df_POOL.loc[:, cols_POOL_norm].median(1, skipna=True) < cutoff_median_reads, 'quality_min_reads'] +=  1\n",
    "    df_POOL.loc[df_POOL.loc[:, cols_POOL_norm].isna().sum(1) > cutoff_NA, 'quality_NA'] += 1\n",
    "    df_POOL.loc[df_POOL.loc[:, cols_POOL_norm].isna().sum(1) == len(profilers) - 1, 'quality_NA'] += 1  # If only one profiler shows the information, we remove it\n",
    "    df_POOL['quality'] = df_POOL.loc[:, ['quality_CV', 'quality_min_reads', 'quality_sum_reads', 'quality_NA']].sum(1)\n",
    "    # df_POOL = df_POOL[~ np.isnan(df_POOL['CV'])]\n",
    "\n",
    "    # We select species with 0 or 1 flag. We are restrictive in that sense to avoid flagging \"low quality\" species\n",
    "    df_POOL_cutoff = df_POOL[df_POOL['quality'] < 2] \n",
    "\n",
    "\n",
    "    df_POOL.to_csv(f'../results_diversity/{today}/{POOL}.diversity_raw_no_correction.tsv', sep='\\t')\n",
    "    df_POOL_cutoff .to_csv(f'../results_diversity/{today}/{POOL}.diversity_cutoff_no_correction.tsv', sep='\\t')\n",
    "\n",
    "    return df_POOL, df_POOL_cutoff\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLOTTING FUNCTIONS\n",
    "\n",
    "def plot_allPOOL_correlations(POOL_list, profilers, corr_method='spearman'):\n",
    "    ncols = int(len(POOL_list) ** 0.5)\n",
    "    nrows = int(len(POOL_list) // ncols) + int(len(POOL_list) % ncols != 0)\n",
    "\n",
    "    \n",
    "\n",
    "    list_mean_heatmaps = []\n",
    "\n",
    "    for type_plot_idx, type_plot in enumerate(['raw', 'cutoff']):\n",
    "        _, axs= plt.subplots(nrows, ncols, figsize=(4 * ncols, 4 * nrows))\n",
    "        corr_mat_list = [] # this is to later sum all correlations and do a big plot\n",
    "        for ax_int, POOL in enumerate(POOL_list):\n",
    "            df_POOL = pd.read_csv(f'../results_diversity/{today}/{POOL}.diversity_{type_plot}_no_correction.tsv', sep='\\t')\n",
    "            df_corr = np.log10(df_POOL.loc[:, profilers].astype(float) + 1)\n",
    "            corr_mat = df_corr.corr(method=corr_method)\n",
    "            corr_mat_list.append(corr_mat)\n",
    "\n",
    "            sns.heatmap(corr_mat, cmap='Blues', annot=True, ax=axs.ravel()[ax_int])\n",
    "            axs.ravel()[ax_int].set_title(POOL)\n",
    "        \n",
    "        plt.suptitle(f'Correlation ({corr_method}, {type_plot})')\n",
    "        plt.tight_layout()\n",
    "\n",
    "        plt.savefig(f'../results_diversity/{today}/correlation_{corr_method}_{type_plot}_no_correction.png', dpi=300)\n",
    "\n",
    "        # Create the mean heatmap\n",
    "\n",
    "        mean_ht = corr_mat_list[0]\n",
    "        NaNs_mat = np.isnan(mean_ht).astype(int)\n",
    "\n",
    "        for corr_mat in corr_mat_list[1:]:\n",
    "            mean_ht = np.nansum(np.dstack((mean_ht,corr_mat)), 2) # This function is used to sum avoiding NaNs\n",
    "            NaNs_mat += np.isnan(corr_mat).astype(int) # To calculate the mean properly, we need to divided buy the number of non NaN elements.\n",
    "        mean_ht /= (len(corr_mat_list) - NaNs_mat)\n",
    "\n",
    "        list_mean_heatmaps.append(mean_ht)\n",
    "\n",
    "    _, axs_all = plt.subplots(1, 2, figsize=(4 * 2, 4 * 1))\n",
    "\n",
    "    for type_plot_idx, type_plot in enumerate(['raw', 'cutoff']):\n",
    "        sns.heatmap(list_mean_heatmaps[type_plot_idx], cmap='Blues', annot=True, ax=axs_all.ravel()[type_plot_idx])\n",
    "        axs_all.ravel()[type_plot_idx].set_title(type_plot)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(f'../results_diversity/{today}/correlation_{corr_method}_mean_no_correction.png', dpi=300)\n",
    "\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dict_FASTQ_len = get_FASTQ_len(POOL_list)\n",
    "# dict_FASTQ_len\n",
    "\n",
    "# dict_FASTQ_len_ctrl = get_FASTQ_len(POOL_list_control)\n",
    "# dict_FASTQ_len_ctrl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_FASTQ_len = {'POOL1': 56548742,\n",
    " 'POOL2': 42205676,\n",
    " 'POOL3': 48792251,\n",
    " 'POOL4': 47145951,\n",
    " 'POOL5': 44694021,\n",
    " 'POOL6': 55640398,\n",
    " 'POOL7': 55697936,\n",
    " 'POOL8': 47826994,\n",
    " 'POOL9': 42852541,\n",
    " 'POOL10': 55067711,\n",
    " 'POOL11': 41656703,\n",
    " 'POOL12': 40367876}\n",
    "\n",
    "\n",
    "\n",
    "dict_FASTQ_len_ctrl = { 'ACIDOLA': 2463392, \n",
    " 'BLACTIS': 2109117}\n",
    "\n",
    "mean_dict_FASTQ_len = np.median(list(dict_FASTQ_len.values()))\n",
    "\n",
    "mean_dict_FASTQ_len_ctrl = np.median(list(dict_FASTQ_len_ctrl.values()))\n",
    "mean_dict_FASTQ_len, mean_dict_FASTQ_len_ctrl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfpool, dfpoolcut = create_POOL_table('ACIDOLA', profilers, dict_FASTQ_len_ctrl, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "display(dfpool)\n",
    "display(dfpoolcut)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfpool, dfpoolcut = create_POOL_table('BLACTIS', profilers, dict_FASTQ_len_ctrl, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "display(dfpool)\n",
    "display(dfpoolcut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfpool, dfpoolcut = create_POOL_table('POOL1', profilers,  dict_FASTQ_len, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "display(dfpool)\n",
    "display(dfpoolcut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfpool, dfpoolcut = create_POOL_table('POOL3', profilers,  dict_FASTQ_len, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "dfpool[dfpool['CV'] > 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfpool, dfpoolcut = create_POOL_table('POOL3', profilers,  dict_FASTQ_len, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "display(dfpool)\n",
    "display(dfpoolcut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfpool, dfpoolcut = create_POOL_table('POOL6', profilers, dict_FASTQ_len, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "display(dfpool)\n",
    "display(dfpoolcut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for pool in POOL_list + ['ACIDOLA', 'BLACTIS']:\n",
    "    dfpool, dfpoolcut = create_POOL_table(pool, profilers, dict_FASTQ_len | dict_FASTQ_len_ctrl, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "    display(pool)\n",
    "    display(dfpoolcut.loc[[True if 'Fungi' in i  else False for i in dfpoolcut['lineage']]])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for pool in POOL_list:\n",
    "    dfpool, dfpoolcut = create_POOL_table(pool, profilers, dict_FASTQ_len, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "    display(pool)\n",
    "    display(dfpoolcut[dfpoolcut['name'] == 'Saccharomycodes'])\n",
    "\n",
    "for pool in ['ACIDOLA', 'BLACTIS']:\n",
    "    dfpool, dfpoolcut = create_POOL_table(pool, profilers, dict_FASTQ_len_ctrl, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "    display(pool)\n",
    "    display(dfpoolcut[dfpoolcut['name'] == 'Saccharomycodes'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for pool in POOL_list:\n",
    "    _, _ = create_POOL_table(pool, profilers, dict_FASTQ_len, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "\n",
    "for pool in ['ACIDOLA', 'BLACTIS']:\n",
    "    _, _ = create_POOL_table(pool, profilers, dict_FASTQ_len_ctrl, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_allPOOL_correlations(POOL_list, profilers, corr_method='spearman')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Export all tables (raw and cut) to an excel samplesheet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_xlsx = pd.ExcelWriter(f\"../results_diversity/{today}/stats_raw_no_correction.xlsx\") \n",
    "cut_xlsx = pd.ExcelWriter(f\"../results_diversity/{today}/stats_cut_no_correction.xlsx\") \n",
    "\n",
    "for pool in POOL_list:\n",
    "    dfpool, dfpoolcut = create_POOL_table(pool, profilers, dict_FASTQ_len, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "    dfpool.to_excel(raw_xlsx, sheet_name=pool, index=False)\n",
    "    dfpoolcut.to_excel(cut_xlsx, sheet_name=pool, index=False)\n",
    "\n",
    "for pool in ['ACIDOLA', 'BLACTIS']:\n",
    "    dfpool, dfpoolcut = create_POOL_table(pool, profilers, dict_FASTQ_len_ctrl, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "    dfpool.to_excel(raw_xlsx, sheet_name=pool, index=False)\n",
    "    dfpoolcut.to_excel(cut_xlsx, sheet_name=pool, index=False)\n",
    "\n",
    "raw_xlsx.close()\n",
    "cut_xlsx.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EVs",
   "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.10.8"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
