notebooks/asc1_white_vs_beige_adipocyte_fate.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary of notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we perform the following analysis:\n",
    "(1) Clustering of the cells of young and old mice, annotation of the clusters using marker genes of immune cells and adipocytes. \n",
    "(2) Separate the preadipocytes from the immune cells and identify the mature subpopulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load packages and set global variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "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",
      "autoreload          NA\n",
      "backcall            0.2.0\n",
      "batchglm            v0.7.4\n",
      "cffi                1.14.5\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",
      "diffxpy             v0.7.4\n",
      "get_version         2.1\n",
      "h5py                3.2.1\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",
      "patsy               0.5.1\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",
      "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",
      "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:08\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": [
    "import numpy as np\n",
    "import scanpy as sc\n",
    "import scipy as sci\n",
    "import scipy.sparse\n",
    "import re\n",
    "import pandas\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import sys\n",
    "\n",
    "from matplotlib import colors\n",
    "\n",
    "import batchglm\n",
    "import diffxpy.api as de\n",
    "\n",
    "#from beakerx import *\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "sc.settings.verbosity = 3 # amount of output\n",
    "dir_in = '/Users/viktorian.miok/Documents/consultation/Altun-Ussar/David/data/'\n",
    "dir_out = '/Users/viktorian.miok/Documents/consultation/Altun-Ussar/David/results/'\n",
    "dir_tables = dir_out+'tables/'\n",
    "sc_settings_figdir = dir_out+'panels/'\n",
    "sc_settings_writedir = dir_out+'anndata/'\n",
    "sc.logging.print_versions()\n",
    "sc.settings.set_figure_params(dpi=80, scanpy=True)\n",
    "print (sys.version)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "v0.7.4\n"
     ]
    }
   ],
   "source": [
    "print(de.__version__)"
   ]
  },
  {
   "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 careful with over-writing 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": 6,
   "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": 7,
   "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": 8,
   "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": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "bool_plot = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp:\n",
    "    # Count matrix:\n",
    "    fn_cnts_young = dir_in+'Data S1 adolescent_filtered_gene_bc_matrices_h5.h5'\n",
    "    adata_young = sc.read_10x_h5(fn_cnts_young)\n",
    "    fn_cnts_old = dir_in+'Data S2 adult_filtered_gene_bc_matrices_h5.h5'\n",
    "    adata_old = sc.read_10x_h5(fn_cnts_old)\n",
    "    adata_raw = adata_young.concatenate([adata_old],\n",
    "                                        batch_key='age',\n",
    "                                        batch_categories=['young', 'old']\n",
    "    )\n",
    "    sc.write(sc_settings_writedir+'adata_raw.h5ad', adata_raw)\n",
    "else:\n",
    "    adata_raw = sc.read(sc_settings_writedir+'adata_raw.h5ad')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Process data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Embeddings and clustering"
   ]
  },
  {
   "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 t-SNE 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": 12,
   "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:\n",
    "    adata_proc = adata_raw.copy()\n",
    "    sc.pp.filter_cells(adata_proc,\n",
    "                       min_counts=500\n",
    "    )\n",
    "    sc.pp.normalize_per_cell(adata_proc)\n",
    "    adata_proc.raw = sc.pp.log1p(adata_proc,\n",
    "                                 copy=True\n",
    "    )\n",
    "    sc.pp.pca(adata_proc,\n",
    "              n_comps=50,\n",
    "              random_state=0, \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",
    "               n_jobs=3\n",
    "    )\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",
    "        pandas.DataFrame(adata_proc.obs).to_csv(\n",
    "            path_or_buf = sc_settings_writedir+\"obs_adata_proc.csv\")\n",
    "    else:\n",
    "        obs = pandas.read_csv(sc_settings_writedir+'obs_adata_proc.csv')\n",
    "        adata_proc.obs['louvain'] = pandas.Series(obs['louvain'].values,\n",
    "                                                  dtype='category'\n",
    "        )\n",
    "    sc.write(sc_settings_writedir+'adata_proc.h5ad', adata_proc)\n",
    "else:\n",
    "    adata_proc = sc.read(sc_settings_writedir+'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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define a nice colour map for gene expression\n",
    "colors2 = plt.cm.Reds(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": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot==True:\n",
    "    sc.pl.tsne(adata_proc,\n",
    "               color=['age'],\n",
    "               size=5,\n",
    "               save=\"_all_age.pdf\"\n",
    "    )\n",
    "    sc.pl.tsne(adata_proc,\n",
    "               color=['louvain'],\n",
    "               size=5,\n",
    "               save=\"_all_louvain.pdf\"\n",
    "    )\n",
    "    sc.pl.tsne(adata_proc,\n",
    "               color=['n_counts'],\n",
    "               size=5, \n",
    "               save=\"_all_n_counts.pdf\"\n",
    "    )\n",
    "    sc.pl.tsne(adata_proc, \n",
    "               color=['Pdgfra'],\n",
    "               size=5,\n",
    "               save=\"_all_Pdgfra.pdf\", \n",
    "               color_map=mymap\n",
    "    )\n",
    "    sc.pl.tsne(adata_proc, \n",
    "               color=['Slc7a10'],\n",
    "               size=5,\n",
    "               save=\"_all_Slc7a10.pdf\", \n",
    "               color_map=mymap\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.paga(adata_proc,\n",
    "               save=\"_all.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of cells in each sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3915\n"
     ]
    }
   ],
   "source": [
    "print(np.sum(adata_proc.obs[\"age\"].values == \"young\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5429\n"
     ]
    }
   ],
   "source": [
    "print(np.sum(adata_proc.obs[\"age\"].values == \"old\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0     1733\n",
       "1     1379\n",
       "2     1222\n",
       "3      948\n",
       "4      892\n",
       "5      816\n",
       "6      799\n",
       "7      518\n",
       "8      514\n",
       "9      421\n",
       "10     102\n",
       "Name: louvain, dtype: int64"
      ]
     },
     "execution_count": 18,
     "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": [
    "## Marker genes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define marker sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define surface marker sets for some of the expected cell types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Leukocyte markers:\n",
    "leukocyte_markers = ['Ptprc']\n",
    "tc_markers = ['Cd3d','Cd3e','Cd3g','Cd4'] \n",
    "nk_markers = ['Nkg7','Il2rb','Ncr1','Klrd1','Klrb1b','Klrb1f']\n",
    "myeloid_markers = ['Cd79a','Itgax','Itgam','Fcgr3','S100a8', 'S100a9']\n",
    "mp_markers = ['Adgre1','Lyz2']\n",
    "dc_markers = ['Cd74','Anpep' , 'Cd33', 'Cd80', 'Cd83', 'Cd86'] \n",
    "bc_markers = ['Cd19']\n",
    "adipocyte_markers = ['Pdgfra','Slc7a10','Pparg','Fermt2','Fbn1','Col4a1','Itgb1','Cd34','Cd24a','Dlk1','Slc7a10']\n",
    "megakaryocyte_markers = ['Ppbp']\n",
    "erythrocyte_markers = ['Gypa']\n",
    "go_adip_dev = ['Aacs','Acat1','Arid5b','Arrdc3','Atf2','Bbs4','Bdh1','Csf1','Dgat2','Dyrk1b','Ebf2','Amer1','Fto',\n",
    "               'Id2','Lep','Lrp5','Nampt','Oxct1','Paxip1','Pik3ca','Ppard','Ppargc1a','Rorc','Sh3pxd2b','Slc25a25',\n",
    "               'Sox8','Spg20','Tbl1xr1','Xbp1']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting routines for marker gene sets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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)+\".pdf\"\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",
    "                         save=\"dasdad\"\n",
    "            )\n",
    "        \n",
    "def plot_tsne_marker(adata, markers, size=5, 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.tsne(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)+\".pdf\"\n",
    "            )\n",
    "        else:\n",
    "            sc.pl.tsne(adata, \n",
    "                       color=markers[(2*i):np.min([2*(i+1), len(markers)])], \n",
    "                       size=size,\n",
    "                       use_raw=use_raw\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Leukocyte markers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc, \n",
    "                       leukocyte_markers, \n",
    "                       save=\"_all_markers_leukocyte\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cluster 1,6,7,8,9 and potentially 0,10 express leukocyte marker Ptprc. These clusters are further validated by leukocyte specific markers below. The remaining clusters 2,3,4 are investigated with non-leukocyte marker sets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Megakaryocyte markers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc,\n",
    "                       megakaryocyte_markers, \n",
    "                       save=\"_all_markers_megakaryocytes\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There do not seem to be many megakaryocytes in this data set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Erythrocyte markers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc, \n",
    "                       erythrocyte_markers, \n",
    "                       save=\"_all_markers_erythrcytes\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There do not seem to be many erythrocytes in this data set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preadipocyte markers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc, \n",
    "                       adipocyte_markers, \n",
    "                       save=\"_all_markers_preadipocytes\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cluster 2,3,4 express adipocyte markers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### T-cell markers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc,\n",
    "                       tc_markers, \n",
    "                       save=\"_all_markers_tcells\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Marker gene expression suggests that cluster 1 and 6 are T-cells, interestingly not Cd4+Cd8+ T-cells it seems as Cd4 expression is low."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Natural killer cell markers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc, \n",
    "                       nk_markers,\n",
    "                       save=\"_all_markers_nk\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cluster 1,6,7, have natural killer cell marker gene expression, cluster 1 also expresses Cd3 so it may contain yd-T-cells-?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Myeloid cell markers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc, \n",
    "                       myeloid_markers,\n",
    "                       save=\"_all_markers_myeloid\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cluster 5,8,9,10 express myeloid cell marker genes. Cluster 9 seems to have bimodal expression in S100a8 and S100a9 so it may need subclustering to subdivide cell types here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Macrophage markers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc,\n",
    "                       mp_markers,\n",
    "                       save=\"_all_markers_macrophages\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cluster 5,7,8,9 express macrophage markers, in line with the myeloid cell marker gene expression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dendritic cell markers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc, \n",
    "                       dc_markers, \n",
    "                       save=\"_all_markers_dendritic\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clusters 0,5,7,8,9 express dendritic cell markers, in line with myeloid marker gene expression, cluster 0,7 could be a non-myeloid dendritic cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### B-cell markers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_proc, \n",
    "                       bc_markers,\n",
    "                       save=\"_all_markers_bcells\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cluster 0 may contain B-cells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary heatmap to characterize cell types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select a few genes to summarize cell type assignments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_leukocyte_markers = ['Ptprc']\n",
    "selected_tc_markers = ['Cd3d','Cd3g']\n",
    "selected_nk_markers = ['Nkg7','Klrd1']\n",
    "selected_myeloid_markers = ['Fcgr3','S100a8']\n",
    "selected_mp_markers = ['Adgre1','Lyz2']\n",
    "selected_dc_markers = ['Cd74','Cd83'] \n",
    "selected_bc_markers = ['Cd19']\n",
    "selected_adipocyte_markers = ['Pdgfra','Fbn1','Col4a1','Cd34','Cd24a','Dlk1','Slc7a10'] \n",
    "selected_megakaryocyte_markers = ['Ppbp']\n",
    "selected_erythrocyte_markers = ['Gypa']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_cell_markers = selected_leukocyte_markers + \\\n",
    "selected_megakaryocyte_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": "markdown",
   "metadata": {},
   "source": [
    "Only keep markers that occur in data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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=True, \n",
    "                  log=False, \n",
    "                  dendrogram=False, \n",
    "                  var_group_rotation=90, \n",
    "                  show_gene_labels=True, \n",
    "                  show=True, \n",
    "                  save=\"_all_markers_celltypes.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preadipocytes only"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Embedding and clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "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)[[x in ['2','3','4'] \n",
    "                                                      for x in np.asarray(adata_proc.obs['louvain'].values)]]\n",
    "    adata_adip = adata_raw[cell_ids_adip,:].copy()\n",
    "    sc.pp.filter_cells(adata_adip,\n",
    "                       min_counts=500\n",
    "    )\n",
    "    sc.pp.normalize_per_cell(adata_adip)\n",
    "    adata_adip.raw = adata_adip.copy()\n",
    "    sc.pp.log1p(adata_adip)\n",
    "    sc.pp.pca(adata_adip,\n",
    "              n_comps=50,\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.tsne(adata_adip,\n",
    "               n_jobs=3\n",
    "    )\n",
    "    if bool_recluster == True:\n",
    "        sc.tl.louvain(adata_adip, \n",
    "                      resolution=1,\n",
    "                      flavor='vtraag', \n",
    "                      random_state=0\n",
    "        )\n",
    "        pandas.DataFrame(adata_adip.obs).to_csv(\n",
    "            path_or_buf = sc_settings_writedir+\"obs_adata_adip.csv\")\n",
    "    else:\n",
    "        obs = pandas.read_csv(sc_settings_writedir+'obs_adata_adip.csv')\n",
    "        adata_adip.obs['louvain'] = pandas.Series(obs['louvain'].values, \n",
    "                                                  dtype='category'\n",
    "        )\n",
    "    sc.write(sc_settings_writedir+'adata_adip.h5ad', adata_adip)\n",
    "else:\n",
    "    adata_adip = sc.read(sc_settings_writedir+'adata_adip.h5ad')\n",
    "sc.tl.paga(adata_adip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot==True:\n",
    "    sc.pl.tsne(adata_adip,\n",
    "               color=['age'],\n",
    "               size=20,\n",
    "               save=\"_preadip_age.pdf\"\n",
    "    )\n",
    "    sc.pl.tsne(adata_adip,\n",
    "               color=['louvain'],\n",
    "               size=20, \n",
    "               save=\"_preadip_louvain.pdf\"\n",
    "    )\n",
    "    sc.pl.tsne(adata_adip,\n",
    "               color=['n_counts'],\n",
    "               size=20,\n",
    "               save=\"_preadip_n_counts.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fairly equal count depth across preadipocytes. Potential batch effect or age effect, indistinguishable in this scenario."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.paga(adata_adip, \n",
    "               save=\"_preadip.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Marker gene sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_adip,\n",
    "                       adipocyte_markers,\n",
    "                       save=\"_preadip_markers_preadipcytes\",\n",
    "                       use_raw=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot==True:\n",
    "    plot_tsne_marker(adata_adip,\n",
    "                     adipocyte_markers, \n",
    "                     size=20, \n",
    "                     save=\"_preadip_markers_preadipcytes\", \n",
    "                     use_raw=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems that Slc7a10+ cells are in clusters 0,2,4, which directly correspond to the samples of the young mice with a few exceptions. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preadipocytes only by age"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Embedding and clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Young mouse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "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_young = np.asarray(adata_proc.obs_names)[\n",
    "        np.logical_and(np.asarray([x in ['2','3','4']\n",
    "                                   for x in np.asarray(adata_proc.obs['louvain'].values)]),\n",
    "                       np.asarray([x=='young'\n",
    "                                   for x in np.asarray(adata_proc.obs['age'].values)]))]\n",
    "    adata_adip_young = adata_raw[cell_ids_adip_young,:].copy()\n",
    "    sc.pp.filter_cells(adata_adip_young,\n",
    "                       min_counts=500\n",
    "    )\n",
    "    sc.pp.normalize_per_cell(adata_adip_young)\n",
    "    adata_adip_young.raw = adata_adip_young.copy()\n",
    "    sc.pp.log1p(adata_adip_young)\n",
    "    sc.pp.pca(adata_adip_young,\n",
    "              n_comps=50, \n",
    "              random_state=0,\n",
    "              svd_solver='arpack'\n",
    "    )\n",
    "    sc.pp.neighbors(adata_adip_young,\n",
    "                    n_neighbors=100,\n",
    "                    knn=True,\n",
    "                    method='umap',\n",
    "                    n_pcs=50,\n",
    "                    random_state=0\n",
    "    )\n",
    "    sc.tl.tsne(adata_adip_young,\n",
    "               n_jobs=3\n",
    "    )\n",
    "    if bool_recluster == True:\n",
    "        sc.tl.louvain(adata_adip_young, \n",
    "                      resolution=1,\n",
    "                      flavor='vtraag',\n",
    "                      random_state=0\n",
    "        )\n",
    "        pandas.DataFrame(adata_adip_young.obs).to_csv(\n",
    "            path_or_buf = sc_settings_writedir+\"obs_adata_adip_young.csv\")\n",
    "    else:\n",
    "        obs = pandas.read_csv(sc_settings_writedir+'obs_adata_adip_young.csv')\n",
    "        adata_adip_young.obs['louvain'] = pandas.Series(obs['louvain'].values, \n",
    "                                                        dtype='category'\n",
    "        )\n",
    "    sc.write(sc_settings_writedir+'adata_adip_young.h5ad', adata_adip_young)\n",
    "else:\n",
    "    adata_adip_young = sc.read(sc_settings_writedir+'adata_adip_young.h5ad')\n",
    "sc.tl.paga(adata_adip_young)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.tsne(adata_adip_young,\n",
    "               color=['louvain'],\n",
    "               size=20,\n",
    "               save=\"_preadipYoung_louvain.pdf\"\n",
    "    )\n",
    "    sc.pl.tsne(adata_adip_young, \n",
    "               color=['n_counts'],\n",
    "               size=20,\n",
    "               save=\"_preadipYoung_n_counts.pdf\"\n",
    "    )\n",
    "    sc.pl.paga(adata_adip_young, \n",
    "               save=\"_preadipYoung.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of preadipocytes observed in young mouse with scRNAseq:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1636\n"
     ]
    }
   ],
   "source": [
    "print(adata_adip_young.X.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Young mouse - coarse clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_adip_young_lowres = adata_adip_young.copy()\n",
    "if bool_recomp == True:\n",
    "    if True:\n",
    "        sc.tl.louvain(adata_adip_young_lowres, \n",
    "                      resolution=0.5, \n",
    "                      flavor='vtraag',\n",
    "                      random_state=0\n",
    "        )\n",
    "        pandas.DataFrame(adata_adip_young_lowres.obs).to_csv(\n",
    "            path_or_buf = sc_settings_writedir+\"obs_adata_adip_young_lowres.csv\")\n",
    "    else:\n",
    "        obs = pandas.read_csv(sc_settings_writedir+'obs_adata_adip_young_lowres.csv')\n",
    "        adata_adip_young_lowres.obs['louvain'] = pandas.Series(obs['louvain'].values, \n",
    "                                                               dtype='category'\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.tsne(adata_adip_young_lowres,\n",
    "               color=['louvain'],\n",
    "               size=20, \n",
    "               save=\"_preadipYoung_louvain_lowres.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    503\n",
       "1    386\n",
       "2    373\n",
       "3    306\n",
       "4     68\n",
       "Name: louvain, dtype: int64"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adata_adip_young_lowres.obs['louvain'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Old mouse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "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_old = np.asarray(adata_proc.obs_names)[\n",
    "        np.logical_and(np.asarray([x in ['2','3','4']\n",
    "                                   for x in np.asarray(adata_proc.obs['louvain'].values)]),\n",
    "                       np.asarray([x=='old'\n",
    "                                   for x in np.asarray(adata_proc.obs['age'].values)]))]\n",
    "    adata_adip_old = adata_raw[cell_ids_adip_old,:].copy()\n",
    "    sc.pp.filter_cells(adata_adip_old,\n",
    "                       min_counts=500\n",
    "    )\n",
    "    sc.pp.normalize_per_cell(adata_adip_old)\n",
    "    adata_adip_old.raw = adata_adip_old.copy()\n",
    "    sc.pp.log1p(adata_adip_old)\n",
    "    sc.pp.pca(adata_adip_old,\n",
    "              n_comps=50,\n",
    "              random_state=0,\n",
    "              svd_solver='arpack'\n",
    "    )\n",
    "    sc.pp.neighbors(adata_adip_old,\n",
    "                    n_neighbors=50,\n",
    "                    knn=True,\n",
    "                    method='umap',\n",
    "                    n_pcs=50,\n",
    "                    random_state=0\n",
    "    )\n",
    "    sc.tl.tsne(adata_adip_old, \n",
    "               n_jobs=3\n",
    "    )\n",
    "    if bool_recluster == True:\n",
    "        sc.tl.louvain(adata_adip_old, \n",
    "                      resolution=1,\n",
    "                      flavor='vtraag',\n",
    "                      random_state=0\n",
    "        )\n",
    "        pandas.DataFrame(adata_adip_old.obs).to_csv(\n",
    "            path_or_buf = sc_settings_writedir+\"obs_adata_adip_old.csv\")\n",
    "    else:\n",
    "        obs = pandas.read_csv(sc_settings_writedir+'obs_adata_adip_old.csv')\n",
    "        adata_adip_old.obs['louvain'] = pandas.Series(obs['louvain'].values, \n",
    "                                                      dtype='category'\n",
    "        )\n",
    "    sc.write(sc_settings_writedir+'adata_adip_old.h5ad', adata_adip_old)\n",
    "else:\n",
    "    adata_adip_old = sc.read(sc_settings_writedir+'adata_adip_old.h5ad')\n",
    "sc.tl.paga(adata_adip_old)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.tsne(adata_adip_old,\n",
    "               color=['louvain'], \n",
    "               size=20,\n",
    "               save=\"_preadipOld_louvain.pdf\"\n",
    "    )\n",
    "    sc.pl.tsne(adata_adip_old, \n",
    "               color=['n_counts'],\n",
    "               size=20, \n",
    "               save=\"_preadipOld_n_counts.pdf\"\n",
    "    )\n",
    "    sc.pl.paga(adata_adip_old,\n",
    "               save=\"_preadipOld.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of preadipocytes observed in old mouse with scRNAseq:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1426\n"
     ]
    }
   ],
   "source": [
    "print(adata_adip_old.X.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Marker gene sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Young mouse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "adipocyte_markers = ['Pdgfra','Slc7a10','Pparg','Fermt2','Fbn1','Col4a1',\n",
    "                     'Itgb1','Cd34','Cd24a','Dlk1','Srr','Fabp4','Dpp4']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_adip_young,\n",
    "                       adipocyte_markers,\n",
    "                       save=\"_preadipYoung_markers_preadipocytes\",\n",
    "                       use_raw=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_adip_young,\n",
    "                       go_adip_dev,\n",
    "                       save=\"_preadipYoung_markers_GO_adipocyte_dev\", \n",
    "                       use_raw=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_tsne_marker(adata_adip_young, \n",
    "                     adipocyte_markers,\n",
    "                     size=20,\n",
    "                     save=\"_preadipYoung_markers_preadipocytes\", \n",
    "                     use_raw=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot==True:\n",
    "    plot_tsne_marker(adata_adip_young,\n",
    "                     go_adip_dev,\n",
    "                     size=20,\n",
    "                     save=\"_preadipYoung_markers_GO_adipocyte_dev\",\n",
    "                     use_raw=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.heatmap(adata=adata_adip_young, \n",
    "                  var_names=adipocyte_markers, \n",
    "                  groupby=\"louvain\", \n",
    "                  use_raw=True, \n",
    "                  log=True, \n",
    "                  cmap=\"viridis\",\n",
    "                  vmin=-1,\n",
    "                  vmax=4,\n",
    "                  dendrogram=False, \n",
    "                  var_group_rotation=90, \n",
    "                  show_gene_labels=True, \n",
    "                  show=True, \n",
    "                  save=\"_preadipYoung_markers_preadipocytes.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "if bool_plot == True:\n",
    "    sc.pl.heatmap(adata=adata_adip_young, \n",
    "                  var_names=go_adip_dev, \n",
    "                  groupby=\"louvain\", \n",
    "                  use_raw=True, \n",
    "                  log=True, \n",
    "                  cmap=\"viridis\",\n",
    "                  vmin=-1,\n",
    "                  vmax=2,\n",
    "                  dendrogram=False, \n",
    "                  var_group_rotation=90, \n",
    "                  show_gene_labels=True, \n",
    "                  show=True, \n",
    "                  save=\"_preadipYoung_markers_GO_adipocyte_dev.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "scrolled": true
   },
   "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"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[{'mean expression of Slc7a10 in louvain group 0': 0.07535989},\n",
       " {'mean expression of Slc7a10 in louvain group 1': 0.13056228},\n",
       " {'mean expression of Slc7a10 in louvain group 2': 0.18283738},\n",
       " {'mean expression of Slc7a10 in louvain group 3': 0.14173967},\n",
       " {'mean expression of Slc7a10 in louvain group 4': 0.00809834}]"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[{'mean expression of Slc7a10 in louvain group '+x:\n",
    "  np.mean(adata_adip_young[adata_adip_young.obs['louvain'].values == x,:][:,'Slc7a10'].X)} \n",
    " for x in adata_adip_young.obs['louvain'].cat.categories.values]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This seems to be a continuum of cell types with two opposing gradients of gene expression: Pdgfra and Slc7a10. This could indicate a developmental lineage."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Young mouse - coarse clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_adip_young_lowres,\n",
    "                       adipocyte_markers,\n",
    "                       save=\"_preadipYoung_lowres_markers_preadipocytes\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_adip_young_lowres, \n",
    "                       go_adip_dev, \n",
    "                       save=\"_preadipYoung_lowres_markers_GO_adipocyte_dev\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.heatmap(adata=adata_adip_young_lowres, \n",
    "                  var_names=adipocyte_markers, \n",
    "                  groupby=\"louvain\", \n",
    "                  use_raw=True, \n",
    "                  log=True, \n",
    "                  dendrogram=False, \n",
    "                  var_group_rotation=90, \n",
    "                  show_gene_labels=True, \n",
    "                  show=True, \n",
    "                  save=\"_preadipYoung_lowres_markers_preadipocytes.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.heatmap(adata=adata_adip_young_lowres, \n",
    "                  var_names=go_adip_dev, \n",
    "                  groupby=\"louvain\", \n",
    "                  use_raw=True, \n",
    "                  log=True, \n",
    "                  dendrogram=False, \n",
    "                  var_group_rotation=90, \n",
    "                  show_gene_labels=True, \n",
    "                  show=True, \n",
    "                  save=\"_preadipYoung_lowres_markers_GO_adipocyte_dev.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Old mouse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_adip_old, \n",
    "                       adipocyte_markers,\n",
    "                       save=\"_preadipOld_markers_preadipocytes\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_violin_marker(adata_adip_old, \n",
    "                       go_adip_dev, \n",
    "                       save=\"_preadipOld_markers_GO_adipocyte_dev\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_tsne_marker(adata_adip_old,\n",
    "                     adipocyte_markers, \n",
    "                     size=10,\n",
    "                     save=\"_preadipOld_markers_preadipocytes\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    plot_tsne_marker(adata_adip_old,\n",
    "                     go_adip_dev, \n",
    "                     size=10,\n",
    "                     save=\"_preadipOld_markers_GO_adipocyte_dev\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.heatmap(adata=adata_adip_old, \n",
    "                  var_names=adipocyte_markers, \n",
    "                  groupby=\"louvain\", \n",
    "                  use_raw=True, \n",
    "                  log=True, \n",
    "                  dendrogram=False, \n",
    "                  var_group_rotation=90, \n",
    "                  show_gene_labels=True, \n",
    "                  show=True, \n",
    "                  save=\"_preadipOld_markers_preadipocytes.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.heatmap(adata=adata_adip_old, \n",
    "                  var_names=go_adip_dev, \n",
    "                  groupby=\"louvain\", \n",
    "                  use_raw=True, \n",
    "                  log=True, \n",
    "                  dendrogram=False, \n",
    "                  var_group_rotation=90, \n",
    "                  show_gene_labels=True, \n",
    "                  show=True, \n",
    "                  save=\"_preadipOld_markers_GO_adipocyte_dev.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks like there are two clusters, one consisting of louvain group 0,4 (Pdgfra-Slc7a10-) and one of group 1,2,3,5 (Pdgfra+Slc7a10-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differential expression analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## By age"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "dets_age = de.test.t_test(data=adata_adip.raw, \n",
    "                          sample_description=adata_adip.obs, \n",
    "                          grouping=\"age\", \n",
    "                          is_logged=False\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The differentially expressed genes with some fold-change and mean expression thresholding are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    dets_age_summary = dets_age.summary(mean_thres=np.log(0.01), \n",
    "                                        qval_thres=0.05, \n",
    "                                        fc_lower_thres=0.5, \n",
    "                                        fc_upper_thres=2\n",
    "    )\n",
    "    dets_age_summary.to_csv(path_or_buf=dir_tables+\"DE_preadip_by_age.tab\", \n",
    "                            sep=\"\\t\"\n",
    "    )\n",
    "    dets_age_summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    dets_age.plot_volcano(alpha=0.01,\n",
    "                          min_fc=1.5,\n",
    "                          size=15, \n",
    "                          log10_p_threshold=-50,\n",
    "                          log2_fc_threshold=7,\n",
    "                          highlight_ids=['Slc7a10'], \n",
    "                          highlight_size=30, \n",
    "                          highlight_col=\"red\",\n",
    "                          save=sc_settings_figdir+\"preadip_DE_age\", \n",
    "                          suffix=\"_volcano.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot==True:\n",
    "    dets_age.plot_ma(alpha=0.01,\n",
    "                     log2_fc_threshold=6,\n",
    "                     highlight_ids=['Slc7a10'], \n",
    "                     highlight_size=30, \n",
    "                     highlight_col=\"red\",\n",
    "                     save=sc_settings_figdir+\"preadip_DE_age\",\n",
    "                     suffix=\"_ma_plot.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Count DE genes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3359"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(np.array([x < 0.01 if x is not np.nan else False for x in dets_age.qval]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## By louvain group"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coarse clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we perform a differential expression test across the coarse-grained louvain clustering which yielded two groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot==True:\n",
    "    dets_young_lowres_louvain = de.test.t_test(data=adata_adip_young_lowres.raw, \n",
    "                                               sample_description=adata_adip_young_lowres.obs, \n",
    "                                               grouping=\"louvain\", \n",
    "                                               is_logged=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The differentially expressed genes with some fold-change and mean expression thresholding are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    dets_young_lowres_louvain_summary = dets_young_lowres_louvain.summary(mean_thres=np.log(0.01), \n",
    "                                                                          qval_thres=0.05, \n",
    "                                                                          fc_lower_thres=0.5,\n",
    "                                                                          fc_upper_thres=2\n",
    "    )\n",
    "    dets_young_lowres_louvain_summary.to_csv(path_or_buf=dir_tables+\"DE_preadipYoung_by_louvain_lowres.csv\",\n",
    "                                             sep=\"\\t\"\n",
    "    )\n",
    "    dets_young_lowres_louvain_summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    dets_young_lowres_louvain.plot_volcano(alpha=0.01,\n",
    "                                           min_fc=1.5,\n",
    "                                           size=15,\n",
    "                                           log10_p_threshold=-15,\n",
    "                                           log2_fc_threshold=5,\n",
    "                                           highlight_ids=['Slc7a10'], \n",
    "                                           highlight_size=30,\n",
    "                                           highlight_col=\"red\",\n",
    "                                           save=sc_settings_figdir+\"preadipYoung_DE_louvain_lowres\",\n",
    "                                           suffix=\"_volcano.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    dets_young_lowres_louvain.plot_ma(alpha=0.01,\n",
    "                                      log2_fc_threshold=5,\n",
    "                                      highlight_ids=['Slc7a10'], \n",
    "                                      highlight_size=30,\n",
    "                                      highlight_col=\"red\",\n",
    "                                      save=sc_settings_figdir+\"preadipYoung_DE_louvain_lowres\",\n",
    "                                      suffix=\"_ma_plot.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Count DE genes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1927"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(np.array([x < 0.01 if x is not np.nan else False for x in dets_young_lowres_louvain.qval]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save DE genes to file for enrichment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "dets_young_lowres_louvain.summary(qval_thres=0.01)[\"gene\"].to_csv(dir_out+\"enrichment/de_genes_young_preadip.csv\",\n",
    "                                                                  index=False\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Heatmaps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select all differentially expressed genes at a corrected p-value threshold of 0.01 and a minimal or maximal log2 fold change of 2 or 0.5 and a minimal mean expression of 0.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dets_young_lowres_louvain_summary_forheatmap = dets_young_lowres_louvain.summary(mean_thres=0.5,\n",
    "                                                                                 qval_thres=0.01, \n",
    "                                                                                 fc_lower_thres=0.5, \n",
    "                                                                                 fc_upper_thres=2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_de_genes_young_coarse = dets_young_lowres_louvain_summary_forheatmap['gene'].values[\n",
    "    np.argsort(dets_young_lowres_louvain_summary_forheatmap['log2fc'].values)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "257\n"
     ]
    }
   ],
   "source": [
    "print(len(all_de_genes_young_coarse))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.pl.heatmap(adata=adata_adip_young_lowres, \n",
    "                  var_names=all_de_genes_young_coarse, \n",
    "                  groupby=\"louvain\", \n",
    "                  use_raw=False, \n",
    "                  log=False, \n",
    "                  cmap=\"viridis\",\n",
    "                  vmin=-1,\n",
    "                  vmax=5,\n",
    "                  dendrogram=False, \n",
    "                  var_group_rotation=90, \n",
    "                  show_gene_labels=True, \n",
    "                  show=True, \n",
    "                  save=\"_preadipYoung_lowres_all_de_genes.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.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": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "295.4px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}