b-mmvec-HPV.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Microbiome-metabolome interaction probabilities via [mmvec](https://github.com/biocore/mmvec) for 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 the application of mmvec to derive the microbiome-metabolome interaction probabilities:             \n",
    "\n",
    "[1) Import packages and set paths](#chap1)      \n",
    "\n",
    "[2) Read datasets](#chap2)      \n",
    "\n",
    "[3) Plot mmvec biplot](#chap3)      \n",
    "\n",
    "[4) Plot mmvec heatmap](#chap4)            \n",
    "\n",
    "      \n",
    "**Beware: This notebook runs in a separate conda environment than the notebook `a-modelling-HPV.ipynb`.** For details on how to setup this environment please refer to the `README.md` file.\n",
    "         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='chap1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Import packages and set paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import os\n",
    "import pandas as pd\n",
    "import qiime2\n",
    "from qiime2.plugins import taxa, mmvec\n",
    "from matplotlib.pylab import plt\n",
    "import skbio\n",
    "from util_mmvec import biplot_from_rhapsody\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "# location to retrieve raw-data\n",
    "path2data = 'data-raw'\n",
    "\n",
    "# location to save outputs\n",
    "output_dir = 'b-mmvec-output'\n",
    "if not os.path.exists(output_dir):\n",
    "    os.makedirs(output_dir)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='chap2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) Read datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read metadata of metabolites\n",
    "metabolite_md = qiime2.Metadata.load(\n",
    "    os.path.join(path2data, 'metabolite_metadata_CC.txt'))\n",
    "\n",
    "# load conditionals created with mmvec\n",
    "df_conditionals = qiime2.Artifact.load(os.path.join(\n",
    "    path2data, 'metab_ranks.qza')).view(pd.DataFrame)\n",
    "\n",
    "# load pcoa created with mmvec\n",
    "pcoa = qiime2.Artifact.load(os.path.join(path2data, 'metab_biplot.qza')).view(\n",
    "    skbio.stats.ordination.OrdinationResults)\n",
    "\n",
    "# load consensus taxonomy\n",
    "path2merged_tax = os.path.join(path2data, 'merged_taxonomy.tsv')\n",
    "merged_taxonomy = pd.read_csv(path2merged_tax, sep='\\t')\n",
    "merged_taxonomy.set_index('Feature ID', inplace=True)\n",
    "taxonomy = merged_taxonomy[['Consensus Taxonomy']].copy(deep=True)\n",
    "taxonomy.columns = ['Taxon']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='chap3'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Plot mmvec biplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metabolite_colors = {\n",
    "    'Lipid': 'crimson',\n",
    "    'Xenobiotics': 'chocolate',\n",
    "    'Peptide': 'orange',\n",
    "    'Cofactors and Vitamins': 'gold',\n",
    "    'Partially Characterized Molecules': 'greenyellow',\n",
    "    'Nucleotide': 'forestgreen',\n",
    "    'Carbohydrate': 'blue',\n",
    "    'Amino Acid': 'indigo',\n",
    "    'Energy': 'violet',\n",
    "}\n",
    "fig = biplot_from_rhapsody(pcoa, taxonomy, metabolite_md.to_dataframe()['SUPER PATHWAY'],\n",
    "                           metabolite_colors, palette=\"tab20\", level=5,\n",
    "                           minradius=.1, label_level=6, feature_labels=None, mincount=3, whitelist=None)\n",
    "axes = plt.gca()\n",
    "plt.xlim(-0.6, 0.65)\n",
    "plt.ylim(-0.6, 0.6)\n",
    "fig.savefig(os.path.join(output_dir, 'rhapsody-biplot.pdf'),\n",
    "            bbox_inches=\"tight\")\n",
    "fig.savefig(os.path.join(output_dir, 'rhapsody-biplot.png'),\n",
    "            bbox_inches=\"tight\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='chap4'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4) Plot mmvec heatmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# process conditionals\n",
    "df_conditionals = df_conditionals[df_conditionals.max(axis=1) >= 1.2]\n",
    "df_conditionals = df_conditionals[df_conditionals.columns[(\n",
    "    df_conditionals >= 2).any()]]\n",
    "\n",
    "# align tax and conditionals\n",
    "taxonomy, df_conditionals = taxonomy[['Taxon']].align(\n",
    "    df_conditionals, join='inner', axis=0)\n",
    "# grab genus + species name (or family if genus is missing)\n",
    "taxa = taxonomy['Taxon'].apply(lambda x: x.split(\n",
    "    ';')[-1] if len(x.split(';')[-1]) > 0 else x.split(';')[-2:][0] + ' unknown')\n",
    "# create new feature IDs\n",
    "taxa = [i[:6] + ' : ' + t for i, t in zip(taxa.index, taxa)]\n",
    "\n",
    "# reset indices\n",
    "taxonomy.index = pd.Index(taxa, name='id')\n",
    "df_conditionals.index = pd.Index(taxa, name='id')\n",
    "\n",
    "# transform tax and conditionals to Q2 types\n",
    "md_taxonomy = qiime2.Metadata(taxonomy)\n",
    "art_conditionals_t = qiime2.Artifact.import_data(\n",
    "    'FeatureData[Conditional]', df_conditionals.T)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# todo adjust Taxon to Genus level\n",
    "hotmap = mmvec.actions.heatmap(ranks=art_conditionals_t,\n",
    "                               microbe_metadata=md_taxonomy.get_column(\n",
    "                                   'Taxon'),\n",
    "                               metabolite_metadata=metabolite_md.get_column(\n",
    "                                   'SUPER PATHWAY'),\n",
    "                               level=6,\n",
    "                               x_labels=True,\n",
    "                               y_labels=True,\n",
    "                               margin_palette='Paired')\n",
    "\n",
    "hotmap.visualization.save(os.path.join(output_dir, 'filtered-heatmap.qzv'))\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "317d0eebc8a73ac18168396754aba1aa263ed9e141510eb0cf15aa05053d1b18"
  },
  "kernelspec": {
   "display_name": "Python 3.6.10 64-bit ('hpv_mmvec': conda)",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": ""
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}