{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import seaborn as sns\n",
    "import os\n",
    "from datetime import datetime\n",
    "from matplotlib.colors import Normalize, ListedColormap\n",
    "from scipy.stats import linregress, pearsonr\n",
    "from scipy.stats import mannwhitneyu, wilcoxon\n",
    "from statsmodels.stats.multitest import multipletests\n",
    "\n",
    "from statannotations.Annotator import Annotator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.options.display.float_format = '{:.3f}'.format\n",
    "pd.set_option('display.max_columns', None)\n",
    "\n",
    "plt.rcParams['figure.dpi']=170\n",
    "\n",
    "plt.rc('axes', linewidth=0.65)  # Adjust the line width of plot frames\n",
    "plt.rc('xtick.major', width=0.0)\n",
    "plt.rc('ytick.major', width=0.0)\n",
    "\n",
    "\n",
    "DPI=250"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from list_vars import LIST_PROFILERS, DIR_FIGURES, RESULTS_DIR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# In silico sample analysis\n",
    "\n",
    "In this notebook we are going to do an analysis on the *in silico* samples, where we are going to study several variables.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How many reads are incorrectly mapped if we do not perfom a host mapping step?\n",
    "\n",
    "It has been reported that not mapping to human databases before profiling increases the number of reads assigned to other organisms. \n",
    "\n",
    "In this case, we are going to do 3 checks with the *in silico* dataset using pass2 (profiling after 2-time host mapping) and pass0 (direct profiling withou host mapping), and we are going to check the influence in parameter sensitivity:\n",
    "-  We are going to see what is the total number of reads mapped to the human dataset, and what is the offset left unmapped which should have been mapped to human.\n",
    "    - We are also going to do the same with the microbial reads, and see if more microbial reads have been assigned to the pass0 dataset.\n",
    "\n",
    "Later in the analysis we are going to do two additional analyses:\n",
    "-  We are going to see the number of species present in total between pass0 and pass2, and their jaccard index.\n",
    "- We are going to calculate the ratio between the number of reads in pass0 and pass2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_host_map_info = pd.read_csv(f'{RESULTS_DIR}/counts/mapping_counts.txt', sep='\\t').set_index('SAMPLE')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "artificial_taxid_counts = pd.read_csv('table_artificial_taxid.csv', sep=';', names=['species', 'taxid', 'reads'])\n",
    "artificial_taxid_counts['reads_true'] = (artificial_taxid_counts['reads'] / 2).astype(int)\n",
    "\n",
    "n_true_human_reads = int(artificial_taxid_counts['reads_true'].iloc[0])\n",
    "n_true_human_reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_mapped_reads_1and2_maps = df_host_map_info.loc['ARTIFICIAL', '1st_mapped'] + df_host_map_info.loc['ARTIFICIAL', '2nd_mapped']\n",
    "\n",
    "print(f'There is a total of {n_mapped_reads_1and2_maps} reads mapped to human during the 1st and 2nd map, which represents around {100 * n_mapped_reads_1and2_maps/n_true_human_reads} % of the total number of reads ({n_true_human_reads}).')\n",
    "print(f'There is a total of {n_true_human_reads - n_mapped_reads_1and2_maps} reads remaining to be mapped.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_host_profile_info = pd.read_csv(f'{RESULTS_DIR}/counts/profiling_counts_ARTIFICIAL.txt', sep='\\t')\n",
    "df_host_profile_info_artificial = df_host_profile_info[df_host_profile_info['SAMPLE'] == 'ARTIFICIAL']\n",
    "\n",
    "df_host_profile_info_artificial['mapped_human_1_2_maps'] = 0\n",
    "df_host_profile_info_artificial.loc[df_host_profile_info_artificial['pass'] == 2, 'mapped_human_1_2_maps'] = n_mapped_reads_1and2_maps\n",
    "\n",
    "df_host_profile_info_artificial['mapped_human_total'] = df_host_profile_info_artificial['mapped_human_1_2_maps'] + df_host_profile_info_artificial['mapped_human']\n",
    "df_host_profile_info_artificial['total_reads'] = df_host_profile_info_artificial['mapped_human_total'] + df_host_profile_info_artificial['mapped_others'] + df_host_profile_info_artificial['unmapped']\n",
    "\n",
    "df_host_profile_info_artificial['observed_human_prop'] = df_host_profile_info_artificial['mapped_human_total'] / df_host_profile_info_artificial['total_reads']\n",
    "df_host_profile_info_artificial['observed_others_prop'] = df_host_profile_info_artificial['mapped_others'] / df_host_profile_info_artificial['total_reads']\n",
    "df_host_profile_info_artificial['observed_unmapped_prop'] = df_host_profile_info_artificial['unmapped'] / df_host_profile_info_artificial['total_reads']\n",
    "\n",
    "df_host_profile_info_artificial['expected_human_prop'] = n_true_human_reads / artificial_taxid_counts['reads_true'].sum() # 0.8\n",
    "df_host_profile_info_artificial['expected_others_prop'] = 1 - n_true_human_reads / artificial_taxid_counts['reads_true'].sum() # 0.8\n",
    "\n",
    "df_host_profile_info_artificial['calculated_unmapped_human_prop'] = df_host_profile_info_artificial['expected_human_prop'] - df_host_profile_info_artificial['observed_human_prop']\n",
    "df_host_profile_info_artificial['calculated_unmapped_others_prop'] = df_host_profile_info_artificial['expected_others_prop'] - df_host_profile_info_artificial['observed_others_prop']\n",
    "\n",
    "df_host_profile_info_artificial['proportion_mapped_other_reads'] = df_host_profile_info_artificial['observed_others_prop'] /  df_host_profile_info_artificial['expected_others_prop']\n",
    "\n",
    "\n",
    "for profiler in LIST_PROFILERS:\n",
    "    display(profiler)\n",
    "    display(df_host_profile_info_artificial[df_host_profile_info_artificial['profiler'] == profiler])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "custom_palette = [\"#648FFF\", \"#785EF0\", \"#DC267F\", \"#FE6100\", \"#FFB000\", \"#848484\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 2, figsize=(8, 3))\n",
    "\n",
    "# 1A) Check if there are differences in human read assignment.\n",
    "\n",
    "# Step 1: Calculate the differences between pass 2 and pass 0 for each profiler and mode\n",
    "\n",
    "pass_diff = (\n",
    "    df_host_profile_info_artificial.pivot_table(\n",
    "        index=[\"profiler\", \"mode\"], columns=\"pass\", values=\"observed_human_prop\"\n",
    "    )\n",
    "    .reset_index()\n",
    ")\n",
    "\n",
    "# Ensure column names are integers\n",
    "pass_diff.columns.name = None  # Remove the columns' name from pivot_table\n",
    "pass_diff.columns = ['profiler', 'mode', 0, 2]  # Explicitly rename columns\n",
    "\n",
    "# Calculate the difference\n",
    "pass_diff[\"difference\"] = 100 * (pass_diff[2] - pass_diff[0])\n",
    "\n",
    "\n",
    "\n",
    "# Step 2: Plot the differences using a lineplot\n",
    "g = sns.lineplot(\n",
    "    data=pass_diff,\n",
    "    x=\"mode\",\n",
    "    y=\"difference\",\n",
    "    hue=\"profiler\",\n",
    "    marker=\"o\",\n",
    "    palette=custom_palette,\n",
    "    ax=axs[0],\n",
    "    legend=False\n",
    ")\n",
    "\n",
    "axs[0].axhline(0, color=\"#848484\", linestyle=\"--\", linewidth=0.8)\n",
    "axs[0].set_title(\"Human assignment\", fontsize=14)\n",
    "axs[0].set_xlabel(\"Mode\", fontsize=12)\n",
    "axs[0].set_ylabel(\"Diff pass2 - pass0 (%)\", fontsize=12)\n",
    "\n",
    "\n",
    "\n",
    "# 1B) Check if there are differences in non-human read assignment.\n",
    "\n",
    "pass_diff = (\n",
    "    df_host_profile_info_artificial.pivot_table(\n",
    "        index=[\"profiler\", \"mode\"], columns=\"pass\", values=\"observed_others_prop\"\n",
    "    )\n",
    "    .reset_index()\n",
    ")\n",
    "\n",
    "pass_diff[\"difference\"] = 100 * (pass_diff[2] - pass_diff[0])\n",
    "\n",
    "\n",
    "\n",
    "# Step 2: Plot the differences using a lineplot\n",
    "\n",
    "h = sns.lineplot(\n",
    "    data=pass_diff,\n",
    "    x=\"mode\",\n",
    "    y=\"difference\",\n",
    "    hue=\"profiler\",\n",
    "    marker=\"o\",\n",
    "    palette=custom_palette,\n",
    "    ax=axs[1],\n",
    ")\n",
    "\n",
    "sns.despine(top=True, right=True)\n",
    "\n",
    "axs[1].axhline(0, color=\"gray\", linestyle=\"--\", linewidth=0.8)\n",
    "axs[1].set_title(\"Non-human assignment\", fontsize=12)\n",
    "axs[1].set_xlabel(\"Mode\", fontsize=12)\n",
    "axs[1].set_ylabel(\"\", fontsize=12)\n",
    "\n",
    "for ax in axs:\n",
    "    ax.set_xticks(range(1, 10))\n",
    "\n",
    "axs[1].legend(title=\"Profiler\", frameon=False, bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=10)\n",
    "plt.tight_layout()\n",
    "plt.rc('axes', linewidth=0.65)\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/figA.{format}', dpi=DPI, )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**What do we see here?**\n",
    "- The number of reads assigned to humans without host mapping is very variable depending on the profiler. Centrifuge, krakenuniq and ganon and map human reads correctly, whereas kaiju, kraken2 fail to map the reads to human. The differences tend to decrease with the sensitivity mode, that is, paradoxically, a more strict read assignment leads to an improved number of human-mapped reads. However, this makes sense because more reads are assigned in general, and thus both human and non-human reads are mapped.\n",
    "- However, this difference does not occur in non-human species. In general, non-human species are assigned equally with or without host mapping. This is interesting because we would expect a higher amount of reads assigned to non-human species originating from a false positive assignment of human reads, but seems not to be the case, even in profilers that have a high ammount of unmapped human reads.\n",
    "    - Still, we have to take into acount that the profiler databases include a host mapping step.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table_artificial_taxcounts = pd.read_csv('../../src/version_2/table_artificial_taxid.csv', sep=';', names=['species', 'taxid', 'count'])\n",
    "table_artificial_taxcounts = table_artificial_taxcounts[table_artificial_taxcounts['taxid'] != 9606]\n",
    "table_artificial_taxcounts['abundance'] = 100 * table_artificial_taxcounts['count'] / table_artificial_taxcounts['count'].sum()\n",
    "table_artificial_taxcounts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing detection stats to aswer the questions\n",
    "\n",
    "One of the parameters used during profiling is the mode of the profilers. Each profiler has a different set of parametters to include reads as valid or not. This may results in the detection of false positives and negatives. \n",
    "\n",
    "Here, we are going to study this effect in *in silico* samples to see if there are major changes. We are can measure the effectivity of several variables: \n",
    "- Categorical values: we can use each of the columns in the flag system to check how well were species assigned. We can use the precision (TP/TP + FP), recall (TP/TP+FN) and F1-score (2 x precision x recall / precision + recall) and Cohen's kappa.\n",
    "$$\\kappa = \\frac{p_0-p_e}{1-p_e} \\quad p_0 = \\frac{TP + TN}{TP + FP + FN + TN} \\quad p_e=\\frac{TP + FP}{TP + FP + FN + TN}\\cdot\\frac{TP + FN }{TP + FP + FN + TN} + \\frac{TN + FP}{TP + FP + FN + TN}\\cdot\\frac{TN + FN}{TP + FP + FN + TN}$$\n",
    "\n",
    "- Numerical values: we can use the normalized value and the abundance to see how well are reads classified. For that we can use the expected number of reads and abundance. With that we will calculate the (1) difference between observed and expected categories and (2) the mean error and the (3) mean absolute error:\n",
    "$$(1) \\qquad DIFF_i = 100\\cdot\\frac{x_{obs,i} - x_{exp,i}}{x_{exp,i}}$$ \n",
    "$$(2) \\qquad ME = E[DIFF_i] = \\frac{100}{N}\\sum\\frac{x_{obs,i} - x_{exp,i}}{x_{exp,i}}$$ \n",
    "$$(3A) \\qquad MAE = E[|DIFF_i|] = \\frac{100}{N}\\sum\\frac{|x_{obs,i} - x_{exp,i}|}{x_{exp,i}}$$ \n",
    "$$(3B) \\qquad MAED = \\sigma[|DIFF_i|]$$ \n",
    "- For numerical values we are also going to calculate the pearson correlation between the observed and expected values, using a log10(1+x) transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Categorical values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_nominal_metrics(df_tax_ground_truth, df_flags_observed, column):\n",
    "    list_expected_taxids = list(df_tax_ground_truth['taxid'].astype(int).values)\n",
    "    list_observed_true_taxids = list(df_flags_observed.loc[df_flags_observed[column] == False, 'taxonomy_id'].astype(int).values)\n",
    "    list_observed_false_taxids = list(df_flags_observed.loc[df_flags_observed[column] == True, 'taxonomy_id'].astype(int).values)\n",
    "\n",
    "    TP = len([i for i in list_expected_taxids if i in list_observed_true_taxids])\n",
    "    FN = len([i for i in list_expected_taxids if i not in list_observed_true_taxids])\n",
    "    FP = len([i for i in list_observed_true_taxids if i not in list_expected_taxids])\n",
    "    TN = len([i for i in list_observed_false_taxids if i not in list_expected_taxids])\n",
    "\n",
    "    assert len(set(list_expected_taxids + list_observed_true_taxids + list_observed_false_taxids)) == TP + FN + FP + TN\n",
    "\n",
    "    precision = TP / (TP + FP)\n",
    "    recall = TP / (TP + FN)\n",
    "\n",
    "    try:\n",
    "        f1 = (2 * precision * recall) / (precision + recall)\n",
    "    except:\n",
    "        f1 = 0\n",
    "    \n",
    "    # Create kappa measures\n",
    "    ALL = TP + FN + FP + TN\n",
    "    p0 = (TP + TN) / (ALL)\n",
    "    pe = (TP + FP)/ALL * (TP + FN)/ALL + (TN + FP)/ALL * (TN + FN)/ALL\n",
    "    kappa = (p0 - pe) / (1 - pe)\n",
    "\n",
    "    return precision, recall, f1, kappa, TP, FN, FP, TN\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_selected = ['centrifuge_norm', 'ganon_norm', 'kaiju_norm', 'kmcp_norm', 'kraken2_norm', 'krakenuniq_norm',\n",
    "                    'centrifuge_relab', 'ganon_relab', 'kaiju_relab', 'kmcp_relab', 'kraken2_relab', 'krakenuniq_relab',\n",
    "                    'mean_norm', 'CV_norm', 'mean_relab', 'CV_relab']\n",
    "df_nominal_stats = {'pass': [], 'mode': [], 'S': [], 'column': [], 'precision': [], 'recall': [], 'f1': [], \n",
    "                    'kappa': [], 'TP|FN|FP|TN': []}\n",
    "\n",
    "for passn in [0, 2]:\n",
    "    for mode in range(1, 10):\n",
    "        for S in [0, 1, 2, 3, 4, 5, 6, 7, 10, 15]:\n",
    "            for column in columns_selected: \n",
    "                summary_table_flags = pd.read_csv(f'{RESULTS_DIR}/summary/ARTIFICIAL_pass{passn}_mode{mode}_taxspecies_S{S}.flags.tsv', sep='\\t')\n",
    "                try:\n",
    "                    precision, recall, f1, kappa, TP, FN, FP, TN = calculate_nominal_metrics(table_artificial_taxcounts, summary_table_flags, column)\n",
    "                except KeyError:\n",
    "                    continue \n",
    "\n",
    "                df_nominal_stats['pass'].append(passn)\n",
    "                df_nominal_stats['mode'].append(mode)\n",
    "                df_nominal_stats['S'].append(S)\n",
    "                df_nominal_stats['column'].append(column)\n",
    "\n",
    "                df_nominal_stats['precision'].append(precision)\n",
    "                df_nominal_stats['recall'].append(recall)\n",
    "                df_nominal_stats['f1'].append(f1)\n",
    "                df_nominal_stats['kappa'].append(kappa)\n",
    "                df_nominal_stats['TP|FN|FP|TN'].append((TP, FN, FP, TN))\n",
    "\n",
    "df_nominal_stats = pd.DataFrame(df_nominal_stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_nominal_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_mad(values):\n",
    "    median = np.median(values)\n",
    "    mad = np.median(np.abs(values - median))\n",
    "    return mad\n",
    "\n",
    "def calculate_numerical_metrics(df_tax_ground_truth, df_counts_observed, profiler, suffix):\n",
    "    df_tax_ground_truth = df_tax_ground_truth.copy().set_index('taxid')\n",
    "    df_counts_observed = df_counts_observed.copy().set_index('taxonomy_id')\n",
    "\n",
    "    list_expected_taxids = df_tax_ground_truth.index.astype(int).values\n",
    "    list_observed_true_taxids = df_counts_observed.index.astype(int).values\n",
    "\n",
    "    combined_taxid = np.intersect1d(list_expected_taxids, list_observed_true_taxids)\n",
    "    species = df_tax_ground_truth.loc[combined_taxid, 'species'].values\n",
    "    observed_counts = df_counts_observed.loc[combined_taxid, f'{profiler}_{suffix}'].fillna(0).values\n",
    "\n",
    "    expected_col = 'count' if suffix == 'norm' else 'abundance'\n",
    "    expected_counts = df_tax_ground_truth.loc[combined_taxid, expected_col].values\n",
    "\n",
    "    diff_counts = 100 * (observed_counts - expected_counts) / expected_counts\n",
    "\n",
    "    MRE_counts = np.mean(diff_counts)\n",
    "    MRED_counts = np.std(diff_counts)\n",
    "    MAE_counts = np.mean(np.abs(diff_counts))\n",
    "    MAED_counts = np.std(np.abs(diff_counts))\n",
    "    MACV_counts = MAED_counts / MAE_counts\n",
    "\n",
    "    if len(combined_taxid) > 35:\n",
    "        x,y = expected_counts, observed_counts \n",
    "\n",
    "        slope, _, r_value, _, _ = linregress(x, y)\n",
    "        r2 = r_value ** 2\n",
    "    else:\n",
    "        slope, r2 = np.nan, np.nan\n",
    "\n",
    "    return diff_counts, MRE_counts, MRED_counts, MAE_counts, MAED_counts, MACV_counts, slope, r2, combined_taxid, species, expected_counts, observed_counts\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for passn in [2]:\n",
    "    for mode in range(5, 6):\n",
    "            for profiler in ['mean']: \n",
    "                try:\n",
    "                    summary_table_flags = pd.read_csv(f'{RESULTS_DIR}/summary/ARTIFICIAL_pass{passn}_mode{mode}_taxspecies_S{S}.flags.tsv', sep='\\t')\n",
    "                    summary_table_counts = pd.read_csv(f'{RESULTS_DIR}/summary/ARTIFICIAL_pass{passn}_mode{mode}_taxspecies_S{S}.diversity.tsv', sep='\\t')\n",
    "                    diff_counts, MRE_counts, MRED_counts, MAE_counts, MAED_counts, MACV_counts, corr_counts, rmse_counts, taxids_counts, species_counts, expected_counts, observed_counts = \\\n",
    "                        calculate_numerical_metrics(table_artificial_taxcounts, summary_table_counts, \\\n",
    "                                                                        profiler, suffix='norm')\n",
    "                    diff_abundance, MRE_abundance, MRED_abundance, MAE_abundance, MAED_abundance, MACV_abundance, corr_abundance, rmse_abundance, taxids_abundance, species_abundance, expected_abundance, observed_abundance = \\\n",
    "                        calculate_numerical_metrics(table_artificial_taxcounts, summary_table_counts, \\\n",
    "                                                                        profiler, suffix='relab')\n",
    "                except KeyError:\n",
    "                    raise\n",
    "\n",
    "df = pd.DataFrame({'species': species_counts, 'expected_counts': expected_counts, 'observed_counts': observed_counts, 'diff_abundance': diff_counts})\n",
    "df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_numerical_stats = {'pass': [], 'mode': [], 'profiler': [], \n",
    "                    'diff_counts': [], 'MRE_counts': [],'MRED_counts': [], 'MAE_counts': [], 'MAED_counts': [], 'MACV_counts': [], \n",
    "                    'corr_counts': [], 'R2_counts': [], 'taxid_counts': [], \n",
    "                    'species_counts': [], 'expected_counts': [], 'observed_counts': [], \n",
    "                    'diff_abundance': [], 'MRE_abundance': [],'MRED_abundance': [], 'MAE_abundance': [], 'MAED_abundance': [], 'MACV_abundance': [], \n",
    "                    'corr_abundance': [], 'R2_abundance': [], 'taxid_abundance': [], \n",
    "                    'species_abundance': [], 'expected_abundance': [], 'observed_abundance': [],  }\n",
    "\n",
    "for passn in [0, 2]:\n",
    "    for mode in range(1, 10):\n",
    "            for profiler in LIST_PROFILERS + ['mean']: \n",
    "                summary_table_flags = pd.read_csv(f'{RESULTS_DIR}/summary/ARTIFICIAL_pass{passn}_mode{mode}_taxspecies_S{S}.flags.tsv', sep='\\t')\n",
    "                summary_table_counts = pd.read_csv(f'{RESULTS_DIR}/summary/ARTIFICIAL_pass{passn}_mode{mode}_taxspecies_S{S}.diversity.tsv', sep='\\t')\n",
    "                try:\n",
    "                    diff_counts, MRE_counts, MRED_counts, MAE_counts, MAED_counts, MACV_counts, corr_counts, rmse_counts, taxids_counts, species_counts, expected_counts, observed_counts = \\\n",
    "                        calculate_numerical_metrics(table_artificial_taxcounts, summary_table_counts, \\\n",
    "                                                                        profiler, suffix='norm')\n",
    "                    diff_abundance, MRE_abundance, MRED_abundance, MAE_abundance, MAED_abundance, MACV_abundance, corr_abundance, rmse_abundance, taxids_abundance, species_abundance, expected_abundance, observed_abundance = \\\n",
    "                        calculate_numerical_metrics(table_artificial_taxcounts, summary_table_counts, \\\n",
    "                                                                        profiler, suffix='relab')\n",
    "                except KeyError:\n",
    "                    raise\n",
    "                    continue \n",
    "                \n",
    "                df_numerical_stats['pass'].append(passn)\n",
    "                df_numerical_stats['mode'].append(mode)\n",
    "                df_numerical_stats['profiler'].append(profiler)\n",
    "\n",
    "                df_numerical_stats['diff_counts'].append(diff_counts)\n",
    "                df_numerical_stats['MRE_counts'].append(MRE_counts)          \n",
    "                df_numerical_stats['MRED_counts'].append(MRED_counts)              \n",
    "                df_numerical_stats['MAE_counts'].append(MAE_counts)                \n",
    "                df_numerical_stats['MAED_counts'].append(MAED_counts)                \n",
    "                df_numerical_stats['MACV_counts'].append(MACV_counts)                \n",
    "                df_numerical_stats['corr_counts'].append(corr_counts)                \n",
    "                df_numerical_stats['R2_counts'].append(rmse_counts) \n",
    "                df_numerical_stats['taxid_counts'].append(taxids_counts)\n",
    "                df_numerical_stats['species_counts'].append(species_counts)\n",
    "                df_numerical_stats['expected_counts'].append(expected_counts)\n",
    "                df_numerical_stats['observed_counts'].append(observed_counts)\n",
    "\n",
    "                df_numerical_stats['diff_abundance'].append(diff_abundance)\n",
    "                df_numerical_stats['MRE_abundance'].append(MRE_abundance)\n",
    "                df_numerical_stats['MRED_abundance'].append(MRED_abundance)\n",
    "\n",
    "                df_numerical_stats['MAE_abundance'].append(MAE_abundance)\n",
    "                df_numerical_stats['MAED_abundance'].append(MAED_abundance)\n",
    "                df_numerical_stats['MACV_abundance'].append(MACV_abundance)\n",
    "                df_numerical_stats['corr_abundance'].append(corr_abundance)                \n",
    "                df_numerical_stats['R2_abundance'].append(rmse_abundance) \n",
    "                df_numerical_stats['taxid_abundance'].append(taxids_abundance)\n",
    "                df_numerical_stats['species_abundance'].append(species_abundance)\n",
    "                df_numerical_stats['expected_abundance'].append(expected_abundance)\n",
    "                df_numerical_stats['observed_abundance'].append(observed_abundance)\n",
    "\n",
    "df_numerical_stats = pd.DataFrame(df_numerical_stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_numerical_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis of kappa/F1\n",
    "\n",
    "F1-score and $\\kappa$ are quite different measures but we observe that they are correlated in this data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute and print correlation (Pearson by default)\n",
    "corr = df_nominal_stats['f1'].corr(df_nominal_stats['kappa'])\n",
    "\n",
    "print(\"Pearson correlation between F1 and Kappa:\", corr)\n",
    "\n",
    "# Scatter plot with the identity line\n",
    "plt.figure(figsize=(4, 4))\n",
    "\n",
    "# Add the identity line (y = x)\n",
    "plt.plot([0, 1], [0, 1], color='red', linestyle='--', alpha=0.5)\n",
    "\n",
    "sns.scatterplot(x='f1', y='kappa', data=df_nominal_stats, s=20)\n",
    "\n",
    "plt.text(0.01, 0.95, f'R$^2$: {corr**2:.5f}')\n",
    "\n",
    "sns.despine(top=True, right=True)\n",
    "\n",
    "# Add title and labels\n",
    "plt.title('F1-score vs Kappa')\n",
    "plt.xlabel('F1-score')\n",
    "plt.ylabel('Kappa')\n",
    "plt.tight_layout()\n",
    "\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/fig0.kappa-f1.{format}', dpi=DPI, bbox_inches='tight' )\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In that sense, we can then use one of the measures to explain the results and don't need the second one. We are going to select the F1 score because it has a more clear interpretability and it is related to precision and recall, which are alrady being used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "def cohen_kappa(tp, fp, fn, tn):\n",
    "    \"\"\"\n",
    "    Compute Cohen's kappa given the confusion matrix counts:\n",
    "      TP = true positives\n",
    "      FP = false positives\n",
    "      FN = false negatives\n",
    "      TN = true negatives\n",
    "    \"\"\"\n",
    "    N = tp + fp + fn + tn  # total\n",
    "    \n",
    "    # Observed agreement\n",
    "    po = (tp + tn) / N\n",
    "    \n",
    "    # Expected agreement\n",
    "    # (actual positives * predicted positives) + (actual negatives * predicted negatives)\n",
    "    p_a = (tp + fn) / N  # actual positives\n",
    "    p_p = (tp + fp) / N  # predicted positives\n",
    "    p_n = (fp + tn) / N  # actual negatives? Actually \"actual positives\" vs \"predicted negatives\" can be spelled out:\n",
    "\n",
    "    p_n_act = (fp + tn) / N   # actual negatives\n",
    "    p_n_pred = (fn + tn) / N  # predicted negatives\n",
    "    \n",
    "    pe = p_a * p_p + p_n_act * p_n_pred\n",
    "    \n",
    "    if np.isclose(pe, 1.0):\n",
    "        # If pe is 1, the denominator (1 - pe) -> 0, can cause issues\n",
    "        return np.nan\n",
    "    \n",
    "    kappa = (po - pe) / (1 - pe)\n",
    "    return kappa\n",
    "\n",
    "def f1_score_simple(tp, fp, fn):\n",
    "    \"\"\"\n",
    "    Compute the F1 score from TP, FP, FN:\n",
    "      F1 = 2 * TP / (2*TP + FP + FN)\n",
    "    \"\"\"\n",
    "    denom = (2 * tp + fp + fn)\n",
    "    if denom == 0:\n",
    "        return np.nan\n",
    "    return 2 * tp / denom\n",
    "\n",
    "# --- Simulation parameters ---\n",
    "N = 3000        # total samples per confusion matrix\n",
    "n_sims = 100000  # number of random simulations\n",
    "\n",
    "f1_values = []\n",
    "kappa_values = []\n",
    "table_values = []\n",
    "\n",
    "for _ in range(n_sims):\n",
    "    # We need TP, FP, FN, TN >= 0 and sum = N.\n",
    "    # One way: sample three random integers and let the fourth be what's left.\n",
    "    # This ensures the sum is N.\n",
    "    \n",
    "    tp = random.randint(0, N)\n",
    "    fp = random.randint(0, N - tp)\n",
    "    fn = random.randint(0, N - tp - fp)\n",
    "    tn = N - tp - fp - fn\n",
    "    \n",
    "    # Compute F1\n",
    "    f1 = f1_score_simple(tp, fp, fn)\n",
    "    # Compute Kappa\n",
    "    kappa = cohen_kappa(tp, fp, fn, tn)\n",
    "    \n",
    "    # We only keep valid calculations (not NaN)\n",
    "    if not np.isnan(f1) and not np.isnan(kappa):\n",
    "        f1_values.append(f1)\n",
    "        kappa_values.append(kappa)\n",
    "        table_values.append((tp, fp, fn, tn))\n",
    "\n",
    "# Convert to numpy arrays for correlation\n",
    "f1_values = np.array(f1_values)\n",
    "kappa_values = np.array(kappa_values)\n",
    "\n",
    "# Compute Pearson correlation\n",
    "corr, p_value = pearsonr(f1_values, kappa_values)\n",
    "\n",
    "print(f\"Number of valid simulations: {len(f1_values)} / {n_sims}\")\n",
    "print(f\"Pearson correlation between F1 and Kappa: {corr:.3f} (p = {p_value:.3e})\")\n",
    "\n",
    "# Plot the scatter of F1 vs Kappa\n",
    "plt.figure(figsize=(7, 5))\n",
    "plt.scatter(f1_values, kappa_values, alpha=0.2, s=10)\n",
    "plt.title(\"Monte Carlo Simulation (TP, FP, FN, TN random)\")\n",
    "plt.xlabel(\"F1 Score\")\n",
    "plt.ylabel(\"Cohen's Kappa\")\n",
    "\n",
    "plt.plot([0, 1], [0, 1])\n",
    "\n",
    "plt.grid(True)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 0.8\n",
    "\n",
    "df = pd.DataFrame({'f1': f1_values, 'kappa': kappa_values, 'table': table_values})\n",
    "\n",
    "df[(df['f1'] > N - 0.03) & (df['f1'] < N + 0.03) & (df['kappa'] < N + 0.03) & (df['kappa'] > N- 0.03)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1\n",
    "\n",
    "df[(df['f1'] > N - 0.03) & (df['f1'] < N + 0.03) & (df['kappa'] < N + 0.03) & (df['kappa'] > N- 0.03)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 0\n",
    "\n",
    "df[(df['f1'] > N - 0.03) & (df['f1'] < N + 0.03) & (df['kappa'] < N + 0.03) & (df['kappa'] > N- 0.03)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How does the S parametter used during curve fitting affect?\n",
    "\n",
    "The S parametter is useful to tweak the detection results, so that we can include more or less species during the flagging step. Since it is a structural parametter, we want to fit it first so that we can answer several other comparisons.\n",
    "\n",
    "To do this we are going to use the nominal variables and their derived statistics.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking recall/precision/F1-score for inclusion/exclusion of species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = [f'{i}_norm' for i in LIST_PROFILERS] + ['mean_norm']\n",
    "modes = range(2, 9)\n",
    "passn = [2]\n",
    "S_values = df_nominal_stats['S'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_df = df_nominal_stats[(df_nominal_stats['pass'].isin(passn)) & \\\n",
    "                             (df_nominal_stats['column'].isin(cols)) & \\\n",
    "                              (df_nominal_stats['mode'].isin(modes))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "melted_df = pd.melt(\n",
    "    subset_df,\n",
    "    id_vars=['mode', 'S', 'column'],\n",
    "    value_vars=['recall', 'precision', 'f1'],\n",
    "    var_name='metric',\n",
    "    value_name='score'\n",
    ")\n",
    "\n",
    "# Create a colormap for 'mode'\n",
    "norm = Normalize(vmin=melted_df['mode'].min(), vmax=melted_df['mode'].max())\n",
    "cmap = plt.cm.viridis  # Choose a colormap (e.g., 'viridis', 'plasma', 'cividis')\n",
    "\n",
    "# Create a FacetGrid: 6x3 grid (row for each profiler, column for each metric)\n",
    "g = sns.FacetGrid(\n",
    "    melted_df, \n",
    "    col='column', \n",
    "    row='metric', \n",
    "    height=2, \n",
    "    sharey=True, \n",
    "    sharex=True\n",
    ")\n",
    "\n",
    "# Map the lineplot to the grid\n",
    "def lineplot_with_cmap(data, **kwargs):\n",
    "    for mode in sorted(data['mode'].unique()):\n",
    "        subset = data[data['mode'] == mode]\n",
    "        plt.plot(subset['S'], subset['score'], label=f\"Mode {mode}\",\n",
    "                 color=cmap(norm(mode)), marker='o')\n",
    "\n",
    "g.map_dataframe(lineplot_with_cmap)\n",
    "\n",
    "# Create a legend for the discrete modes\n",
    "handles = [\n",
    "    plt.Line2D([0], [0], color=cmap(norm(mode)), marker='o', linestyle='', label=f\"Mode {mode}\")\n",
    "    for mode in sorted(melted_df['mode'].unique())\n",
    "]\n",
    "plt.legend(\n",
    "    handles=handles, \n",
    "    title=\"\", \n",
    "    bbox_to_anchor=(1.05, 3), \n",
    "    loc='center left', \n",
    "    frameon=False\n",
    ")\n",
    "\n",
    "# Set x-axis ticks (if you have specific S values)\n",
    "g.set(xticks=subset_df['S'].unique())\n",
    "\n",
    "for ax in g.axes.ravel():\n",
    "    ax.set_title('')\n",
    "\n",
    "# Add axis labels and titles\n",
    "for ax, profiler in zip(g.axes[0, :], melted_df['column'].unique()):\n",
    "    ax.set_title(profiler.replace('_norm', ''))\n",
    "\n",
    "for ax, score in zip(g.axes[:, 0], ['recall', 'precision', 'F1-score']):\n",
    "    ax.set_ylabel(score)\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/figE.{format}', dpi=DPI, bbox_inches='tight')\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.FacetGrid(\n",
    "    melted_df, \n",
    "    col='column', \n",
    "    row='metric', \n",
    "    height=3, \n",
    "    sharey=True, \n",
    "    sharex=True\n",
    ")\n",
    "\n",
    "g.map(sns.boxplot, 'S', 'score')\n",
    "\n",
    "# Set x-axis ticks (if you have specific S values)\n",
    "for ax in g.axes.ravel():\n",
    "    ax.set_title('')\n",
    "\n",
    "# Add axis labels and titles\n",
    "for ax, profiler in zip(g.axes[0, :], melted_df['column'].unique()):\n",
    "    ax.set_title(profiler.replace('_norm', ''))\n",
    "\n",
    "for ax, score in zip(g.axes[:, 0], ['recall', 'precision', 'F1-score']):\n",
    "    ax.set_ylabel(score)\n",
    "\n",
    "plt.subplots_adjust(top=0.9)\n",
    "g.fig.suptitle(\"Metrics by Profiler and Mode\", fontsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The aim of this part of the analysis was to select the \"optimal\" `S` to then make other comparisons and extract proper conclusions. \n",
    "If we look at individual profilers, the aim is not the select the `S` with best F1 score, but to select the smallest `S` that provides a sufficiently high recall, ensuring that we don't lose TP species. This threshold depends on the profiler. For CEN it is 6-7, GAN is 7-10, KAI is 1-2, KR2 is 4-5, KRU is 5-6. We see that at these values the precision drops (expectedly), but it remains stable afterwards for most profilers. Therefore, a value of `S=2` should be sufficient to ensure that the results are correct. \n",
    "\n",
    "**IMPORTANT**: conceptually, if we are choosng the number of the species based on the mean number of reads, using a low `S` is still good, because the number of reads reported by profilers that are not flagged are still considered. That is, the # of reads assigned to one profiler is the same regardless of S. \n",
    "\n",
    "The advantage of using the mean value instead of the individual profilers is that it tends to retrieve a better stability on the precision throughout the S values and modes. In fact, each profiler has an individual dinamic, and the mean value averages them all. Therefore, using the averaged value for `S` allows us to choose the parameter with better predictability, ensuring that we don't choose a very high value, while at the same time keeping the correct number of reads.\n",
    "\n",
    "Therefore, we are going to choose `S=2` and `S=7` for comparisons of robustness with biological samples. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Does pass0/pass2 (no host pre-mapping vs host pre-mapping) affect the detection of the species?\n",
    "\n",
    "For this part we are going to run run analyses:\n",
    "- Retrieve the raw detection of species with the passes, and calculate their jaccard index.\n",
    "- Calculate $\\chi^2$ for each case and see if there are significative differences.\n",
    "- Calculate the Pearson correlation + RMSE for several mode values.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Total number of species and Jaccard index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot Jaccard index between the different species\n",
    "\n",
    "\"\"\"\n",
    "For this part we are going to read all the .standardised.species reports and simply read the number of species and compute a table with the total number of species and the jaccard index\n",
    "\"\"\"\n",
    "\n",
    "def get_n_species(sample, mode, profiler, S=15):\n",
    "    pass0_df = pd.read_csv(f'{RESULTS_DIR}/summary/{sample}_pass0_mode{mode}_taxspecies_S{S}.diversity.tsv', sep='\\t').set_index('taxonomy_id')\n",
    "    pass2_df = pd.read_csv(f'{RESULTS_DIR}/summary/{sample}_pass2_mode{mode}_taxspecies_S{S}.diversity.tsv', sep='\\t').set_index('taxonomy_id')\n",
    "\n",
    "    pass0_profiler = pass0_df[profiler].dropna()\n",
    "    pass2_profiler = pass2_df[profiler].dropna()\n",
    "\n",
    "    pass0_taxids = pass0_profiler.index.values\n",
    "    pass2_taxids = pass2_profiler.index.values\n",
    "\n",
    "    jaccard_index = len(np.intersect1d(pass0_taxids, pass2_taxids)) / len(np.union1d(pass0_taxids, pass2_taxids))\n",
    "\n",
    "    return len(pass0_taxids), len(pass2_taxids), jaccard_index\n",
    "\n",
    "\n",
    "dict_n_species = {'profiler': [],\n",
    "                  'mode': [],\n",
    "                  'pass 0 species': [],\n",
    "                  'pass 2 species': [],\n",
    "                  'jaccard': []}\n",
    "\n",
    "for profiler in LIST_PROFILERS + ['mean']:\n",
    "    for mode in range(1,10):\n",
    "        try:\n",
    "            pass0_n_species, pass2_n_species, jaccard_index = get_n_species('ARTIFICIAL', mode, profiler + '_norm')\n",
    "            dict_n_species['profiler'].append(profiler)\n",
    "            dict_n_species['mode'].append(mode)\n",
    "            dict_n_species['pass 0 species'].append(pass0_n_species)\n",
    "            dict_n_species['pass 2 species'].append(pass2_n_species)\n",
    "            dict_n_species['jaccard'].append(jaccard_index)\n",
    "        except:\n",
    "            print(f'No entry added for profiler {profiler} and mode {mode}')\n",
    "\n",
    "df_n_species = pd.DataFrame(dict_n_species)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc('axes', linewidth=0.65)\n",
    "\n",
    "# Initialize a 3x2 grid for subplots\n",
    "fig, axes = plt.subplots(2, 3, figsize=(11, 6))\n",
    "axes = axes.flatten()\n",
    "\n",
    "# Set Seaborn style\n",
    "sns.set_theme(style=\"white\")\n",
    "\n",
    "# Define custom labels for the shared legend\n",
    "legend_elements = [\n",
    "    plt.Line2D([0], [0], marker='o', color='#648FFF', label='Pass 0', markersize=8, linestyle='None'),\n",
    "    plt.Line2D([0], [0], marker='o', color='#785EF0', label='Pass 2', markersize=8, linestyle='None'),\n",
    "    plt.Line2D([0], [0], marker='o', color='#DC267F', label='Jaccard Index', markersize=8, linestyle='-')\n",
    "]\n",
    "\n",
    "# Plot for each profiler\n",
    "for i, profiler in enumerate(LIST_PROFILERS + ['mean']):\n",
    "    df_profiler = df_n_species[df_n_species['profiler'] == profiler]\n",
    "    ax1 = axes[i]\n",
    "    \n",
    "    # Scatter plots for `pass 0 species` and `pass 2 species`\n",
    "    sns.scatterplot(\n",
    "        data=df_profiler, x='mode', y='pass 0 species', ax=ax1, color='#648FFF', s=50\n",
    "    )\n",
    "    sns.scatterplot(\n",
    "        data=df_profiler, x='mode', y='pass 2 species', ax=ax1, color='#785EF0', s=50\n",
    "    )\n",
    "    \n",
    "    sns.despine(top=True, right=True)\n",
    "    \n",
    "    sns.set_theme(style=\"white\")\n",
    "    \n",
    "    # Configure the primary y-axis\n",
    "    ax1.set_xticks(range(1, 10))\n",
    "    ax1.set_xlabel('Mode')\n",
    "    ax1.set_ylabel('Species Count', color='black')\n",
    "    #ax1.set_yscale('log')  # Optional: Log scale\n",
    "    ax1.set_title(f'{profiler.capitalize()}')\n",
    "    \n",
    "    # Create a secondary y-axis for Jaccard index\n",
    "    ax2 = ax1.twinx()\n",
    "    sns.lineplot(\n",
    "        data=df_profiler, x='mode', y='jaccard', ax=ax2, color='#DC267F', marker='o'\n",
    "    )\n",
    "    ax2.set_ylabel('Jaccard Index', color='#DC267F')\n",
    "\n",
    "    sns.despine(top=True, right=True)\n",
    "# Hide unused axes\n",
    "for j in range(len(LIST_PROFILERS) + 1, len(axes)):\n",
    "    axes[j].axis('off')\n",
    "\n",
    "# Add a shared legend\n",
    "fig.legend(\n",
    "    handles=legend_elements, loc='center right', frameon=False, title=\"\", bbox_to_anchor=(1.15, 0.5)\n",
    ")\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/figB.{format}', dpi=DPI, bbox_inches='tight' )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.lineplot(df_n_species, x='mode', y='jaccard', hue='profiler', marker='o', palette=custom_palette)\n",
    "plt.legend(\n",
    "    loc='center right', frameon=False, bbox_to_anchor=(1.35, 0.5),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of detected species tends to go up with the mode, which is expected (with the exception of krakenuniq). For most profilers the increase is linear, but for centrifuge the increase is exponential, although it may reach the top of detected species at some point.\n",
    "\n",
    "We see that in general Jaccard indexes are very high. We should consider that this plot is not \"strictly\" relevant because many, many species are false positive, so we don't really care about the number of falsely assigned species. Still, it is important to note that pass does not have, grosso modo, a relevant impact. It is interesting to find that for kaiju and centrifuge the mode has an impact. For Kaiju it decreases; this it might be because the number of spurious species has increased; while for centrifuge the jaccard index may increase because their number of species is so high (12000!!) that we may be reaching the maximum of discoverable amount of species, and therefore it is obvious that the Jaccard index will increase in that case. For the rest of profilers the mode does not seem to affect the Jaccard index, regardless of the total number of species detected. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking statistical differences in the truth table.\n",
    "\n",
    "We are going to use the truth tables to compute statistically differential capture of species. We are going to use `S=2` and `S=7` throughout the different modes and profilers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import chi2_contingency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "SMALL_VAL = 0.1\n",
    "\n",
    "df_sub = df_nominal_stats[(df_nominal_stats['column'].isin([f'{i}_norm' for i in LIST_PROFILERS] + ['mean_norm'])) & \n",
    "                          (df_nominal_stats['S'].isin([0, 1, 2, 3, 4, 5, 6, 7, 10, 15]))].copy()\n",
    "\n",
    "# Initialize a list to store chi-squared results\n",
    "chi2_results = []\n",
    "\n",
    "# Iterate over unique combinations of mode, S, and column\n",
    "for (mode, S, column), group in df_sub.groupby(['mode', 'S', 'column']):\n",
    "    try:\n",
    "        # Extract contingency tables for full (TP, FN, FP, TN) and partial (TP, FP, FN)\n",
    "        full_contingency_table = []\n",
    "        partial_contingency_table = []\n",
    "        precision_diff, recall_diff, f1_diff = None, None, None  # Initialize differences\n",
    "        \n",
    "        for p in [0, 2]:\n",
    "            data = group[group['pass'] == p]\n",
    "            if len(data):\n",
    "                (tp, fn, fp, tn) = data.iloc[0]['TP|FN|FP|TN']\n",
    "                full_contingency_table.append([tp, fn, fp, tn])  # Full table\n",
    "                partial_contingency_table.append([tp, fp, fn])   # Partial table\n",
    "\n",
    "                # Calculate differences for precision, recall, and f1\n",
    "                if p == 0:\n",
    "                    precision_0 = data.iloc[0]['precision']\n",
    "                    recall_0 = data.iloc[0]['recall']\n",
    "                    f1_0 = data.iloc[0]['f1']\n",
    "                elif p == 2:\n",
    "                    precision_2 = data.iloc[0]['precision']\n",
    "                    recall_2 = data.iloc[0]['recall']\n",
    "                    f1_2 = data.iloc[0]['f1']\n",
    "\n",
    "        # Compute differences\n",
    "        if 'precision_0' in locals() and 'precision_2' in locals():\n",
    "            precision_diff = precision_2 - precision_0\n",
    "            recall_diff = recall_2 - recall_0\n",
    "            f1_diff = f1_2 - f1_0\n",
    "\n",
    "        # Perform chi-squared test for full contingency table\n",
    "        if len(full_contingency_table) == 2:\n",
    "            chi2_full, p_value_full, _, _ = chi2_contingency(np.array(full_contingency_table) + SMALL_VAL)\n",
    "        \n",
    "        # Perform chi-squared test for partial contingency table\n",
    "        if len(partial_contingency_table) == 2:\n",
    "            chi2_partial, p_value_partial, _, _ = chi2_contingency(np.array(partial_contingency_table) + SMALL_VAL)\n",
    "        \n",
    "        # Store results\n",
    "        chi2_results.append({\n",
    "            'mode': mode,\n",
    "            'S': S,\n",
    "            'column': column,\n",
    "            'chi2_full': chi2_full if len(full_contingency_table) == 2 else None,\n",
    "            'p_value_full': p_value_full if len(full_contingency_table) == 2 else None,\n",
    "            'chi2_partial': chi2_partial if len(partial_contingency_table) == 2 else None,\n",
    "            'p_value_partial': p_value_partial if len(partial_contingency_table) == 2 else None,\n",
    "            'precision_diff': precision_diff,\n",
    "            'recall_diff': recall_diff,\n",
    "            'f1_diff': f1_diff,\n",
    "            'stats_pass0': full_contingency_table[0],\n",
    "            'stats_pass2': full_contingency_table[1],\n",
    "        })\n",
    "    except Exception as e:\n",
    "        print(f\"Error processing {mode}, {S}, {column}: {e}\")\n",
    "\n",
    "# Convert results to a DataFrame\n",
    "df_pass_chi2_stats = pd.DataFrame(chi2_results)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add significance categories based on p-value thresholds\n",
    "alpha = 0.1 \n",
    "\n",
    "df_pass_chi2_stats['significance'] = df_pass_chi2_stats.apply(\n",
    "    lambda row: (\n",
    "        'Both Significant' if row['p_value_full'] < alpha and row['p_value_partial'] < alpha else\n",
    "        'Only Full Significant' if row['p_value_full'] < alpha else\n",
    "        'Only Partial Significant' if row['p_value_partial'] < alpha else\n",
    "        'Neither Significant'\n",
    "    ),\n",
    "    axis=1\n",
    ")\n",
    "\n",
    "# Count points in each quadrant\n",
    "quadrant_counts = {\n",
    "    'Both Significant': len(df_pass_chi2_stats[(df_pass_chi2_stats['p_value_full'] < alpha) & (df_pass_chi2_stats['p_value_partial'] < alpha)]),\n",
    "    'Only Full Significant': len(df_pass_chi2_stats[(df_pass_chi2_stats['p_value_full'] < alpha) & (df_pass_chi2_stats['p_value_partial'] >= alpha)]),\n",
    "    'Only Partial Significant': len(df_pass_chi2_stats[(df_pass_chi2_stats['p_value_full'] >= alpha) & (df_pass_chi2_stats['p_value_partial'] < alpha)]),\n",
    "    'Neither Significant': len(df_pass_chi2_stats[(df_pass_chi2_stats['p_value_full'] >= alpha) & (df_pass_chi2_stats['p_value_partial'] >= alpha)])\n",
    "}\n",
    "\n",
    "# Set up the plots\n",
    "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n",
    "\n",
    "# First plot: chi2_full vs chi2_partial\n",
    "sns.scatterplot(data=df_pass_chi2_stats, x='chi2_full', y='chi2_partial', ax=axs[0])\n",
    "axs[0].set_title('Chi2 Full vs Partial')\n",
    "axs[0].set_xlabel('Chi2 Full')\n",
    "axs[0].set_ylabel('Chi2 Partial')\n",
    "\n",
    "# Second plot: p_value_full vs p_value_partial with significance categories\n",
    "sns.scatterplot(\n",
    "    data=df_pass_chi2_stats,\n",
    "    x='p_value_full',\n",
    "    y='p_value_partial',\n",
    "    hue='significance',\n",
    "    palette={\n",
    "        'Both Significant': '#bc0000',\n",
    "        'Only Full Significant': '#0000bc',\n",
    "        'Only Partial Significant': '#00bc00',\n",
    "        'Neither Significant': '#aaaaaa'\n",
    "    },\n",
    "    ax=axs[1]\n",
    ")\n",
    "\n",
    "# Add thresholds\n",
    "axs[1].axhline(alpha, color='black', linestyle='--', linewidth=1)\n",
    "axs[1].axvline(alpha, color='black', linestyle='--', linewidth=1)\n",
    "\n",
    "# Annotate quadrant counts\n",
    "axs[1].text(0.25, 0.75, f\"Both Significant\\n{quadrant_counts['Both Significant']}\", \n",
    "            ha='center', color='#bc0000', fontsize=12, bbox=dict(facecolor='white', edgecolor='none'))\n",
    "axs[1].text(0.75, 0.75, f\"Only Full Significant\\n{quadrant_counts['Only Full Significant']}\", \n",
    "            ha='center', color='#0000bc', fontsize=12, bbox=dict(facecolor='white', edgecolor='none'))\n",
    "axs[1].text(0.25, 0.25, f\"Only Partial Significant\\n{quadrant_counts['Only Partial Significant']}\", \n",
    "            ha='center', color='#00bc00', fontsize=12, bbox=dict(facecolor='white', edgecolor='none'))\n",
    "axs[1].text(0.75, 0.25, f\"Neither Significant\\n{quadrant_counts['Neither Significant']}\", \n",
    "            ha='center', color='#aaaaaa', fontsize=12, bbox=dict(facecolor='white', edgecolor='none'))\n",
    "\n",
    "# Customize the plot\n",
    "axs[1].set_title('P-Value Comparison with Significance')\n",
    "axs[1].set_xlabel('P-Value Full')\n",
    "axs[1].set_ylabel('P-Value Partial')\n",
    "axs[1].set_xlim(0, 1)\n",
    "axs[1].set_ylim(0, 1)\n",
    "\n",
    "plt.legend(loc='center right', frameon=False, bbox_to_anchor=(1.65, 0.5))\n",
    "# Adjust layout\n",
    "plt.tight_layout()\n",
    "\n",
    "# Show the plots\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pass_chi2_stats[df_pass_chi2_stats['significance'] != 'Neither Significant'].sort_values(by=['mode', 'S'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Only using the full table (for publication)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc('axes', linewidth=0.65)\n",
    "\n",
    "\n",
    "df_pass_chi2_stats['are_diff'] = [True if i < alpha else False for i in df_pass_chi2_stats['p_value_full']]\n",
    "\n",
    "\n",
    "print((df_pass_chi2_stats['are_diff'] == True).sum(), (df_pass_chi2_stats['are_diff'] == False).sum())\n",
    "\n",
    "display(df_pass_chi2_stats)\n",
    "display(df_pass_chi2_stats[df_pass_chi2_stats['are_diff'] == True])\n",
    "\n",
    "\n",
    "\n",
    "# Create a grid of three axes for the plots\n",
    "fig, axs = plt.subplots(1, 3, figsize=(7.5, 2.7), sharey=True)\n",
    "\n",
    "# Define the metrics to plot\n",
    "metrics = ['f1_diff', 'precision_diff', 'recall_diff']\n",
    "titles = ['F1 Score', 'Precision', 'Recall']\n",
    "\n",
    "\n",
    "# Iterate through the metrics and create a boxplot for each\n",
    "for ax, metric, title in zip(axs, metrics, titles):\n",
    "    sns.boxplot(\n",
    "        data=df_pass_chi2_stats, \n",
    "        x=\"are_diff\", \n",
    "        y=metric, \n",
    "        palette=sns.color_palette(['#4E79A7', '#A0CBE8']),\n",
    "        ax=ax\n",
    "    )\n",
    "    sns.despine(top=True, right=True, ax=ax)\n",
    "\n",
    "    ax.set_title(title, fontsize=14, pad=40)\n",
    "    ax.set_xlabel(\"Group\", fontsize=12)\n",
    "    if metric == 'f1_diff':  # Label y-axis only for the first plot\n",
    "        ax.set_ylabel(\"Difference\", fontsize=12)\n",
    "    else:\n",
    "        ax.set_ylabel(\"\")\n",
    "\n",
    "    # Highlight medians\n",
    "    medians = df_pass_chi2_stats.groupby(\"are_diff\")[metric].median()\n",
    "    for i, median in enumerate(medians):\n",
    "        ax.text(i, median, f\"{median:.2f}\", ha='center', va='bottom', fontsize=10)\n",
    "    \n",
    "    pairs=[(False, True)]\n",
    "\n",
    "    annotator = Annotator(ax, pairs, data=df_pass_chi2_stats, x='are_diff', y=metric)\n",
    "    annotator.configure(test='Mann-Whitney', text_format='simple', loc='outside', line_width=1)\n",
    "    annotator.apply_and_annotate()\n",
    "\n",
    "\n",
    "\n",
    "# Adjust layout and save the plot\n",
    "plt.tight_layout()\n",
    "\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/figC.{format}', dpi=DPI)\n",
    "\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pass_chi2_stats[df_pass_chi2_stats['are_diff'] == True].groupby('S').count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pass_chi2_stats[(df_pass_chi2_stats['are_diff'] == True) & (df_pass_chi2_stats['S'] == 0)][['precision_diff',\t'recall_diff',\t'f1_diff']].median()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pass_chi2_stats[(df_pass_chi2_stats['are_diff'] == True) & (df_pass_chi2_stats['S'] == 3)][['precision_diff',\t'recall_diff',\t'f1_diff']].median()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pass_chi2_stats[df_pass_chi2_stats['are_diff'] == True].groupby('column').count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pass_chi2_stats[(df_pass_chi2_stats['are_diff'] == True) & (df_pass_chi2_stats['column'] == 'ganon_norm')][['precision_diff',\t'recall_diff',\t'f1_diff']].median()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pass_chi2_stats[(df_pass_chi2_stats['are_diff'] == True) & (df_pass_chi2_stats['column'] == 'mean_norm')][['precision_diff',\t'recall_diff',\t'f1_diff']].median()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pass_chi2_stats[df_pass_chi2_stats['are_diff'] == True].groupby('mode').count()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we see, using a second pass generally does not significantly change the assignment of TP/FP/FN | TN.\n",
    "\n",
    "In the cases where it changes, the f1-score generally improves, either by increasing the number of TP (by reducing the amount of FN) or reduces the amount of false positives.\n",
    "\n",
    "Therefore, looking at nominal info **choosing pass2 is the best option**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking at the correlation between the read counts\n",
    "\n",
    "Now we are going to check the proportion of reads that are correctly assigned to gold truth species. For that we are goin two show: \n",
    "- A correlation plot with pass0 and pass2 reads. We are also going to calculate the correlation between pass0 / pass2 and the expected counts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_numerical_stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter for pass=0 and pass=2\n",
    "df_pass_0 = df_numerical_stats[df_numerical_stats['pass'] == 0].set_index(['mode', 'profiler'])\n",
    "df_pass_2 = df_numerical_stats[df_numerical_stats['pass'] == 2].set_index(['mode', 'profiler'])\n",
    "\n",
    "# Align the two DataFrames to ensure consistent indexing\n",
    "aligned_df = pd.concat([df_pass_0, df_pass_2], axis=1, keys=['pass0', 'pass2'])\n",
    "\n",
    "# Extract relevant columns and observed counts for pass=0 and pass=2\n",
    "comparison_df = aligned_df[['pass0', 'pass2']].apply(\n",
    "    lambda x: pd.Series({\n",
    "        'passn': 'comparison',  # Mark as comparison for clarity\n",
    "        'mode': x.name[0],\n",
    "        'profiler': x.name[1],\n",
    "        'taxid_counts': x['pass0']['taxid_counts'],  # Keep consistent metadata\n",
    "        'expected_counts': x['pass0']['expected_counts'],\n",
    "        'pass0_counts': x['pass0']['observed_counts'],\n",
    "        'pass2_counts': x['pass2']['observed_counts'],\n",
    "        'pass0_MAE': x['pass0']['MAE_counts'],\n",
    "        'pass0_diff': x['pass0']['diff_counts'],\n",
    "        'pass2_MAE': x['pass2']['MAE_counts'],\n",
    "        'pass2_diff': x['pass2']['diff_counts'],\n",
    "    }),\n",
    "    axis=1\n",
    ")\n",
    "\n",
    "# Reset index for the final DataFrame\n",
    "comparison_df = comparison_df.reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc('axes', linewidth=0.65)\n",
    "\n",
    "# Filter the DataFrame for the specific mode (5 in this case)\n",
    "mode_df = comparison_df[comparison_df['mode'] == 5]\n",
    "\n",
    "# Get unique profilers for layout\n",
    "profilers = mode_df['profiler'].unique()\n",
    "num_profilers = len(profilers)\n",
    "\n",
    "# Set up the figure with GridSpec for custom layout\n",
    "fig, axs = plt.subplots(3, num_profilers, figsize=(2 * num_profilers, 6), sharex=True, sharey=True)\n",
    "\n",
    "# Define plotting functions with regression line\n",
    "def plot_pass0_vs_pass2(ax, data, profiler, i):\n",
    "    x = data['pass0_counts'].values[0]\n",
    "    y = data['pass2_counts'].values[0]\n",
    "    \n",
    "    # Scatter plot\n",
    "    ax.scatter(x, y, alpha=0.7)\n",
    "    \n",
    "    # Plot y=x line\n",
    "    min_val, max_val = min(x.min(), y.min()), max(x.max(), y.max())\n",
    "    ax.plot([0, 350000], [0, 350000], '--', color='#DC267F')\n",
    "    \n",
    "    # Pearson correlation and R^2\n",
    "    slope, intercept, r_value, _, _ = linregress(x, y)\n",
    "    r2 = r_value ** 2\n",
    "    ax.annotate(f\"Slope: {slope:.2f}\\nR$^2$: {r2:.2f}\", xy=(0.05, 0.95),\n",
    "                xycoords='axes fraction', fontsize=12, verticalalignment='top')\n",
    "    if i == 0:\n",
    "        ax.set_ylabel(\"Pass0 vs Pass2\", fontsize=12)\n",
    "    # Set profiler title on top\n",
    "    ax.set_title(profiler, fontsize=14)\n",
    "\n",
    "    sns.despine(top=True, right=True, ax=ax)\n",
    "\n",
    "    \n",
    "def plot_with_regression(ax, x, y, xlabel, i):\n",
    "    # Scatter plot\n",
    "    ax.scatter(x, y, alpha=0.7)\n",
    "\n",
    "    # Regression line\n",
    "    slope, intercept, r_value, _, _ = linregress(x, y)\n",
    "    ax.plot(x, slope * x + intercept, '-', color='#DC267F')\n",
    "\n",
    "    # Pearson correlation and R^2\n",
    "    r2 = r_value ** 2\n",
    "    ax.annotate(f\"Slope: {slope:.2f}\\nR$^2$: {r2:.2f}\", xy=(0.05, 0.95),\n",
    "                xycoords='axes fraction', fontsize=12, verticalalignment='top')\n",
    "    if i == 0:\n",
    "        ax.set_ylabel(xlabel, fontsize=12)\n",
    "    sns.despine(top=True, right=True, ax=ax)\n",
    "\n",
    "# Iterate through profilers and create subplots for each\n",
    "for i, profiler in enumerate(profilers):\n",
    "    profiler_data = mode_df[mode_df['profiler'] == profiler]\n",
    "    x_pass0 = profiler_data['pass0_counts'].values[0]\n",
    "    x_expected = profiler_data['expected_counts'].values[0]\n",
    "    y_pass2 = profiler_data['pass2_counts'].values[0]\n",
    "    \n",
    "    # Pass0 vs Pass2 (Row 1)\n",
    "    ax = fig.add_subplot(axs[0, i])\n",
    "    plot_pass0_vs_pass2(ax, profiler_data, profiler, i)\n",
    "    \n",
    "    # Pass2 vs Expected (Row 2)\n",
    "    ax = fig.add_subplot(axs[1, i])\n",
    "    plot_with_regression(ax, x_expected, y_pass2, \"Pass2 vs Expected\", i)\n",
    "    \n",
    "    # Pass0 vs Expected (Row 3)\n",
    "    ax = fig.add_subplot(axs[2, i])\n",
    "    plot_with_regression(ax, x_expected, x_pass0, \"Pass0 vs Expected\", i)\n",
    "\n",
    "    ax.set_xticks([0, 300000])\n",
    "    ax.set_yticks([0, 300000])\n",
    "\n",
    "\n",
    "\n",
    "# Aesthetic adjustments\n",
    "for ax in axs[-1, :]:\n",
    "    ax.set_xlabel(\"Counts\", fontsize=12)  # Only bottom row has x-labels\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/figD.{format}', dpi=DPI)\n",
    "\n",
    "\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the correlation in number of counts for the ground truth species is almost one. Additionally, if there are differences, they improve the correlation to the expected amount of counts, so it has a positive impact."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking at MAE values\n",
    "\n",
    "We are going to see if the MAE values are equal or higher in pass2 vs pass0. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "comparison_df_exploded = comparison_df.melt(\n",
    "    id_vars=['mode', 'profiler'],\n",
    "    value_vars=['pass0_diff', 'pass2_diff'],\n",
    "    var_name='pass_type',\n",
    "    value_name='diff_value'\n",
    ")\n",
    "\n",
    "# Convert `diff_value` from a list to separate rows\n",
    "comparison_df_exploded = comparison_df_exploded.explode('diff_value')\n",
    "\n",
    "# Ensure `diff_value` is numeric after exploding\n",
    "comparison_df_exploded['diff_value'] = pd.to_numeric(comparison_df_exploded['diff_value'])\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Get unique profilers\n",
    "profilers = comparison_df_exploded['profiler'].unique()\n",
    "num_profilers = len(profilers)\n",
    "\n",
    "# Set up the figure\n",
    "fig, axs = plt.subplots(1, num_profilers, figsize=(6 * num_profilers, 6), sharey=True)\n",
    "\n",
    "# Iterate through profilers and plot the data\n",
    "for i, profiler in enumerate(profilers):\n",
    "    ax = axs[i]\n",
    "    profiler_data = comparison_df_exploded[comparison_df_exploded['profiler'] == profiler]\n",
    "\n",
    "    # Create the plot\n",
    "    sns.boxplot(\n",
    "        data=profiler_data,\n",
    "        x='mode',\n",
    "        y='diff_value',\n",
    "        hue='pass_type',\n",
    "        ax=ax,\n",
    "        palette='viridis'\n",
    "    )\n",
    "    \n",
    "    # Set plot title and labels\n",
    "    ax.set_title(f'Profiler: {profiler}', fontsize=14)\n",
    "    ax.set_xlabel('Mode', fontsize=12)\n",
    "    if i == 0:\n",
    "        ax.set_ylabel('Diff Values', fontsize=12)\n",
    "    else:\n",
    "        ax.set_ylabel('')\n",
    "\n",
    "# Adjust layout and display\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that, again, the differences in MAE are not due to the pass, with the exception of ganon (and mean by extension), where running the data with pass2 improves the error rate. \n",
    "\n",
    "Thus, **both based on nominal and numerical criteria, running with pass2 is the best option**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Which mode is best for the data?\n",
    "\n",
    "We have used several modes to study how it affects the detection of species, and the number of reads it detects. We are going to use the same methods as before to check for the answer.\n",
    "- We are going to check which mode (using S=2 and S=7) retains the best F1 scores.\n",
    "- We are going to check which mode keeps the best correlation / MAE with expected counts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking F1-scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_nominal_stats_sub = df_nominal_stats[(df_nominal_stats['pass'] == 2) & (df_nominal_stats['S'].isin([0, 2, 7, 15])) & (df_nominal_stats['column'].isin([f'{i}_norm' for i in LIST_PROFILERS + ['mean']]))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_nominal_stats_sub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc('axes', linewidth=0.65)\n",
    "\n",
    "melted_df = pd.melt(\n",
    "    df_nominal_stats_sub,\n",
    "    id_vars=['mode', 'S', 'column'],\n",
    "    value_vars=['recall', 'precision', 'f1'],\n",
    "    var_name='metric',\n",
    "    value_name='score'\n",
    ")\n",
    "\n",
    "# Create a colormap for 'mode'\n",
    "norm = Normalize(vmin=melted_df['mode'].min(), vmax=melted_df['mode'].max())\n",
    "cmap = plt.cm.viridis  # Choose a colormap (e.g., 'viridis', 'plasma', 'cividis')\n",
    "\n",
    "# Create a FacetGrid: 6x3 grid (row for each profiler, column for each metric)\n",
    "g = sns.FacetGrid(\n",
    "    melted_df, \n",
    "    col='column', \n",
    "    row='metric', \n",
    "    height=2, \n",
    "    sharey=True, \n",
    "    sharex=True\n",
    ")\n",
    "\n",
    "# Map the lineplot to the grid\n",
    "def lineplot_with_cmap(data, **kwargs):\n",
    "    for S in sorted(data['S'].unique()):\n",
    "        subset = data[data['S'] == S]\n",
    "        plt.plot(subset['mode'], subset['score'], label=f\"Mode {mode}\",\n",
    "                 color=cmap(norm(S)), marker='o')\n",
    "\n",
    "g.map_dataframe(lineplot_with_cmap)\n",
    "\n",
    "# Create a legend for the discrete modes\n",
    "handles = [\n",
    "    plt.Line2D([0], [0], color=cmap(norm(mode)), marker='o', linestyle='', label=f\"S {mode}\")\n",
    "    for mode in sorted(melted_df['S'].unique())\n",
    "]\n",
    "plt.legend(\n",
    "    handles=handles, \n",
    "    title=\"\", \n",
    "    bbox_to_anchor=(1.05, 3), \n",
    "    loc='center left', \n",
    "    frameon=False\n",
    ")\n",
    "\n",
    "# Set x-axis ticks (if you have specific S values)\n",
    "g.set(xticks=df_nominal_stats_sub['mode'].unique())\n",
    "\n",
    "for ax in g.axes.ravel():\n",
    "    ax.set_title('')\n",
    "\n",
    "# Add axis labels and titles\n",
    "for ax, profiler in zip(g.axes[0, :], melted_df['column'].unique()):\n",
    "    ax.set_title(profiler.replace('_norm', ''))\n",
    "\n",
    "for ax, score in zip(g.axes[:, 0], ['recall', 'precision', 'F1-score']):\n",
    "    ax.set_ylabel(score)\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/figF.{format}', dpi=DPI, bbox_inches='tight')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the table of truth, the precision, recall and f1-score are more S-dependent than mode-dependent; and their patterns are more profiler dependent than anything else. Therefore, the mode is not completely relevant for a proper species detection and, if so, since higher mode values lead to a lower precision (more FP), we could consider using a lower mode if we want to minimize the FPs values. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking correlation and MAE values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_numerical_stats_sub = df_numerical_stats[(df_numerical_stats['pass'] == 2) & (df_numerical_stats['mode'].isin([1, 3, 5, 7, 9]))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get unique profilers for layout\n",
    "profilers = df_numerical_stats_sub['profiler'].unique()\n",
    "modes = df_numerical_stats_sub['mode'].unique()\n",
    "\n",
    "\n",
    "# Set up the figure with GridSpec for custom layout\n",
    "fig, axs = plt.subplots(len(profilers), len(modes), figsize=(2 * len(modes), 2 * len(profilers)), sharex=True, sharey=True)\n",
    "\n",
    "# Define plotting functions with regression line\n",
    "def plot_regression(ax, data, profiler, i):\n",
    "    x = data['expected_counts'].values[0]\n",
    "    y = data['observed_counts'].values[0]\n",
    "    \n",
    "    # Scatter plot\n",
    "    ax.scatter(x, y, alpha=0.7)\n",
    "    \n",
    "    # Plot y=x line\n",
    "    min_val, max_val = min(x.min(), y.min()), max(x.max(), y.max())\n",
    "    ax.plot([0, 350000], [0, 350000], 'r--')\n",
    "    \n",
    "    # Pearson correlation and R^2\n",
    "    slope, intercept, r_value, _, _ = linregress(x, y)\n",
    "    r2 = r_value ** 2\n",
    "    ax.annotate(f\"Slope: {slope:.2f}\\nR2: {r2:.2f}\", xy=(0.05, 0.95),\n",
    "                xycoords='axes fraction', fontsize=10, verticalalignment='top')\n",
    "    \n",
    "\n",
    "# Iterate through profilers and create subplots for each\n",
    "for i, profiler in enumerate(profilers):\n",
    "    for j, mode in enumerate(modes):\n",
    "        profiler_data = df_numerical_stats_sub[(df_numerical_stats_sub['profiler'] == profiler) & (df_numerical_stats_sub['mode'] == mode)]\n",
    "        \n",
    "        # Pass0 vs Pass2 (Row 1)\n",
    "        ax = fig.add_subplot(axs[i, j])\n",
    "        plot_regression(ax, profiler_data, profiler, i)\n",
    "\n",
    "        if j == 0:\n",
    "            ax.set_ylabel(profiler, fontsize=12)\n",
    "            # Set profiler title on top\n",
    "        if i == 0:\n",
    "            ax.set_title(f'Mode {mode}', fontsize=14)\n",
    "\n",
    "# Aesthetic adjustments\n",
    "for ax in axs[-1, :]:\n",
    "    ax.set_xlabel(\"Counts\", fontsize=12)  # Only bottom row has x-labels\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare the data for plotting\n",
    "expanded_data = []\n",
    "for _, row in df_numerical_stats_sub.iterrows():\n",
    "    profiler = row['profiler']\n",
    "    mode = row['mode']\n",
    "    for value in row['diff_counts']:\n",
    "        expanded_data.append({'profiler': profiler, 'mode': mode, 'diff_counts': value})\n",
    "\n",
    "# Convert to a new DataFrame\n",
    "plot_data = pd.DataFrame(expanded_data)\n",
    "\n",
    "# Initialize the grid of plots with seaborn\n",
    "g = sns.FacetGrid(plot_data, col=\"profiler\", height=4, aspect=1)\n",
    "\n",
    "# Plot the data on each facet\n",
    "g.map_dataframe(\n",
    "    sns.boxplot,\n",
    "    x=\"mode\",\n",
    "    y=\"diff_counts\",\n",
    "    palette=\"viridis\"\n",
    ")\n",
    "\n",
    "# Adjust the legend and titles\n",
    "g.set_axis_labels(\"Mode\", \"Difference Counts\")\n",
    "g.set_titles(col_template=\"{col_name} Profiler\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.ticker import MaxNLocator\n",
    "\n",
    "plt.rc('axes', linewidth=0.65)\n",
    "\n",
    "# Initialize the grid of plots with seaborn\n",
    "g = sns.FacetGrid(df_numerical_stats_sub, col=\"profiler\", height=2.3, aspect=0.85, sharex=False, sharey=False)\n",
    "\n",
    "# Plot the data on each facet\n",
    "g.map_dataframe(\n",
    "    sns.scatterplot,\n",
    "    x=\"corr_counts\",\n",
    "    y=\"R2_counts\",\n",
    "    hue=\"mode\",\n",
    "    palette=\"viridis\"\n",
    ")\n",
    "\n",
    "# Adjust the legend and titles\n",
    "g.add_legend(title=\"Mode\", bbox_to_anchor=(1.05, 0.5))\n",
    "\n",
    "g.set_axis_labels(\"Correlation Counts\", \"R² Counts\")\n",
    "g.set_titles(col_template=\"{col_name}\")\n",
    "\n",
    "for ax in g.axes.flat:\n",
    "    ax.xaxis.set_major_locator(MaxNLocator(nbins=2)) \n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/figGA.{format}', dpi=DPI, bbox_inches='tight')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.FacetGrid(df_numerical_stats_sub, col=\"profiler\", height=2.3, aspect=0.85, sharex=False, sharey=False)\n",
    "\n",
    "# Plot the data on each facet\n",
    "g.map_dataframe(\n",
    "    sns.scatterplot,\n",
    "    x=\"MAE_counts\",\n",
    "    y=\"MAED_counts\",\n",
    "    hue=\"mode\",\n",
    "    palette=\"viridis\"\n",
    ")\n",
    "\n",
    "# Adjust the legend and titles\n",
    "g.add_legend(title=\"Mode\", bbox_to_anchor=(1.05, 0.5))\n",
    "\n",
    "g.set_axis_labels(\"MAE\", \"MAED\")\n",
    "g.set_titles(col_template=\"{col_name}\")\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/figGB.{format}', dpi=DPI, bbox_inches='tight')\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.FacetGrid(df_numerical_stats_sub, col=\"profiler\", height=2.3, aspect=0.85, sharex=False, sharey=False)\n",
    "\n",
    "# Plot the data on each facet\n",
    "g.map_dataframe(\n",
    "    sns.scatterplot,\n",
    "    x=\"MRE_counts\",\n",
    "    y=\"MRED_counts\",\n",
    "    hue=\"mode\",\n",
    "    palette=\"viridis\"\n",
    ")\n",
    "\n",
    "# Adjust the legend and titles\n",
    "g.add_legend(title=\"Mode\", bbox_to_anchor=(1.05, 0.5))\n",
    "\n",
    "g.set_axis_labels(\"MRE\", \"MRED\")\n",
    "g.set_titles(col_template=\"{col_name}\")\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/figGC.{format}', dpi=DPI, bbox_inches='tight')\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the DataFrame\n",
    "df = df_numerical_stats[df_numerical_stats['pass'] == 2].copy()\n",
    "\n",
    "\n",
    "\n",
    "# Calculate the differences from the mean for each mode\n",
    "df_mean = df[df['profiler'] == 'mean'][['mode', 'MRED_counts']].rename(columns={'MRED_counts': 'mean_MRED'})\n",
    "df = df.merge(df_mean, on='mode')\n",
    "df['difference'] = df['MRED_counts'] - df['mean_MRED']\n",
    "\n",
    "\n",
    "# Filter out the mean profiler itself\n",
    "df_filtered = df[df['profiler'] != 'mean'].sort_values('profiler')\n",
    "\n",
    "\n",
    "# Perform Mann-Whitney U tests\n",
    "profiler_groups = df_filtered['profiler'].unique()\n",
    "p_values = []\n",
    "results = []\n",
    "\n",
    "\n",
    "for profiler in profiler_groups:\n",
    "    group_data = df_filtered[df_filtered['profiler'] == profiler]\n",
    "    mean_data = group_data['MRED_counts'].values\n",
    "    prof_data = group_data['mean_MRED'].values\n",
    "    stat, p = wilcoxon(mean_data, prof_data, alternative='greater')\n",
    "    p_values.append(p)\n",
    "    results.append({'profiler': profiler, 'statistic': stat, 'p_value': p})\n",
    "\n",
    "# Correct p-values for multiple testing\n",
    "p_values_corrected = multipletests(p_values, method='fdr_bh')[1]\n",
    "\n",
    "for i, result in enumerate(results):\n",
    "    result['p_value_corrected'] = p_values_corrected[i]\n",
    "\n",
    "\n",
    "# Convert results to DataFrame\n",
    "results_df = pd.DataFrame(results)\n",
    "\n",
    "display(results_df)\n",
    "\n",
    "# Plot the differences\n",
    "fig, ax = plt.subplots(1, 1, figsize=(5, 3))\n",
    "\n",
    "sns.boxplot(data=df_filtered, x='profiler', y='difference', palette=custom_palette, ax=ax)\n",
    "# sns.stripplot(data=df_filtered, x='profiler', y='difference', color='black', alpha=0.5, jitter=True)\n",
    "medians = df_filtered.groupby(\"profiler\")['difference'].median()\n",
    "for i, median in enumerate(medians):\n",
    "        ax.text(i, median , f\"{median:.2f}\", ha='center', va='bottom', fontsize=8)\n",
    "\n",
    "\n",
    "plt.axhline(0, linestyle='--', color='#848484', linewidth=0.5)\n",
    "plt.ylabel('MRED difference from mean')\n",
    "plt.xlabel('')\n",
    "\n",
    "sns.despine(top=True, right=True)\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/figJ.{format}', dpi=DPI, bbox_inches='tight')\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the DataFrame\n",
    "df = df_numerical_stats[df_numerical_stats['pass'] == 2].copy()\n",
    "\n",
    "\n",
    "\n",
    "# Calculate the differences from the mean for each mode\n",
    "df_mean = df[df['profiler'] == 'mean'][['mode', 'MAED_counts']].rename(columns={'MAED_counts': 'mean_MAED'})\n",
    "df = df.merge(df_mean, on='mode')\n",
    "df['difference'] = df['MAED_counts'] - df['mean_MAED']\n",
    "\n",
    "\n",
    "# Filter out the mean profiler itself\n",
    "df_filtered = df[df['profiler'] != 'mean'].sort_values('profiler')\n",
    "\n",
    "\n",
    "# Perform Mann-Whitney U tests\n",
    "profiler_groups = df_filtered['profiler'].unique()\n",
    "p_values = []\n",
    "results = []\n",
    "\n",
    "\n",
    "for profiler in profiler_groups:\n",
    "    group_data = df_filtered[df_filtered['profiler'] == profiler]\n",
    "    mean_data = group_data['MAED_counts'].values\n",
    "    prof_data = group_data['mean_MAED'].values\n",
    "    stat, p = wilcoxon(mean_data, prof_data, alternative='greater')\n",
    "    p_values.append(p)\n",
    "    results.append({'profiler': profiler, 'statistic': stat, 'p_value': p})\n",
    "\n",
    "# Correct p-values for multiple testing\n",
    "p_values_corrected = multipletests(p_values, method='fdr_bh')[1]\n",
    "\n",
    "for i, result in enumerate(results):\n",
    "    result['p_value_corrected'] = p_values_corrected[i]\n",
    "\n",
    "\n",
    "# Convert results to DataFrame\n",
    "results_df = pd.DataFrame(results)\n",
    "\n",
    "display(results_df)\n",
    "\n",
    "# Plot the differences\n",
    "fig, ax = plt.subplots(1, 1, figsize=(5, 3))\n",
    "\n",
    "sns.boxplot(data=df_filtered, x='profiler', y='difference', palette=custom_palette, ax=ax)\n",
    "# sns.stripplot(data=df_filtered, x='profiler', y='difference', color='black', alpha=0.5, jitter=True)\n",
    "medians = df_filtered.groupby(\"profiler\")['difference'].median()\n",
    "for i, median in enumerate(medians):\n",
    "        ax.text(i, median , f\"{median:.2f}\", ha='center', va='bottom', fontsize=8)\n",
    "\n",
    "\n",
    "\n",
    "plt.axhline(0, linestyle='--', color='#848484', linewidth=0.5)\n",
    "plt.ylabel('MAED difference from mean')\n",
    "plt.xlabel('')\n",
    "\n",
    "sns.despine(top=True, right=True)\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/figJ.{format}', dpi=DPI, bbox_inches='tight')\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_numerical_stats_sub[df_numerical_stats_sub['mode'].isin([3, 7])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_prof_couns(df_combo):\n",
    "    # Assuming your data is in a DataFrame called `df`\n",
    "    # Expand the diff_counts column into individual rows for plotting\n",
    "    df_expanded = df_combo.explode('diff_counts')\n",
    "    df_expanded['diff_counts'] = pd.to_numeric(df_expanded['diff_counts'])\n",
    "\n",
    "    # Create the plot\n",
    "    plt.figure(figsize=(10, 2.5))\n",
    "    \n",
    "\n",
    "    for profiler, countsdiff in zip(df_combo['profiler'], df_combo['diff_counts']):\n",
    "        me = np.median(countsdiff)\n",
    "        med = np.std(countsdiff)\n",
    "        plt.scatter(me, profiler, color='#DC267F', label='_nolegend_', s=100, zorder=3, marker = '|')\n",
    "        plt.plot([me-med, me+med], [profiler, profiler], color='#DC267F', label='_nolegend_',)\n",
    "\n",
    "    sns.stripplot(\n",
    "        data=df_expanded,\n",
    "        x='diff_counts',\n",
    "        y='profiler',\n",
    "        jitter=True,  # Adds jitter for better visibility of points\n",
    "        size=5,  # Adjust point size\n",
    "        alpha=0.7,  # Slight transparency\n",
    "        color=\"#648FFF\"\n",
    "    )\n",
    "\n",
    "    # Add a vertical line at x=0\n",
    "    plt.axvline(0, color='#848484', linestyle='--', linewidth=1)\n",
    "\n",
    "    sns.despine(top=True, right=True)\n",
    "\n",
    "    # Add labels and title\n",
    "    # plt.title('Diff Counts Across Profilers', fontsize=14)\n",
    "    plt.xlabel('PRE', fontsize=12)\n",
    "    plt.ylabel('', fontsize=12)\n",
    "    plt.grid(True, axis='x', linestyle='--', alpha=0.6)\n",
    "    plt.tight_layout()\n",
    "\n",
    "    for format in ['png', 'tiff']: \n",
    "        plt.savefig(f'{RESULTS_DIR}/figures/paper/figH{mode}.{format}', dpi=DPI, bbox_inches='tight')\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for mode in [3,5,7]:\n",
    "    print(mode)\n",
    "    df_combo = df_numerical_stats[(df_numerical_stats['pass'] == 2) & \\\n",
    "                              (df_numerical_stats['mode'] == mode)] # We choose a large number because S in not relevant here (but with small S we may select few datasets)\n",
    "\n",
    "    plot_prof_couns(df_combo)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on these results, we see that a higher mode increases the correlation (slope) and R2; and decreases the MAE (although it increases the MAED). However, again, we see that this effect is very profiler-dependent. It is interesting to note, however, that there are some species that have a very bad detection rate throughout the profiler/mode/S values.\n",
    "\n",
    "Based on that, considering that higher modes (1) lower the MAE / increase the correlation but (2) decrease the precision; but considering that the effect on (1) is more pronounced that in (2), we are going to choose a high mode value, but not extreme. For instance, `mode = 7`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Which species have the worst assignment?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "species_to_kingdom = {\n",
    "    # Bacteria\n",
    "    \"Cutibacterium acnes\": \"Bacteria\",\n",
    "    \"Lactobacillus acidophilus\": \"Bacteria\",\n",
    "    \"Bifidobacterium bifidum\": \"Bacteria\",\n",
    "    \"Akkermansia muciniphila\": \"Bacteria\",\n",
    "    \"Blautia coccoides\": \"Bacteria\",\n",
    "    \"Blautia luti\": \"Bacteria\",\n",
    "    \"Bacteroides ovatus\": \"Bacteria\",\n",
    "    \"Bacteroides intestinalis\": \"Bacteria\",\n",
    "    \"Bacteroides fragilis\": \"Bacteria\",\n",
    "    \"Escherichia coli\": \"Bacteria\",\n",
    "    \"Dietzia lutea\": \"Bacteria\",\n",
    "    \"Ruthenibacterium lactatiformans\": \"Bacteria\",\n",
    "    \"Faecalibacterium prausnitzii\": \"Bacteria\",\n",
    "    \"Parabacteroides distasonis\": \"Bacteria\",\n",
    "    \"Parabacteroides merdae\": \"Bacteria\",\n",
    "    \"Fusicatenibacter saccharivorans\": \"Bacteria\",\n",
    "    \"Erysipelatoclostridium ramosum\": \"Bacteria\",\n",
    "    \"Streptococcus salivarius\": \"Bacteria\",\n",
    "    \"Hungatella hathewayi\": \"Bacteria\",\n",
    "    \"Eisenbergiella porci\": \"Bacteria\",\n",
    "    \"Butyricimonas faecalis\": \"Bacteria\",\n",
    "    \"Alistipes indistinctus\": \"Bacteria\",\n",
    "    \"Alistipes finegoldii\": \"Bacteria\",\n",
    "    \"Eubacterium callanderi\": \"Bacteria\",\n",
    "    \"Acidaminococcus intestini\": \"Bacteria\",\n",
    "    \n",
    "    # Fungi\n",
    "    \"Aspergillus chevalieri\": \"Fungi\",\n",
    "    \"Aspergillus flavus\": \"Fungi\",\n",
    "    \"Saccharomyces cerevisiae\": \"Fungi\",\n",
    "    \"Saccharomyces kudriavzevii\": \"Fungi\",\n",
    "    \"Saccharomyces mikatae\": \"Fungi\",\n",
    "    \"Candida albicans\": \"Fungi\",\n",
    "    \"Candida dubliniensis\": \"Fungi\",\n",
    "    \"Candida orthopsilosis\": \"Fungi\",\n",
    "    \"Malassezia restricta\": \"Fungi\",\n",
    "    \"Alternaria dauci\": \"Fungi\",\n",
    "    \"Kazachstania africana\": \"Fungi\",\n",
    "    \"Penicillium digitatum\": \"Fungi\",\n",
    "    \"Pichia kudriavzevii\": \"Fungi\",\n",
    "    \"Trichoderma asperellum\": \"Fungi\",\n",
    "    \"Akanthomyces muscarius\": \"Fungi\",\n",
    "    \"Fusarium falciforme\": \"Fungi\",\n",
    "    \"Eremothecium sinecaudum\": \"Fungi\",\n",
    "    \"Cryptococcus decagattii\": \"Fungi\",\n",
    "    \"Kwoniella shandongensis\": \"Fungi\",\n",
    "    \"Puccinia triticina\": \"Fungi\",\n",
    "    \n",
    "    # Virus\n",
    "    \"Tobacco mosaic virus\": \"Virus\",\n",
    "    \"Rotavirus A\": \"Virus\",\n",
    "    \"Rotavirus B\": \"Virus\",\n",
    "    \"Rotavirus C\": \"Virus\",\n",
    "    \"Bacteriophage P2\": \"Virus\",\n",
    "    \"Escherichia phage T4\": \"Virus\",\n",
    "    \"Human immunodeficiency virus 1\": \"Virus\",\n",
    "    \"Human adenovirus 7\": \"Virus\",\n",
    "    \"Hepatitis C virus\": \"Virus\",\n",
    "    \"Bovine alphaherpesvirus 2\": \"Virus\",\n",
    "    \"Human herpesvirus 4 type 2\": \"Virus\",\n",
    "    \"Mimivirus terra2\": \"Virus\",\n",
    "    \"Dengue virus\": \"Virus\",\n",
    "    \"Norovirus GI\": \"Virus\",\n",
    "    \"Zaire ebolavirus\": \"Virus\"\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "mode = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_combo = df_numerical_stats[(df_numerical_stats['pass'] == 2) & \\\n",
    "                              (df_numerical_stats['mode'] == mode)].set_index('profiler')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_combo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "# RUN FOR MEAN IN MODE 1 3 5 7 9\n",
    "# RUN FOR MODE 5 AND THE REST OF PROFILERS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prof = 'mean'\n",
    "\n",
    "df_diffab = pd.DataFrame({'species': df_combo.loc[prof, 'species_counts'],  \n",
    "                          'expected_counts': df_combo.loc[prof, 'expected_counts'], \\\n",
    "                          'observed_counts': df_combo.loc[prof, 'observed_counts'], \\\n",
    "                          'diff_abundance': df_combo.loc[prof, 'diff_counts']})\n",
    "df_diffab['kingdom'] = [species_to_kingdom[i] for i in df_diffab['species']]\n",
    "df_diffab['isfamily'] = False\n",
    "\n",
    "family_species = ['Blautia coccoides', 'Blautia luti', 'Bacteroides ovatus', 'Bacteroides intestinalis','Bacteroides fragilis', \n",
    "                  'Parabacteroides distasonis', 'Parabacteroides merdae', 'Alistipes indistinctus', 'Alistipes finegoldii', 'Aspergillus chevalieri', \n",
    "                  'Aspergillus flavus', 'Saccharomyces cerevisiae', 'Saccharomyces kudriavzevii', 'Saccharomyces mikatae', 'Candida albicans', 'Candida dubliniensis', \n",
    "                  'Candida orthopsilosis', 'Rotavirus A', 'Rotavirus B', 'Rotavirus C']\n",
    "df_diffab.loc[df_diffab['species'].isin(family_species), 'isfamily'] = True\n",
    "df_diffab.sort_values(by=['expected_counts', 'diff_abundance'], ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a grid of three axes for the plots\n",
    "fig, axs = plt.subplots(1, 3, figsize=(9, 3.5), sharey=True)\n",
    "\n",
    "pairs=[(True, False)]\n",
    "\n",
    "\n",
    "ax1 = sns.boxplot(df_diffab, x='isfamily', y='diff_abundance', ax=axs[0], palette=sns.color_palette(['#4E79A7', '#A0CBE8']))\n",
    "medians = df_diffab.groupby(\"isfamily\")['diff_abundance'].median()\n",
    "for i, median in enumerate(medians):\n",
    "        ax1.text(i, median - 8, f\"{median:.2f}\", ha='center', va='bottom', fontsize=8)\n",
    "\n",
    "axs[0].set_ylabel(\"PRE\", fontsize=12)\n",
    "annotator = Annotator(ax1, pairs, data=df_diffab, x='isfamily', y='diff_abundance')\n",
    "annotator.configure(test='Mann-Whitney', text_format='simple', loc='outside', line_width=1)\n",
    "annotator.apply_and_annotate()\n",
    "axs[0].set_xlabel(\"Species belongs to\\na shared genus\", fontsize=12)\n",
    "\n",
    "\n",
    "\n",
    "ax2 = sns.boxplot(df_diffab, x='kingdom', y='diff_abundance', ax=axs[1], palette=sns.color_palette(['#4E79A7', '#77a2c8', '#A0CBE8']))\n",
    "pairs=[(\"Virus\", 'Bacteria'), ('Fungi', 'Bacteria'), ('Fungi', 'Virus')]\n",
    "annotator = Annotator(ax2, pairs, data=df_diffab, x='kingdom', y='diff_abundance')\n",
    "annotator.configure(test='Mann-Whitney', text_format='simple', loc='outside', line_width=1)\n",
    "annotator.apply_and_annotate()\n",
    "medians = df_diffab.groupby(\"kingdom\")['diff_abundance'].median()\n",
    "for i, median in enumerate(medians):\n",
    "        ax2.text(i, median - 8, f\"{median:.2f}\", ha='center', va='bottom', fontsize=8)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "df_bac = df_diffab[df_diffab['kingdom'] == 'Bacteria']\n",
    "\n",
    "ax3 = sns.boxplot(df_bac, x='isfamily', y='diff_abundance', ax=axs[2], palette=sns.color_palette(['#4E79A7', '#A0CBE8']))\n",
    "pairs=[(True, False)]\n",
    "annotator = Annotator(ax3, pairs, data=df_bac, x='isfamily', y='diff_abundance')\n",
    "annotator.configure(test='Mann-Whitney', text_format='simple', loc='outside', line_width=1)\n",
    "annotator.apply_and_annotate()\n",
    "axs[2].set_xlabel(\"Bacterial species belongs\\nto a shared genus\", fontsize=12)\n",
    "medians = df_bac.groupby(\"isfamily\")['diff_abundance'].median()\n",
    "for i, median in enumerate(medians):\n",
    "        ax3.text(i, median - 8, f\"{median:.2f}\", ha='center', va='bottom', fontsize=8)\n",
    "\n",
    "\n",
    "for ax in axs:\n",
    "    sns.despine(top=True, right=True, ax=ax)\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "for format in ['png', 'tiff']: \n",
    "    plt.savefig(f'{RESULTS_DIR}/figures/paper/figI-{mode}-{prof}.{format}', dpi=DPI, bbox_inches='tight')\n",
    "\n",
    "plt.show()\n"
   ]
  }
 ],
 "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.13.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
