notebooks/brown_adipocytes_subtypes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of Mouse Adipose Data Set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of contents:\n",
    "\n",
    "* <a href=#Load>Load Packages and Set Global Variables</a>\n",
    "    * <a href=#Imports>Imports and Settings</a>\n",
    "    * <a href=#Global>Global Variables</a> \n",
    "* <a href=#Dataloading>Loading Data, Quality Control and Preprocessing</a>\n",
    "* <a href=#Define>Define Cell Types</a>\n",
    "* <a href=#adipocytes>Adipocytes Only</a>\n",
    "    * <a href=#Embedding>Embeddings and Clustering</a>\n",
    "    * <a href=#adipmarkers>Adipocyte Marker Analysis</a>\n",
    "    * <a href=#topde>Top ranking DE Genes</a>\n",
    "    * <a href=#count_dist>Count distribution for Tcf25, Bin1 and Eif5</a>\n",
    "* <a href=#traject>Preadipocytes only</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"Load\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Packages and Set Global Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"imports\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports and Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scanpy as sc\n",
    "import sys\n",
    "import pandas as pd\n",
    "import seaborn as sb\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import colors\n",
    "from gprofiler import GProfiler\n",
    "import matplotlib_venn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: If you miss a compact list, please try `print_header`!\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----\n",
      "anndata     0.7.5\n",
      "scanpy      1.7.1\n",
      "sinfo       0.3.1\n",
      "-----\n",
      "PIL                 8.1.2\n",
      "PyObjCTools         NA\n",
      "anndata             0.7.5\n",
      "appnope             0.1.2\n",
      "backcall            0.2.0\n",
      "certifi             2020.12.05\n",
      "cffi                1.14.5\n",
      "chardet             4.0.0\n",
      "cloudpickle         1.6.0\n",
      "colorama            0.4.4\n",
      "cycler              0.10.0\n",
      "cython_runtime      NA\n",
      "dask                2021.03.0\n",
      "dateutil            2.8.1\n",
      "decorator           4.4.2\n",
      "get_version         2.1\n",
      "gprofiler           1.0.0\n",
      "h5py                3.2.1\n",
      "idna                2.10\n",
      "igraph              0.9.0\n",
      "ipykernel           5.4.3\n",
      "ipython_genutils    0.2.0\n",
      "ipywidgets          7.6.3\n",
      "jedi                0.17.2\n",
      "joblib              1.0.1\n",
      "kiwisolver          1.3.1\n",
      "legacy_api_wrap     1.2\n",
      "leidenalg           0.8.3\n",
      "llvmlite            0.36.0\n",
      "matplotlib          3.3.4\n",
      "mpl_toolkits        NA\n",
      "natsort             7.1.1\n",
      "numba               0.53.0\n",
      "numexpr             2.7.3\n",
      "numpy               1.20.1\n",
      "packaging           20.9\n",
      "pandas              1.2.3\n",
      "parso               0.7.0\n",
      "pexpect             4.8.0\n",
      "pickleshare         0.7.5\n",
      "pkg_resources       NA\n",
      "prompt_toolkit      3.0.14\n",
      "psutil              5.8.0\n",
      "ptyprocess          0.7.0\n",
      "pyexpat             NA\n",
      "pygments            2.7.4\n",
      "pyparsing           2.4.7\n",
      "pytz                2021.1\n",
      "requests            2.25.1\n",
      "scanpy              1.7.1\n",
      "scipy               1.6.1\n",
      "seaborn             0.11.1\n",
      "sinfo               0.3.1\n",
      "six                 1.15.0\n",
      "sklearn             0.24.1\n",
      "sparse              0.12.0\n",
      "sphinxcontrib       NA\n",
      "statsmodels         0.12.2\n",
      "storemagic          NA\n",
      "tables              3.6.1\n",
      "texttable           1.6.3\n",
      "tlz                 0.11.1\n",
      "toolz               0.11.1\n",
      "tornado             6.1\n",
      "traitlets           5.0.5\n",
      "typing_extensions   NA\n",
      "urllib3             1.26.3\n",
      "wcwidth             0.2.5\n",
      "yaml                5.4.1\n",
      "zmq                 22.0.2\n",
      "zope                NA\n",
      "-----\n",
      "IPython             7.20.0\n",
      "jupyter_client      6.1.11\n",
      "jupyter_core        4.7.1\n",
      "notebook            6.2.0\n",
      "-----\n",
      "Python 3.9.2 (v3.9.2:1a79785e3e, Feb 19 2021, 09:06:10) [Clang 6.0 (clang-600.0.57)]\n",
      "macOS-10.16-x86_64-i386-64bit\n",
      "12 logical CPU cores, i386\n",
      "-----\n",
      "Session information updated at 2021-06-07 06:22\n",
      "\n",
      "3.9.2 (v3.9.2:1a79785e3e, Feb 19 2021, 09:06:10) \n",
      "[Clang 6.0 (clang-600.0.57)]\n"
     ]
    }
   ],
   "source": [
    "base_dir = '/Users/viktorian.miok/Documents/consultation/Ruth/David/data/bat/MUC8387/'\n",
    "dir_out = '/Users/viktorian.miok/Documents/consultation/Ruth/David/results/'\n",
    "\n",
    "sc.settings.verbosity = 3 # amount of output\n",
    "dir_in = base_dir+'mm10_ensrel94/'\n",
    "dir_tables = dir_out+'tables/'\n",
    "sc.settings.figdir = dir_out+'panels/'\n",
    "dir_adata = dir_out+'anndata/'\n",
    "sc.logging.print_versions()\n",
    "sc.settings.set_figure_params(dpi=80, \n",
    "                              scanpy=True\n",
    ")\n",
    "print (sys.version)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define a nice colour map for gene expression\n",
    "colors2 = plt.cm.Blues(np.linspace(0, 1, 128))\n",
    "colors3 = plt.cm.Greys_r(np.linspace(0.7,0.8,20))\n",
    "colorsComb = np.vstack([colors3, colors2])\n",
    "mymap = colors.LinearSegmentedColormap.from_list('my_colormap',\n",
    "                                                 colorsComb\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Global Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All embeddings and clusterings can be saved and loaded into this script. Be carful with overwriting cluster caches as soon as cell type annotation has started as cluster labels may be shuffled."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set whether anndata objects are recomputed or loaded from cache."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "bool_recomp = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set whether clustering is recomputed or loaded from saved .obs file. Loading makes sense if the clustering changes due to a change in scanpy or one of its dependencies and the number of clusters or the cluster labels change accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "bool_recluster = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set whether cluster cache is overwritten. Note that the cache exists for reproducibility of clustering, see above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "bool_write_cluster_cache = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set whether to produce plots, set to False for test runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "bool_plot = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set whether observations should be calculated. If false, it is necessary to read cacheed file that contains the necssary information. It then shows the the distributions of counts and genes, as well as mt_frac after filtering. \n",
    "Set to true in order to see the data before filtering and follow the decisions for cutoffs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "bool_create_observations = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"Dataloading\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loading Data, Quality Control and Preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the data in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_raw = sc.read(dir_in+'raw_gene_bc_matrices.h5ad')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"QC\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quality Control"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Summary of steps performed here: Only cells with at least 500 UMIs are kept. Counts per cell are cell library depth normalized. The gene (feature) space is reduced with PCA to 50 PCs. A nearest neighbour graph and umap are computed based on the PC space. Cell are clustered with louvain clustering based on the nearest neighbour graph. Graph abstraction is computed based on the louvain clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(737280, 31125)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adata_raw.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data contains 737280 observations with 31125 different genes. Due to dropouts, some of the observations might not show any counts and genes. In order to calculate the fraction of mitochondrial RNA in the next steps, each observations without counts must be filtered out to prevent NaN from emerging. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "filtered out 352632 cells that have less than 1 counts\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of cells: 737280\n",
      "Number of cells after min count filter: 384648\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(384648, 31125)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('Total number of cells: {:d}'.format(adata_raw.n_obs))\n",
    "\n",
    "sc.pp.filter_cells(adata_raw, min_counts = 1)\n",
    "print('Number of cells after min count filter: {:d}'.format(adata_raw.n_obs))\n",
    "\n",
    "adata_raw.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the 737280 observations, only 384648 are left with counts greater than 0, which is ~52%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"n_gen_counts\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gene numbers and counts with and without mitochondrial RNA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create necessary obs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_qc = adata_raw.copy()\n",
    "adata_qc.obs['n_genes'] = (adata_qc.X > 0).sum(1)\n",
    "mt_gene_mask = [gene.startswith('mt-') for gene in adata_qc.var_names]\n",
    "temp_mt_sum = adata_qc[:,mt_gene_mask].X.sum(1)\n",
    "temp_mt_sum = np.squeeze(np.asarray(temp_mt_sum))\n",
    "temp_n_counts = adata_qc.obs['n_counts']\n",
    "adata_qc.obs['mt_frac'] = temp_mt_sum/adata_qc.obs['n_counts']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot n_counts and mt_frac:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    t1 = sc.pl.violin(adata_qc,\n",
    "                      'n_counts',\n",
    "                      size =2, \n",
    "                      log=True, \n",
    "                      cut = 0\n",
    "    )\n",
    "    t2 = sc.pl.violin(adata_qc, \n",
    "                      'mt_frac'\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Overall, the data contains a lot of observations with high fractions of mitochondrial RNA. Additionally, most observations show counts below 100, suggesting poor data quality. To further investigate the distributions counts over genes per observations, scatterplots are created:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Number of Genes versus Number of Counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    p1 = sc.pl.scatter(adata_qc, \n",
    "                       'n_counts',\n",
    "                       'n_genes',\n",
    "                       color='mt_frac',\n",
    "                       size = 5\n",
    "    )\n",
    "    p2 = sc.pl.scatter(adata_qc[adata_qc.obs['n_counts'] < 10000],\n",
    "                       'n_counts',\n",
    "                       'n_genes',\n",
    "                       color='mt_frac',\n",
    "                       size=5\n",
    "    )\n",
    "    p3 = sc.pl.scatter(adata_qc[adata_qc.obs['n_counts'] < 1000],\n",
    "                       'n_counts',\n",
    "                       'n_genes',\n",
    "                       color='mt_frac', \n",
    "                       size=5\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distribution of Counts and Genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    p6 = sb.distplot(adata_qc.obs['n_counts'], \n",
    "                     kde = False\n",
    "    )\n",
    "    plt.show()\n",
    "    p7 = sb.distplot(adata_qc.obs['n_counts'][adata_qc.obs['n_counts'] < 1000], \n",
    "                     kde = False\n",
    "    )\n",
    "    plt.show()\n",
    "    p8 = sb.distplot(adata_qc.obs['n_counts'][adata_qc.obs['n_counts'] < 100],\n",
    "                     kde = False\n",
    "    )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    p9 = sb.distplot(adata_qc.obs['n_genes'],\n",
    "                     kde=False, \n",
    "                     bins=60\n",
    "    )\n",
    "    plt.show()\n",
    "    p10 = sb.distplot(adata_qc.obs['n_genes'][adata_qc.obs['n_genes'] < 100],\n",
    "                      kde = False, \n",
    "                      bins=60\n",
    "    )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "filtered out 380040 cells that have less than 200 counts\n",
      "filtered out 32 cells that have more than 15000 counts\n",
      "filtered out 150 cells that have less than 100 genes expressed\n",
      "Trying to set attribute `.obs` of view, copying.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of cells: 384648\n",
      "Number of cells after min count filter: 4608\n",
      "Number of cells after max count filter: 4576\n",
      "Number of cells after MT filter: 4163\n",
      "Number of cells after gene filter: 4013\n"
     ]
    }
   ],
   "source": [
    "# Filter cells according to identified QC thresholds:\n",
    "print('Total number of cells: {:d}'.format(adata_qc.n_obs))\n",
    "\n",
    "sc.pp.filter_cells(adata_qc, min_counts = 200)\n",
    "print('Number of cells after min count filter: {:d}'.format(adata_qc.n_obs))\n",
    "\n",
    "sc.pp.filter_cells(adata_qc, max_counts = 15000)\n",
    "print('Number of cells after max count filter: {:d}'.format(adata_qc.n_obs))\n",
    "\n",
    "adata_qc = adata_qc[adata_qc.obs['mt_frac'] < 0.2]\n",
    "print('Number of cells after MT filter: {:d}'.format(adata_qc.n_obs))\n",
    "\n",
    "sc.pp.filter_cells(adata_qc, min_genes = 100)\n",
    "print('Number of cells after gene filter: {:d}'.format(adata_qc.n_obs))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only 4013 observations are left that pass the filters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    p12 = sc.pl.scatter(adata_qc,\n",
    "                        'n_counts',\n",
    "                        'n_genes',\n",
    "                        color='mt_frac', \n",
    "                        size=5\n",
    "    )\n",
    "    p22 = sc.pl.scatter(adata_qc[adata_qc.obs['n_counts'] < 10000],\n",
    "                        'n_counts',\n",
    "                        'n_genes',\n",
    "                        color='mt_frac', \n",
    "                        size=5\n",
    "    )\n",
    "    p52 = sc.pl.scatter(adata_qc[adata_qc.obs['n_counts'] < 1000],\n",
    "                        'n_counts',\n",
    "                        'n_genes',\n",
    "                        color='mt_frac',\n",
    "                        size=5\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.violin(adata_qc, \n",
    "                 'n_counts',\n",
    "                 size=2,\n",
    "                 log=True,\n",
    "                 cut=0\n",
    "    )\n",
    "    sc.pl.violin(adata_qc,\n",
    "                 'mt_frac'\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the remaining observations, the fraction of mitochondrial RNA is generally very low and at most 20%"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    p62 = sb.distplot(adata_qc.obs['n_counts'],\n",
    "                      kde=False\n",
    "    )\n",
    "    plt.show()\n",
    "    p72 = sb.distplot(adata_qc.obs['n_counts'][adata_qc.obs['n_counts'] < 1000],\n",
    "                      kde=False\n",
    "    )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    p92= sb.distplot(adata_qc.obs['n_genes'],\n",
    "                     kde=False,\n",
    "                     bins=60\n",
    "    )\n",
    "    plt.show()\n",
    "    p100 = sb.distplot(adata_qc.obs['n_genes'][adata_qc.obs['n_genes'] < 1000],\n",
    "                       kde=False,\n",
    "                       bins=60\n",
    "    )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "filtered out 18501 genes that are detected in less than 20 cells\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of genes: 31125\n",
      "Number of genes after cell filter: 12624\n"
     ]
    }
   ],
   "source": [
    "#Filter genes:\n",
    "print('Total number of genes: {:d}'.format(adata_qc.n_vars))\n",
    "\n",
    "# Min 20 cells - filters out 0 count genes\n",
    "sc.pp.filter_genes(adata_qc,\n",
    "                   min_cells=20\n",
    ")\n",
    "print('Number of genes after cell filter: {:d}'.format(adata_qc.n_vars))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"normalization_and_clustering\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normalization and Clustering with highly variable genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/anndata/compat/__init__.py:179: FutureWarning: Moving element from .uns['neighbors']['distances'] to .obsp['distances'].\n",
      "\n",
      "This is where adjacency matrices should go now.\n",
      "  warn(\n",
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/anndata/compat/__init__.py:179: FutureWarning: Moving element from .uns['neighbors']['connectivities'] to .obsp['connectivities'].\n",
      "\n",
      "This is where adjacency matrices should go now.\n",
      "  warn(\n",
      "running PAGA\n",
      "    finished: added\n",
      "    'paga/connectivities', connectivities adjacency (adata.uns)\n",
      "    'paga/connectivities_tree', connectivities subtree (adata.uns) (0:00:00)\n"
     ]
    }
   ],
   "source": [
    "if bool_recomp == True:\n",
    "        \n",
    "    adata_proc = adata_qc.copy()\n",
    "    adata_proc.raw = adata_proc\n",
    "    sc.pp.normalize_per_cell(adata_proc)\n",
    "    sc.pp.log1p(adata_proc)\n",
    "    sc.pp.highly_variable_genes(adata_proc,\n",
    "                                flavor='cell_ranger',\n",
    "                                n_top_genes=4000\n",
    "    )\n",
    "    sc.pl.highly_variable_genes(adata_proc)\n",
    "   \n",
    "    sc.pp.pca(adata_proc,\n",
    "              n_comps=50, \n",
    "              random_state=0, \n",
    "              use_highly_variable=True, \n",
    "              svd_solver='arpack'\n",
    "    )\n",
    "    sc.pp.neighbors(adata_proc,\n",
    "                    n_neighbors=100, \n",
    "                    knn=True,\n",
    "                    method='umap',\n",
    "                    n_pcs=50,\n",
    "                    random_state=0\n",
    "    )\n",
    "    #sc.tl.tsne(adata_proc, n_jobs=3)\n",
    "    sc.tl.umap(adata_proc)\n",
    "    if bool_recluster == True:\n",
    "        sc.tl.louvain(adata_proc, \n",
    "                      resolution=1,\n",
    "                      flavor='vtraag',\n",
    "                      random_state=0\n",
    "        )\n",
    "        pd.DataFrame(adata_proc.obs).to_csv(path_or_buf=dir_adata+\"obs_adata_proc.csv\")\n",
    "    else:\n",
    "        obs = pd.read_csv(dir_adata+'obs_adata_proc.csv')\n",
    "        adata_proc.obs['louvain'] = pd.Series(obs['louvain'].values,\n",
    "                                              dtype='category')\n",
    "    sc.write(dir_adata+'adata_proc.h5ad', adata_proc)\n",
    "else:\n",
    "    adata_proc = sc.read(dir_adata+'adata_proc.h5ad') \n",
    "sc.tl.paga(adata_proc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Produce some summarizing plots that show the global characteristics of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.umap(adata_proc,\n",
    "               color=['louvain', 'mt_frac'],\n",
    "               size=5,\n",
    "               save=\"_all_louvain.svg\",\n",
    "               use_raw=False, \n",
    "               color_map=mymap\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A high fraction of mitochondrial RNA is in cluster 7 and around the central cluster 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.umap(adata_proc, \n",
    "               color=['n_genes', 'n_counts'],\n",
    "               size=5, \n",
    "               save=\"_all_n_counts.svg\", \n",
    "               use_raw=False,\n",
    "               color_map=mymap\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the number of genes and number of counts, cluster 0 and 8 could contain preadipocytes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.paga(adata_proc,\n",
    "               save=\"_all.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of cells in each cluster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0     614\n",
       "1     473\n",
       "2     422\n",
       "3     366\n",
       "4     328\n",
       "5     304\n",
       "6     299\n",
       "7     272\n",
       "8     270\n",
       "9     254\n",
       "10    207\n",
       "11    204\n",
       "Name: louvain, dtype: int64"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adata_proc.obs[\"louvain\"].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define Cell Types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"DE\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DE Genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.tl.rank_genes_groups(adata_proc, \n",
    "                            groupby='louvain',\n",
    "                            key_added='rank_genes'\n",
    "    )\n",
    "    sc.pl.rank_genes_groups(adata_proc,\n",
    "                            key='rank_genes', \n",
    "                            groups=['0','1','2'],\n",
    "                            fontsize=12,\n",
    "                            save=\"_adata_proc_genes_1.svg\"\n",
    "    )\n",
    "    sc.pl.rank_genes_groups(adata_proc,\n",
    "                            key='rank_genes', \n",
    "                            groups=['3','4','5'], \n",
    "                            fontsize=12,\n",
    "                            save=\"_adata_proc_genes_2.svg\"\n",
    "    )\n",
    "    sc.pl.rank_genes_groups(adata_proc,\n",
    "                            key='rank_genes',\n",
    "                            groups=['6','7','8'], \n",
    "                            fontsize=12,\n",
    "                            save=\"_adata_proc_genes_2.svg\"\n",
    "    )\n",
    "    sc.pl.rank_genes_groups(adata_proc, \n",
    "                            key='rank_genes', \n",
    "                            groups=['9','10'], \n",
    "                            fontsize=12,\n",
    "                            save=\"_adata_proc_genes_2.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Marker Sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define marker sets for some of the expected cell types and add DE genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Leukocyte markers:\n",
    "astrocyte_markers = ['Gfap','Aldh1l1','Slc1a2','Slc1a3','Slc27a1']\n",
    "leukocyte_markers = ['Ptprc','Tcf25']\n",
    "tc_markers = ['Cd3d','Cd3e','Cd3g','Cd4','Cd8a1','Cd8b','Ms4a4b','Trbc2']\n",
    "nk_markers = ['Nkg7','Il2rb', 'Ncr1', 'Klrd1', 'Klrb1b', 'Klrb1f']\n",
    "myeloid_markers = ['Cd79a', 'Itgax', 'Itgam', 'Fcgr3','S100a8', 'S100a9']\n",
    "mp_markers = ['Adgre1', 'Lyz2','C1qa','Pf4']\n",
    "dc_markers = ['Cd74', 'Anpep' , 'Cd33', 'Cd80', 'Cd83', 'Cd86'] \n",
    "bc_markers = ['Cd19','Igkc','Ighm','Cd79a']\n",
    "adipocyte_markers = ['Fgf10','Bmp4','Zfp423','Psmb8','Pdgfra', 'Slc7a10', 'Pparg', 'Cd34', 'Dlk1', 'Dpp4', 'Fabp4',\n",
    "                     'Ucp1', 'Slc36a2', 'Tmem26','Adipoq','Prdm16','Ppargc1a','Tcf25','Bin1','Eif5','Slc2a4',\n",
    "                     'Cebpa','Lpl','Pnpla2','Plin1'] # 'Slc7a10' is 'Asc1'; Dlk1 is Pref1; Cd29 is Itgb1  ,'P2rd5'\n",
    "differentialtion_markers = ['Pparg','Cebpa','Cebpb','Cebpd','Dlk1','Fabp4']\n",
    "megakaryocyte_markers = ['Ppbp']\n",
    "erythrocyte_markers = ['Gypa','Hba-a2','Hbb-bt','Alas2']\n",
    "go_adip_dev = ['Aacs', 'Acat1', 'Arid5b','Arrdc3', 'Atf2', 'Bbs4', 'Bdh1', 'Csf1', 'Dgat2','Dyrk1b', 'Ebf2',\n",
    "               'Fam123b', 'Fto','Hgmcs2', 'Id2', 'Lep', 'Lrp5', 'Nampt', 'Oxct1', 'Paxip1', 'Pik3ca', 'Ppard',\n",
    "               'Ppargc1a', 'Rorc', 'Sh3pxd2b', 'Slc25a25',  'Sox8', 'Spg20', 'Tbl1xr1', 'Umold1', 'Xbp1', 'Znf516']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only keep markers occurring in data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "astrocyte_markers = np.array([x for x in astrocyte_markers if x in adata_proc.var_names])\n",
    "leukocyte_markers = np.array([x for x in leukocyte_markers if x in adata_proc.var_names])\n",
    "tc_markers = np.array([x for x in tc_markers if x in adata_proc.var_names])\n",
    "nk_markers = np.array([x for x in nk_markers if x in adata_proc.var_names])\n",
    "myeloid_markers = np.array([x for x in myeloid_markers if x in adata_proc.var_names])\n",
    "mp_markers = np.array([x for x in mp_markers if x in adata_proc.var_names])\n",
    "dc_markers = np.array([x for x in dc_markers if x in adata_proc.var_names])\n",
    "bc_markers = np.array([x for x in bc_markers if x in adata_proc.var_names])\n",
    "megakaryocyte_markers = np.array([x for x in megakaryocyte_markers if x in adata_proc.var_names])\n",
    "erythrocyte_markers = np.array([x for x in erythrocyte_markers if x in adata_proc.var_names])\n",
    "adipocyte_markers = np.array([x for x in adipocyte_markers if x in adata_proc.var_names])\n",
    "differentialtion_markers = np.array([x for x in differentialtion_markers if x in adata_proc.var_names])\n",
    "go_adip_dev = np.array([x for x in go_adip_dev if x in adata_proc.var_names])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"Plotting\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting Routines for Markers Gene Sets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_violin_marker(adata, markers, save=None, use_raw=True):\n",
    "    for i in range(len(markers) // 2 + len(markers) % 2):\n",
    "        if save is not None:\n",
    "            sc.pl.violin(adata, \n",
    "                         groupby='louvain', \n",
    "                         keys=markers[(2*i):np.min([2*(i+1), len(markers)])], \n",
    "                         use_raw=use_raw, \n",
    "                         rotation=90,\n",
    "                         size=5,\n",
    "                         save=save+\"_\"+str(i)+\".svg\"\n",
    "            )\n",
    "        else:\n",
    "            sc.pl.violin(adata, \n",
    "                         groupby='louvain', \n",
    "                         keys=markers[(2*i):np.min([2*(i+1), len(markers)])], \n",
    "                         use_raw=use_raw, \n",
    "                         rotation=90, \n",
    "                         size=5\n",
    "            )\n",
    "def plot_umap_marker(adata, markers, size=3, save=None, use_raw=True):\n",
    "    for i in range(len(markers) // 2 + len(markers) % 2):\n",
    "        print(markers[(2*i):np.min([2*(i+1), len(markers)])])\n",
    "        if save is not None:\n",
    "            sc.pl.umap(adata, \n",
    "                       color=markers[(2*i):np.min([2*(i+1), len(markers)])], \n",
    "                       size=size,\n",
    "                       use_raw=use_raw,\n",
    "                       color_map=mymap,\n",
    "                       save=save+\"_\"+str(i)+\".svg\"\n",
    "            )\n",
    "        else:\n",
    "            sc.pl.umap(adata, \n",
    "                       color=markers[(2*i):np.min([2*(i+1), len(markers)])], \n",
    "                       size=size,\n",
    "                       use_raw=use_raw,\n",
    "                       color_map=mymap\n",
    "            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc,\n",
    "                       astrocyte_markers.tolist()\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Leukocyte Markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc, \n",
    "                       leukocyte_markers.tolist(),\n",
    "                       save=\"_all_markers_leukocyte\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cluster 1,9 and 10 show expression for the leukocyte marker Ptprc. Tcf25 found in many clusters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Erythrocyte Markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc,\n",
    "                       erythrocyte_markers.tolist(),\n",
    "                       save=\"_all_markers_erythrcytes\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The haemoglobin markers and Alas2 show high counts in cluster 5 indicating that cluster 5 corresponds to erythrocytes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"Preadipocyte\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preadipocyte Markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc,\n",
    "                       adipocyte_markers.tolist(),\n",
    "                       save=\"_all_markers_preadipocytes\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## T-Cell Markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc,\n",
    "                       tc_markers.tolist(),\n",
    "                       save=\"_all_markers_tcells\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "T-cell markers are strongly expressed in cluster 6."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"Natural\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Natural Killer Cells Markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc, \n",
    "                       nk_markers.tolist(),\n",
    "                       save=\"_all_markers_nk\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NK cell markers are expressed in cluster 6, which is inline with their transcriptomic similarity to T cells which are also found in cluster 6."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Myeloid Cell Markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc, \n",
    "                       myeloid_markers.tolist(),\n",
    "                       save=\"_all_markers_myeloid\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Macrophage Markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc,\n",
    "                       mp_markers.tolist(),\n",
    "                       save=\"_all_markers_macrophages\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The macrophage markers are strongly expressed in cluster 3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"Dendridic\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dendritic cell Markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc, \n",
    "                       dc_markers.tolist(),\n",
    "                       save=\"_all_markers_dendritic\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dendritic cells markers are stongly expressed in cluster 3 and 4."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"B-Cell\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## B-Cell Markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc, \n",
    "                       bc_markers.tolist(),\n",
    "                       save=\"_all_markers_bcells\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "B cell markers are highly expressed in cluster 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=heatmap></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary Heatmap for cluster assignments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_leukocyte_markers = ['Ptprc']\n",
    "selected_tc_markers = ['Cd3d', 'Cd3g','Cd3e']\n",
    "selected_nk_markers = ['Nkg7','Klrd1','Ncr1','Il2rb']\n",
    "selected_myeloid_markers = ['Fcgr3','S100a8']\n",
    "selected_mp_markers = ['Adgre1', 'C1qa','Pf4']\n",
    "selected_dc_markers = ['Cd74', 'Cd83','Cd86'] \n",
    "selected_bc_markers = ['Cd19','Igkc','Ighm','Cd79a']\n",
    "selected_adipocyte_markers = ['Pdgfra','Cd34', 'Dlk1','Fabp4','Fgf10','Bmp4','Pparg','Ucp1','Prdm16','Lpl'] # 'Slc7a10' is 'Asc1'; Dlk1 is Pref1; Cd29 is Itgb1\n",
    "#selected_megakaryocyte_markers = ['Ppbp']\n",
    "selected_erythrocyte_markers = ['Gypa','Hba-a2','Hbb-bt','Alas2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_cell_markers = selected_leukocyte_markers + \\\n",
    "selected_erythrocyte_markers + \\\n",
    "selected_myeloid_markers + \\\n",
    "selected_mp_markers + \\\n",
    "selected_dc_markers + \\\n",
    "selected_bc_markers + \\\n",
    "selected_tc_markers + \\\n",
    "selected_nk_markers + \\\n",
    "selected_adipocyte_markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot==True:\n",
    "    sc.pl.heatmap(adata=adata_proc, \n",
    "                  var_names=selected_cell_markers, \n",
    "                  groupby=\"louvain\", \n",
    "                  use_raw=False, \n",
    "                  log=False, \n",
    "                  dendrogram=True, \n",
    "                  var_group_rotation=90, \n",
    "                  show_gene_labels=True, \n",
    "                  show=True, \n",
    "                  save=\"_all_markers_celltypes.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## UMAP with assigned cell types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_cluster_names = {\n",
    "    '0': \"preadipocytes\",\n",
    "    '1': \"low expression\",\n",
    "    '2': \"lymphoid (B cell)\",\n",
    "    '3': \"myeloid (DC, macrophage)\",\n",
    "    '4': \"myeloid (DC)\",\n",
    "    '5': \"erythrocytes\",\n",
    "    '6': \"lymphoid (TC, NK)\",\n",
    "    '7': \"adipocytes\",\n",
    "    '8': \"preadipocytes\",\n",
    "    '9': \"lymphoid (TC)\",\n",
    "    '10': \"myeloid (macrophage)\",\n",
    "    '11': \"myeloid (s100a8 high)\"\n",
    "}\n",
    "adata_proc.obs['celltypes'] = [new_cluster_names[x] for x in  adata_proc.obs['louvain']]\n",
    "\n",
    "if bool_plot == True:\n",
    "    sc.pl.umap(adata_proc,\n",
    "               size=5,\n",
    "               color=['celltypes'], \n",
    "               save=\"_all_louvain_named.svg\"\n",
    "    )\n",
    "    sc.pl.paga(adata_proc,\n",
    "               color=['celltypes'], \n",
    "               save=\"_all_named.svg\"\n",
    "    )\n",
    "    #sc.pl.umap(adata_proc, size=5, color=adipocyte_markers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"adipocytes\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adipocytes only"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"Embedding\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Embedding and Clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/anndata/compat/__init__.py:179: FutureWarning: Moving element from .uns['neighbors']['distances'] to .obsp['distances'].\n",
      "\n",
      "This is where adjacency matrices should go now.\n",
      "  warn(\n",
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/anndata/compat/__init__.py:179: FutureWarning: Moving element from .uns['neighbors']['connectivities'] to .obsp['connectivities'].\n",
      "\n",
      "This is where adjacency matrices should go now.\n",
      "  warn(\n",
      "running PAGA\n",
      "    finished: added\n",
      "    'paga/connectivities', connectivities adjacency (adata.uns)\n",
      "    'paga/connectivities_tree', connectivities subtree (adata.uns) (0:00:00)\n"
     ]
    }
   ],
   "source": [
    "if bool_recomp == True:  \n",
    "    cell_ids_adip = np.asarray(adata_proc.obs_names)[\n",
    "        [x in ['preadipocytes','adipocytes'] \n",
    "         for x in np.asarray(adata_proc.obs['celltypes'].values)]\n",
    "    ]\n",
    "    adata_adip = adata_raw[cell_ids_adip,:].copy()\n",
    "    adata_adip.obs['n_genes'] = (adata_adip.X > 0).sum(1)\n",
    "    mt_gene_mask = [gene.startswith('mt-') for gene in adata_adip.var_names]\n",
    "    temp_mt_sum = adata_adip[:,mt_gene_mask].X.sum(1)\n",
    "    temp_mt_sum = np.squeeze(np.asarray(temp_mt_sum))\n",
    "    temp_n_counts = adata_adip.obs['n_counts']\n",
    "    adata_adip.obs['mt_frac'] = temp_mt_sum/adata_adip.obs['n_counts']\n",
    "    adata_adip.raw = adata_adip\n",
    "    sc.pp.normalize_per_cell(adata_adip)\n",
    "    sc.pp.log1p(adata_adip)\n",
    "    sc.pp.highly_variable_genes(adata_adip,\n",
    "                                n_top_genes=4000\n",
    "    )\n",
    "    sc.pl.highly_variable_genes(adata_adip)\n",
    "    sc.pp.pca(adata_adip,\n",
    "              n_comps=50, \n",
    "              use_highly_variable=True, \n",
    "              random_state=0,\n",
    "              svd_solver='arpack'\n",
    "    )\n",
    "    sc.pp.neighbors(adata_adip, \n",
    "                    n_neighbors=100,\n",
    "                    knn=True,\n",
    "                    method='umap',\n",
    "                    n_pcs=50,\n",
    "                    random_state=0\n",
    "    )\n",
    "    sc.tl.umap(adata_adip)\n",
    "    if bool_recluster == True:\n",
    "        sc.tl.louvain(adata_adip, \n",
    "                      resolution=1,\n",
    "                      flavor='vtraag', \n",
    "                      random_state=0\n",
    "        )\n",
    "        pd.DataFrame(adata_adip.obs).to_csv(path_or_buf=dir_adata+'obs_adata_adip.csv')\n",
    "    else:\n",
    "        obs = pd.read_csv(dir_adata+'obs_adata_adip.csv')\n",
    "        adata_adip.obs['louvain']=pd.Series(obs['louvain'].values, \n",
    "                                            dtype='category'\n",
    "        )\n",
    "    sc.write(dir_adata+'adata_adip.h5ad', adata_adip)\n",
    "else:\n",
    "    adata_adip = sc.read(dir_adata+'adata_adip.h5ad') \n",
    "sc.tl.paga(adata_adip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.umap(adata_adip, \n",
    "               color=['louvain','mt_frac'],\n",
    "               size=30,\n",
    "               save=\"_adipocytes_louvain_0.svg\", \n",
    "               color_map=mymap\n",
    "    )\n",
    "    sc.pl.umap(adata_adip, \n",
    "               color=['n_genes','n_counts'],\n",
    "               size=30,\n",
    "               save=\"_adipocytes_louvain_1.svg\",\n",
    "               color_map=mymap\n",
    "    )\n",
    "    sc.pl.paga(adata_adip, \n",
    "               color=['louvain','n_counts'],\n",
    "               save=\"_adipocytes_louvain.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cluster 0,1,3 (top) shows many more counts than cluster 2 (bottom)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"adipmarkers\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adipocyte Marker Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_adip,\n",
    "                       adipocyte_markers.tolist(),\n",
    "                       save=\"_adipo_markers_adipocytes\",\n",
    "                       use_raw=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "FABP4 is present only in cluster 2. Together with Pparg expression in cluster 2, it can be concluded that cluster 2 represents mature adipocytes instead of preadipocytes. Clusters 0,1 and 3 express multiple preadpocyte markers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.heatmap(adata=adata_adip, \n",
    "                  var_names=adipocyte_markers, \n",
    "                  groupby=\"louvain\", \n",
    "                  use_raw=False, \n",
    "                  log=False, \n",
    "                  dendrogram=True, \n",
    "                  var_group_rotation=90, \n",
    "                  show_gene_labels=True, \n",
    "                  show=True, \n",
    "                  save=\"_adipocyte_markers_celltypes_in_adipocytes.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"topde\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Top ranking DE Genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Default of the method has been changed to 't-test' from 't-test_overestim_var'\n",
      "ranking genes\n",
      "    finished: added to `.uns['rank_genes']`\n",
      "    'names', sorted np.recarray to be indexed by group ids\n",
      "    'scores', sorted np.recarray to be indexed by group ids\n",
      "    'logfoldchanges', sorted np.recarray to be indexed by group ids\n",
      "    'pvals', sorted np.recarray to be indexed by group ids\n",
      "    'pvals_adj', sorted np.recarray to be indexed by group ids (0:00:00)\n"
     ]
    }
   ],
   "source": [
    "sc.tl.rank_genes_groups(adata_adip,\n",
    "                        groupby='louvain', \n",
    "                        n_genes=4000, \n",
    "                        key_added='rank_genes'\n",
    ")\n",
    "if bool_plot == True:\n",
    "    sc.pl.rank_genes_groups(adata_adip, \n",
    "                            key='rank_genes',\n",
    "                            groups=['0','1','2','3'],\n",
    "                            fontsize=12, \n",
    "                            save=\"_adipocytes.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = adata_adip.uns['rank_genes']\n",
    "groups = result['names'].dtype.names\n",
    "diffExp = pd.DataFrame(\n",
    "    {group + '_' + key[:1]: result[key][group]\n",
    "    for group in groups for key in ['names', 'pvals_adj','logfoldchanges']}).head(4000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "diffExp.columns = ['0_name','0_pval_ajd','0_logFC','1_name','1_pval_ajd','1_logFC',\n",
    "                   '2_name','2_pval_ajd','2_logFC','3_name','3_pval_ajd','3_logFC']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save a DE table\n",
    "diffExp.to_csv(path_or_buf=dir_tables+\"DE_adipocytes_by_cluster.csv\", sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"marker_counts\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Marker Counts Raw and Normalized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_umap_marker(adata_adip,\n",
    "                     adipocyte_markers,\n",
    "                     size=50,\n",
    "                     save=\"_preadip_markers_adipocytes_raw\",\n",
    "                     use_raw=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_umap_marker(adata_adip, \n",
    "                     adipocyte_markers, \n",
    "                     size=50,\n",
    "                     save=\"_preadip_markers_adipocytes_norm\", \n",
    "                     use_raw=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"count_dist\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Count distribution for Tcf25, Bin1 and Eif5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"countspg\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Single Double and Triple Positive Counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define booleans for single, double and triple positive counts of bin1, eif5 and tcf25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp == True:  \n",
    "    non_boolean_int = np.array((adata_adip[:,'Bin1'].X <= 0) & (adata_adip[:,'Eif5'].X <= 0) &\n",
    "                               (adata_adip[:,'Tcf25'].X <= 0), dtype=int)\n",
    "\n",
    "    bin1_single_boolean = (adata_adip[:,'Bin1'].X > 0) & (adata_adip[:,'Eif5'].X <= 0) & (adata_adip[:,'Tcf25'].X <= 0)\n",
    "    eif5_single_boolean = (adata_adip[:,'Eif5'].X > 0) & (adata_adip[:,'Bin1'].X <= 0) & (adata_adip[:,'Tcf25'].X <= 0)\n",
    "    tcf25_single_boolean = (adata_adip[:,'Tcf25'].X > 0) & (adata_adip[:,'Eif5'].X <= 0) & (adata_adip[:,'Bin1'].X <= 0)\n",
    "    single_boolean_int = np.array((bin1_single_boolean | eif5_single_boolean | tcf25_single_boolean), dtype=int)*1\n",
    "\n",
    "    bin1_eif5_double_boolean = (adata_adip[:,'Bin1'].X > 0) & (adata_adip[:,'Eif5'].X > 0) & (adata_adip[:,'Tcf25'].X <= 0)\n",
    "    bin1_tcf25_double_boolean = (adata_adip[:,'Eif5'].X <= 0) & (adata_adip[:,'Bin1'].X > 0) & (adata_adip[:,'Tcf25'].X > 0)\n",
    "    tcf25_eif5_double_boolean = (adata_adip[:,'Tcf25'].X > 0) & (adata_adip[:,'Eif5'].X > 0) & (adata_adip[:,'Bin1'].X <= 0)\n",
    "    double_boolean_int = np.array((bin1_eif5_double_boolean | bin1_tcf25_double_boolean | tcf25_eif5_double_boolean), dtype=int)*2\n",
    "\n",
    "    triple_boolean = (adata_adip[:,'Tcf25'].X > 0) & (adata_adip[:,'Eif5'].X > 0) & (adata_adip[:,'Bin1'].X > 0)\n",
    "    triple_boolean_int = np.array(triple_boolean, dtype=int)*3\n",
    "    non_boolean_int *= 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp == True: \n",
    "    bin1_single_pos = adata_adip[:,'Bin1'].X[bin1_single_boolean]\n",
    "    eif5_single_pos = adata_adip[:,'Eif5'].X[eif5_single_boolean]\n",
    "    tcf25_single_pos = adata_adip[:,'Tcf25'].X[tcf25_single_boolean]\n",
    "\n",
    "    print('Bin1 Single Positive ', len(bin1_single_pos))\n",
    "    print('eif5 Single Positive ', len(eif5_single_pos))\n",
    "    print('tcf25 Single Positive ', len(tcf25_single_pos))\n",
    "\n",
    "    bin1_eif5_double_pos = adata_adip[:,'Bin1'].X[bin1_eif5_double_boolean]\n",
    "    bin1_tcf25_double_pos = adata_adip[:,'Eif5'].X[bin1_tcf25_double_boolean]\n",
    "    tcf25_eif5_double_pos = adata_adip[:,'Tcf25'].X[tcf25_eif5_double_boolean]\n",
    "\n",
    "    print('Bin1/Eif5 Double Positive ', len(bin1_eif5_double_pos))\n",
    "    print('Bin1/Tcf25 Double Positive ', len(bin1_tcf25_double_pos))\n",
    "    print('Tcf25/Eif5 Double Positive ', len(tcf25_eif5_double_pos))\n",
    "\n",
    "    triple_pos = adata_adip[:,'Bin1'].X[triple_boolean]\n",
    "    print('Triple Positive ', len(triple_pos))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define observations for single, double and triple positive counts of bin1, eif5 and tcf25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_recomp == True:\n",
    "    adata_adip.obs['sdt_pos'] = np.array((non_boolean_int + single_boolean_int + double_boolean_int + triple_boolean_int), dtype=str)\n",
    "    adata_adip.obs['s_pos'] = np.array((np.array(bin1_single_boolean, dtype=int)*1)+(np.array(eif5_single_boolean, dtype=int)*2)+(np.array(tcf25_single_boolean, dtype=int)*3), dtype=str)\n",
    "    adata_adip.obs['d_pos'] = np.array((np.array(bin1_eif5_double_boolean, dtype=int)*1)+(np.array(bin1_tcf25_double_boolean, dtype=int)*2)+(np.array(tcf25_eif5_double_boolean, dtype=int)*3), dtype=str)\n",
    "    adata_adip.obs['t_pos'] = np.array(triple_boolean, dtype=str)\n",
    "    # make them categorical\n",
    "    adata_adip.obs['sdt_pos'] = pd.Series(adata_adip.obs['sdt_pos'], dtype=\"category\")\n",
    "    adata_adip.obs['s_pos'] = pd.Series(adata_adip.obs['s_pos'], dtype=\"category\")\n",
    "    adata_adip.obs['d_pos'] = pd.Series(adata_adip.obs['d_pos'], dtype=\"category\")\n",
    "    adata_adip.obs['t_pos'] = pd.Series(adata_adip.obs['t_pos'], dtype=\"category\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    new_cluster_names = ['Non Positives','Single Positives','Double Positives','Triple Positives']\n",
    "    adata_adip.rename_categories('sdt_pos',\n",
    "                                 new_cluster_names\n",
    "    )\n",
    "    sc.pl.umap(adata_adip,\n",
    "               color=['sdt_pos'],\n",
    "               size=20, \n",
    "               save=\"_non-single-double-triple_positive.svg\"\n",
    "    )\n",
    "    new_cluster_names = ['Non Single Positives','Bin1 Single Positives','Eif5 Single Positives','Tcf25 Single Positives']\n",
    "    adata_adip.rename_categories('s_pos',\n",
    "                                 new_cluster_names\n",
    "    )\n",
    "    sc.pl.umap(adata_adip[adata_adip.obs['s_pos']!='Non Single Positives'],\n",
    "               color=['s_pos'],\n",
    "               size=20,\n",
    "               save=\"_Bin1_Eif5_Tcf25_single_positive.svg\"\n",
    "    )\n",
    "    new_cluster_names = ['Non Double Positives','Bin1/Eif5 Double Positive','Bin1/Tcf25 Double Positive','Tcf25/Eif5 Double Positive']\n",
    "    adata_adip.rename_categories('d_pos', \n",
    "                                 new_cluster_names\n",
    "    )\n",
    "    sc.pl.umap(adata_adip[adata_adip.obs['d_pos']!='Non Double Positives'], \n",
    "               color=['d_pos'],\n",
    "               size=20,\n",
    "               save=\"_Bin1Eif5-Bin1Tcf25-Tcf25Eif5_double_positive.svg\"\n",
    "    )\n",
    "    new_cluster_names = ['Non Triple Positives','Triple Positives']\n",
    "    adata_adip.rename_categories('t_pos', \n",
    "                                 new_cluster_names\n",
    "    )\n",
    "    sc.pl.umap(adata_adip[adata_adip.obs['t_pos']!='Non Triple Positives'],\n",
    "               color=['t_pos'],\n",
    "               size=20,\n",
    "               save=\"_triple_positive.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Single positives are found in all 4 clusters. Double positives are mostly in the clusters 0,1 and 3 (preadipocyte clusters) and far less present in cluster 2 (mature adipocyte cluster). Triple positives are found only in clusters 0,1 and 3 with only one exception."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis separately for preadipocytes and mature adipocytes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/pandas/core/arrays/categorical.py:2487: FutureWarning: The `inplace` parameter in pandas.Categorical.remove_unused_categories is deprecated and will be removed in a future version.\n",
      "  res = method(*args, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "adata_adip.obs['louvain'].head()\n",
    "adata_preadip = adata_adip[adata_adip.obs['louvain'] != '2']\n",
    "adata_matadip = adata_adip[adata_adip.obs['louvain'] == '2']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Preadipocyte clusters (0,1,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.umap(adata_preadip,\n",
    "               color=['sdt_pos'],\n",
    "               size=20\n",
    "    )\n",
    "    sc.pl.umap(adata_preadip[adata_preadip.obs['s_pos'] != 'Non Single Positives'],\n",
    "               color=['s_pos'],\n",
    "               size=20\n",
    "    )\n",
    "    sc.pl.umap(adata_preadip[adata_preadip.obs['d_pos'] != 'Non Double Positives'],\n",
    "               color=['d_pos'],\n",
    "               size=20\n",
    "    )\n",
    "    sc.pl.umap(adata_preadip[adata_preadip.obs['t_pos'] != 'Non Triple Positives'],\n",
    "               color=['t_pos'], \n",
    "               size=20\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    import matplotlib_venn\n",
    "\n",
    "    bin1_pos = adata_preadip.obs_names[np.asarray(adata_preadip[:,'Bin1'].X).flatten() > 0]\n",
    "    eif5_pos = adata_preadip.obs_names[np.asarray(adata_preadip[:,'Eif5'].X).flatten() > 0]\n",
    "    tcf25_pos = adata_preadip.obs_names[np.asarray(adata_preadip[:,'Tcf25'].X).flatten() > 0]\n",
    "\n",
    "    matplotlib_venn.venn3([set(bin1_pos),\n",
    "                           set(eif5_pos),\n",
    "                           set(tcf25_pos)],\n",
    "                          set_labels=(\"Bin1\", \"Eif5\", \"Tcf25\"))\n",
    "    plt.savefig(dir_out + \"panels/venn_target_positives_adipocytes.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mature adipocyte cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.umap(adata_matadip, \n",
    "               color=['sdt_pos'],\n",
    "               size=20\n",
    "    )\n",
    "    sc.pl.umap(adata_matadip[adata_matadip.obs['s_pos'] != 'Non Single Positives'],\n",
    "               color=['s_pos'], \n",
    "               size=20\n",
    "    )\n",
    "    sc.pl.umap(adata_matadip[adata_matadip.obs['d_pos'] != 'Non Double Positives'],\n",
    "               color=['d_pos'],\n",
    "               size=20\n",
    "    )\n",
    "    sc.pl.umap(adata_matadip[adata_matadip.obs['t_pos'] != 'Non Triple Positives'], \n",
    "               color=['t_pos'],\n",
    "               size=20\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"means\"></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    \n",
    "\n",
    "    bin1_pos = adata_matadip.obs_names[np.asarray(adata_matadip[:,'Bin1'].X).flatten() > 0]\n",
    "    eif5_pos = adata_matadip.obs_names[np.asarray(adata_matadip[:,'Eif5'].X).flatten() > 0]\n",
    "    tcf25_pos = adata_matadip.obs_names[np.asarray(adata_matadip[:,'Tcf25'].X).flatten() > 0]\n",
    "\n",
    "    matplotlib_venn.venn3([\n",
    "        set(bin1_pos),\n",
    "        set(eif5_pos),\n",
    "        set(tcf25_pos)\n",
    "    ], set_labels = (\"Bin1\", \"Eif5\", \"Tcf25\"))\n",
    "    plt.savefig(dir_out + \"panels/venn_target_positives_mature_adipocytes.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean Counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normalized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bin1 mean counts:  0.23998265\n",
      "Eif5 mean counts:  0.6047478\n",
      "Tcf25 mean counts:  0.37402338\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/pandas/core/arrays/categorical.py:2487: FutureWarning: The `inplace` parameter in pandas.Categorical.remove_unused_categories is deprecated and will be removed in a future version.\n",
      "  res = method(*args, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "print('Bin1 mean counts: ',np.mean(adata_adip[:,'Bin1'].X))\n",
    "print('Eif5 mean counts: ',np.mean(adata_adip[:,'Eif5'].X))\n",
    "print('Tcf25 mean counts: ',np.mean(adata_adip[:,'Tcf25'].X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"traject\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preadipocytes only"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/anndata/compat/__init__.py:179: FutureWarning: Moving element from .uns['neighbors']['distances'] to .obsp['distances'].\n",
      "\n",
      "This is where adjacency matrices should go now.\n",
      "  warn(\n",
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/anndata/compat/__init__.py:179: FutureWarning: Moving element from .uns['neighbors']['connectivities'] to .obsp['connectivities'].\n",
      "\n",
      "This is where adjacency matrices should go now.\n",
      "  warn(\n",
      "running PAGA\n",
      "    finished: added\n",
      "    'paga/connectivities', connectivities adjacency (adata.uns)\n",
      "    'paga/connectivities_tree', connectivities subtree (adata.uns) (0:00:00)\n"
     ]
    }
   ],
   "source": [
    "if bool_recomp == True:\n",
    "    cell_ids_adip2 = np.asarray(adata_proc.obs_names)[\n",
    "        [x in ['preadipocytes'] \n",
    "         for x in np.asarray(adata_proc.obs['celltypes'].values)]\n",
    "    ]\n",
    "    adata_adip2 = adata_raw[cell_ids_adip2,:].copy()\n",
    "    adata_adip2.obs['n_genes'] = (adata_adip2.X > 0).sum(1)\n",
    "    temp_n_counts2 = adata_adip2.obs['n_counts']\n",
    "    adata_adip2.raw = adata_adip2\n",
    "    adata_adip2.obs['louvain_old'] = adata_proc[cell_ids_adip2, :].obs['louvain']\n",
    "    sc.pp.normalize_per_cell(adata_adip2)\n",
    "    sc.pp.log1p(adata_adip2)\n",
    "    sc.pp.highly_variable_genes(adata_adip2,\n",
    "                                n_top_genes=4000\n",
    "    )\n",
    "    sc.pl.highly_variable_genes(adata_adip2)\n",
    "    sc.pp.pca(adata_adip2,\n",
    "              n_comps=50,\n",
    "              use_highly_variable=True,\n",
    "              random_state=0, \n",
    "              svd_solver='arpack'\n",
    "    )\n",
    "    sc.pp.neighbors(adata_adip2,\n",
    "                    n_neighbors=100,\n",
    "                    knn=True,\n",
    "                    method='umap',\n",
    "                    n_pcs=50, \n",
    "                    random_state=0\n",
    "    )\n",
    "    sc.tl.umap(adata_adip2)\n",
    "    if bool_recluster == True:\n",
    "        sc.tl.louvain(adata_adip2,\n",
    "                      resolution=1, \n",
    "                      flavor='vtraag',\n",
    "                      random_state=0\n",
    "        )\n",
    "        pd.DataFrame(adata_adip2.obs).to_csv(path_or_buf=dir_adata+'obs_adata_adip2.csv')\n",
    "    else:\n",
    "        obs2 = pd.read_csv(dir_adata+'obs_adata_adip2.csv')\n",
    "        adata_adip2.obs['louvain'] = pd.Series(obs2['louvain'].values,\n",
    "                                               dtype='category'\n",
    "        )\n",
    "    sc.write(dir_adata+'adata_adip2.h5ad',adata_adip2)\n",
    "else:\n",
    "    adata_adip2 = sc.read(dir_adata+'adata_adip2.h5ad') \n",
    "sc.tl.paga(adata_adip2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.umap(adata_adip2, \n",
    "               color=['louvain_old', 'louvain'],\n",
    "               size=30,\n",
    "               save=\"_preadipocytes_louvain_0.svg\"\n",
    "    )\n",
    "    sc.pl.paga(adata_adip2, \n",
    "               color=['louvain'],\n",
    "               save=\"_preadipocytes_louvain.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_umap_marker(adata_adip2, \n",
    "                     differentialtion_markers, \n",
    "                     size=50,\n",
    "                     save=\"_differentialtion_markers_preadipcytes\",\n",
    "                     use_raw=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_umap_marker(adata_adip2, \n",
    "                     adipocyte_markers,\n",
    "                     size=50,\n",
    "                     save=\"_preadip_markers_preadipcytes\", \n",
    "                     use_raw=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_adip,\n",
    "                       differentialtion_markers.tolist(),\n",
    "                       save=\"_differentialtion_markers_preadipcytes\",\n",
    "                       use_raw=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.heatmap(adata=adata_adip2, \n",
    "                  var_names=differentialtion_markers, \n",
    "                  groupby=\"louvain\", \n",
    "                  use_raw=False, \n",
    "                  log=False, \n",
    "                  dendrogram=True, \n",
    "                  var_group_rotation=90, \n",
    "                  show_gene_labels=True, \n",
    "                  show=True, \n",
    "                  save=\"_differentialtion_markers_celltypes_preadipocytes.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.heatmap(adata=adata_adip2, \n",
    "                  var_names=adipocyte_markers, \n",
    "                  groupby=\"louvain\", \n",
    "                  use_raw=False, \n",
    "                  log=False, \n",
    "                  dendrogram=True, \n",
    "                  var_group_rotation=90, \n",
    "                  show_gene_labels=True, \n",
    "                  show=True, \n",
    "                  save=\"_preadipocyte_markers_celltypes_preadipocytes.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Default of the method has been changed to 't-test' from 't-test_overestim_var'\n",
      "ranking genes\n",
      "    finished: added to `.uns['rank_genes']`\n",
      "    'names', sorted np.recarray to be indexed by group ids\n",
      "    'scores', sorted np.recarray to be indexed by group ids\n",
      "    'logfoldchanges', sorted np.recarray to be indexed by group ids\n",
      "    'pvals', sorted np.recarray to be indexed by group ids\n",
      "    'pvals_adj', sorted np.recarray to be indexed by group ids (0:00:00)\n"
     ]
    }
   ],
   "source": [
    "sc.tl.rank_genes_groups(adata_adip2, \n",
    "                        groupby='louvain', \n",
    "                        n_genes=4000,\n",
    "                        key_added='rank_genes'\n",
    ")\n",
    "if bool_plot == True:\n",
    "    sc.pl.rank_genes_groups(adata_adip2,\n",
    "                            key='rank_genes',\n",
    "                            groups=['0','1','2','3'],\n",
    "                            fontsize=12,\n",
    "                            save=\"_preadipocytes.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = adata_adip2.uns['rank_genes']\n",
    "groups = result['names'].dtype.names\n",
    "diffExp = pd.DataFrame(\n",
    "    {group + '_' + key[:1]: result[key][group]\n",
    "    for group in groups for key in ['names', 'pvals_adj','logfoldchanges']}).head(4000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "diffExp.columns = ['0_name','0_pval_ajd','0_logFC','1_name','1_pval_ajd','1_logFC',\n",
    "                   '2_name','2_pval_ajd','2_logFC','3_name','3_pval_ajd','3_logFC']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save a DE table\n",
    "diffExp.to_csv(path_or_buf=dir_tables+\"DE_preadipocytes_by_cluster.csv\", sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp = GProfiler(return_dataframe=True)\n",
    "cluster_0 = gp.profile(organism='mmusculus',\n",
    "                       query=diffExp['0_name'][1:676].tolist()\n",
    ")\n",
    "cluster_1 = gp.profile(organism='mmusculus', \n",
    "                       query=diffExp['1_name'][1:93].tolist()\n",
    ")\n",
    "cluster_2 = gp.profile(organism='mmusculus', \n",
    "                       query=diffExp['2_name'][1:467].tolist()\n",
    ")\n",
    "cluster_3 = gp.profile(organism='mmusculus', \n",
    "                       query=diffExp['3_name'][1:129].tolist()\n",
    ")\n",
    "#cluster_0.sort_values(by='recall', ascending=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_0.to_csv(path_or_buf=dir_tables+\"enrichment_DE_genes_cluster_0.csv\", sep=\"\\t\")\n",
    "cluster_1.to_csv(path_or_buf=dir_tables+\"enrichment_DE_genes_cluster_1.csv\", sep=\"\\t\")\n",
    "cluster_2.to_csv(path_or_buf=dir_tables+\"enrichment_DE_genes_cluster_2.csv\", sep=\"\\t\")\n",
    "cluster_3.to_csv(path_or_buf=dir_tables+\"enrichment_DE_genes_cluster_3.csv\", sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diffusion Pseudotime Adipocytes¶"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "computing neighbors\n",
      "    using 'X_pca' with n_pcs = 50\n",
      "    finished: added to `.uns['neighbors']`\n",
      "    `.obsp['distances']`, distances for each pair of neighbors\n",
      "    `.obsp['connectivities']`, weighted adjacency matrix (0:00:11)\n",
      "computing Diffusion Maps using n_comps=15(=n_dcs)\n",
      "computing transitions\n",
      "    finished (0:00:00)\n",
      "    eigenvalues of transition matrix\n",
      "    [1.         0.989646   0.98074245 0.95198727 0.9318196  0.9060198\n",
      "     0.8811513  0.877778   0.86690754 0.8443195  0.84319466 0.81539434\n",
      "     0.80603105 0.779319   0.7653546 ]\n",
      "    finished: added\n",
      "    'X_diffmap', diffmap coordinates (adata.obsm)\n",
      "    'diffmap_evals', eigenvalues of transition matrix (adata.uns) (0:00:00)\n",
      "computing neighbors\n",
      "    using 'X_pca' with n_pcs = 50\n",
      "    finished: added to `.uns['neighbors']`\n",
      "    `.obsp['distances']`, distances for each pair of neighbors\n",
      "    `.obsp['connectivities']`, weighted adjacency matrix (0:00:00)\n",
      "computing Diffusion Maps using n_comps=15(=n_dcs)\n",
      "computing transitions\n",
      "    finished (0:00:00)\n",
      "    eigenvalues of transition matrix\n",
      "    [1.         0.989646   0.98074245 0.95198727 0.9318196  0.9060198\n",
      "     0.8811513  0.877778   0.86690754 0.8443195  0.84319466 0.81539434\n",
      "     0.80603105 0.779319   0.7653546 ]\n",
      "    finished: added\n",
      "    'X_diffmap', diffmap coordinates (adata.obsm)\n",
      "    'diffmap_evals', eigenvalues of transition matrix (adata.uns) (0:00:00)\n"
     ]
    }
   ],
   "source": [
    "adata_adip_raw = adata_adip.copy()\n",
    "adata_adip_raw.X = adata_adip.raw.X.copy()\n",
    "sc.pp.log1p(adata_adip_raw)\n",
    "sc.pp.neighbors(adata_adip_raw)\n",
    "sc.tl.diffmap(adata_adip_raw)\n",
    "\n",
    "sc.pp.neighbors(adata_adip)\n",
    "sc.tl.diffmap(adata_adip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_diffmap_marker(adata, markers, save=None, use_raw=True):\n",
    "    for i in range(len(markers) // 2 + len(markers) % 2):\n",
    "        print(markers[(2*i):np.min([2*(i+1), len(markers)])])\n",
    "        if save is not None:\n",
    "            sc.pl.diffmap(adata, \n",
    "                          components = '1,2',\n",
    "                          color=markers[(2*i):np.min([2*(i+1), len(markers)])], \n",
    "                          use_raw=use_raw,\n",
    "                          color_map=mymap,\n",
    "                          save=save+\"_\"+str(i)+\".svg\"\n",
    "            )\n",
    "        else:\n",
    "            sc.pl.diffmap(adata, \n",
    "                          components,\n",
    "                          color=markers[(2*i):np.min([2*(i+1), len(markers)])], \n",
    "                          use_raw=use_raw,\n",
    "                          color_map=mymap\n",
    "            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.diffmap(adata_adip,\n",
    "                  components='1,2', \n",
    "                  color='louvain', \n",
    "                  save=\"_adipocytes_pseudotime_louvain_DC_1-2_norm.svg\"\n",
    "    )\n",
    "    sc.pl.diffmap(adata_adip,\n",
    "                  components='1,3',\n",
    "                  color='louvain',\n",
    "                  save=\"_adipocytes_pseudotime_louvain_DC_1-3_norm.svg\"\n",
    "    )\n",
    "    sc.pl.diffmap(adata_adip_raw, \n",
    "                  components='1,2', \n",
    "                  color='louvain', \n",
    "                  save=\"_adipocytes_pseudotime_louvain_DC_1-2.svg\"\n",
    "    )\n",
    "    sc.pl.diffmap(adata_adip_raw,\n",
    "                  components='1,3',\n",
    "                  color='louvain',\n",
    "                  save=\"_adipocytes_pseudotime_louvain_DC_1-3.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "computing Diffusion Pseudotime using n_dcs=10\n",
      "    finished: added\n",
      "    'dpt_pseudotime', the pseudotime (adata.obs) (0:00:00)\n",
      "computing Diffusion Pseudotime using n_dcs=10\n",
      "    finished: added\n",
      "    'dpt_pseudotime', the pseudotime (adata.obs) (0:00:00)\n"
     ]
    }
   ],
   "source": [
    "adipocyte_mask21 = np.isin(adata_adip.obs['louvain'], '0')\n",
    "max_adipocyte_id21 = np.argmax(adata_adip.obsm['X_diffmap'][adipocyte_mask21, 2])\n",
    "max_adipocyte_id2_raw1 = np.argmax(adata_adip.obsm['X_diffmap'][adipocyte_mask21, 2])\n",
    "root_id21 = np.arange(len(adipocyte_mask21))[adipocyte_mask21][max_adipocyte_id21]\n",
    "adata_adip.uns['iroot'] = root_id21\n",
    "adata_adip_raw.uns['iroot'] = root_id21\n",
    "\n",
    "#Compute dpt\n",
    "sc.tl.dpt(adata_adip,\n",
    "          n_branchings=0\n",
    ")\n",
    "sc.tl.dpt(adata_adip_raw, \n",
    "          n_branchings=0\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_diffmap_marker(adata_adip,\n",
    "                        adipocyte_markers, \n",
    "                        save=\"_adipocytes_pseudotime_adipocyte_markers_DC_1-2_norm\"\n",
    "    )\n",
    "    plot_diffmap_marker(adata_adip_raw,\n",
    "                        adipocyte_markers,\n",
    "                        save=\"_adipocytes_pseudotime_adipocyte_markers_DC_1-2\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Visualize pseudotime over differentiation\n",
    "if bool_plot == True:\n",
    "    sc.pl.diffmap(adata_adip, \n",
    "                  components='1,2',\n",
    "                  color='dpt_pseudotime',\n",
    "                  color_map=mymap, \n",
    "                  save=\"_adipocytes_dpt_pseudotime_DC_1-2_norm.svg\"\n",
    "    )\n",
    "    sc.pl.diffmap(adata_adip,\n",
    "                  components='1,3', \n",
    "                  color='dpt_pseudotime', \n",
    "                  color_map=mymap, \n",
    "                  save=\"_adipocytes_dpt_pseudotime_DC_1-3_norm.svg\"\n",
    "    )\n",
    "    sc.pl.diffmap(adata_adip_raw,\n",
    "                  components='1,2',\n",
    "                  color='dpt_pseudotime',\n",
    "                  color_map=mymap,\n",
    "                  save=\"_adipocytes_dpt_pseudotime_DC_1-2.svg\"\n",
    "    )\n",
    "    sc.pl.diffmap(adata_adip_raw, \n",
    "                  components='1,3',\n",
    "                  color='dpt_pseudotime', \n",
    "                  color_map=mymap,\n",
    "                  save=\"_adipocytes_dpt_pseudotime_DC_1-3.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diffusion Pseudotime Preadipocytes¶"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "computing neighbors\n",
      "    using 'X_pca' with n_pcs = 50\n",
      "    finished: added to `.uns['neighbors']`\n",
      "    `.obsp['distances']`, distances for each pair of neighbors\n",
      "    `.obsp['connectivities']`, weighted adjacency matrix (0:00:00)\n",
      "computing Diffusion Maps using n_comps=15(=n_dcs)\n",
      "computing transitions\n",
      "    finished (0:00:00)\n",
      "    eigenvalues of transition matrix\n",
      "    [1.         0.981316   0.9460798  0.91849977 0.89427274 0.86765224\n",
      "     0.8456241  0.82347286 0.80268925 0.78898954 0.76468396 0.7550992\n",
      "     0.7451544  0.7367611  0.7321832 ]\n",
      "    finished: added\n",
      "    'X_diffmap', diffmap coordinates (adata.obsm)\n",
      "    'diffmap_evals', eigenvalues of transition matrix (adata.uns) (0:00:00)\n",
      "computing neighbors\n",
      "    using 'X_pca' with n_pcs = 50\n",
      "    finished: added to `.uns['neighbors']`\n",
      "    `.obsp['distances']`, distances for each pair of neighbors\n",
      "    `.obsp['connectivities']`, weighted adjacency matrix (0:00:00)\n",
      "computing Diffusion Maps using n_comps=15(=n_dcs)\n",
      "computing transitions\n",
      "    finished (0:00:00)\n",
      "    eigenvalues of transition matrix\n",
      "    [1.         0.981316   0.9460798  0.91849977 0.89427274 0.86765224\n",
      "     0.8456241  0.82347286 0.80268925 0.78898954 0.76468396 0.7550992\n",
      "     0.7451544  0.7367611  0.7321832 ]\n",
      "    finished: added\n",
      "    'X_diffmap', diffmap coordinates (adata.obsm)\n",
      "    'diffmap_evals', eigenvalues of transition matrix (adata.uns) (0:00:00)\n"
     ]
    }
   ],
   "source": [
    "adata_adip2_raw = adata_adip2.copy()\n",
    "adata_adip2_raw.X = adata_adip2.raw.X.copy()\n",
    "sc.pp.log1p(adata_adip2_raw)\n",
    "sc.pp.neighbors(adata_adip2_raw)\n",
    "sc.tl.diffmap(adata_adip2_raw)\n",
    "\n",
    "sc.pp.neighbors(adata_adip2)\n",
    "sc.tl.diffmap(adata_adip2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.diffmap(adata_adip2, \n",
    "                  components='1,2',\n",
    "                  color='louvain',\n",
    "                  save=\"_preadipocytes_pseudotime_louvain_DC_1-2_norm.svg\"\n",
    "    )\n",
    "    sc.pl.diffmap(adata_adip2,\n",
    "                  components='1,3', \n",
    "                  color='louvain', \n",
    "                  save=\"_preadipocytes_pseudotime_louvain_DC_1-3_norm.svg\"\n",
    "    )\n",
    "    sc.pl.diffmap(adata_adip2_raw,\n",
    "                  components='1,2',\n",
    "                  color='louvain', \n",
    "                  save=\"_preadipocytes_pseudotime_louvain_DC_1-2.svg\"\n",
    "    )\n",
    "    sc.pl.diffmap(adata_adip2_raw, \n",
    "                  components='1,3', \n",
    "                  color='louvain',\n",
    "                  save=\"_preadipocytes_pseudotime_louvain_DC_1-3.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "computing Diffusion Pseudotime using n_dcs=10\n",
      "    finished: added\n",
      "    'dpt_pseudotime', the pseudotime (adata.obs) (0:00:00)\n",
      "computing Diffusion Pseudotime using n_dcs=10\n",
      "    finished: added\n",
      "    'dpt_pseudotime', the pseudotime (adata.obs) (0:00:00)\n"
     ]
    }
   ],
   "source": [
    "adipocyte_mask2 = np.isin(adata_adip2.obs['louvain'], '0')\n",
    "max_adipocyte_id2 = np.argmax(adata_adip2.obsm['X_diffmap'][adipocyte_mask2, 2])\n",
    "max_adipocyte_id2_raw = np.argmax(adata_adip2.obsm['X_diffmap'][adipocyte_mask2, 2])\n",
    "root_id2 = np.arange(len(adipocyte_mask2))[adipocyte_mask2][max_adipocyte_id2]\n",
    "adata_adip2.uns['iroot'] = root_id2\n",
    "adata_adip2_raw.uns['iroot'] = root_id2\n",
    "\n",
    "#Compute dpt\n",
    "sc.tl.dpt(adata_adip2,\n",
    "          n_branchings=0\n",
    ")\n",
    "sc.tl.dpt(adata_adip2_raw, \n",
    "          n_branchings=0\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_diffmap_marker(adata_adip2, \n",
    "                        adipocyte_markers, \n",
    "                        save=\"_preadipocytes_pseudotime_adipocyte_markers_DC_1-2_norm\"\n",
    "    )\n",
    "    plot_diffmap_marker(adata_adip2_raw,\n",
    "                        adipocyte_markers, \n",
    "                        save=\"_preadipocytes_pseudotime_adipocyte_markers_DC_1-2\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Visualize pseudotime over differentiation\n",
    "if bool_plot == True:\n",
    "    sc.pl.diffmap(adata_adip2,\n",
    "                  components='1,2',\n",
    "                  color='dpt_pseudotime', \n",
    "                  color_map=mymap,\n",
    "                  save=\"_preadipocytes_dpt_pseudotime_DC_1-2_norm.svg\"\n",
    "    )\n",
    "    sc.pl.diffmap(adata_adip2, \n",
    "                  components='1,3', \n",
    "                  color='dpt_pseudotime',\n",
    "                  color_map=mymap, \n",
    "                  save=\"_preadipocytes_dpt_pseudotime_DC_1-3_norm.svg\"\n",
    "    )\n",
    "    sc.pl.diffmap(adata_adip2_raw,\n",
    "                  components='1,2',\n",
    "                  color='dpt_pseudotime', \n",
    "                  color_map=mymap, \n",
    "                  save=\"_preadipocytes_dpt_pseudotime_DC_1-2.svg\"\n",
    "    )\n",
    "    sc.pl.diffmap(adata_adip2_raw, \n",
    "                  components='1,3', \n",
    "                  color='dpt_pseudotime', \n",
    "                  color_map=mymap,\n",
    "                  save=\"_preadipocytes_dpt_pseudotime_DC_1-3.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Visualize pseudotime over differentiation\n",
    "if bool_plot == True:\n",
    "    sc.pl.diffmap(adata_adip2_raw, \n",
    "                  components='1,2',\n",
    "                  color=['Pparg','Cebpa','Cebpb','Cebpd','Dlk1','Fabp4','Cd34','louvain'],\n",
    "                  color_map=mymap\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_adip2, \n",
    "                       ['Pparg','Cebpa','Cebpb','Cebpd','Dlk1','Fabp4','Cd34'],\n",
    "                       use_raw=False,\n",
    "                       save=\"_Nature-rebuttal_differentialtion_markers_preadipcytes.svg\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_adip, \n",
    "                       ['Pparg','Cebpa','Cebpb','Cebpd','Dlk1','Fabp4'],\n",
    "                       use_raw=False, \n",
    "                       save=\"_Nature-rebuttal_differentialtion_markers_adipcytes.svg\"\n",
    "    )"
   ]
  }
 ],
 "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.9.2"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}