incubator/microbiome-case-study.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import war\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis\n",
    "\n",
    "Microbiome composition can be shifted by diet. The experiments that have backed this conclusion involve randomly assigning mice into control and case groups, feeding them different diets, and then quantifying the composition of the microbiota in the gut. I want to know whether we are able to use a Bayesian Dirichlet-Multinomial model to quantify the uncertainty surrounding the measured proportions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "microbiome = pd.read_csv(\"../datasets/MicrobiomeWithMetadata.csv\")\n",
    "microbiome.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The metadata file that's associated with this CSV file has to be re-coded from a CSV file to a YAML file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "\n",
    "md = defaultdict(dict)  # \"md\" stands for \"metadata dictionary\"\n",
    "\n",
    "sex = [\"Male\", \"Female\"]\n",
    "donor = [\"HMouseLFPP\", \"CONVR\", \"Human\", \"Fresh\", \"Frozen\", \"HMouseWestern\", \"CONVD\"]\n",
    "diet = [\"LFPP\", \"Western\", \"CARBR\", \"FATR\", \"Suckling\", \"Human\"]\n",
    "source = [\n",
    "    \"Cecum1\",\n",
    "    \"Cecum2\",\n",
    "    \"Colon1\",\n",
    "    \"Colon2\",\n",
    "    \"Feces\",\n",
    "    \"SI1\",\n",
    "    \"SI13\",\n",
    "    \"SI15\",\n",
    "    \"SI2\",\n",
    "    \"SI5\",\n",
    "    \"SI9\",\n",
    "    \"Stomach\",\n",
    "    \"Cecum\",\n",
    "]\n",
    "collection_met = [\"Contents\", \"Scraping\"]\n",
    "\n",
    "for i, s in enumerate(sex):\n",
    "    md[\"sex\"][i] = s\n",
    "\n",
    "for i, d in enumerate(donor):\n",
    "    md[\"donor\"][i] = d\n",
    "\n",
    "for i, d in enumerate(diet):\n",
    "    md[\"diet\"][i] = d\n",
    "\n",
    "for i, s in enumerate(source):\n",
    "    md[\"source\"][i] = s\n",
    "\n",
    "for i, c in enumerate(collection_met):\n",
    "    md[\"collection_met\"][i] = c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import yaml\n",
    "\n",
    "print(yaml.dump(md))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"datasets/MicrobiomeMetadataDictionary.yml\", \"w+\") as f:\n",
    "    f.write(yaml.dump(md))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"datasets/MicrobiomeMetadataDictionary.yml\", \"r+\") as f:\n",
    "    metadata = yaml.load(f)\n",
    "\n",
    "metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(microbiome[\"Diet\"].values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(microbiome[\"Source\"].values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "otu_cols = [c for c in microbiome.columns if \"OTU\" in c]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as dirichlet_model:\n",
    "    mu = pm.HalfNormal(\"mu\", sd=100 ** 2)\n",
    "    n_seq_reads = pm.Poisson(\"n_seq_reads\", mu=mu, observed=healthy_reads.sum(axis=1))\n",
    "    proportions = pm.Dirichlet(\"proportions\", a=np.ones(3), shape=(3,))\n",
    "    for i in range(healthy_reads.shape[0]):\n",
    "        draws = pm.Multinomial(\n",
    "            f\"draws_{i}\",\n",
    "            n=healthy_reads[i].sum(),\n",
    "            p=proportions,\n",
    "            observed=healthy_reads[i, :],\n",
    "        )\n",
    "    dirichlet_trace = pm.sample(draws=2000)\n",
    "    pm.traceplot(dirichlet_trace)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bayesian",
   "language": "python",
   "name": "bayesian"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}