07_compute_retina.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import numpy as np\n",
    "\n",
    "from tools import kobak_tsne, PCA_sklearn, pearson_residuals, prepare_largest_batch, run_glmpca, sqrt_lazy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute normalization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Of 19252 total genes, returning 17973 genes that are detected in 5 or more cells.\n",
      "Output shape: (24769, 17973)\n",
      "Of 18149 total genes, returning 16520 genes that are detected in 5 or more cells.\n",
      "Output shape: (13987, 16520)\n",
      "Of 18137 total genes, returning 17685 genes that are detected in 5 or more cells.\n",
      "Output shape: (15750, 17685)\n"
     ]
    }
   ],
   "source": [
    "files =['datasets/retina/macosko_all/GSE63472_P14Retina_merged_digital_expression.txt_preprocessed.pickle',\n",
    "        'datasets/retina/shekhar_bipolar/GSE81904_BipolarUMICounts_Cell2016.txt.gz_preprocessed.pickle',\n",
    "        'datasets/retina/tran_ganglion/GSE133382_AtlasRGCs_CountMatrix.csv_preprocessed.pickle']\n",
    "\n",
    "labels = ['macosko',\n",
    "          'shekhar_bipolar',\n",
    "          'tran_ganglion']\n",
    "\n",
    "datasets=[]\n",
    "for file,label in zip(files,labels):\n",
    "    with open(file, \"rb\") as f:\n",
    "        dataset = pickle.load(f)\n",
    "        dataset['labelshort'] = label\n",
    "    \n",
    "    \n",
    "    if label in ['macosko','shekhar_bipolar']:\n",
    "        prepare_largest_batch(dataset,extra_keys=['replicates'])\n",
    "    else:\n",
    "        prepare_largest_batch(dataset)\n",
    "    \n",
    "    n_clusters = len(np.unique(dataset['clusters']))\n",
    "    dataset['gene_means'] = np.array(np.mean(dataset['counts'],axis=0)).flatten()\n",
    "    dataset['depths'] = np.array(np.sum(dataset['counts'],axis=1)).flatten()\n",
    "    \n",
    "    datasets.append(dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analytical Pearson residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 20.5 s, sys: 4.32 s, total: 24.8 s\n",
      "Wall time: 24.8 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for dataset in datasets:\n",
    "    dataset['residuals_theta100'] = pearson_residuals(dataset['counts'].toarray(),100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sqrt transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.82 s, sys: 1.07 s, total: 7.89 s\n",
      "Wall time: 7.89 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for dataset in datasets:\n",
    "    dataset['sqrt_lazy'] = sqrt_lazy(dataset['counts'].toarray())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GLM PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "for dataset in datasets:\n",
    "    run_glmpca(dataset['counts'].toarray(),fam='nb',optimize_nb_theta=False,theta=100,dataset_label=dataset['labelshort'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "for dataset in datasets:\n",
    "    with open('glmpca_results/%sglmpca-py_nb_fixedTheta100_penalty1.pickle' % (dataset['labelshort']),\"rb\") as f:\n",
    "        dataset['glmpca-py_nb_thetaFixed100'] = pickle.load(f)['factors']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute tSNEs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "for dataset in datasets:\n",
    "\n",
    "    _,sqrt_lazy_init = PCA_sklearn(dataset['sqrt_lazy'],50,42)\n",
    "    _ = kobak_tsne(dataset['sqrt_lazy'],name=dataset['labelshort']+'_sqrt_lazy_initsqrtLazy',do_pca=True,init=sqrt_lazy_init)\n",
    "    _ = kobak_tsne(dataset['residuals_theta100'],name=dataset['labelshort']+'_residuals_theta100_initsqrtLazy',do_pca=True,init=sqrt_lazy_init)\n",
    "    _ = kobak_tsne(dataset['glmpca-py_nb_thetaFixed100'],name=dataset['labelshort']+'_glmpca-py_nb_thetaFixed100_initsqrtLazy',do_pca=False,init=sqrt_lazy_init)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext watermark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Last updated: 2021-07-31T19:14:03.340668+02:00\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.8.0\n",
      "IPython version      : 7.21.0\n",
      "\n",
      "Compiler    : GCC 8.3.0\n",
      "OS          : Linux\n",
      "Release     : 3.10.0-957.el7.x86_64\n",
      "Machine     : x86_64\n",
      "Processor   : x86_64\n",
      "CPU cores   : 40\n",
      "Architecture: 64bit\n",
      "\n"
     ]
    }
   ],
   "source": [
    "watermark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "numpy: 1.20.1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "watermark --iversions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}