ppcs_pbmc3k.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%matplotlib inline\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.chdir(\"/home/romain/PycharmProjects/scVI-harmo/\")\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scvi\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "mkdir: cannot create directory ‘data’: File exists\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "\r",
      "  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0\r",
      "  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0\r",
      "  0 11.0M    0     0    0     0      0      0 --:--:--  0:00:01 --:--:--     0\r",
      "  2 11.0M    2  330k    0     0   155k      0  0:01:13  0:00:02  0:01:11  155k\r",
      " 22 11.0M   22 2515k    0     0   804k      0  0:00:14  0:00:03  0:00:11  804k\r",
      " 52 11.0M   52 5978k    0     0  1448k      0  0:00:07  0:00:04  0:00:03 1448k\r",
      " 99 11.0M   99 11.0M    0     0  2203k      0  0:00:05  0:00:05 --:--:-- 2305k\r",
      "100 11.0M  100 11.0M    0     0  2207k      0  0:00:05  0:00:05 --:--:-- 2877k\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "\r",
      "  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0\r",
      "100 11071  100 11071    0     0  45862      0 --:--:-- --:--:-- --:--:-- 45937\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "\r",
      "  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0\r",
      " 14 20.4M   14 3111k    0     0  4704k      0  0:00:04 --:--:--  0:00:04 4700k\r",
      "100 20.4M  100 20.4M    0     0  13.6M      0  0:00:01  0:00:01 --:--:-- 13.6M\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "\r",
      "  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0\r",
      "  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0\r",
      " 73 15.0M   73 11.0M    0     0   9.7M      0  0:00:01  0:00:01 --:--:--  9.7M\r",
      "100 15.0M  100 15.0M    0     0  11.4M      0  0:00:01  0:00:01 --:--:-- 11.4M\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "\r",
      "  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0\r",
      " 41 14.9M   41 6375k    0     0  8206k      0  0:00:01 --:--:--  0:00:01 8205k\r",
      "100 14.9M  100 14.9M    0     0  12.5M      0  0:00:01  0:00:01 --:--:-- 12.5M\n"
     ]
    }
   ],
   "source": [
    "%%bash \n",
    "mkdir data\n",
    "curl https://ndownloader.figshare.com/files/18789473?private_link=6e12f5f56bb46842f5ee --output data/pbmc3k.h5ad \n",
    "curl https://raw.githubusercontent.com/normjam/ppcs/master/genes_1000_pbmc3k.csv --output data/genes_1000_pbmc3k.csv\n",
    "curl https://raw.githubusercontent.com/normjam/ppcs/master/pbmc3k_generated_sctransform.h5 --output data/pbmc3k_generated_sctransform.h5\n",
    "curl https://raw.githubusercontent.com/normjam/ppcs/master/pbmc3k_1kgenes_glmpca10_poi_samples.npz --output data/pbmc3k_1kgenes_glmpca10_poi_samples.npz\n",
    "curl https://raw.githubusercontent.com/normjam/ppcs/master/pbmc3k_1kgenes_glmpca10_nb_samples.npz --output data/pbmc3k_1kgenes_glmpca10_nb_samples.npz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import anndata\n",
    "anndataset = anndata.read(\"data/pbmc3k.h5ad\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2019-11-20 10:22:27,487] INFO - scvi.dataset.anndataset | Dense size under 1Gb, casting to dense format (np.ndarray).\n",
      "[2019-11-20 10:22:27,633] INFO - scvi.dataset.dataset | Remapping batch_indices to [0,N]\n",
      "[2019-11-20 10:22:27,635] INFO - scvi.dataset.dataset | Remapping labels to [0,N]\n",
      "[2019-11-20 10:22:27,704] INFO - scvi.dataset.dataset | Computing the library size for the new data\n",
      "[2019-11-20 10:22:27,779] INFO - scvi.dataset.dataset | Downsampled from 2695 to 2695 cells\n"
     ]
    }
   ],
   "source": [
    "from scvi.dataset import AnnDatasetFromAnnData\n",
    "dataset = AnnDatasetFromAnnData(ad = anndataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2019-11-20 10:22:27,880] INFO - scvi.dataset.dataset | Computing the library size for the new data\n"
     ]
    }
   ],
   "source": [
    "select_genes = np.array([x.split(\",\")[1].upper() for x in np.loadtxt(\"data/genes_1000_pbmc3k.csv\", dtype=np.str)])\n",
    "dataset.reorder_genes(select_genes, drop_omitted_genes=True)\n",
    "assert(np.all(select_genes == dataset.gene_names))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "193296f2507c49538c63b1a4f2e2b08e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='training', max=200, style=ProgressStyle(description_width='in…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a7d72b63881e4a29ae101a887f90899d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='training', max=200, style=ProgressStyle(description_width='in…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bd41b3bd634f4781888d78cbedf36cc6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='training', max=200, style=ProgressStyle(description_width='in…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# Model training\n",
    "from scvi.models import VAE\n",
    "from scvi.inference import UnsupervisedTrainer\n",
    "\n",
    "n_epochs = 200\n",
    "warmup = 100\n",
    "n_latent = 10\n",
    "dict_post = {}\n",
    "\n",
    "# ZINB scVI\n",
    "for reconstruction_loss in [\"zinb\", \"nb\", \"poisson\"]:\n",
    "    model = VAE(n_input=dataset.nb_genes, n_latent=n_latent, n_batch=0, reconstruction_loss=reconstruction_loss)\n",
    "    trainer = UnsupervisedTrainer(\n",
    "        model=model, \n",
    "        gene_dataset=dataset, \n",
    "        n_epochs_kl_warmup=warmup,\n",
    "        frequency=1,\n",
    "        train_size=1.0,\n",
    "    )    \n",
    "    trainer.train(n_epochs=n_epochs, lr=1e-2)\n",
    "    dict_post[reconstruction_loss] = trainer.create_posterior(trainer.model, dataset, \n",
    "                                                              indices=np.arange(len(dataset))).sequential()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PPC for the entire dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scvi.inference import PosteriorPredictiveCheck"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppc = PosteriorPredictiveCheck(dict_post, n_samples=50)\n",
    "# add factor analysis\n",
    "ppc.store_fa_samples(n_components=30, normalization=\"log\")\n",
    "ppc.store_fa_samples(n_components=30, normalization=\"log_rate\")\n",
    "ppc.store_fa_samples(n_components=30, normalization=\"rate\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/romain/PycharmProjects/scVI-harmo/scvi/inference/ppc.py:59: RuntimeWarning: invalid value encountered in true_divide\n",
      "  np.std(samples, axis=axis) / np.mean(samples, axis=axis), axis=-1\n",
      "/home/romain/PycharmProjects/scVI-harmo/scvi/inference/ppc.py:65: RuntimeWarning: invalid value encountered in true_divide\n",
      "  self.raw_counts, axis=axis\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b832a08e4bba49a49d7d4cc135d9be07",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=6), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ppc.coeff_of_variation(cell_wise=True)\n",
    "ppc.coeff_of_variation(cell_wise=False)\n",
    "ppc.mean(cell_wise=False)\n",
    "ppc.variance(cell_wise=False)\n",
    "ppc.dropout_ratio()\n",
    "ppc.gene_gene_correlation(gene_indices=np.arange(1000))\n",
    "ppc.calibration_error([0, 100])\n",
    "ppc.mean_squared_log_error()\n",
    "ppc.median_absolute_error()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "# open a file, where you ant to store the data\n",
    "file = open('results_scVI_and_FA.pickle', 'wb')\n",
    "\n",
    "# dump information to that file\n",
    "pickle.dump(ppc.metrics, file)\n",
    "\n",
    "# close the file\n",
    "file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppc = PosteriorPredictiveCheck(dict_post, n_samples=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import numpy as np \n",
    "\n",
    "scTransform = h5py.File(\"data/pbmc3k_generated_sctransform.h5\", 'r')  \n",
    "scTransform = np.array(scTransform['x'])\n",
    "scTransform = scTransform.swapaxes(0, 2)\n",
    "ppc.store_external_samples(scTransform, \"sctransform\")\n",
    "\n",
    "glm_pca_poisson = np.load(\"data/pbmc3k_1kgenes_glmpca10_poi_samples.npz\")\n",
    "glm_pca_poisson = list(glm_pca_poisson.values())[0].astype(np.float64)\n",
    "ppc.store_external_samples(glm_pca_poisson, \"glm_pca_poisson\")\n",
    "\n",
    "glm_pca_nb = np.load(\"data/pbmc3k_1kgenes_glmpca10_nb_samples.npz\")\n",
    "glm_pca_nb = list(glm_pca_nb.values())[0].astype(np.float64)\n",
    "ppc.store_external_samples(glm_pca_nb, \"glm_pca_nb\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/romain/PycharmProjects/scVI-harmo/scvi/inference/ppc.py:59: RuntimeWarning: invalid value encountered in true_divide\n",
      "  np.std(samples, axis=axis) / np.mean(samples, axis=axis), axis=-1\n",
      "/home/romain/PycharmProjects/scVI-harmo/scvi/inference/ppc.py:59: RuntimeWarning: Mean of empty slice\n",
      "  np.std(samples, axis=axis) / np.mean(samples, axis=axis), axis=-1\n",
      "/home/romain/PycharmProjects/scVI-harmo/scvi/inference/ppc.py:65: RuntimeWarning: invalid value encountered in true_divide\n",
      "  self.raw_counts, axis=axis\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9c0781a782354b8d823418f640cf33eb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=6), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ppc.coeff_of_variation(cell_wise=True)\n",
    "ppc.coeff_of_variation(cell_wise=False)\n",
    "ppc.mean(cell_wise=False)\n",
    "ppc.variance(cell_wise=False)\n",
    "ppc.dropout_ratio()\n",
    "ppc.gene_gene_correlation(gene_indices=np.arange(1000))\n",
    "ppc.calibration_error([0, 100])\n",
    "ppc.mean_squared_log_error()\n",
    "ppc.median_absolute_error()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "# open a file, where you ant to store the data\n",
    "file = open('results_scVI_and_others.pickle', 'wb')\n",
    "\n",
    "# dump information to that file\n",
    "pickle.dump(ppc.metrics, file)\n",
    "\n",
    "# close the file\n",
    "file.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py37",
   "language": "python",
   "name": "py37"
  },
  "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}