{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of in silico samples"
   ]
  },
  {
   "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/artificial_reads/'\n",
    "fastq_basename = 'POOL_R1.fastq.gz'\n",
    "dir_results_profiling = '../results_profiling/artificial_reads'\n",
    "pools_file=\"../data/artificial_reads/samples_profiling.txt\""
   ]
  },
  {
   "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",
    "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.tsv', sep='\\t')\n",
    "    df_POOL_cutoff .to_csv(f'../results_diversity/{today}/{POOL}.diversity_cutoff.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}.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}.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.png', dpi=300)\n",
    "\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_FASTQ_len = {'ARTIFICIAL': 200000128/4}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfpool, dfpoolcut = create_POOL_table('ARTIFICIAL', profilers, dict_FASTQ_len, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "\n",
    "display(dfpool)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(dfpool.sum()[['kaiju', 'kraken_2', 'krakenuniq', 'centrifuge']])\n",
    "display(dfpool[dfpool['name'] != 'Homo'].sum()[['kaiju', 'centrifuge', 'kraken_2', 'krakenuniq',]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfpool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list_taxids = [1912216, 37914, 1678, 1578, 816, 841, 40323, 13687, 29574, 683756, 5052, 55193, 4930, 118259, \n",
    "               5598, 2948645, 1984786, 11646, 2560188, 10375]\n",
    "list_n_reads = [375000, 250000, 250000, 125000, 125000, 125000, 75000, 75000, 50000, 25000, 250000, 125000, \n",
    "                125000, 75000, 75000, 125000, 75000, 75000, 75000, 25000, ]\n",
    "\n",
    "# Some species may not make the cut to dfpoolcut, so we select the index from the\n",
    "dfpoolcut = dfpool.loc[list_taxids + [i for i in dfpoolcut.index if i not in list_taxids]]\n",
    "\n",
    "dfpoolcut = dfpoolcut.fillna(0)\n",
    "dfpoolcut.loc[:, 'reads'] = 0\n",
    "\n",
    "dfpoolcut.loc[list_taxids, 'reads'] = list_n_reads\n",
    "\n",
    "\n",
    "dfpoolcut = dfpoolcut[['name', 'reads', 'kaiju', 'kraken_2', 'krakenuniq', 'centrifuge', 'mean']]\n",
    "dfpoolcut = dfpoolcut[dfpoolcut['name'] != 'Homo']\n",
    "display(dfpoolcut)\n",
    "\n",
    "len(dfpoolcut)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(3.5, 3))\n",
    "\n",
    "g = sns.swarmplot([np.log10(dfpoolcut[dfpoolcut['reads'] > 0]['mean'].values), np.log10(dfpoolcut[dfpoolcut['reads'] == 0]['mean'].values)],)\n",
    "g.set_xticklabels(['Existing genera', 'Non-existing genera'])\n",
    "g.set_ylabel(\"log$_{10}$(Mean number of reads \\nacross profilers)\")\n",
    "\n",
    "display(dfpoolcut[dfpoolcut['reads'] > 0]['mean'].values.mean(), dfpoolcut[dfpoolcut['reads'] == 0]['mean'].values.mean())\n",
    "plt.savefig(f\"../results_diversity/{today}/mean_reads_artificial.png\", dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysing reads assigned to existing species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfpoolcut_with_reads = dfpoolcut[dfpoolcut['reads'] > 0]\n",
    "\n",
    "dfpoolcut_with_reads.loc[:, 'kaiju_norm'] = 1 + (dfpoolcut_with_reads['kaiju'] - dfpoolcut_with_reads['reads']) / dfpoolcut_with_reads['reads']\n",
    "dfpoolcut_with_reads.loc[:, 'kraken_2_norm'] = 1 + (dfpoolcut_with_reads['kraken_2'] - dfpoolcut_with_reads['reads']) / dfpoolcut_with_reads['reads']\n",
    "dfpoolcut_with_reads.loc[:, 'krakenuniq_norm'] = 1 + (dfpoolcut_with_reads['krakenuniq'] - dfpoolcut_with_reads['reads']) / dfpoolcut_with_reads['reads']\n",
    "dfpoolcut_with_reads.loc[:, 'centrifuge_norm'] = 1 + (dfpoolcut_with_reads['centrifuge'] - dfpoolcut_with_reads['reads']) / dfpoolcut_with_reads['reads']\n",
    "dfpoolcut_with_reads.loc[:, 'mean_norm'] = 1 + (dfpoolcut_with_reads['mean'] - dfpoolcut_with_reads['reads']) / dfpoolcut_with_reads['reads']\n",
    "\n",
    "display(dfpoolcut_with_reads)\n",
    "\n",
    "df_stats = pd.DataFrame({'MEDIAN': dfpoolcut_with_reads[['kaiju_norm', 'kraken_2_norm', 'krakenuniq_norm', 'centrifuge_norm', 'mean_norm']].median(), \n",
    "                         'CV_MAD': dfpoolcut_with_reads[['kaiju_norm', 'kraken_2_norm', 'krakenuniq_norm', 'centrifuge_norm', 'mean_norm']].mad() / dfpoolcut_with_reads[['kaiju_norm', 'kraken_2_norm', 'krakenuniq_norm', 'centrifuge_norm', 'mean_norm']].median(),\n",
    "                         'MEAN': dfpoolcut_with_reads[['kaiju_norm', 'kraken_2_norm', 'krakenuniq_norm', 'centrifuge_norm', 'mean_norm']].mean(), \n",
    "                         'CV': dfpoolcut_with_reads[['kaiju_norm', 'kraken_2_norm', 'krakenuniq_norm', 'centrifuge_norm', 'mean_norm']].std() / dfpoolcut_with_reads[['kaiju_norm', 'kraken_2_norm', 'krakenuniq_norm', 'centrifuge_norm', 'mean_norm']].mean()})\n",
    "display(df_stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,3))\n",
    "plt.plot([1,1],[0,3.5], c='#bcbcbc')\n",
    "bw_adjust = 0.25\n",
    "sns.kdeplot(dfpoolcut_with_reads.loc[:, 'kaiju_norm'], bw_adjust=bw_adjust, label='Kaiju')\n",
    "sns.kdeplot(dfpoolcut_with_reads.loc[:, 'kraken_2_norm'], bw_adjust=bw_adjust, label='Kraken2')\n",
    "sns.kdeplot(dfpoolcut_with_reads.loc[:, 'krakenuniq_norm'], bw_adjust=bw_adjust, label='KrakenUniq')\n",
    "sns.kdeplot(dfpoolcut_with_reads.loc[:, 'centrifuge_norm'], bw_adjust=bw_adjust, label='Centrifuge')\n",
    "sns.kdeplot(dfpoolcut_with_reads.loc[:, 'mean_norm'], bw_adjust=bw_adjust, label='Mean')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(7,2.5), dpi=300)\n",
    "plt.plot([1, 1], [0.9, 5.1], c='#bcbcbc')\n",
    "random_noise = np.random.normal(loc=0.0, scale=0.02, size=len(dfpoolcut_with_reads))\n",
    "\n",
    "plt.scatter(dfpoolcut_with_reads.loc[:, 'kaiju_norm'] + random_noise, [1] * len(dfpoolcut_with_reads), label='Kaiju', alpha=0.8)\n",
    "plt.plot([np.mean(dfpoolcut_with_reads.loc[:, 'kaiju_norm'])] * 2, [1 - 0.25, 1 + 0.25], alpha=0.8, c='#454545')\n",
    "\n",
    "plt.scatter(dfpoolcut_with_reads.loc[:, 'kraken_2_norm'] + random_noise, [2] * len(dfpoolcut_with_reads), label='Kraken2', alpha=0.8)\n",
    "plt.plot([np.mean(dfpoolcut_with_reads.loc[:, 'kraken_2_norm'])] * 2, [2 - 0.25, 2 + 0.25], alpha=0.8, c='#454545')\n",
    "\n",
    "plt.scatter(dfpoolcut_with_reads.loc[:, 'krakenuniq_norm'] + random_noise, [3] * len(dfpoolcut_with_reads), label='KrakenUniq', alpha=0.8)\n",
    "plt.plot([np.mean(dfpoolcut_with_reads.loc[:, 'krakenuniq_norm'])] * 2, [3 - 0.25, 3 + 0.25], alpha=0.8, c='#454545')\n",
    "\n",
    "plt.scatter(dfpoolcut_with_reads.loc[:, 'centrifuge_norm'] + random_noise, [4] * len(dfpoolcut_with_reads), label='Centrifuge', alpha=0.8)\n",
    "plt.plot([np.mean(dfpoolcut_with_reads.loc[:, 'centrifuge_norm'])] * 2, [4 - 0.25, 4 + 0.25], alpha=0.8, c='#454545')\n",
    "\n",
    "plt.scatter(dfpoolcut_with_reads.loc[:, 'mean_norm'] + random_noise, [5] * len(dfpoolcut_with_reads), label='Mean')\n",
    "plt.plot([np.mean(dfpoolcut_with_reads.loc[:, 'mean_norm'])] * 2, [5 - 0.25, 5 + 0.25], alpha=0.8, c='#454545')\n",
    "\n",
    "plt.yticks([])\n",
    "plt.legend()\n",
    "plt.xlabel('Efficiency')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each of the profilers outputs a different pattern:\n",
    "- Kaiju and kraken2 seem to work similarly, assigning, in mean, half of the reads to the species. \n",
    "- Kraquenuniq: It shows the highest similarity to the original number of reads, except for especies that does not recognize, probably because they are not integrated in the database.\n",
    "- Centrifuge: In general it detects reads in a much higher proportion (2-3 times) than the expected value, although some species (some of them \"common\" like Dietzia) fail to be assigned.\n",
    "\n",
    "Ironically, despite the underrepresentation of reads in 2/4 profiles, the mean value approaches the expected value more."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysing reads assigned to non-existing species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfpoolcut_non_reads = dfpoolcut[dfpoolcut['reads'] == 0]\n",
    "dfpoolcut_non_reads.loc[:, 'kaiju_norm'] = np.log10(dfpoolcut_non_reads['kaiju'] + 1)\n",
    "dfpoolcut_non_reads.loc[:, 'kraken_2_norm'] = np.log10(dfpoolcut_non_reads['kraken_2'] + 1)\n",
    "dfpoolcut_non_reads.loc[:, 'krakenuniq_norm'] = np.log10(dfpoolcut_non_reads['krakenuniq'] + 1)\n",
    "dfpoolcut_non_reads.loc[:, 'centrifuge_norm'] = np.log10(dfpoolcut_non_reads['centrifuge'] + 1)\n",
    "dfpoolcut_non_reads.loc[:, 'mean_norm'] = np.log10(dfpoolcut_non_reads['mean'] + 1)\n",
    "\n",
    "dfpoolcut_non_reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,3))\n",
    "plt.plot([0, 0],[0,3.5], c='#bcbcbc')\n",
    "bw_adjust = 0.3\n",
    "sns.kdeplot(dfpoolcut_non_reads.loc[:, 'kaiju_norm'], bw_adjust=bw_adjust, label='Kaiju')\n",
    "sns.kdeplot(dfpoolcut_non_reads.loc[:, 'kraken_2_norm'], bw_adjust=bw_adjust, label='Kraken2')\n",
    "sns.kdeplot(dfpoolcut_non_reads.loc[:, 'krakenuniq_norm'], bw_adjust=bw_adjust, label='KrakenUniq')\n",
    "sns.kdeplot(dfpoolcut_non_reads.loc[:, 'centrifuge_norm'], bw_adjust=bw_adjust, label='Centrifuge')\n",
    "sns.kdeplot(dfpoolcut_non_reads.loc[:, 'mean_norm'], bw_adjust=bw_adjust, label='Mean')\n",
    "plt.legend()\n",
    "plt.xticks([0, 1, 2, 3, 4], [0, 10, 100, 1000, 10000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(7,2.5), dpi=300)\n",
    "plt.plot([0, 0], [0.9, 5.1], c='#bcbcbc')\n",
    "random_noise = np.log10(np.random.normal(loc=0.0, scale=0.05, size=len(dfpoolcut_non_reads)) + 1)\n",
    "\n",
    "\n",
    "plt.scatter(dfpoolcut_non_reads.loc[:, 'kaiju_norm'] + random_noise, [1] * len(dfpoolcut_non_reads), label='Kaiju', alpha=0.8)\n",
    "plt.plot([np.mean(dfpoolcut_non_reads.loc[:, 'kaiju_norm'])] * 2, [1 - 0.25, 1 + 0.25], alpha=0.8, c='#454545')\n",
    "\n",
    "plt.scatter(dfpoolcut_non_reads.loc[:, 'kraken_2_norm'] + random_noise, [2] * len(dfpoolcut_non_reads), label='Kraken2', alpha=0.8)\n",
    "plt.plot([np.mean(dfpoolcut_non_reads.loc[:, 'kraken_2_norm'])] * 2, [2 - 0.25, 2 + 0.25], alpha=0.8, c='#454545')\n",
    "\n",
    "plt.scatter(dfpoolcut_non_reads.loc[:, 'krakenuniq_norm'] + random_noise, [3] * len(dfpoolcut_non_reads), label='KrakenUniq', alpha=0.8)\n",
    "plt.plot([np.mean(dfpoolcut_non_reads.loc[:, 'krakenuniq_norm'])] * 2, [3 - 0.25, 3 + 0.25], alpha=0.8, c='#454545')\n",
    "\n",
    "plt.scatter(dfpoolcut_non_reads.loc[:, 'centrifuge_norm'] + random_noise, [4] * len(dfpoolcut_non_reads), label='Centrifuge', alpha=0.8)\n",
    "plt.plot([np.mean(dfpoolcut_non_reads.loc[:, 'centrifuge_norm'])] * 2, [4 - 0.25, 4 + 0.25], alpha=0.8, c='#454545')\n",
    "\n",
    "plt.scatter(dfpoolcut_non_reads.loc[:, 'mean_norm'] + random_noise, [5] * len(dfpoolcut_non_reads), label='Mean')\n",
    "plt.plot([np.mean(dfpoolcut_non_reads.loc[:, 'mean_norm'])] * 2, [5 - 0.25, 5 + 0.25], alpha=0.8, c='#454545')\n",
    "\n",
    "plt.xticks([0, 1, 2, 3, 4], [0, 10, 100, 1000, 10000])\n",
    "plt.yticks([])\n",
    "\n",
    "plt.xlabel('Number of assigned reads')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Export all tables (raw and cut) to an excel samplesheet\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_xlsx = pd.ExcelWriter(f\"../results_diversity/{today}/stats_raw.xlsx\", mode='a') \n",
    "cut_xlsx = pd.ExcelWriter(f\"../results_diversity/{today}/stats_cut.xlsx\", mode='a') \n",
    "\n",
    "\n",
    "dfpool, dfpoolcut = create_POOL_table('ARTIFICIAL', profilers, dict_FASTQ_len, cutoff_NA=1, cutoff_CV=1.0, cutoff_min_reads=100)\n",
    "dfpool.to_excel(raw_xlsx, sheet_name='ARTIFICIAL', index=False)\n",
    "dfpoolcut.to_excel(cut_xlsx, sheet_name='ARTIFICIAL', index=False)\n",
    "\n",
    "\n",
    "\n",
    "raw_xlsx.close()\n",
    "cut_xlsx.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regarding the non-assigned reads, it looks like kaiju and kraken2 (kraken2 more than kaiju) has lower spurious asignment rates; whereas krakenuniq and centrifuge have high assignment rates. However, this assignment rates are, on average, 1-2 orders of magnitude lower than the number of reads of the other species in general."
   ]
  }
 ],
 "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
}
