Notebooks/Auxiliary/compute_performance_per_sample_test_set.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performance per sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('../../Src/')\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Read otus original, reconstructed with autoencoder and predicted from environmental features in test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_otus_original = pd.read_csv('../../Results/ReconstructionAndPredictionMicrobialComposition/otus_original_test.tsv', index_col=0, sep='\\t').T\n",
    "df_otus_reconstructed = pd.read_csv('../../Results/ReconstructionAndPredictionMicrobialComposition/otus_reconstAEfromBiome.tsv', index_col=0, sep='\\t').T\n",
    "df_otus_predicted = pd.read_csv('../../Results/ReconstructionAndPredictionMicrobialComposition/otus_predFromDomain.tsv', index_col=0, sep='\\t').T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_otus_original = df_otus_original.to_numpy(dtype=np.float32)\n",
    "data_otus_reconstructed = df_otus_reconstructed.to_numpy(dtype=np.float32)\n",
    "data_otus_predicted = df_otus_predicted.to_numpy(dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df_otus_original.shape)\n",
    "print(df_otus_reconstructed.shape)\n",
    "print(df_otus_predicted.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Absolute abundance transformed to TSS (with epsilon=1E-6)\n",
    "def transform_to_rel_abundance(dataset):\n",
    "    epsilon=1E-6\n",
    "    sum_per_sample = dataset.sum(axis=1)\n",
    "    num_samples = sum_per_sample.shape\n",
    "    num_OTUs = np.shape(dataset)[-1] \n",
    "    sum_per_sample = sum_per_sample + (num_OTUs * epsilon)\n",
    "    dividend=dataset+epsilon\n",
    "    dataset_rel_abund = np.divide(dividend,sum_per_sample[:,None])\n",
    "    #display(Markdown(\"{}</p>\".format(np.array2string(actual_array,precision=6,floatmode='fixed'))))\n",
    "    #actual_array.sum(axis=1)\n",
    "    return dataset_rel_abund"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_otus_original_rel = transform_to_rel_abundance(data_otus_original)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pearson correlation per sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "def compute_pearson_correlation_per_sample(original,predicted):\n",
    "    tot = 0.0\n",
    "    count = 0\n",
    "    corr_sample=np.zeros(original.shape[0],dtype=np.float32)\n",
    "    for i, (actual,pred) in enumerate(zip(original,predicted)): \n",
    "        r, _ = stats.pearsonr(actual,pred)\n",
    "        if not np.isnan(r):\n",
    "            count += 1\n",
    "            tot += r\n",
    "            #print(i,\": \",r)\n",
    "            corr_sample[i]=r\n",
    "    r_final = tot/count\n",
    "    print('Average Pearson: ',r_final)\n",
    "    return corr_sample   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pearson correlation: original vs reconstructed (AE)\n",
    "pearson_per_sample_reconst = compute_pearson_correlation_per_sample(data_otus_original_rel,data_otus_reconstructed)\n",
    "# Pearson correlation: original vs predicted (from env.features)\n",
    "pearson_per_sample_predicted = compute_pearson_correlation_per_sample(data_otus_original_rel,data_otus_predicted)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bray-Curtis per sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skbio.diversity import beta_diversity\n",
    "\n",
    "def compute_braycurtis_distance_per_sample(original,predicted):\n",
    "    tot = 0.0\n",
    "    count = 0\n",
    "    bc_sample=np.zeros(original.shape[0],dtype=np.float32)\n",
    "    for i, (actual,pred) in enumerate(zip(original,predicted)): \n",
    "        bc_dm = beta_diversity(\"braycurtis\", [actual,pred]) # Source: http://scikit-bio.org/docs/0.4.2/diversity.html\n",
    "        bc = bc_dm[0,1]\n",
    "        count += 1\n",
    "        tot += bc\n",
    "        #print(i,\": \",bc)\n",
    "        bc_sample[i]=bc\n",
    "    bc_final = tot/count\n",
    "    print('Average Bray-Curtis: ',bc_final)\n",
    "    return bc_sample   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pearson correlation: original vs reconstructed (AE)\n",
    "braycurtis_per_sample_reconst = compute_braycurtis_distance_per_sample(data_otus_original_rel,data_otus_reconstructed)\n",
    "# Pearson correlation: original vs predicted (from env.features)\n",
    "braycurtis_per_sample_predicted = compute_braycurtis_distance_per_sample(data_otus_original_rel,data_otus_predicted)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Save performance per sample in a file\n",
    "To plot a graph in R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_performance_per_sample(perf1,perf2,perf3,perf4,sample_names,suffix=''):\n",
    "    df_performance = pd.DataFrame(perf1, index=sample_names, columns=['pearson_reconstructed'])\n",
    "    df_temp = pd.DataFrame(perf2, index=sample_names, columns=['pearson_predicted'])\n",
    "    df_performance = df_performance.join(df_temp)\n",
    "    df_temp = pd.DataFrame(perf3, index=sample_names, columns=['braycurtis_reconstructed'])\n",
    "    df_performance = df_performance.join(df_temp)\n",
    "    df_temp = pd.DataFrame(perf4, index=sample_names, columns=['braycurtis_predicted'])\n",
    "    df_performance = df_performance.join(df_temp)\n",
    "    df_performance.to_csv(suffix+'.tsv', index=True, header=True, sep='\\t')\n",
    "    return df_performance\n",
    "\n",
    "df_perf = save_performance_per_sample( \\\n",
    "            pearson_per_sample_reconst,pearson_per_sample_predicted,\\\n",
    "            braycurtis_per_sample_reconst,braycurtis_per_sample_predicted,\\\n",
    "            df_otus_original.index,'../../Results/ReconstructionAndPredictionMicrobialComposition/performance_per_sample')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}