scEU_seq_organoid_analysis_kinetic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# scEU-seq organoid\n",
    "\n",
    "This tutorial uses the intestine organoid data from Battich, et al (2020). This tutorial is the second one of the two tutorials for demonstrating how dynamo can use used to analyze the scEU-seq data. Please refer the [cell cycle](https://dynamo-release.readthedocs.io/en/latest/scEU_seq_rpe1_analysis_kinetic.html) tutorial for details on how to analyze the cell cycle dataset.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>package</th>\n",
       "      <th>dynamo-release</th>\n",
       "      <th>pre-commit</th>\n",
       "      <th>cvxopt</th>\n",
       "      <th>trimap</th>\n",
       "      <th>numdifftools</th>\n",
       "      <th>colorcet</th>\n",
       "      <th>python-igraph</th>\n",
       "      <th>pynndescent</th>\n",
       "      <th>hdbscan</th>\n",
       "      <th>loompy</th>\n",
       "      <th>matplotlib</th>\n",
       "      <th>networkx</th>\n",
       "      <th>numba</th>\n",
       "      <th>numpy</th>\n",
       "      <th>pandas</th>\n",
       "      <th>scikit-learn</th>\n",
       "      <th>scipy</th>\n",
       "      <th>seaborn</th>\n",
       "      <th>setuptools</th>\n",
       "      <th>statsmodels</th>\n",
       "      <th>tqdm</th>\n",
       "      <th>umap-learn</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>version</th>\n",
       "      <td>1.0.0</td>\n",
       "      <td>2.11.1</td>\n",
       "      <td>1.2.6</td>\n",
       "      <td>1.4.3.dev1</td>\n",
       "      <td>0.9.39</td>\n",
       "      <td>2.0.6</td>\n",
       "      <td>0.9.0</td>\n",
       "      <td>0.5.2</td>\n",
       "      <td>0.8.27</td>\n",
       "      <td>3.0.6</td>\n",
       "      <td>3.4.3</td>\n",
       "      <td>2.5.1</td>\n",
       "      <td>0.53.1</td>\n",
       "      <td>1.19.5</td>\n",
       "      <td>1.3.3</td>\n",
       "      <td>0.23.2</td>\n",
       "      <td>1.7.0</td>\n",
       "      <td>0.11.1</td>\n",
       "      <td>54.2.0</td>\n",
       "      <td>0.12.2</td>\n",
       "      <td>4.58.0</td>\n",
       "      <td>0.5.1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "package dynamo-release pre-commit cvxopt      trimap numdifftools colorcet  \\\n",
       "version          1.0.0     2.11.1  1.2.6  1.4.3.dev1       0.9.39    2.0.6   \n",
       "\n",
       "package python-igraph pynndescent hdbscan loompy matplotlib networkx   numba  \\\n",
       "version         0.9.0       0.5.2  0.8.27  3.0.6      3.4.3    2.5.1  0.53.1   \n",
       "\n",
       "package   numpy pandas scikit-learn  scipy seaborn setuptools statsmodels  \\\n",
       "version  1.19.5  1.3.3       0.23.2  1.7.0  0.11.1     54.2.0      0.12.2   \n",
       "\n",
       "package    tqdm umap-learn  \n",
       "version  4.58.0      0.5.1  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import dynamo as dyn\n",
    "import anndata\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.sparse\n",
    "\n",
    "from anndata import AnnData\n",
    "from scipy.sparse import csr_matrix\n",
    "\n",
    "dyn.get_all_dependencies_version()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "|-----> Downloading scEU_seq data\n",
      "|-----> Downloading data to ./data/organoid.h5ad\n",
      "|-----> in progress: 5.5272%IOPub message rate exceeded.\n",
      "The notebook server will temporarily stop sending output\n",
      "to the client in order to avoid crashing it.\n",
      "To change this limit, set the config variable\n",
      "`--NotebookApp.iopub_msg_rate_limit`.\n",
      "\n",
      "Current values:\n",
      "NotebookApp.iopub_msg_rate_limit=1000.0 (msgs/sec)\n",
      "NotebookApp.rate_limit_window=3.0 (secs)\n",
      "\n",
      "|-----> in progress: 11.5010%IOPub message rate exceeded.\n",
      "The notebook server will temporarily stop sending output\n",
      "to the client in order to avoid crashing it.\n",
      "To change this limit, set the config variable\n",
      "`--NotebookApp.iopub_msg_rate_limit`.\n",
      "\n",
      "Current values:\n",
      "NotebookApp.iopub_msg_rate_limit=1000.0 (msgs/sec)\n",
      "NotebookApp.rate_limit_window=3.0 (secs)\n",
      "\n",
      "|-----> in progress: 25.0421%IOPub message rate exceeded.\n",
      "The notebook server will temporarily stop sending output\n",
      "to the client in order to avoid crashing it.\n",
      "To change this limit, set the config variable\n",
      "`--NotebookApp.iopub_msg_rate_limit`.\n",
      "\n",
      "Current values:\n",
      "NotebookApp.iopub_msg_rate_limit=1000.0 (msgs/sec)\n",
      "NotebookApp.rate_limit_window=3.0 (secs)\n",
      "\n",
      "|-----> in progress: 31.1626%IOPub message rate exceeded.\n",
      "The notebook server will temporarily stop sending output\n",
      "to the client in order to avoid crashing it.\n",
      "To change this limit, set the config variable\n",
      "`--NotebookApp.iopub_msg_rate_limit`.\n",
      "\n",
      "Current values:\n",
      "NotebookApp.iopub_msg_rate_limit=1000.0 (msgs/sec)\n",
      "NotebookApp.rate_limit_window=3.0 (secs)\n",
      "\n",
      "|-----> in progress: 37.4440%IOPub message rate exceeded.\n",
      "The notebook server will temporarily stop sending output\n",
      "to the client in order to avoid crashing it.\n",
      "To change this limit, set the config variable\n",
      "`--NotebookApp.iopub_msg_rate_limit`.\n",
      "\n",
      "Current values:\n",
      "NotebookApp.iopub_msg_rate_limit=1000.0 (msgs/sec)\n",
      "NotebookApp.rate_limit_window=3.0 (secs)\n",
      "\n",
      "|-----> in progress: 43.6039%IOPub message rate exceeded.\n",
      "The notebook server will temporarily stop sending output\n",
      "to the client in order to avoid crashing it.\n",
      "To change this limit, set the config variable\n",
      "`--NotebookApp.iopub_msg_rate_limit`.\n",
      "\n",
      "Current values:\n",
      "NotebookApp.iopub_msg_rate_limit=1000.0 (msgs/sec)\n",
      "NotebookApp.rate_limit_window=3.0 (secs)\n",
      "\n",
      "|-----> in progress: 49.6913%IOPub message rate exceeded.\n",
      "The notebook server will temporarily stop sending output\n",
      "to the client in order to avoid crashing it.\n",
      "To change this limit, set the config variable\n",
      "`--NotebookApp.iopub_msg_rate_limit`.\n",
      "\n",
      "Current values:\n",
      "NotebookApp.iopub_msg_rate_limit=1000.0 (msgs/sec)\n",
      "NotebookApp.rate_limit_window=3.0 (secs)\n",
      "\n",
      "|-----> in progress: 65.7143%IOPub message rate exceeded.\n",
      "The notebook server will temporarily stop sending output\n",
      "to the client in order to avoid crashing it.\n",
      "To change this limit, set the config variable\n",
      "`--NotebookApp.iopub_msg_rate_limit`.\n",
      "\n",
      "Current values:\n",
      "NotebookApp.iopub_msg_rate_limit=1000.0 (msgs/sec)\n",
      "NotebookApp.rate_limit_window=3.0 (secs)\n",
      "\n",
      "|-----> in progress: 71.9058%IOPub message rate exceeded.\n",
      "The notebook server will temporarily stop sending output\n",
      "to the client in order to avoid crashing it.\n",
      "To change this limit, set the config variable\n",
      "`--NotebookApp.iopub_msg_rate_limit`.\n",
      "\n",
      "Current values:\n",
      "NotebookApp.iopub_msg_rate_limit=1000.0 (msgs/sec)\n",
      "NotebookApp.rate_limit_window=3.0 (secs)\n",
      "\n",
      "|-----> in progress: 78.0137%IOPub message rate exceeded.\n",
      "The notebook server will temporarily stop sending output\n",
      "to the client in order to avoid crashing it.\n",
      "To change this limit, set the config variable\n",
      "`--NotebookApp.iopub_msg_rate_limit`.\n",
      "\n",
      "Current values:\n",
      "NotebookApp.iopub_msg_rate_limit=1000.0 (msgs/sec)\n",
      "NotebookApp.rate_limit_window=3.0 (secs)\n",
      "\n",
      "|-----> in progress: 84.1499%IOPub message rate exceeded.\n",
      "The notebook server will temporarily stop sending output\n",
      "to the client in order to avoid crashing it.\n",
      "To change this limit, set the config variable\n",
      "`--NotebookApp.iopub_msg_rate_limit`.\n",
      "\n",
      "Current values:\n",
      "NotebookApp.iopub_msg_rate_limit=1000.0 (msgs/sec)\n",
      "NotebookApp.rate_limit_window=3.0 (secs)\n",
      "\n",
      "|-----> [download] in progress: 100.0000%\n",
      "|-----> [download] finished [11566.5517s]\n"
     ]
    }
   ],
   "source": [
    "organoid = dyn.sample_data.scEU_seq_organoid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AnnData object with n_obs × n_vars = 5422 × 11848\n",
       "    obs: 'Plate_Id', 'Condition_Id', 'Well_Id', 'RFP_log10_corrected', 'GFP_log10_corrected', 'Cell_cycle_possition', 'Cell_cycle_relativePos', 'exp_type', 'time'\n",
       "    var: 'Gene_Id'\n",
       "    layers: 'sl', 'su', 'ul', 'uu'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "organoid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'DataFrame' object has no attribute 'som_cluster_id'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-7-eb09c2f6d4d8>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     14\u001b[0m  }\n\u001b[1;32m     15\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 16\u001b[0;31m \u001b[0morganoid\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mobs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'cell_type'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0morganoid\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mobs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msom_cluster_id\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcell_mapper\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'str'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/python37test-env/lib/python3.7/site-packages/pandas/core/generic.py\u001b[0m in \u001b[0;36m__getattr__\u001b[0;34m(self, name)\u001b[0m\n\u001b[1;32m   5485\u001b[0m         ):\n\u001b[1;32m   5486\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 5487\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mobject\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__getattribute__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   5488\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   5489\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__setattr__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m->\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'DataFrame' object has no attribute 'som_cluster_id'"
     ]
    }
   ],
   "source": [
    "# mapping:\n",
    "cell_mapper = {\n",
    "    '1': 'Enterocytes',\n",
    "    '2': 'Enterocytes',\n",
    "    '3': 'Enteroendocrine',\n",
    "    '4': 'Enteroendocrine progenitor',\n",
    "    '5': 'Tuft cells',\n",
    "    '6': 'TA cells',\n",
    "    '7': 'TA cells',\n",
    "    '8': 'Stem cells',\n",
    "    '9': 'Paneth cells',\n",
    "    '10': 'Goblet cells',\n",
    "    '11': 'Stem cells',\n",
    " }\n",
    "\n",
    "organoid.obs['cell_type'] = organoid.obs.som_cluster_id.map(cell_mapper).astype('str')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## typical dynamo analysis workflow\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dyn.pl.basic_stats(organoid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "organoid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "organoid.obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "organoid.obs.groupby(['exp_type', 'time']).agg('count')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = organoid.copy()\n",
    "adata.obs.time = adata.obs.time.astype('str')\n",
    "adata.obs.loc[adata.obs['time'] == 'dmso', 'time'] = -1\n",
    "adata.obs['time'] = adata.obs['time'].astype(float)\n",
    "adata = adata[adata.obs.time != -1, :]\n",
    "adata = adata[adata.obs.exp_type == 'Pulse', :]\n",
    "adata.layers['new'], adata.layers['total'] = adata.layers['ul'] + adata.layers['sl'], adata.layers['su'] + adata.layers['sl'] + adata.layers['uu'] + adata.layers['ul']\n",
    "del adata.layers['uu'], adata.layers['ul'], adata.layers['su'], adata.layers['sl']\n",
    "dyn.pp.recipe_monocle(adata, n_top_genes=1000, total_layers=False)\n",
    "dyn.pl.basic_stats(adata)\n",
    "dyn.pl.show_fraction(organoid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs.time = adata.obs.time/60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs.time  = adata.obs.time.astype('float')\n",
    "dyn.tl.dynamics(adata, model='deterministic', tkey='time', assumption_mRNA='ss')\n",
    "\n",
    "dyn.tl.reduceDimension(adata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dyn.tl.cell_velocities(adata, ekey='M_t', vkey='velocity_T', enforce=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obsm['X_umap_ori'] = adata.obs.loc[:, ['rotated_umap1', 'rotated_umap2']].values.astype(float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize time-resolved vector flow learned with dynamo "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dyn.tl.cell_velocities(adata, basis='umap_ori')\n",
    "\n",
    "dyn.pl.streamline_plot(adata, color='cell_type', basis='umap_ori')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dyn.pl.streamline_plot(adata, color='cell_cycle_phase', basis='umap_ori')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.var_names[adata.var.use_for_transition][:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dyn.pl.phase_portraits(adata, genes=['Brat1', 'Ccnd2', 'Ckmt1', 'Pdgfb', 'Gpa33'],\n",
    "                       color='som_cluster_id', basis='umap_ori')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Animate intestine organoid differentiation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dyn.vf.VectorField(adata, basis='umap_ori')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "progenitor = adata.obs_names[adata.obs.cell_type == 'Stem cells']\n",
    "len(progenitor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(19491001)\n",
    "\n",
    "from matplotlib import animation\n",
    "info_genes = adata.var_names[adata.var.use_for_transition]\n",
    "dyn.pd.fate(adata, basis='umap_ori', init_cells=progenitor[:100], interpolation_num=100,  direction='forward',\n",
    "   inverse_transform=False, average=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax = dyn.pl.topography(adata, basis='umap_ori', color='cell_type', ax=ax, save_show_or_return='return',  figsize=(24, 24))\n",
    "ax.set_aspect(0.8)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "adata.obs['time'] = adata.obs.time.astype('float')\n",
    "instance = dyn.mv.StreamFuncAnim(adata=adata, basis='umap_ori', color='cell_type', ax=ax)"
   ]
  }
 ],
 "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.9.13"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}