{
 "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 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\n",
    "from statsmodels.stats.multitest import multipletests\n",
    "import json\n",
    "from kneed import KneeLocator\n"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from list_vars import LIST_PROFILERS, DIR_FIGURES, RESULTS_DIR, POOLS, CONTROLS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Biological sample analysis\n",
    "\n",
    "In this notebook we are going to do an analysis on the *biological* samples (POOL samples + controls).\n",
    "\n",
    "There are two main variables that we are going to consider:\n",
    "- The importance of including the biological control samples to ensure that false positives are not considered.\n",
    "- The importance of normalizing the reads considering the biooogical samples.\n",
    "\n",
    "These two concepts are intertwined, so what we are going to do is the following:\n",
    "- Load all pool + control tables.\n",
    "- Get the species that are selected using the flags. With that we are going to generate a cut-off table with the species. The way to merge the table is \"outer\", that is, we are going to include any species that appears in any sample. We can later discard them.\n",
    "- We are going to discard the species that have less than X times more expression in control samples than in pools. X is determined dynamically using the kneed method.\n",
    "- Then we are going to compare which species have been discarded using the | normalization (normalizing controls and pools separately) and the + normalization (normalizng controls and pools jointly)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading the tables\n",
    "def process_samples(pass_num, mode, S, NORM, samples, verbose, min_sample_flag='dynamic'):\n",
    "    # Initialize empty DataFrames for counts and flags\n",
    "    joined_counts = pd.DataFrame()\n",
    "    taxid_list = []\n",
    "\n",
    "    for sample in samples:\n",
    "        # Define file paths for counts and flags\n",
    "        counts_file = f'{RESULTS_DIR}/summary/{sample}_pass{pass_num}_mode{mode}_taxgenus_S{S}_{NORM}.diversity.tsv'\n",
    "        flags_file =  f'{RESULTS_DIR}/summary/{sample}_pass{pass_num}_mode{mode}_taxgenus_S{S}_{NORM}.flags.tsv'\n",
    "\n",
    "        # Load the data\n",
    "        df_counts = pd.read_csv(counts_file, sep='\\t').set_index('taxonomy_id')[['name', 'lineage', 'mean_norm']]\n",
    "        df_flags = pd.read_csv(flags_file, sep='\\t').set_index('taxonomy_id')[['name', 'lineage', 'mean_norm']]\n",
    "\n",
    "        # Ensure name and lineage columns are retained correctly\n",
    "        if joined_counts.empty:\n",
    "            joined_counts = df_counts.rename(columns={'mean_norm': sample})\n",
    "        else:\n",
    "            df_counts = df_counts.rename(columns={'mean_norm': sample})\n",
    "            joined_counts = joined_counts.join(df_counts, how='outer', rsuffix=f'_{sample}')\n",
    "\n",
    "        # Add tax_ids where mean_norm is False in flags dataframe\n",
    "        taxid_list += df_flags[df_flags['mean_norm'] == False].index.tolist()\n",
    "\n",
    "    # Consolidate name and lineage columns to avoid suffix issues\n",
    "    if not joined_counts.empty:\n",
    "        joined_counts['name'] = joined_counts.filter(like='name').bfill(axis=1).iloc[:, 0]\n",
    "        joined_counts['lineage'] = joined_counts.filter(like='lineage').bfill(axis=1).iloc[:, 0]\n",
    "        joined_counts = joined_counts.drop(columns=joined_counts.filter(like='name_').columns)\n",
    "        joined_counts = joined_counts.drop(columns=joined_counts.filter(like='lineage_').columns)\n",
    "\n",
    "    # Rename pool columns based on their ranges\n",
    "    rename_mapping = {\n",
    "        'POOL1': 'RR1', 'POOL2': 'RR2', 'POOL3': 'RR3', 'POOL4': 'RR4',\n",
    "        'POOL5': 'SP1', 'POOL6': 'SP2', 'POOL7': 'SP3', 'POOL8': 'SP4',\n",
    "        'POOL9': 'HC1', 'POOL10': 'HC2', 'POOL11': 'HC3', 'POOL12': 'HC4'\n",
    "    }\n",
    "    joined_counts = joined_counts.rename(columns=rename_mapping)\n",
    "\n",
    "    # Create a \"cut\" DataFrame containing only the tax_ids in the taxid_list\n",
    "    taxidvalues, samplecounts = np.unique(taxid_list, return_counts=True)\n",
    "    joined_counts['n_samples_flag'] = 0\n",
    "    joined_counts.loc[taxidvalues, 'n_samples_flag'] = samplecounts\n",
    "\n",
    "\n",
    "    n_samples, counts_ntaxids = np.unique(joined_counts['n_samples_flag'].values, return_counts=True)\n",
    "            \n",
    "    if min_sample_flag == 'dynamic':\n",
    "        kneedle = KneeLocator(n_samples, np.cumsum(counts_ntaxids), curve='concave', direction='increasing', S=0)\n",
    "        min_sample_flag = kneedle.knee\n",
    "    \n",
    "\n",
    "    joined_counts['selected_flag'] = joined_counts['n_samples_flag'] >= min_sample_flag\n",
    "    cut_df = joined_counts[joined_counts['selected_flag'] == True]\n",
    "\n",
    "    if verbose: \n",
    "        print('TaxIDs species count:', n_samples, counts_ntaxids)\n",
    "        print(f'Flag threshold: {min_sample_flag} | Number of species: {len(joined_counts)} | Species selected: {len(cut_df)} ({100 * len(cut_df) / len(joined_counts):.2f}%)')\n",
    "\n",
    "    \n",
    "\n",
    "    # Reset index and sort by mean counts (descending order)\n",
    "    joined_counts = joined_counts.reset_index().sort_values(by=list(rename_mapping.values()), ascending=False)\n",
    "    cut_df = cut_df.reset_index().sort_values(by=list(rename_mapping.values()), ascending=False)\n",
    "\n",
    "    return joined_counts, taxid_list, cut_df\n",
    "# Example usage\n",
    "samples = [\n",
    "    'POOL1', 'POOL2', 'POOL3', 'POOL4', 'POOL5', 'POOL6',\n",
    "    'POOL7', 'POOL8', 'POOL9', 'POOL10', 'POOL11', 'POOL12',\n",
    "    'ACIDOLA', 'BLACTIS'\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def filter_by_nan_percentage(df, per_cutoff=0.35):\n",
    "    # Identify sample columns (excluding taxonomy_id, name, lineage, and controls)\n",
    "    control_cols = ['ACIDOLA', 'BLACTIS']\n",
    "    sample_cols = [col for col in df.columns if col not in ['taxonomy_id', 'name', 'lineage'] + control_cols]\n",
    "\n",
    "    # Calculate the percentage of NaNs in sample columns\n",
    "    nan_percentage = df[sample_cols].isna().mean(axis=1)\n",
    "\n",
    "    # Retain species with less than 35% NaNs\n",
    "    filtered_df = df[nan_percentage < per_cutoff]\n",
    "\n",
    "    return filtered_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_retained_discarded(df, threshold, verbose):\n",
    "        df = df.copy()\n",
    "        \n",
    "        # Separate ACIDOLA and BLACTIS columns\n",
    "        control_cols = ['ACIDOLA', 'BLACTIS']\n",
    "        sample_cols = [col for col in df.columns if col not in ['taxonomy_id', 'name', 'lineage'] + control_cols]\n",
    "\n",
    "        # Calculate mean across samples\n",
    "        df['median_across_samples'] = df[sample_cols].median(axis=1, skipna=True)\n",
    "\n",
    "        # Calculate max of ACIDOLA and BLACTIS\n",
    "        df['max_control'] = df[control_cols].max(axis=1, skipna=True).fillna(0)\n",
    "\n",
    "        if threshold == 'dynamic':\n",
    "            list_len_discarded = []\n",
    "\n",
    "            for threshold in range(1, 500): # In theory the dataframe len is not related but it is just a number to add, which should be bigger the bigger the dataframe \n",
    "                discarded = df[~((df['median_across_samples'] > (df['max_control'] * threshold)))]\n",
    "                list_len_discarded.append(len(discarded) / len(df))\n",
    "\n",
    "            kneedle = KneeLocator(np.arange(1,500), list_len_discarded, curve='concave', direction='increasing', S=0)\n",
    "            threshold = kneedle.knee\n",
    "\n",
    "\n",
    "        # Define retention logic\n",
    "        retained = df[(df['median_across_samples'] > (df['max_control'] * threshold)) | (df['max_control'].isna())]\n",
    "        discarded = df[~((df['median_across_samples'] > (df['max_control'] * threshold)) | (df['max_control'].isna()))]\n",
    "\n",
    "        if verbose: \n",
    "                print(f'Threshold: {threshold} | Number of species: {len(df)} | Species discarded: {len(discarded)} ({100 * len(discarded) / len(df):.2f}%)')\n",
    "\n",
    "        return retained, discarded\n",
    "\n",
    "def filter_species_ids(joined_counts_norm_plus, joined_counts_norm_pipe, threshold='dynamic', verbose=True):\n",
    "    retained_norm_plus, discarded_norm_plus = calculate_retained_discarded(joined_counts_norm_plus, threshold, verbose)\n",
    "    retained_norm_pipe, discarded_norm_pipe = calculate_retained_discarded(joined_counts_norm_pipe, threshold, verbose)\n",
    "\n",
    "    # Extract taxonomy IDs\n",
    "    retained_ids_norm_plus = retained_norm_plus['taxonomy_id'].tolist()\n",
    "    discarded_ids_norm_plus = discarded_norm_plus['taxonomy_id'].tolist()\n",
    "    retained_ids_norm_pipe = retained_norm_pipe['taxonomy_id'].tolist()\n",
    "    discarded_ids_norm_pipe = discarded_norm_pipe['taxonomy_id'].tolist()\n",
    "\n",
    "\n",
    "\n",
    "    discarded_common = np.intersect1d(discarded_ids_norm_plus, discarded_ids_norm_pipe).tolist()\n",
    "    discarded_exclusive_norm_plus = [i for i in discarded_ids_norm_plus if i not in discarded_ids_norm_pipe]\n",
    "    discarded_exclusive_norm_pipe = [i for i in discarded_ids_norm_pipe if i not in discarded_ids_norm_plus]    \n",
    "\n",
    "    return {\n",
    "        'discarded_common': discarded_common,\n",
    "        'discarded_exclusive_norm_plus': discarded_exclusive_norm_plus,\n",
    "        'discarded_exclusive_norm_pipe': discarded_exclusive_norm_pipe,\n",
    "        'retained_norm_plus': retained_ids_norm_plus,\n",
    "        'retained_norm_pipe': retained_ids_norm_pipe\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def differential_abundance_analysis(df, condition_cols, reference_cols):\n",
    "    \"\"\"\n",
    "    Perform differential abundance analysis between condition and reference groups.\n",
    "\n",
    "    Parameters:\n",
    "        df (pd.DataFrame): Input dataframe containing species counts and metadata.\n",
    "        condition_cols (list): Column names for the condition group.\n",
    "        reference_cols (list): Column names for the reference group.\n",
    "        output_file (str, optional): Path to save results to an Excel file. Default is None.\n",
    "        sheet_name (str): Sheet name for Excel output. Default is 'Results'.\n",
    "\n",
    "    Returns:\n",
    "        pd.DataFrame: Dataframe containing p-values, log2 fold changes, and sorted results.\n",
    "    \"\"\"\n",
    "    list_pvals_mannwhitney = []\n",
    "    L2FC = []\n",
    "\n",
    "    for row in range(len(df)):\n",
    "        # Extract condition and reference values\n",
    "        condition_vals = df.iloc[row][condition_cols].astype(float).dropna().values\n",
    "        reference_vals = df.iloc[row][reference_cols].astype(float).dropna().values\n",
    "\n",
    "        condition_vals, reference_vals = condition_vals, reference_vals\n",
    "\n",
    "        # Mann-Whitney U test\n",
    "        res_mw = mannwhitneyu(condition_vals, reference_vals, alternative='two-sided')\n",
    "        list_pvals_mannwhitney.append(res_mw.pvalue)\n",
    "\n",
    "        # Log2 fold change\n",
    "        L2FC.append(np.log2(np.nanmedian(condition_vals) / np.nanmedian(reference_vals)))\n",
    "\n",
    "    # Compile results\n",
    "    df_pval = df.copy()\n",
    "    df_pval['log2FC'] = L2FC\n",
    "    df_pval['pval_MW'] = list_pvals_mannwhitney\n",
    "\n",
    "    # Add the corrected p-values to the dataset\n",
    "    _, pvals_corrected, _, _ = multipletests(df_pval['pval_MW'].values, alpha=0.05, method='fdr_bh')\n",
    "    df_pval['pval_MW_corrected'] = pvals_corrected\n",
    "\n",
    "\n",
    "    # Sort by p-values\n",
    "    df_pval = df_pval.sort_values(by=['pval_MW'])\n",
    "\n",
    "    return df_pval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.makedirs(f'{RESULTS_DIR}/merged_counts', exist_ok=True)\n",
    "os.makedirs(f'{RESULTS_DIR}/differential_abundance', exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for mode in [3, 5, 7]:\n",
    "    for S in [0, 1, 2, 3, 4, 5, 6, 7, 10, 15]:\n",
    "        print(f'MODE: {mode} | S {S}')\n",
    "        df_all_normplus, taxid_list, df_cut_normplus = process_samples(pass_num=2, mode=mode, S=S, NORM='NORM+', samples=samples, verbose=True, min_sample_flag=2)\n",
    "        df_cut_nan_percentage_normplus = filter_by_nan_percentage(df_cut_normplus, per_cutoff=0.35)\n",
    "\n",
    "        df_all_normpipe, taxid_list, df_cut_normpipe = process_samples(pass_num=2, mode=mode, S=S, NORM='NORMx', samples=samples, verbose=True, min_sample_flag=2)\n",
    "        df_cut_nan_percentage_normpipe = filter_by_nan_percentage(df_cut_normpipe, per_cutoff=0.35)\n",
    "\n",
    "        dict_filternorm_cut = filter_species_ids(df_cut_nan_percentage_normplus, df_cut_nan_percentage_normpipe)\n",
    "        with open(f'{RESULTS_DIR}/merged_counts/mode{mode}_S{S}_dict_norm+|_species.tsv', \"w\") as file:\n",
    "            json.dump(dict_filternorm_cut, file)\n",
    "        \n",
    "        print([(i, len(dict_filternorm_cut[i])) for i in dict_filternorm_cut.keys()])\n",
    "\n",
    "        df_cut_nan_percentage_normplus_discarded_common = df_cut_nan_percentage_normplus[df_cut_nan_percentage_normplus['taxonomy_id'].isin(dict_filternorm_cut['discarded_common'])]\n",
    "        df_cut_nan_percentage_normplus_discarded_common.to_csv(f'{RESULTS_DIR}/merged_counts/mode{mode}_S{S}_NORM+_discarded_common.tsv', sep='\\t', index=None)\n",
    "\n",
    "        df_cut_nan_percentage_normplus_discarded_normplus = df_cut_nan_percentage_normplus[df_cut_nan_percentage_normplus['taxonomy_id'].isin(dict_filternorm_cut['discarded_exclusive_norm_plus'])]\n",
    "        df_cut_nan_percentage_normplus_discarded_normplus.to_csv(f'{RESULTS_DIR}/merged_counts/mode{mode}_S{S}_NORM+_discarded_norm+.tsv', sep='\\t', index=None)\n",
    "        \n",
    "        df_cut_nan_percentage_normplus_retained = df_cut_nan_percentage_normplus[df_cut_nan_percentage_normplus['taxonomy_id'].isin(dict_filternorm_cut['retained_norm_plus'])]\n",
    "        df_cut_nan_percentage_normplus_retained.to_csv(f'{RESULTS_DIR}/merged_counts/mode{mode}_S{S}_NORM+_retained.tsv', sep='\\t', index=None)\n",
    "\n",
    "        df_cut_nan_percentage_normpipe_discarded_common = df_cut_nan_percentage_normpipe[df_cut_nan_percentage_normpipe['taxonomy_id'].isin(dict_filternorm_cut['discarded_common'])]\n",
    "        df_cut_nan_percentage_normpipe_discarded_common.to_csv(f'{RESULTS_DIR}/merged_counts/mode{mode}_S{S}_NORMx_discarded_common.tsv', sep='\\t', index=None)\n",
    "\n",
    "        df_cut_nan_percentage_normpipe_discarded_normplus = df_cut_nan_percentage_normpipe[df_cut_nan_percentage_normpipe['taxonomy_id'].isin(dict_filternorm_cut['discarded_exclusive_norm_plus'])]\n",
    "        df_cut_nan_percentage_normpipe_discarded_normplus.to_csv(f'{RESULTS_DIR}/merged_counts/mode{mode}_S{S}_NORMx_discarded_norm+.tsv', sep='\\t', index=None)\n",
    "\n",
    "        df_cut_nan_percentage_normpipe_retained = df_cut_nan_percentage_normpipe[df_cut_nan_percentage_normpipe['taxonomy_id'].isin(dict_filternorm_cut['retained_norm_pipe'])]\n",
    "        df_cut_nan_percentage_normpipe_retained.to_csv(f'{RESULTS_DIR}/merged_counts/mode{mode}_S{S}_NORMx_retained.tsv', sep='\\t', index=None)\n",
    "\n",
    "\n",
    "\n",
    "        df_pval_HCvsRR = differential_abundance_analysis(df_cut_nan_percentage_normplus_retained, ['HC1', 'HC2', 'HC3', 'HC4'], ['RR1', 'RR2', 'RR3', 'RR4'])\n",
    "        df_pval_HCvsRR.to_csv(f'{RESULTS_DIR}/differential_abundance/mode{mode}_S{S}_HCvsRR.tsv', sep='\\t', index=None)\n",
    "\n",
    "        df_pval_HCvsSP = differential_abundance_analysis(df_cut_nan_percentage_normplus_retained, ['HC1', 'HC2', 'HC3', 'HC4'], ['SP1', 'SP2', 'SP3', 'SP4'])\n",
    "        df_pval_HCvsSP.to_csv(f'{RESULTS_DIR}/differential_abundance/mode{mode}_S{S}_HCvsSP.tsv', sep='\\t', index=None)\n",
    "\n",
    "        df_pval_RRvsSP = differential_abundance_analysis(df_cut_nan_percentage_normplus_retained, ['RR1', 'RR2', 'RR3', 'RR4'], ['SP1', 'SP2', 'SP3', 'SP4'])\n",
    "        df_pval_RRvsSP.to_csv(f'{RESULTS_DIR}/differential_abundance/mode{mode}_S{S}_RRvsSP.tsv', sep='\\t', index=None)\n",
    "\n",
    "        df_pval_sex = differential_abundance_analysis(df_cut_nan_percentage_normplus_retained, ['HC1', 'HC2', 'RR1', 'RR2', 'SP1', 'SP2'], ['HC3', 'HC4', 'RR3', 'RR4', 'SP3', 'SP4'])\n",
    "        df_pval_sex.to_csv(f'{RESULTS_DIR}/differential_abundance/mode{mode}_S{S}_sex.tsv', sep='\\t', index=None)\n",
    "\n",
    "        print('\\n\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}
