a-modelling-HPV.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Multi-Omics Analysis & Modelling of Cervicovaginal Microenvironment\n",
    "\n",
    "Notebook accompanies publication: \"Multi-omics data integration reveals metabolome as the top predictor of the cervicovaginal microenvironment\".\n",
    "\n",
    "It includes:\n",
    "\n",
    "[1) Import packages and set paths](#chap1)      \n",
    "\n",
    "[2) Read & process datasets](#chap2)\n",
    "\n",
    "[3) Exploratory analysis of omics datasets](#chap3)\n",
    "\n",
    "[4) Training and evaluation of classifiers](#chap4)\n",
    "\n",
    "[5) Training and evaluation of regressors](#chap5)\n",
    "\n",
    "**To run this notebook, please setup the conda environment as described in `README.md` section \"Setup for `a-modelling-HPV.ipynb`\".**\n",
    "         "
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "<a id='chap1'></a>"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 1) Import packages and set paths"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "import warnings\n",
    "import os\n",
    "import pandas as pd\n",
    "import skbio\n",
    "from itertools import compress\n",
    "from skbio.stats.ordination import pcoa\n",
    "\n",
    "import qiime2\n",
    "from qiime2.plugins import diversity\n",
    "\n",
    "import matplotlib as mpl\n",
    "from matplotlib.pylab import plt\n",
    "from IPython.display import Image\n",
    "\n",
    "# custom functions\n",
    "from util_data_proc import (add_microbiome_data,\n",
    "                            add_metabolome_data,\n",
    "                            perform_taxonomic_classification, add_targets,\n",
    "                            extract_microbiome_artifact)\n",
    "from util_eda import (plot_data_avail_per_target,\n",
    "                      return_pcoa_metrics_microbiome,\n",
    "                      calc_pca_metrics_metabolome,\n",
    "                      calc_pca_metrics_proteome,\n",
    "                      merge_all_pca_metrics, run_manova)\n",
    "from util_classification import train_n_eval_classifier\n",
    "from util_regression import train_n_eval_regressors\n",
    "\n",
    "# R settings\n",
    "import rpy2\n",
    "\n",
    "plt.rcParams['axes.grid'] = True\n",
    "plt.style.use('seaborn-pastel')\n",
    "mpl.rcParams['figure.dpi'] = 250\n",
    "\n",
    "# jupyter magic commands\n",
    "%load_ext rpy2.ipython\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# location to retrieve raw-data\n",
    "path2data = 'data-raw'\n",
    "\n",
    "# location to save outputs\n",
    "output_dir = 'a-modelling-output'\n",
    "if not os.path.exists(output_dir):\n",
    "    os.makedirs(output_dir)\n",
    "\n",
    "# set random seed across all experiments\n",
    "seed = 12\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "%%R -w 800 --type=cairo\n",
    "\n",
    "# importing all required R packages required for omics visualisation\n",
    "source(\"util_R_plotter.R\")"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "<a id='chap2'></a>"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 2) Read & process datasets"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Get all features"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# # read patient data (prefix 'F_pcov_') + immunoproteo data (prefix 'F_proteo_')\n",
    "df_dataset = pd.read_csv(os.path.join(path2data, 'sample_metadata_n_proteome.tsv'), sep='\\t',\n",
    "                         index_col='sample-id')\n",
    "print('Number of patient covariate features (tagged with F_pcov_): {}'.format(\n",
    "    len([x for x in df_dataset.columns if x.startswith('F_pcov_')])))\n",
    "print('Number of immunoproteome features (tagged with F_proteo_): {}'.format(\n",
    "    len([x for x in df_dataset.columns if x.startswith('F_proteo_')])))\n",
    "\n",
    "# read cancer biomarkers - needed for cancer biomarker regression\n",
    "df_cancer_biomarkers = pd.read_csv(os.path.join(\n",
    "    'data-raw', 'cancer_biomarkers.tsv'), sep='\\t')\n",
    "ls_cancer_biomarkers = df_cancer_biomarkers['cancer_biomarkers'].values.tolist(\n",
    ")\n",
    "\n",
    "# add microbiome data\n",
    "df_dataset = add_microbiome_data(df_dataset, path2data)\n",
    "\n",
    "# add metabolome data\n",
    "df_dataset = add_metabolome_data(df_dataset, path2data)\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Perform taxonomic classification (`taxa`)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# perform taxonomic classification\n",
    "merged_taxonomy, taxonomy_qza = perform_taxonomic_classification(path2data)\n",
    "\n",
    "# extract consensus taxonomy\n",
    "taxonomy = merged_taxonomy[['Consensus Taxonomy']].copy(deep=True)\n",
    "taxonomy.columns = ['Taxon']\n",
    "\n",
    "# extract lowest available level of taxonomic classification - required for eval of classifiers\n",
    "taxa = taxonomy['Taxon'].apply(lambda x: x.split(';')[-1] if len(x.split(';')[-1])>0 else x.split(';')[-2:][0] + ' unknown')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Get all Targets (`df_data_w_targets`)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# add targets ('T_')\n",
    "df_data_w_targets = add_targets(df_dataset, path2data, taxonomy_qza)\n",
    "\n",
    "# Targets to consider\n",
    "ls_targets = [col for col in df_data_w_targets.columns if col.startswith(\n",
    "    'T_') and col != 'T_infl_score_flt']\n",
    "ls_targets\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Select samples with all omics data available"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# get omics availability in original dataset per target class\n",
    "\n",
    "plot_data_avail_per_target(df_data_w_targets, ls_targets)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# remove all samples that don't have all omics and patient covariate features available\n",
    "ls_microbiome_cols = [\n",
    "    x for x in df_data_w_targets.columns if x.startswith('F_micro')]\n",
    "ls_metabolome_cols = [\n",
    "    x for x in df_data_w_targets.columns if x.startswith('F_metabo')]\n",
    "ls_proteome_cols = [\n",
    "    x for x in df_data_w_targets.columns if x.startswith('F_proteo')]\n",
    "ls_pcov = [x for x in df_data_w_targets.columns if x.startswith('F_pcov_')]\n",
    "\n",
    "nb_samples_before = df_data_w_targets.shape[0]\n",
    "print('Total number of samples before: {}'.format(nb_samples_before))\n",
    "ls_all_omics = ls_microbiome_cols + \\\n",
    "    ls_metabolome_cols + ls_proteome_cols + ls_pcov\n",
    "\n",
    "df_data_w_targets_sel = df_data_w_targets[\n",
    "    ~df_data_w_targets.loc[:, ls_all_omics+ls_targets]\n",
    "    .isnull().any(axis=1)]\n",
    "print('Total number of samples with all info: {}'.format(\n",
    "    df_data_w_targets_sel.shape[0]))\n",
    "\n",
    "print('Class count in all info available:')\n",
    "for target in ls_targets:\n",
    "    print(df_data_w_targets_sel[target].value_counts(dropna=False))\n",
    "\n",
    "# remove other dataframe - ensuring only selected is used below\n",
    "del df_data_w_targets\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Remove features with same value in all samples  (`df_data_w_targets_sel`)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# remove micro features that don't occur in any of the\n",
    "# 72 samples (value == 0.0 everywhere) - 68 features\n",
    "\n",
    "print('Before: {}'.format(df_data_w_targets_sel.shape))\n",
    "\n",
    "ls_microbiome_features = [\n",
    "    x for x in df_data_w_targets_sel.columns if x.startswith('F_micro_')]\n",
    "bool_missing_micro = (\n",
    "    df_data_w_targets_sel[ls_microbiome_features] == 0.0).all().values\n",
    "ls_missing_micro = list(compress(ls_microbiome_features, bool_missing_micro))\n",
    "ls_features2select = [\n",
    "    x for x in df_data_w_targets_sel.columns if x not in ls_missing_micro]\n",
    "\n",
    "df_data_w_targets_sel = df_data_w_targets_sel[ls_features2select].copy(\n",
    "    deep=True)\n",
    "print('After: {}'.format(df_data_w_targets_sel.shape))\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# remove metabolites that all have the same value (stddev = 0.0) (8 features affected)\n",
    "print('Before: {}'.format(df_data_w_targets_sel.shape))\n",
    "ls_metabolite_features = [\n",
    "    x for x in df_data_w_targets_sel.columns if x.startswith('F_metabo_')]\n",
    "\n",
    "std_metabo = df_data_w_targets_sel[ls_metabolite_features].describe(\n",
    ").loc['std']\n",
    "ls_zero_stddev = std_metabo[std_metabo == 0.0].index.to_list()\n",
    "feat2select = [\n",
    "    x for x in df_data_w_targets_sel.columns if x not in ls_zero_stddev]\n",
    "\n",
    "df_data_w_targets_sel = df_data_w_targets_sel[feat2select].copy(deep=True)\n",
    "print('After: {}'.format(df_data_w_targets_sel.shape))\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# proteome does not have any features where all values are same\n",
    "ls_proteome = [\n",
    "    x for x in df_data_w_targets_sel.columns if x.startswith('F_proteo_')]\n",
    "std_proteo = df_data_w_targets_sel[ls_proteome].describe().loc['std']\n",
    "ls_zero_stddev = std_proteo[std_proteo == 0.0].index.to_list()\n",
    "print(len(ls_zero_stddev))\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# extract updated feature names\n",
    "ls_microbiome_cols = [x for x in df_data_w_targets_sel.columns if x.startswith('F_micro')]\n",
    "ls_metabolome_cols = [x for x in df_data_w_targets_sel.columns if x.startswith('F_metabo')]\n",
    "ls_proteome_cols = [x for x in df_data_w_targets_sel.columns if x.startswith('F_proteo')]"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "<a id='chap3'></a>"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 3) Exploratory analysis"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# create subfolder for exploratory data analysis results\n",
    "output_dir_eda = os.path.join(output_dir, 'eda')\n",
    "if not os.path.exists(output_dir_eda):\n",
    "    os.makedirs(output_dir_eda)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Microbiome: alpha diversity analysis"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "microbiome_table = extract_microbiome_artifact(df_data_w_targets_sel)\n",
    "\n",
    "# calculating alpha diversity\n",
    "art_alpha_diversity,  = diversity.actions.alpha(\n",
    "    table=microbiome_table, metric='shannon')\n",
    "df_alpha_div = art_alpha_diversity.view(pd.Series).to_frame()\n",
    "df_alpha_div.index.name = 'SampleID'\n",
    "\n",
    "# compare alpha group significance\n",
    "alpha_div_viz, = diversity.actions.alpha_group_significance(\n",
    "    alpha_diversity=art_alpha_diversity,\n",
    "    metadata=qiime2.Metadata(df_data_w_targets_sel[ls_targets]))\n",
    "\n",
    "print('Alpha Div Boxplots per cat class and Kruskal-Wallis '\n",
    "      'test results are saved as qiime viz file in:')\n",
    "alpha_div_viz.save(os.path.join(output_dir_eda, 'alpha-div-viz.qzv'))\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Microbiome: beta diversity analysis"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# calculate beta diversity (bray-curtis, jaccard) & perform PCoA from distance matrix\n",
    "dict_beta_results = {}\n",
    "\n",
    "for beta_metric in ['braycurtis', 'jaccard']:\n",
    "    art_beta_div_dis_matrix, = diversity.actions.beta(\n",
    "        table=microbiome_table, metric=beta_metric)\n",
    "\n",
    "    beta_pcoa_results = pcoa(art_beta_div_dis_matrix.view(\n",
    "        skbio.DistanceMatrix), number_of_dimensions=2)\n",
    "\n",
    "    dict_beta_results[beta_metric+'_dis_matrix'] = art_beta_div_dis_matrix\n",
    "    dict_beta_results[beta_metric+'_pcoa_res'] = beta_pcoa_results\n",
    "\n",
    "# analyse stat. trends of beta diversity with PERMANOVA\n",
    "# = tests the hypothesis that samples within a group are more\n",
    "# similar to each other than they are to samples in another group\n",
    "for metric in ['braycurtis', 'jaccard']:\n",
    "    for target in ls_targets:\n",
    "        beta_div_viz, = diversity.actions.beta_group_significance(\n",
    "            distance_matrix=dict_beta_results[metric+'_dis_matrix'],\n",
    "            metadata=qiime2.Metadata(\n",
    "                df_data_w_targets_sel[ls_targets]).get_column(target),\n",
    "            pairwise=False)\n",
    "\n",
    "        print('Beta Div {} PERMANOVA results for {}:'.format(metric, target))\n",
    "        path2save = os.path.join(\n",
    "            output_dir_eda, 'beta-div-permanova-{}-{}.qzv'.format(metric, target))\n",
    "        print(path2save)\n",
    "        beta_div_viz.save(path2save)\n",
    "        print()\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Omics: PCoA beta-div of micro- and PCA of meta-/proteome"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# microbiome pcoa metrics\n",
    "dict_explained_var = {}\n",
    "df_pcoa_micro, dict_explained_var = return_pcoa_metrics_microbiome(dict_beta_results,\n",
    "                                                                   dict_explained_var,\n",
    "                                                                   df_data_w_targets_sel,\n",
    "                                                                   ls_targets)\n",
    "# metabolome pca metrics\n",
    "df_pca_metab, dict_explained_var = calc_pca_metrics_metabolome(dict_explained_var,\n",
    "                                                               df_data_w_targets_sel,\n",
    "                                                               ls_targets)\n",
    "\n",
    "# proteome pca metrics\n",
    "df_pca_proteo, dict_explained_var = calc_pca_metrics_proteome(dict_explained_var,\n",
    "                                                              df_data_w_targets_sel,\n",
    "                                                              ls_targets)\n",
    "                                                              \n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# merge derived metrics to one df for explained variance\n",
    "# and pca coordinates with targets\n",
    "beta_div2_choose = 'jaccard'\n",
    "df_explained_var = pd.DataFrame(dict_explained_var)\n",
    "df_pca_all = merge_all_pca_metrics(df_pcoa_micro, beta_div2_choose,\n",
    "                                   df_pca_metab, df_pca_proteo)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Calculate Manova pvalue for metabo/proteome features"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "ls_omics = ['Metabolome', 'Immunoproteome']\n",
    "df_manova_results = pd.DataFrame(columns=ls_targets, index=ls_omics)\n",
    "\n",
    "for target in ls_targets:\n",
    "    for omics in ls_omics:\n",
    "        df_pca_sel = df_pca_all[df_pca_all['Omics']==omics].copy(deep=True)\n",
    "        ls_features = ['PC1', 'PC2']\n",
    "        pvalue = run_manova(df_pca_sel, ls_features, target)\n",
    "        df_manova_results.loc[omics,target] = round(pvalue,3)\n",
    "\n",
    "df_manova_results"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "visualise in ggplot2"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "%%R -w 800 --type=cairo --input beta_div2_choose,df_pca_all,df_explained_var,ls_targets,output_dir_eda --res 100 -u px 9000\n",
    "\n",
    "plot_pca_scatters(beta_div2_choose,df_pca_all,df_explained_var,ls_targets,output_dir_eda)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Notes on PC(o)A plots:           \n",
    "`stat_ellipse` calculates confidence interval ellipse assuming a multivariate normal distribution (`type='norm'`) drawn at a 95 % confidence level (`level=0.95`).\n"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "print(\"Plotting image saved in R function `plot_pca_scatters`:\")\n",
    "Image(filename=os.path.join(output_dir_eda, \"omics-pca-ggplot-log.png\"))\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "<a id='chap4'></a>"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 4) Train and evaluate classifiers"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict `disease_status`"
   ],
   "metadata": {
    "heading_collapsed": true
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "ls_features_disease = [\n",
    "    x for x in df_data_w_targets_sel.columns if x.startswith('F_')]\n",
    "\n",
    "ls_class_order = ['Ctrl_HPV_neg', 'Ctrl_HPV_pos', 'LGD', 'HGD', 'ICC']\n",
    "dic_color_palette = {'Ctrl_HPV_neg': 'blue',\n",
    "                     'Ctrl_HPV_pos': 'green',\n",
    "                     'LGD': 'yellow',\n",
    "                     'HGD': 'orange',\n",
    "                     'ICC': 'red'}\n",
    "\n",
    "train_n_eval_classifier('T_disease_state',\n",
    "                        ls_features_disease,\n",
    "                        df_data_w_targets_sel,\n",
    "                        taxa,\n",
    "                        ls_class_order,\n",
    "                        dic_color_palette,\n",
    "                        seed,\n",
    "                        output_dir)\n"
   ],
   "outputs": [],
   "metadata": {
    "hidden": true
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict `genital_inflammation`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "ls_inflammatory_markers = ['CCL20', 'IL-1alpha',\n",
    "                           'IL-1beta', 'IL-8',\n",
    "                           'MIP-1beta', 'RAES', 'TNFalpha']\n",
    "ls_inflammatory_markers = [\"F_proteo_\" + x for x in ls_inflammatory_markers]\n",
    "ls_features4inflammation = [x for x in df_data_w_targets_sel.columns if (x.startswith('F_')\n",
    "                                                                         and x not in ls_inflammatory_markers\n",
    "                                                                         and not x.startswith('F_pcov_'))]\n",
    "\n",
    "ls_class_order = ['None', 'Low', 'High']\n",
    "dic_color_palette = {'None': 'grey', 'Low': 'blue', 'High': 'red'}\n",
    "\n",
    "train_n_eval_classifier('T_inflammation',\n",
    "                        ls_features4inflammation,\n",
    "                        df_data_w_targets_sel,\n",
    "                        taxa,\n",
    "                        ls_class_order,\n",
    "                        dic_color_palette,\n",
    "                        seed,\n",
    "                        output_dir)\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict `lactobacillus_dominance`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "ls_features4lacto = [x for x in df_data_w_targets_sel.columns if\n",
    "                     (x.startswith('F_metabo_') or x.startswith('F_proteo_'))]\n",
    "\n",
    "ls_class_order=['LD', 'NLD']\n",
    "dic_color_palette={'LD': 'blue', 'NLD': 'red'}\n",
    "\n",
    "train_n_eval_classifier('T_lactobacillus_dominance',\n",
    "                        ls_features4lacto,\n",
    "                        df_data_w_targets_sel,\n",
    "                        taxa,\n",
    "                        ls_class_order,\n",
    "                        dic_color_palette,\n",
    "                        seed,\n",
    "                        output_dir)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict `vaginal_pH`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "ls_features4ph = [x for x in df_data_w_targets_sel.columns if\n",
    "                  ((x.startswith('F_') and not x.startswith('F_pcov_')))]\n",
    "\n",
    "ls_class_order=['Low', 'High']\n",
    "dic_color_palette={'Low': 'blue', 'High': 'red'}\n",
    "\n",
    "train_n_eval_classifier('T_pH',\n",
    "                        ls_features4ph,\n",
    "                        df_data_w_targets_sel,\n",
    "                        taxa,\n",
    "                        ls_class_order,\n",
    "                        dic_color_palette,\n",
    "                        seed,\n",
    "                        output_dir)\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "<a id='chap5'></a>"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 5) Train and evaluate regressors"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict `metabolites`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# features\n",
    "ls_features_metabolites = [x for x in df_data_w_targets_sel.columns if\n",
    "                           ((x.startswith('F_micro_') or x.startswith('F_proteo_')))]\n",
    "\n",
    "# targets: 94 selected metabolites\n",
    "ls_targets_metabolites = list(set([\n",
    "    '1-(1-enyl-palmitoyl)-2-arachidonoyl-GPC (P-16:0/20:4)',\n",
    "    '1-(1-enyl-palmitoyl)-2-arachidonoyl-GPE (P-16:0/20:4)',\n",
    "    '1-(1-enyl-palmitoyl)-2-linoleoyl-GPE (P-16:0/18:2)',\n",
    "    '1-(1-enyl-palmitoyl)-2-oleoyl-GPC (P-16:0/18:1)',\n",
    "    '1-(1-enyl-palmitoyl)-2-oleoyl-GPE (P-16:0/18:1)',\n",
    "    '1-(1-enyl-palmitoyl)-2-palmitoyl-GPC (P-16:0/16:0)',\n",
    "    '1-(1-enyl-stearoyl)-2-arachidonoyl-GPE (P-18:0/20:4)',\n",
    "    '1-(1-enyl-stearoyl)-2-linoleoyl-GPE (P-18:0/18:2)',\n",
    "    '1-(1-enyl-stearoyl)-2-oleoyl-GPE (P-18:0/18:1)',\n",
    "    '3-indoxyl sulfate',\n",
    "    '4-acetamidobutanoate',\n",
    "    '5-methylthioadenosine (MTA)',\n",
    "    'agmatine',\n",
    "    'arachidonate (20:4n6)',\n",
    "    'C-glycosyltryptophan',\n",
    "    'diacetylspermidine',\n",
    "    'dihomo-linolenate (20:3n3 or n6)',\n",
    "    'glycochenodeoxycholate',\n",
    "    'indolelactate',\n",
    "    'kynurenate',\n",
    "    'kynurenine',\n",
    "    'linoleate (18:2n6)',\n",
    "    'linolenate [alpha or gamma; (18:3n3 or 6)]',\n",
    "    \"N('1)-acetylspermidine\",\n",
    "    'N(1)-acetylspermine',\n",
    "    'N1,N12-diacetylspermine',\n",
    "    'N-acetyl-isoputreanine',\n",
    "    'N-acetylputrescine',\n",
    "    'putrescine',\n",
    "    'spermidine',\n",
    "    'spermine',\n",
    "    'tryptamine',\n",
    "    'tryptophan',\n",
    "    'tryptophan betaine',\n",
    "    '2-hydroxyhippurate (salicylurate)',\n",
    "    '3-hydroxyhippurate',\n",
    "    '4-hydroxyhippurate',\n",
    "    'hippurate',\n",
    "    'oleate/vaccenate (18:1)',\n",
    "    'histamine',\n",
    "    '1-methylhistamine',\n",
    "    'tryptamine',\n",
    "    'tyramine',\n",
    "    'agmatine',\n",
    "    'cadaverine',\n",
    "    'N-acetyl-cadaverine',\n",
    "    'N-acetylputrescine',\n",
    "    'putrescine',\n",
    "    'diacetylspermidine',\n",
    "    \"N('1)-acetylspermidine\",\n",
    "    'N(1)-acetylspermine',\n",
    "    'N1,N12-diacetylspermine',\n",
    "    'spermidine',\n",
    "    'spermine',\n",
    "    '1-(1-enyl-palmitoyl)-2-oleoyl-GPC (P-16:0/18:1)',\n",
    "    '1,2-dilinoleoyl-GPC (18:2/18:2)', '1,2-dipalmitoyl-GPC (16:0/16:0)',\n",
    "    '1-linoleoyl-2-arachidonoyl-GPC (18:2/20:4n6)',\n",
    "    '1-linoleoyl-GPC (18:2)', '1-oleoyl-2-linoleoyl-GPC (18:1/18:2)',\n",
    "    '1-oleoyl-GPC (18:1)', '1-palmitoyl-2-arachidonoyl-GPC (16:0/20:4n6)',\n",
    "    '1-palmitoyl-2-arachidonoyl-GPE (16:0/20:4)',\n",
    "    '1-palmitoyl-2-docosahexaenoyl-GPC (16:0/22:6)',\n",
    "    '1-palmitoyl-2-linoleoyl-GPC (16:0/18:2)',\n",
    "    '1-palmitoyl-2-oleoyl-GPC (16:0/18:1)',\n",
    "    '1-palmitoyl-2-palmitoleoyl-GPC (16:0/16:1)', '1-palmitoyl-GPC (16:0)',\n",
    "    '1-stearoyl-2-arachidonoyl-GPC (18:0/20:4)',\n",
    "    '1-stearoyl-2-docosahexaenoyl-GPC (18:0/22:6)',\n",
    "    '1-stearoyl-2-oleoyl-GPC (18:0/18:1)', '1-stearoyl-GPC (18:0)',\n",
    "    '2-hydroxyhippurate (salicylurate)', '3-hydroxybutyrate (BHBA)',\n",
    "    'behenoyl sphingomyelin (d18:1/22:0)', 'cefazolin', 'cholesterol',\n",
    "    'cytosine', 'deoxycarnitine', 'glutathione, oxidized (GSSG)',\n",
    "    'glycerol', 'glycerophosphoglycerol', 'imidazole propionate',\n",
    "    \"inosine 5'-monophosphate (IMP)\",\n",
    "    # 'lamotrigine',  # excluded during feature processing (all values were the same: stddev = 0.0)\n",
    "    'linolenate [alpha or gamma; (18:3n3 or 6)]', 'maltopentaose',\n",
    "    'mannonate', 'meglumine', 'N-acetyl-cadaverine',\n",
    "    'N-acetylmethionine sulfoxide', 'N-acetylserine',\n",
    "    'N-alpha-acetylornithine', 'oleate/vaccenate (18:1)',\n",
    "    'phenylalanylglycine', 'pipecolate', 'pyroglutamine',\n",
    "    'sphingomyelin (d18:1/14:0, d16:1/16:0)',\n",
    "    'sphingomyelin (d18:1/22:1, d18:2/22:0, d16:1/24:1)',\n",
    "    'sphingomyelin (d18:2/16:0, d18:1/16:1)',\n",
    "    'tricosanoyl sphingomyelin (d18:1/23:0)', 'tryptamine',\n",
    "    'eicosenoate (20:1)', 'salicylate',\n",
    "    'N6-methyllysine',\n",
    "    'dihomo-linolenate (20:3n3 or n6)',\n",
    "    'uridine',\n",
    "    'arachidonate (20:4n6)',\n",
    "    'dehydroepiandrosterone sulfate (DHEA-S)',\n",
    "    '2-hydroxypalmitate',\n",
    "    '10-nonadecenoate (19:1n9)',\n",
    "    '3-carboxy-4-methyl-5-propyl-2-furanpropanoate (CMPF)',\n",
    "    'phenylpyruvate']))\n",
    "\n",
    "# train and eval\n",
    "train_n_eval_regressors(ls_targets_metabolites,\n",
    "                        'metabolites',\n",
    "                        transform_target2log=True,\n",
    "                        ls_features=ls_features_metabolites,\n",
    "                        df_data=df_data_w_targets_sel,\n",
    "                        taxa=taxa,\n",
    "                        seed=seed,\n",
    "                        output_dir=output_dir)\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict `metabolites` removing cancer cases"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "df_data_wo_cancer = df_data_w_targets_sel[df_data_w_targets_sel['T_disease_state'] != 'ICC'].copy(\n",
    "    deep=True)\n",
    "\n",
    "# same features as before\n",
    "# train and eval\n",
    "train_n_eval_regressors(ls_targets_metabolites,\n",
    "                        'metabolites_no_cancer',\n",
    "                        transform_target2log=True,\n",
    "                        ls_features=ls_features_metabolites,\n",
    "                        df_data=df_data_wo_cancer,\n",
    "                        taxa=taxa,\n",
    "                        seed=seed,\n",
    "                        output_dir=output_dir)\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict `biomarkers`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# features\n",
    "ls_features_biomarkers = [x for x in df_data_w_targets_sel.columns if\n",
    "                           ((x.startswith('F_micro_') or x.startswith('F_metabo_')))]\n",
    "\n",
    "# targets were defined at the start when reading in the list ls_cancer_biomarkers\n",
    "\n",
    "# train and eval\n",
    "train_n_eval_regressors(ls_cancer_biomarkers,\n",
    "                        'biomarkers',\n",
    "                        transform_target2log=True,\n",
    "                        ls_features=ls_features_biomarkers,\n",
    "                        df_data=df_data_w_targets_sel,\n",
    "                        taxa=taxa,\n",
    "                        seed=seed,\n",
    "                        output_dir=output_dir)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict `biomarkers` removing cancer cases"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# print()\n",
    "train_n_eval_regressors(ls_cancer_biomarkers,\n",
    "                        'biomarkers_no_cancer',\n",
    "                        transform_target2log=True,\n",
    "                        ls_features=ls_features_biomarkers,\n",
    "                        df_data=df_data_wo_cancer,\n",
    "                        taxa=taxa,\n",
    "                        seed=seed,\n",
    "                        output_dir=output_dir)\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict actual value of Vaginal pH (`F_pcov_pH`)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "ls_features4ph = [x for x in df_data_w_targets_sel.columns if (\n",
    "    x.startswith('F_') and not x.startswith('F_pcov_'))]\n",
    "ls_target_num_pH = ['F_pcov_pH']\n",
    "\n",
    "train_n_eval_regressors(ls_target_num_pH,\n",
    "                        'num_pH',\n",
    "                        True,\n",
    "                        ls_features4ph,\n",
    "                        df_data_w_targets_sel,\n",
    "                        taxa,\n",
    "                        seed,\n",
    "                        output_dir)\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# Publication: Supplementary figure 10: numeric vaginal pH distribution\n",
    "fig, ax = plt.subplots(1, figsize=(2, 3))\n",
    "g = df_data_w_targets_sel['F_pcov_pH'].hist(ax=ax, bins=7)\n",
    "ax.set_xlabel('pH')\n",
    "ax.set_ylabel('Sample Count')\n",
    "#ax.set_xticklabels([4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5])\n",
    "fig.savefig(os.path.join(output_dir, 'regressors',\n",
    "            'num_pH_histogram.jpeg'), bbox_inches='tight')\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict numeric inflammation score (`T_infl_score_flt`)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "ls_target_num_inf = ['T_infl_score_flt']\n",
    "train_n_eval_regressors(ls_target_num_inf,\n",
    "                        'num_inflammation', False,\n",
    "                        ls_features4inflammation,\n",
    "                        df_data_w_targets_sel,\n",
    "                        taxa,\n",
    "                        seed,\n",
    "                        output_dir)\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# Publication: Supplementary Figure 12: Inflammation score distribution\n",
    "\n",
    "fig, ax = plt.subplots(1, figsize=(2, 3))\n",
    "g = df_data_w_targets_sel['T_infl_score_flt'].hist(ax=ax, bins=8)\n",
    "ax.set_xlabel('Inflammation Score')\n",
    "ax.set_ylabel('Sample Count')\n",
    "fig.savefig(os.path.join(output_dir, 'regressors',\n",
    "            'num_inflammation_histogram.jpeg'), bbox_inches='tight')\n",
    "plt.show()\n"
   ],
   "outputs": [],
   "metadata": {}
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "762ae4eb0dbf910a4d6be6e2c7a1c874bbda83ad0aae04f1ceabd5b2d7623b1a"
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.8.8 64-bit ('hpv_modelling': conda)"
  },
  "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.8.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "491.5px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}