demo_scVI_cortex.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Application of scVI on Cortex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import torch\n",
    "import os\n",
    "import time\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from scvi.dataset import CortexDataset, RetinaDataset\n",
    "from scvi.models import VAE\n",
    "from scvi.inference import UnsupervisedTrainer, load_posterior\n",
    "from scvi import set_seed\n",
    "from scvi.dataset.csv import CsvDataset\n",
    "from utils import *\n",
    "np.random.seed(1234)\n",
    "torch.manual_seed(1234)\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = '2'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data and fit scVI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = \"/home/jzhaoaz/jiazhao/scPI_v2/package/datasets/cortex/\"\n",
    "expression_train = np.loadtxt(data_path + \"data_train\", dtype='float32')\n",
    "label_train = np.loadtxt(data_path + \"label_train\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "save expression_train as .csv file named \"expression_train.csv\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2020-07-15 13:24:49,689] INFO - scvi.dataset.csv | Preprocessing dataset\n",
      "[2020-07-15 13:24:49,916] INFO - scvi.dataset.csv | Finished preprocessing dataset\n",
      "[2020-07-15 13:24:49,922] INFO - scvi.dataset.dataset | Remapping labels to [0,N]\n",
      "[2020-07-15 13:24:49,923] INFO - scvi.dataset.dataset | Remapping batch_indices to [0,N]\n",
      "[2020-07-15 13:24:49,951] INFO - scvi.dataset.dataset | Computing the library size for the new data\n",
      "[2020-07-15 13:24:49,953] INFO - scvi.dataset.dataset | Downsampled from 2253 to 2253 cells\n"
     ]
    }
   ],
   "source": [
    "exp_train = CsvDataset(filename = data_path + \"expression_train.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2020-07-15 13:24:52,904] INFO - scvi.inference.inference | KL warmup for 400 epochs\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "23785f21310c45c9a7c158f8c89ef317",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='training', max=400.0, style=ProgressStyle(description_wid…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Computational time = 1.9 minutes.\n"
     ]
    }
   ],
   "source": [
    "n_epochs = 400\n",
    "lr = 1e-3\n",
    "use_cuda = True\n",
    "n_celltype = 7\n",
    "t_start = time.time()\n",
    "vae = VAE(exp_train.nb_genes)\n",
    "trainer = UnsupervisedTrainer(\n",
    "    vae,\n",
    "    exp_train,\n",
    "    train_size=0.90,\n",
    "    use_cuda=use_cuda,\n",
    "    frequency=5,\n",
    ")\n",
    "trainer.train(n_epochs=n_epochs, lr=lr)\n",
    "full = trainer.create_posterior(trainer.model, exp_train, indices=np.arange(len(exp_train)))\n",
    "latent, batch_indices, labels = full.sequential().get_latent()\n",
    "t_end = time.time()\n",
    "print('Computational time = %1.1f minutes.' % ((t_end-t_start)/60.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate clustering scores (ASW, NMI, ARI) and visualize latent distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Clustering scores [ASW, NMI, ARI] = [0.236863, 0.767987, 0.743195].\n"
     ]
    }
   ],
   "source": [
    "scores = cluster_scores(latent, 7, label_train)\n",
    "print('Clustering scores [ASW, NMI, ARI] = [%f, %f, %f].' % (scores[0], scores[1], scores[2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imputation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = \"/home/jzhaoaz/jiazhao/scPI_v2/package/datasets/cortex/\"\n",
    "X_true = np.loadtxt(data_path + \"data_train\", dtype='float32')\n",
    "X_zero = np.loadtxt(data_path + \"X_zero.txt\", dtype='float32')\n",
    "i = np.load(data_path + \"i.npy\")\n",
    "j = np.load(data_path + \"j.npy\")\n",
    "ix = np.load(data_path + \"ix.npy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "save X_zero as .csv file named \"X_zero.csv\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2020-07-15 13:26:45,877] INFO - scvi.dataset.csv | Preprocessing dataset\n",
      "[2020-07-15 13:26:46,119] INFO - scvi.dataset.csv | Finished preprocessing dataset\n",
      "[2020-07-15 13:26:46,123] INFO - scvi.dataset.dataset | Remapping labels to [0,N]\n",
      "[2020-07-15 13:26:46,125] INFO - scvi.dataset.dataset | Remapping batch_indices to [0,N]\n",
      "[2020-07-15 13:26:46,153] INFO - scvi.dataset.dataset | Computing the library size for the new data\n",
      "[2020-07-15 13:26:46,156] INFO - scvi.dataset.dataset | Downsampled from 2253 to 2253 cells\n"
     ]
    }
   ],
   "source": [
    "data_path = \"/home/jzhaoaz/jiazhao/scPI_v2/package/datasets/cortex/\"\n",
    "exp_train = CsvDataset(filename = data_path + \"X_zero.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2020-07-15 13:26:46,225] INFO - scvi.inference.inference | KL warmup for 400 epochs\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a8deb71ed6d8402bb61b621c13adb59e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='training', max=400.0, style=ProgressStyle(description_wid…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Computational time = 1.9 minutes.\n"
     ]
    }
   ],
   "source": [
    "t_start = time.time()\n",
    "vae = VAE(exp_train.nb_genes)\n",
    "trainer = UnsupervisedTrainer(\n",
    "    vae,\n",
    "    exp_train,\n",
    "    train_size=0.90,\n",
    "    use_cuda=use_cuda,\n",
    "    frequency=5,\n",
    ")\n",
    "trainer.train(n_epochs=n_epochs, lr=lr)\n",
    "full = trainer.create_posterior(trainer.model, exp_train, indices=np.arange(len(exp_train)))\n",
    "latent, batch_indices, labels = full.sequential().get_latent()\n",
    "t_end = time.time()\n",
    "print('Computational time = %1.1f minutes.' % ((t_end-t_start)/60.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate imputation error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imputation error on randomly dropout data = 2.328291.\n"
     ]
    }
   ],
   "source": [
    "imputed_values = full.sequential().imputation()\n",
    "impu_err = imputation_error(imputed_values, X_true, X_zero, i, j, ix)\n",
    "print('Imputation error on randomly dropout data = %f.' % impu_err)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}