Notebooks/reference_model_predictions_and_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('../Src/')\n",
    "from data import *\n",
    "from transfer_learning import *\n",
    "from test_functions import *\n",
    "from layers import *\n",
    "from utils import *\n",
    "from loss import *\n",
    "from metric import *\n",
    "from results import *\n",
    "import tensorflow as tf\n",
    "import tensorflow.keras as keras\n",
    "from tensorflow.keras import layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "np.set_printoptions(precision=4)\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reference model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_microbioma_train, df_microbioma_test, _, _, \\\n",
    "df_domain_train, df_domain_test, _, _, otu_columns, domain_columns = \\\n",
    "    read_df_with_transfer_learning_subset_fewerDomainFeatures( \\\n",
    "        metadata_names=['age','Temperature','Precipitation3Days'], \\\n",
    "        otu_filename='../Datasets/otu_table_all_80.csv', \\\n",
    "        metadata_filename='../Datasets/metadata_table_all_80.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_microbioma_train = df_microbioma_train.to_numpy(dtype=np.float32)\n",
    "data_microbioma_test = df_microbioma_test.to_numpy(dtype=np.float32)\n",
    "data_domain_train = df_domain_train.to_numpy(dtype=np.float32)\n",
    "data_domain_test = df_domain_test.to_numpy(dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_microbioma_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train the selected model (the best one from those with the smallest latent space (10)): no.351\n",
    "experiment_metrics, models, results = perform_experiment_2(cv_folds=0, \n",
    "                        epochs=100, \n",
    "                        batch_size=64, \n",
    "                        learning_rate=0.001, \n",
    "                        optimizer=optimizers.Adam,\n",
    "                        learning_rate_scheduler=None,\n",
    "                        input_transform=Percentage,\n",
    "                        output_transform=tf.keras.layers.Softmax,\n",
    "                        reconstruction_loss=MakeLoss(LossBrayCurtis, Percentage, None), \n",
    "                        latent_space=10, \n",
    "                        layers=[512,256],\n",
    "                        activation='tanh', \n",
    "                        activation_latent='tanh', \n",
    "                        data_microbioma_train=data_microbioma_train,\n",
    "                        data_domain_train=None,\n",
    "                        show_results=True, \n",
    "                        device='/CPU:0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To get encoders and decoders to use in domain->latent model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model, encoder, _ ,decoder = models[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_domain_train.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To predict latent space for samples in domain->latent model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latent_train = encoder.predict(data_microbioma_train)\n",
    "latent_test = encoder.predict(data_microbioma_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To build model to predict latent space "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_fn_latent():\n",
    "    in_layer = layers.Input(shape=(3,))\n",
    "    net = layers.Dense(128, activation='tanh')(in_layer)\n",
    "    net = layers.Dense(64, activation='tanh')(net)\n",
    "    net = layers.Dense(32, activation='tanh')(net)\n",
    "    net = layers.Dense(16, activation='tanh')(net)\n",
    "    out_layer = layers.Dense(latent_train.shape[1], activation=None)(net) # 'tanh already'\n",
    "    model = keras.Model(inputs=[in_layer], outputs=[out_layer], name='model')\n",
    "    model.compile(optimizer=optimizers.Adam(lr=0.001), loss=tf.keras.losses.MeanSquaredError(),\n",
    "                 metrics=[tf.keras.metrics.MeanSquaredError()])\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_latent, model_latent = train_tl_noEnsemble(model_fn_latent,\n",
    "                            latent_train,\n",
    "                            latent_train,\n",
    "                            data_domain_train,\n",
    "                            data_domain_train,\n",
    "                            epochs=100,\n",
    "                            batch_size=16,\n",
    "                            verbose=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "predictions = test_model_tl_noEnsemble(model_latent, decoder, Percentage, tf.keras.layers.Softmax, otu_columns, data_microbioma_test, data_domain_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save_predictions(predictions, 'experiment_testSet_AE_OTUlatent_3var.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoder.save('../Results/Models/encoder_biome.h5')\n",
    "decoder.save('../Results/Models/decoder.h5')\n",
    "model_latent.save('../Results/Models/encoder_domain_model_latent.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoder_biome = encoder\n",
    "encoder_domain = model_latent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Predictions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save predictions in files for external use (R plots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_microbioma_test.T.to_csv('../Results/otus_original_test.tsv', index=True, header=True, sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_predicted_otu_table_and_latent(pred,pred_latent,sample_names,otu_names,suffix=''):\n",
    "    df_otu = pd.DataFrame(pred, index=sample_names, columns=otu_names)\n",
    "    df_otu.T.to_csv('../Results/otus_'+suffix+'.tsv', index=True, header=True, sep='\\t')\n",
    "\n",
    "    df_latent = pd.DataFrame(pred_latent, index=sample_names)\n",
    "    df_latent.T.to_csv('../Results/latent_'+suffix+'.tsv', index=True, sep='\\t')\n",
    "    \n",
    "    return df_otu, df_latent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run prediction test set from microbiome, i.e. reconstructed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input only domain (i.e. environmental features)\n",
    "pred_latent_biome = encoder_biome.predict(data_microbioma_test)\n",
    "pred_biome = decoder.predict(pred_latent_biome)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, _ = save_predicted_otu_table_and_latent(pred_biome,pred_latent_biome,df_microbioma_test.index,df_microbioma_test.columns,'reconstAEfromBiome')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run prediction test set from domain (i.e. environmental features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input only domain (i.e. environmental features)\n",
    "pred_latent = encoder_domain.predict(data_domain_test)\n",
    "pred_domain = decoder.predict(pred_latent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pred_otu, df_pred_latent = save_predicted_otu_table_and_latent(pred_domain,pred_latent,df_microbioma_test.index,df_microbioma_test.columns,'predFromDomain')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute metrics by OTU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compute Relative Squared Error per OTU"
   ]
  },
  {
   "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": [
    "def compute_relative_squared_error(actual,pred):\n",
    "    rse_otu=np.zeros(actual.shape[1],dtype=np.float32)\n",
    "    actual=actual_array.transpose()\n",
    "    pred=pred_domain.transpose()\n",
    "    # for each OTU\n",
    "    for i, (act_otu,pred_otu) in enumerate(zip(actual,pred)):\n",
    "        mean_otu = act_otu.mean()\n",
    "        #display(Markdown(\"{} (mean(act)): {}</p>\".format(i,mean_otu)))\n",
    "        #display(Markdown(\"{} (act): {}</p>\".format(i,np.array2string(act_otu,precision=8,floatmode='fixed'))))\n",
    "        #display(Markdown(\"{} (pred): {}</p>\".format(i,np.array2string(pred_otu,precision=8,floatmode='fixed'))))\n",
    "        div_up=((pred_otu-act_otu)**2).sum()\n",
    "        div_down=((pred_otu-mean_otu)**2).sum()\n",
    "        rse_otu[i]=div_up/div_down\n",
    "    return rse_otu, np.sqrt(rse_otu)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "actual_array = transform_to_rel_abundance(data_microbioma_test)\n",
    "RSE_perOTU, RRSE_perOTU  = compute_relative_squared_error(actual_array,pred_domain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RRSE_perOTU.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RRSE_perOTU.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RRSE_perOTU.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Save errors per OTU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_errors_per_OTU(RSE,RRSE,otu_names,suffix=''):\n",
    "    df_RSE = pd.DataFrame(RSE, index=otu_names, columns=['RSE'])\n",
    "    df_RRSE = pd.DataFrame(RRSE, index=otu_names, columns=['RRSE'])\n",
    "    df_error = df_RSE.join(df_RRSE)\n",
    "    df_error.to_csv(suffix+'.tsv', index=True, header=True, sep='\\t')\n",
    "    \n",
    "    return df_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_error_perOTU = save_errors_per_OTU(RSE_perOTU,RRSE_perOTU,df_microbioma_test.columns,'../Results/errors_perOTU')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_error_perOTU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Novel samples: prediction of novel ecosystems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metadata = pd.read_csv('../Datasets/NovelEcosystems/metadata_novel_samples_only3envFeatures.csv', sep='\\t')\n",
    "metadata = metadata.set_index('X.SampleID')\n",
    "domain = metadata[['age',\n",
    "                   'Temperature',\n",
    "                   'Precipitation3Days']]\n",
    "print(domain.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "domain.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "domain_novel_samples = domain.to_numpy(dtype=np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decoder = tf.keras.models.load_model('../Results/Models/decoder.h5')\n",
    "encoder_domain = tf.keras.models.load_model('../Results/Models/encoder_domain_model_latent.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compute new predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input only domain (i.e. environmental features)\n",
    "pred_latent_novel_samples = encoder_domain.predict(domain_novel_samples)\n",
    "pred_domain_novel_samples = decoder.predict(pred_latent_novel_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pred_otu, df_pred_latent = save_predicted_otu_table_and_latent(pred_domain_novel_samples,pred_latent_novel_samples,domain.index,df_microbioma_test.columns,'predFromDomain_novelSamples')"
   ]
  }
 ],
 "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": 1
}