astrocyte_scRNAseq_cluster_approach.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of the scRNA-seq data from astrocytes in  the ARC"
   ]
  },
  {
   "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=#Counts>Gene numbers and counts with and without mitochondrial RNA</a>\n",
    "* <a href=#Allcells>All cells - normalization, projection and clustering</a>\n",
    "* <a href=#Define>Define Cell Types</a>\n",
    "* <a href=#astrocytes>Astrocytes Only</a>\n",
    "    * <a href=#Embedding>Embeddings and Clustering</a>\n",
    "    * <a href=#adipmarkers>Astrocyte Marker Analysis</a>\n",
    "    * <a href=#topde>Top ranking DE Genes</a>\n",
    "    * <a href=#count_dist>Count distribution for Gfap, Aldh1l1 and Slc1a3</a>\n",
    "* <a href=#traject>Gfap and Aldh1l1 only</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": 1,
   "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",
      "appdirs             1.4.4\n",
      "appnope             0.1.2\n",
      "autoreload          NA\n",
      "backcall            0.2.0\n",
      "batchglm            v0.7.4\n",
      "bioservices         1.7.11\n",
      "bs4                 4.9.3\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",
      "colorlog            NA\n",
      "custom_functions    NA\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",
      "docutils            0.16\n",
      "easydev             0.11.0\n",
      "get_version         2.1\n",
      "gprofiler           1.0.0\n",
      "gseapy              0.10.4\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",
      "littleutils         NA\n",
      "llvmlite            0.36.0\n",
      "lxml                4.6.3\n",
      "matplotlib          3.3.4\n",
      "matplotlib_venn     0.11.6\n",
      "mpl_toolkits        NA\n",
      "natsort             7.1.1\n",
      "numba               0.53.1\n",
      "numexpr             2.7.3\n",
      "numpy               1.20.1\n",
      "outdated            0.2.1\n",
      "packaging           20.9\n",
      "pandas              1.2.3\n",
      "pandas_flavor       NA\n",
      "parso               0.7.0\n",
      "patsy               0.5.1\n",
      "pexpect             4.8.0\n",
      "pickleshare         0.7.5\n",
      "pingouin            0.3.12\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",
      "pylab               NA\n",
      "pyparsing           2.4.7\n",
      "pytz                2021.1\n",
      "requests            2.25.1\n",
      "requests_cache      0.5.2\n",
      "scanpy              1.7.1\n",
      "scipy               1.6.1\n",
      "scvelo              0.2.3\n",
      "seaborn             0.11.1\n",
      "sinfo               0.3.1\n",
      "six                 1.15.0\n",
      "sklearn             0.24.1\n",
      "soupsieve           2.2.1\n",
      "sparse              0.12.0\n",
      "sphinxcontrib       NA\n",
      "statsmodels         0.12.2\n",
      "storemagic          NA\n",
      "tables              3.6.1\n",
      "tabulate            0.8.9\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",
      "wrapt               1.12.1\n",
      "xarray              0.18.2\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-10-07 10:50\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 pandas as pd\n",
    "import seaborn as sb\n",
    "import scvelo as scv\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import colors\n",
    "from gprofiler import GProfiler\n",
    "import custom_functions as cf\n",
    "from matplotlib_venn import venn3_unweighted\n",
    "from scipy import stats\n",
    "import pingouin as pg\n",
    "import matplotlib_venn\n",
    "import statistics\n",
    "import gseapy\n",
    "import sys\n",
    "import re\n",
    "import os\n",
    "\n",
    "import batchglm\n",
    "import diffxpy.api as de\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "sc.settings.verbosity=3 # amount of output\n",
    "\n",
    "base_dir = '/Users/viktorian.miok/Documents/consultation/Luiza/single_cell/data/scanpy_AnnData/'\n",
    "dir_out = '/Users/viktorian.miok/Documents/consultation/Luiza/single_cell/results/'\n",
    "dir_scv = '/Users/viktorian.miok/Documents/consultation/Luiza/single_cell/data/velocyto/'\n",
    "dir_tables = dir_out+'tables/'\n",
    "sc_settings_figdir = dir_out+'figures/'\n",
    "sc_settings_writedir = dir_out+'anndata/'\n",
    "sc.logging.print_versions()\n",
    "os.chdir(dir_out)\n",
    "sc.settings.set_figure_params(dpi=80, scanpy=True)\n",
    "print (sys.version)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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",
    ")\n",
    "sc.set_figure_params(scanpy=True, \n",
    "                     fontsize=17\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": 3,
   "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": 4,
   "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": 5,
   "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": 6,
   "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": 7,
   "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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp:\n",
    "    adata_raw1 = sc.read(base_dir + 'MUC26030/filtered_feature_bc_matrix.h5ad')\n",
    "    adata_raw2 = sc.read(base_dir + 'MUC26031/filtered_feature_bc_matrix.h5ad')\n",
    "    adata_raw3 = sc.read(base_dir + 'MUC26032/filtered_feature_bc_matrix.h5ad')\n",
    "    adata_raw = adata_raw1.concatenate([adata_raw2, adata_raw3],\n",
    "                                       batch_key='diet', \n",
    "                                       batch_categories=['chow', 'hfd_5', 'hfd_15']\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": [
    "<a id=\"QC\"></a>"
   ]
  },
  {
   "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": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.filter_cells(adata_raw,\n",
    "                   min_counts=1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data contains 21143 observations with 31253 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": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of cells: 21143\n",
      "Number of genes: 31253\n",
      "Number of cells per diet:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "chow      7116\n",
       "hfd_5     6204\n",
       "hfd_15    7823\n",
       "Name: diet, dtype: int64"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('Number of cells: {:d}'.format(adata_raw.n_obs))\n",
    "print('Number of genes: {:d}'.format(adata_raw.shape[1]))\n",
    "print('Number of cells per diet:')\n",
    "adata_raw.obs['diet'].value_counts().sort_index()"
   ]
  },
  {
   "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": 11,
   "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",
    "adata_qc.obs['n_counts'] = adata_qc.X.sum(1)\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": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    t1 = sc.pl.violin(adata_qc, \n",
    "                      ['n_counts', 'n_genes', 'mt_frac'],\n",
    "                      size=1, \n",
    "                      log=False,\n",
    "                      jitter=3,\n",
    "                      multi_panel=True\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot==True:\n",
    "    sc.pl.highest_expr_genes(adata_qc,\n",
    "                             n_top=20\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": 14,
   "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']<5000], \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": "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": 15,
   "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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']<500],\n",
    "                      kde=False,\n",
    "                      bins=60\n",
    "    )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of cells: 21143\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "filtered out 35 cells that have more than 100000 counts\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of cells after min count filter: 21143\n",
      "Number of cells after max count filter: 21108\n",
      "Number of cells after MT filter: 20491\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "filtered out 496 cells that have less than 350 genes expressed\n",
      "Trying to set attribute `.obs` of view, copying.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of cells after gene filter: 19995\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,\n",
    "                   min_counts=200\n",
    ")\n",
    "print('Number of cells after min count filter: {:d}'.format(adata_qc.n_obs))\n",
    "\n",
    "sc.pp.filter_cells(adata_qc,\n",
    "                   max_counts=100000\n",
    ")\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.5]\n",
    "print('Number of cells after MT filter: {:d}'.format(adata_qc.n_obs))\n",
    "\n",
    "sc.pp.filter_cells(adata_qc,\n",
    "                   min_genes=350\n",
    ")\n",
    "print('Number of cells after gene filter: {:d}'.format(adata_qc.n_obs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of genes: 31253\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "filtered out 13164 genes that are detected in less than 20 cells\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of genes after cell filter: 18089\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": "code",
   "execution_count": 19,
   "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",
    "    p3 = sc.pl.scatter(adata_qc[adata_qc.obs['n_counts']<5000],\n",
    "                       'n_counts',\n",
    "                       'n_genes',\n",
    "                       color='mt_frac',\n",
    "                       size=5\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of cells: 19995\n",
      "Number of genes: 18089\n",
      "Number of cells per diet:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "chow      6741\n",
       "hfd_5     5886\n",
       "hfd_15    7368\n",
       "Name: diet, dtype: int64"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('Number of cells: {:d}'.format(adata_qc.n_obs))\n",
    "print('Number of genes: {:d}'.format(adata_qc.shape[1]))\n",
    "print('Number of cells per diet:')\n",
    "adata_qc.obs['diet'].value_counts().sort_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## All cells - normalization, projection and clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "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:02)\n"
     ]
    }
   ],
   "source": [
    "if bool_recomp == True:\n",
    "        \n",
    "    adata_proc = adata_qc.copy()\n",
    "    adata_proc.raw = adata_qc\n",
    "    sc.pp.normalize_per_cell(adata_proc)\n",
    "    sc.pp.log1p(adata_proc)\n",
    "    sc.pp.combat(adata_proc,\n",
    "                 key='diet'\n",
    "    )\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",
    "    #adata_proc.X=adata_proc.X.toarray()\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.umap(adata_proc)\n",
    "    if bool_recluster == True:\n",
    "        #sc.tl.louvain(adata_proc, resolution=0.5, flavor='vtraag', random_state=0)\n",
    "        sc.tl.leiden(adata_proc,\n",
    "                     resolution=0.3\n",
    "        )\n",
    "        pd.DataFrame(adata_proc.obs).to_csv(path_or_buf =sc_settings_writedir+\"obs_adata_proc.csv\")\n",
    "    else:\n",
    "        obs = pd.read_csv(sc_settings_writedir+'obs_adata_proc.csv')\n",
    "        adata_proc.obs['leiden'] = pd.Series(obs['leiden'].values,\n",
    "                                             dtype = 'category'\n",
    "        )\n",
    "    sc.write(sc_settings_writedir+'adata_proc.h5ad',\n",
    "             adata_proc\n",
    "    )\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": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    cf.plot_umap_marker(adata_proc,\n",
    "                        ['leiden'],\n",
    "                        save=\"_all_cells_leiden\",\n",
    "                        use_raw=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "#####################################################################################################################\n",
    "if bool_plot == True:\n",
    "    plt.rcParams['figure.figsize']=[5,5]\n",
    "    cf.plot_umap_marker(adata_proc, \n",
    "                        ['leiden'],\n",
    "                        save=\"_all_cells_leiden_ondata\",\n",
    "                        use_raw=False, \n",
    "                        legend_loc='on data',\n",
    "                        frameon=False,\n",
    "                        title='',\n",
    "                        size=20\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    cf.plot_umap_marker(adata_proc,\n",
    "                        ['Nnat'],\n",
    "                        color_map=mymap,\n",
    "                        use_raw=False,\n",
    "                        frameon=False,\n",
    "                        size=20\n",
    "    )\n",
    "    cf.plot_umap_marker(adata_proc, \n",
    "                        ['Peg3'],\n",
    "                        color_map=mymap,\n",
    "                        use_raw=False,\n",
    "                        frameon=False,\n",
    "                        size=20\n",
    "    )\n",
    "    cf.plot_umap_marker(adata_proc,\n",
    "                        ['Plagl1'],\n",
    "                        color_map=mymap,\n",
    "                        use_raw=False,\n",
    "                        frameon=False, \n",
    "                        size=20\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sd=adata_proc[adata_proc.obs['diet'] == 'chow']\n",
    "    h5=adata_proc[adata_proc.obs['diet'] == 'hfd_5']\n",
    "    h15=adata_proc[adata_proc.obs['diet'] == 'hfd_15']\n",
    "    cf.plot_umap_marker(sd, \n",
    "                        ['Ucp2', 'Cpt1a'],\n",
    "                        color_map=mymap,\n",
    "                        use_raw=False,\n",
    "                        frameon=False, \n",
    "                        title='',\n",
    "                        size=20\n",
    "    )\n",
    "    cf.plot_umap_marker(h5, \n",
    "                        ['Ucp2', 'Cpt1a'],\n",
    "                        color_map=mymap,\n",
    "                        use_raw=False,\n",
    "                        frameon=False, \n",
    "                        size=20\n",
    "    )\n",
    "    cf.plot_umap_marker(h15,\n",
    "                        ['Ucp2', 'Cpt1a'],\n",
    "                        color_map=mymap,\n",
    "                        use_raw=False, \n",
    "                        frameon=False, \n",
    "                        size=20\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of cells in each cluster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0     3921\n",
       "1     2547\n",
       "2     2386\n",
       "3     2247\n",
       "4     1466\n",
       "5     1403\n",
       "6     1233\n",
       "7     1080\n",
       "8     1075\n",
       "9     1056\n",
       "10     877\n",
       "11     515\n",
       "12     109\n",
       "13      80\n",
       "Name: leiden, dtype: int64"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adata_proc.obs[\"leiden\"].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define Cell Types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"DE\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary heatmap, dotplot and stacked_violin for cluster assignments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "#####################################################################################################################\n",
    "if bool_plot==True:\n",
    "    marker_genes_dict = {'Astrocytes': ['Slc1a2', 'Slc1a3', 'Aqp4', 'Gfap', 'Aldh1l1', 'Gja1', 'Gjb6', 'Atp1b2'],\n",
    "                         'Endothelial cells': ['Cldn5', 'Pecam1', 'Slco1c1'],\n",
    "                         'Ependymal cells': ['Ccdc153', 'Rarres2', 'Hdc', 'Tm4sf1'],\n",
    "                         'Microglia': ['Itgam', 'Tmem119', 'Cx3cr1','Csf1r', 'Aif1', 'P2ry12'],\n",
    "                         'Mural cells': ['Mustn1', 'Pdgfrb', 'Des'],\n",
    "                         'Neurons': ['Rbfox3', 'Syp', 'Tubb3', 'Snap25', 'Syt1'],\n",
    "                         'Oligodendrocytes': ['Olig1','Mog', 'Mag'],\n",
    "                         'Mural cells': ['Mustn1', 'Pdgfrb','Des'],\n",
    "                         'Tanycytes': ['Rax', 'Lhx2', 'Col23a1', 'Slc16a2', 'Crym', 'Adm'],\n",
    "                         'VLMCs': ['Lum', 'Col1a1', 'Col3a1']}\n",
    "    \n",
    "    for i in list(marker_genes_dict.keys()):\n",
    "        marker_genes_dict[i].sort()\n",
    "        \n",
    "    sc.pl.dotplot(adata=adata_proc,\n",
    "                  var_names=marker_genes_dict, \n",
    "                  groupby='leiden',\n",
    "                  use_raw=False, \n",
    "                  log=False, \n",
    "                  dendrogram=True, \n",
    "                  var_group_rotation=90, \n",
    "                  show=True, \n",
    "                  #size_title=5,\n",
    "                  save=\"all_cells_celltypes_markers.pdf\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot==True:\n",
    "    adata_proc.obs['Astrocytes'] = (adata_proc[:,'Gfap'].X > 0.01) & (adata_proc[:,'Aldh1l1'].X > 0.01) & \n",
    "                                    (adata_proc[:,'Slc1a3'].X > 0.01)*1\n",
    "    adata_proc.obs['Endothelial'] = (adata_proc[:,'Cldn5'].X > 0.01) & (adata_proc[:,'Pecam1'].X > 0.01) & \n",
    "                                    (adata_proc[:,'Slco1c1'].X > 0.01)*1\n",
    "    adata_proc.obs['Ependymal'] = (adata_proc[:,'Ccdc153'].X > 0.01) & (adata_proc[:,'Rarres2'].X > 0.01) & \n",
    "                                  (adata_proc[:,'Hdc'].X > 0.01) & (adata_proc[:,'Tm4sf1'].X > 0.01)*1\n",
    "    adata_proc.obs['Microglia'] = (adata_proc[:,'Itgam'].X > 0.01) & (adata_proc[:,'Tmem119'].X > 0.01) & \n",
    "                                  (adata_proc[:,'Cx3cr1'].X > 0.01) & (adata_proc[:,'Csf1r'].X > 0.01) & \n",
    "                                  (adata_proc[:,'Aif1'].X > 0.01) & (adata_proc[:,'P2ry12'].X > 0.01)*1\n",
    "    adata_proc.obs['Mural'] = (adata_proc[:,'Mustn1'].X > 0.01) & (adata_proc[:,'Pdgfrb'].X > 0.01) & \n",
    "                               (adata_proc[:,'Des'].X > 0.01)*1\n",
    "    adata_proc.obs['Neurons'] = (adata_proc[:,'Rbfox3'].X > 0.01) & (adata_proc[:,'Syp'].X > 0.01) & \n",
    "                                (adata_proc[:,'Tubb3'].X > 0.01) & (adata_proc[:,'Snap25'].X > 0.01) & \n",
    "                                (adata_proc[:,'Syt1'].X > 0.01)*1\n",
    "    adata_proc.obs['Oligodendro'] = (adata_proc[:,'Olig1'].X > 0.01) & (adata_proc[:,'Mog'].X > 0.01) &\n",
    "                                    (adata_proc[:,'Mag'].X > 0.01)*1\n",
    "    adata_proc.obs['Tanycytes'] = (adata_proc[:,'Rax'].X > 0.01) & (adata_proc[:,'Lhx2'].X > 0.01) & \n",
    "                                  (adata_proc[:,'Col23a1'].X > 0.01) & (adata_proc[:,'Slc16a2'].X > 0.01)*1\n",
    "    adata_proc.obs['VLMCs'] = (adata_proc[:,'Lum'].X > 0.01) & (adata_proc[:,'Col1a1'].X > 0.01) & \n",
    "                              (adata_proc[:,'Col3a1'].X > 0.01)*1\n",
    "\n",
    "\n",
    "    adata_proc.obs['gfap_lepr'] = (adata_proc[:,'Gfap'].X > 0.01) & (adata_proc[:,'Lepr'].X > 0.01) &\n",
    "                                  (adata_proc[:,'Aldh1l1'].X == 0)*1\n",
    "    adata_proc.obs['gfap_insr'] = (adata_proc[:,'Gfap'].X > 0.01) & (adata_proc[:,'Insr'].X > 0.01) & \n",
    "    (adata_proc[:,'Aldh1l1'].X ==0)*1\n",
    "    adata_proc.obs['gfap_igf1r'] = (adata_proc[:,'Gfap'].X > 0.01) & (adata_proc[:,'Igf1r'].X > 0.01) &\n",
    "    (adata_proc[:,'Aldh1l1'].X == 0)*1\n",
    "    adata_proc.obs['gfap_glp1r'] = (adata_proc[:,'Gfap'].X > 0.01) & (adata_proc[:,'Glp1r'].X > 0.01) & \n",
    "    (adata_proc[:,'Aldh1l1'].X == 0)*1\n",
    "    adata_proc.obs['gfap_lpl'] = (adata_proc[:,'Gfap'].X > 0.01) & (adata_proc[:,'Lpl'].X > 0.01) & \n",
    "    (adata_proc[:,'Aldh1l1'].X == 0)*1\n",
    "    adata_proc.obs['gfap_slc2a1'] = (adata_proc[:,'Gfap'].X > 0.01) & (adata_proc[:,'Slc2a1'].X > 0.01) & \n",
    "    (adata_proc[:,'Aldh1l1'].X == 0)*1\n",
    "\n",
    "    adata_proc.obs['aldh_lepr'] = (adata_proc[:,'Aldh1l1'].X > 0.01) & (adata_proc[:,'Lepr'].X > 0.01) &\n",
    "                                  (adata_proc[:,'Gfap'].X == 0)*1\n",
    "    adata_proc.obs['aldh_insr'] = (adata_proc[:,'Aldh1l1'].X > 0.01) & (adata_proc[:,'Insr'].X > 0.01) & \n",
    "                                  (adata_proc[:,'Gfap'].X == 0)*1\n",
    "    adata_proc.obs['aldh_igf1r'] = (adata_proc[:,'Aldh1l1'].X > 0.01) & (adata_proc[:,'Igf1r'].X > 0.01) & \n",
    "                                   (adata_proc[:,'Gfap'].X == 0)*1\n",
    "    adata_proc.obs['aldh_glp1r'] = (adata_proc[:,'Aldh1l1'].X > 0.01) & (adata_proc[:,'Glp1r'].X > 0.01) &\n",
    "                                   (adata_proc[:,'Gfap'].X == 0)*1\n",
    "    adata_proc.obs['aldh_lpl'] = (adata_proc[:,'Aldh1l1'].X > 0.01) & (adata_proc[:,'Lpl'].X > 0.01) & \n",
    "                                 (adata_proc[:,'Gfap'].X == 0)*1\n",
    "    adata_proc.obs['aldh_slc2a1'] = (adata_proc[:,'Aldh1l1'].X > 0.01) & (adata_proc[:,'Slc2a1'].X > 0.01) &\n",
    "                                    (adata_proc[:,'Gfap'].X == 0)*1\n",
    "\n",
    "    adata_proc.obs['both_lepr'] = (adata_proc[:,'Aldh1l1'].X > 0.01) & (adata_proc[:,'Lepr'].X > 0.01) & \n",
    "                                  (adata_proc[:,'Gfap'].X > 0.01)*1\n",
    "    adata_proc.obs['both_insr'] = (adata_proc[:,'Aldh1l1'].X > 0.01) & (adata_proc[:,'Insr'].X > 0.01) & \n",
    "                                  (adata_proc[:,'Gfap'].X > 0.01)*1\n",
    "    adata_proc.obs['both_igf1r'] = (adata_proc[:,'Aldh1l1'].X > 0.01) & (adata_proc[:,'Igf1r'].X > 0.01) & \n",
    "                                   (adata_proc[:,'Gfap'].X > 0.01)*1\n",
    "    adata_proc.obs['both_glp1r'] = (adata_proc[:,'Aldh1l1'].X > 0.01) & (adata_proc[:,'Glp1r'].X > 0.01) & \n",
    "                                   (adata_proc[:,'Gfap'].X > 0.01)*1\n",
    "    adata_proc.obs['both_lpl'] = (adata_proc[:,'Aldh1l1'].X > 0.01) & (adata_proc[:,'Lpl'].X > 0.01) &\n",
    "                                 (adata_proc[:,'Gfap'].X > 0.01)*1\n",
    "    adata_proc.obs['both_slc2a1'] = (adata_proc[:,'Aldh1l1'].X > 0.01) & (adata_proc[:,'Slc2a1'].X > 0.01) &\n",
    "                                    (adata_proc[:,'Gfap'].X > 0.01)*1\n",
    "    \n",
    "    adata_proc.obs.to_csv(\"adata_proc_upset.obs.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## UMAP with assigned cell types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "#####################################################################################################################\n",
    "if bool_plot == True:\n",
    "    new_cluster_names = {'0': \"Astrocytes\",\n",
    "                         '1': \"Tanycytes\",\n",
    "                         '2': \"Ependymal cells\",\n",
    "                         '3': \"Neurons\",\n",
    "                         '4': \"Endothelial cells\",\n",
    "                         '5': \"Oligodendrocytes\",\n",
    "                         '6': \"Microglia\",\n",
    "                         '7': \"Oligodendrocytes\",\n",
    "                         '8': \"Neurons\",\n",
    "                         '9': \"Microglia\",\n",
    "                         '10': \"Oligodendrocytes\",\n",
    "                         '11': \"Mural cells\",\n",
    "                         '12': \"VLMCs\",\n",
    "                         '13': \"Oligodendrocytes\"\n",
    "    }\n",
    "    adata_proc.obs['celltypes'] = [new_cluster_names[x] for x in  adata_proc.obs['leiden']]\n",
    "    \n",
    "    cf.plot_umap_marker(adata_proc,\n",
    "                        ['celltypes'],\n",
    "                        save=\"_all_cells_celltypes.png\",\n",
    "                        use_raw=False, \n",
    "                        frameon=False,\n",
    "                        size=20,\n",
    "                        title='', \n",
    "                        palette=['#1f77b4','#aa40fc','#279e68','#aec7e8','#98df8a',\n",
    "                                 '#d62728','#ffbb78','#ff7f0e','#ff9896'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Melanie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    adata_mel = adata_proc[(adata_proc.obs['celltypes']=='Neurons'),:] # Astrocytes Neurons\n",
    "    adata_mel.obs['diet_celltype'] = adata_mel.obs['celltypes'].astype(str) + '_' + adata_mel.obs['diet'].astype(str)\n",
    "\n",
    "    melanie=['Opa1','Dnm1l','Mfn1','Ngb','Apoe','Cpt1a','Lpl','Pdk4','Acadm','Acsl6','Acacb','Mlycd','Gja1']\n",
    "\n",
    "    adata_mel.layers['scaled'] = sc.pp.scale(adata_mel, copy=True).X\n",
    "\n",
    "    sc.pl.matrixplot(adata_mel,\n",
    "                     melanie,\n",
    "                     'diet_celltype',\n",
    "                     dendrogram=False,\n",
    "                     swap_axes=True,\n",
    "                     colorbar_title='mean z-score', \n",
    "                     layer='scaled',\n",
    "                     cmap='RdBu_r' #, vmin=-0.2, vmax=0.2\n",
    "    )\n",
    "    sc.pl.stacked_violin(adata_mel, melanie,\n",
    "                         'diet_celltype',\n",
    "                         dendrogram=False,\n",
    "                         swap_axes=True,\n",
    "                         colorbar_title='mean z-score', \n",
    "                         layer='scaled',\n",
    "                         cmap='RdBu_r'\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#####################################################################################################################\n",
    "if bool_plot == True:\n",
    "    adata_proc_hv = adata_proc[:, adata_proc.var.highly_variable]\n",
    "    adata_proc_hv.raw = adata_qc[:, adata_proc.var.highly_variable]\n",
    "    adata_proc_hv.obs['diet_leiden'] = adata_proc_hv.obs['diet'].str.cat(adata_proc_hv.obs['leiden'], sep='_')\n",
    "    adata_proc_hv.obs['diet_celltypes'] = adata_proc_hv.obs['diet'].str.cat(adata_proc_hv.obs['celltypes'], sep='_')\n",
    "    \n",
    "    sc.tl.rank_genes_groups(adata_proc_hv,\n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_5_Astrocytes'], \n",
    "                            reference='chow_Astrocytes',\n",
    "                            key_added=\"ct5_ast\",\n",
    "                            method='t-test'\n",
    "    ) #  wilcoxon\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_5_Endothelial cells'], \n",
    "                            reference='chow_Endothelial cells',\n",
    "                            key_added=\"ct5_end\", \n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_5_Ependymal cells'],\n",
    "                            reference='chow_Ependymal cells',\n",
    "                            key_added=\"ct5_epe\",\n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_5_Microglia'], \n",
    "                            reference='chow_Microglia',\n",
    "                            key_added=\"ct5_mic\",\n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_5_Mural cells'],\n",
    "                            reference='chow_Mural cells',\n",
    "                            key_added=\"ct5_mur\",\n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_5_Neurons'],\n",
    "                            reference='chow_Neurons',\n",
    "                            key_added=\"ct5_neu\",\n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_5_Oligodendrocytes'],\n",
    "                            reference='chow_Oligodendrocytes',\n",
    "                            key_added=\"ct5_oli\", \n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_5_Tanycytes'],\n",
    "                            reference='chow_Tanycytes',\n",
    "                            key_added=\"ct5_tan\",\n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_5_VLMCs'],\n",
    "                            reference='chow_VLMCs', \n",
    "                            key_added=\"ct5_vlmc\", \n",
    "                            method='t-test'\n",
    "    )\n",
    "    ##############\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv,\n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_15_Astrocytes'],\n",
    "                            reference='chow_Astrocytes',\n",
    "                            key_added=\"ct15_ast\",\n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_15_Endothelial cells'],\n",
    "                            reference='chow_Endothelial cells',\n",
    "                            key_added=\"ct15_end\",\n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_15_Ependymal cells'], \n",
    "                            reference='chow_Ependymal cells', \n",
    "                            key_added=\"ct15_epe\",\n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv,\n",
    "                            'diet_celltypes', \n",
    "                            groups=['hfd_15_Microglia'],\n",
    "                            reference='chow_Microglia',\n",
    "                            key_added=\"ct15_mic\", \n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes', \n",
    "                            groups=['hfd_15_Mural cells'],\n",
    "                            reference='chow_Mural cells', \n",
    "                            key_added=\"ct15_mur\", \n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv,\n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_15_Neurons'],\n",
    "                            reference='chow_Neurons', \n",
    "                            key_added=\"ct15_neu\", \n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv,\n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_15_Oligodendrocytes'],\n",
    "                            reference='chow_Oligodendrocytes',\n",
    "                            key_added=\"ct15_oli\", \n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_15_Tanycytes'],\n",
    "                            reference='chow_Tanycytes', \n",
    "                            key_added=\"ct15_tan\",\n",
    "                            method='t-test'\n",
    "    )\n",
    "    sc.tl.rank_genes_groups(adata_proc_hv, \n",
    "                            'diet_celltypes',\n",
    "                            groups=['hfd_15_VLMCs'],\n",
    "                            reference='chow_VLMCs', \n",
    "                            key_added=\"ct15_vlmc\",\n",
    "                            method='t-test'\n",
    "    )\n",
    "\n",
    "    ct0_5=[]\n",
    "    ct0_5up=[]\n",
    "    ct0_5down=[]\n",
    "    df0_5=[]\n",
    "    for i in ['ct5_ast','ct5_end','ct5_epe','ct5_mic','ct5_mur','ct5_neu','ct5_oli','ct5_tan','ct5_vlmc']:\n",
    "        result = adata_proc_hv.uns[i]\n",
    "        groups = result['names'].dtype.names\n",
    "        df5 = pd.DataFrame({group + '_' + key[:1]: result[key][group]\n",
    "                            for group in groups for key in ['pvals_adj', 'logfoldchanges']}) \n",
    "        #print(sum(df5.iloc[:, 0] < 0.05))\n",
    "        print(sum((df5.iloc[:, 0] < 0.05) & (df5.iloc[:, 1] > 0)))\n",
    "        print(-sum((df5.iloc[:, 0] < 0.05) & (df5.iloc[:, 1] <= 0)))\n",
    "        ct0_5.append(sum(df5.iloc[:, 0] < 0.05))  \n",
    "        ct0_5up.append(sum((df5.iloc[:, 0] < 0.05) & (df5.iloc[:, 1] > 0)))  \n",
    "        ct0_5down.append(sum((df5.iloc[:, 0] < 0.05) & (df5.iloc[:, 1] <= 0)))  \n",
    "        df0_5.append(df5) \n",
    "        \n",
    "    ct0_15=[]\n",
    "    ct0_15up=[]\n",
    "    ct0_15down=[]\n",
    "    df0_15=[]\n",
    "    for i in ['ct15_ast','ct15_end','ct15_epe','ct15_mic','ct15_mur',\n",
    "              'ct15_neu','ct15_oli','ct15_tan','ct15_vlmc']:\n",
    "        result = adata_proc_hv.uns[i]\n",
    "        groups = result['names'].dtype.names\n",
    "        df15 = pd.DataFrame({group + '_' + key[:1]: result[key][group]\n",
    "                             for group in groups for key in ['pvals_adj', 'logfoldchanges']}) \n",
    "        #print(sum(df15.iloc[:, 0] < 0.05))\n",
    "        #print(sum((df15.iloc[:, 0] < 0.05) & (df15.iloc[:, 1] > 0)))\n",
    "        #print(-sum((df15.iloc[:, 0] < 0.05) & (df15.iloc[:, 1] <= 0)))\n",
    "        ct0_15.append(sum(df15.iloc[:, 0] < 0.05))\n",
    "        df0_15.append(df15)\n",
    "    \n",
    "    d = {'0_5':ct0_5,'0_15':ct0_15}\n",
    "    dfin = pd.DataFrame(d)\n",
    "    dfin['0_0'] = [0,0,0,0,0,0,0,0,0]\n",
    "    dfin=dfin.reindex(['0_0','0_5','0_15'],\n",
    "                      axis=1\n",
    "    )\n",
    "    dfin.index=['Astrocytes','Endothelial cells','Ependymal cells','Microglia',\n",
    "                'Mural cells','Neurons','Oligodendrocytes','Tanycytes','VLMCs']\n",
    "    \n",
    "    plt.rcParams[\"figure.figsize\"] = (7,7)\n",
    "\n",
    "    line_chart1 = plt.plot(['SD', '5d HFHS diet', '15d HFHS diet'],\n",
    "                           dfin.iloc[0],\n",
    "                           linewidth=3,\n",
    "                           c='#1f77b4'\n",
    "    )\n",
    "    line_chart1 = plt.plot(['SD', '5d HFHS diet', '15d HFHS diet'], \n",
    "                           dfin.iloc[1],\n",
    "                           linewidth=3,\n",
    "                           c='#aa40fc'\n",
    "    )\n",
    "    line_chart1 = plt.plot(['SD', '5d HFHS diet', '15d HFHS diet'],\n",
    "                           dfin.iloc[2],\n",
    "                           linewidth=3, \n",
    "                           c='#279e68'\n",
    "    )\n",
    "    line_chart1 = plt.plot(['SD', '5d HFHS diet', '15d HFHS diet'], \n",
    "                           dfin.iloc[3], \n",
    "                           linewidth=3,\n",
    "                           c='#aec7e8'\n",
    "    )\n",
    "    line_chart1 = plt.plot(['SD', '5d HFHS diet', '15d HFHS diet'], \n",
    "                           dfin.iloc[4], \n",
    "                           linewidth=3, \n",
    "                           c='#98df8a'\n",
    "    )\n",
    "    line_chart1 = plt.plot(['SD', '5d HFHS diet', '15d HFHS diet'], \n",
    "                           dfin.iloc[5], \n",
    "                           linewidth=3,\n",
    "                           c='#d62728'\n",
    "    )\n",
    "    line_chart1 = plt.plot(['SD', '5d HFHS diet', '15d HFHS diet'], \n",
    "                           dfin.iloc[6],\n",
    "                           linewidth=3,\n",
    "                           c='#ffbb78'\n",
    "    )\n",
    "    line_chart1 = plt.plot(['SD', '5d HFHS diet', '15d HFHS diet'],\n",
    "                           dfin.iloc[7], \n",
    "                           linewidth=3,\n",
    "                           c='#ff7f0e'\n",
    "    )\n",
    "    line_chart1 = plt.plot(['SD', '5d HFHS diet', '15d HFHS diet'], \n",
    "                           dfin.iloc[8], \n",
    "                           linewidth=3,\n",
    "                           c='#ff9896'\n",
    "    )\n",
    "\n",
    "    #plt.xlabel(\"X axis label\")\n",
    "    plt.ylabel(\"Number of DEGs\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    result = adata_proc_hv.uns['ct5_ast']\n",
    "    groups = result['names'].dtype.names\n",
    "    dfa5 = pd.DataFrame({group + '_' + key[:1]: result[key][group]\n",
    "                        for group in groups for key in ['names', 'pvals_adj', 'logfoldchanges']}) \n",
    "\n",
    "    result = adata_proc_hv.uns['ct5_neu']\n",
    "    groups = result['names'].dtype.names\n",
    "    dfn5 = pd.DataFrame({group + '_' + key[:1]: result[key][group]\n",
    "                        for group in groups for key in ['names', 'pvals_adj', 'logfoldchanges']}) \n",
    "\n",
    "    result = adata_proc_hv.uns['ct5_mic']\n",
    "    groups = result['names'].dtype.names\n",
    "    dfm5 = pd.DataFrame({group + '_' + key[:1]: result[key][group]\n",
    "                        for group in groups for key in ['names', 'pvals_adj', 'logfoldchanges']}) \n",
    "    \n",
    "    result = adata_proc_hv.uns['ct15_ast']\n",
    "    groups = result['names'].dtype.names\n",
    "    dfa15 = pd.DataFrame({group + '_' + key[:1]: result[key][group]\n",
    "                        for group in groups for key in ['names', 'pvals_adj', 'logfoldchanges']}) \n",
    "\n",
    "    result = adata_proc_hv.uns['ct15_neu']\n",
    "    groups = result['names'].dtype.names\n",
    "    dfn15 = pd.DataFrame({group + '_' + key[:1]: result[key][group]\n",
    "                        for group in groups for key in ['names', 'pvals_adj', 'logfoldchanges']}) \n",
    "    \n",
    "    result = adata_proc_hv.uns['ct15_mic']\n",
    "    groups = result['names'].dtype.names\n",
    "    dfm15 = pd.DataFrame({group + '_' + key[:1]: result[key][group]\n",
    "                        for group in groups for key in ['names', 'pvals_adj', 'logfoldchanges']}) \n",
    "    \n",
    "    #dfa5.to_csv(\"DGE0_5astro.csv\")\n",
    "    #dfn5.to_csv(\"DGE0_5neuro.csv\")\n",
    "    #dfa15.to_csv(\"DGE0_15astro.csv\")\n",
    "    #dfn15.to_csv(\"DGE0_15neuro.csv\")\n",
    "    \n",
    "    #df0_5.to_csv(\"DGE0_5celltypes.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True: \n",
    "    dfa5[dfa5['hfd_5_Astrocytes_n'].isin(melanie)][['hfd_5_Astrocytes_n', 'hfd_5_Astrocytes_p']]\n",
    "    dfa15[dfa15['hfd_15_Astrocytes_n'].isin(melanie)][['hfd_15_Astrocytes_n', 'hfd_15_Astrocytes_p']]\n",
    "    dfn5[dfn5['hfd_5_Neurons_n'].isin(melanie)][['hfd_5_Neurons_n', 'hfd_5_Neurons_p']]\n",
    "    dfn15[dfn15['hfd_15_Neurons_n'].isin(melanie)][['hfd_15_Neurons_n', 'hfd_15_Neurons_p']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:  \n",
    "    dfa5[dfa5['hfd_5_Astrocytes_n'].isin(np.concatenate((marker_pan_a1_a2['PAN_reactive'],\n",
    "                                                         marker_pan_a1_a2['A1_specific'],\n",
    "                                                         marker_pan_a1_a2['A2_specific']), \n",
    "                                                         axis=None))][['hfd_5_Astrocytes_n', 'hfd_5_Astrocytes_p']]\n",
    "    dfa15[dfa15['hfd_15_Astrocytes_n'].isin(np.concatenate((marker_pan_a1_a2['PAN_reactive'],\n",
    "                                                            marker_pan_a1_a2['A1_specific'],\n",
    "                                                            marker_pan_a1_a2['A2_specific']),\n",
    "                                                            axis=None))][['hfd_15_Astrocytes_n', 'hfd_15_Astrocytes_p']]\n",
    "\n",
    "    dfm5[dfm5['hfd_5_Microglia_n'].isin(marker_genes_mglia['activated microglia'])][['hfd_5_Microglia_n', 'hfd_5_Microglia_p']]\n",
    "    dfm15[dfm15['hfd_15_Microglia_n'].isin(marker_genes_mglia['activated microglia'])][['hfd_15_Microglia_n', 'hfd_15_Microglia_p']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### microglia inflamation markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "#####################################################################################################################\n",
    "if bool_plot == True:  #  (adata_proc.obs['celltypes']=='Astrocytes') |\n",
    "    adata_inf=adata_proc[(adata_proc.obs['celltypes'] == 'Microglia'),:]\n",
    "    adata_inf.obs['diet_celltype'] = adata_inf.obs['celltypes'].astype(str) + '_' + adata_inf.obs['diet'].astype(str)\n",
    "\n",
    "    marker_genes_mglia = {#'inflam_mark': ['Il1b','Il2','Il6','Tnf','Il1r1','Il6ra','Tlr4','Lcn2','Mmp2','Ccl5','Ccl12',\n",
    "                          #                'Vegfa','Tpo','Plaur','Axin2'], # ,'Il5','Mmp13',\n",
    "                          'activated microglia': ['Itgam','Cd14','Fcgr3','Fcgr2b','Cd40','Ptprc','Cd68','Cd80',\n",
    "                                                  'Cd86','Cx3cr1','Adgre1','Fcer1g']} \n",
    "                     # ,'H2''C3H2-T23','Serping1','H2-D1','Ggta1','Iigp1','Gbp2','Fbln5','Fkbp5','Psmb8','Srgn','Amigo2'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "#####################################################################################################################\n",
    "if bool_plot == True:\n",
    "    adata_inf.layers['scaled'] = sc.pp.scale(adata_inf, copy=True).X\n",
    "\n",
    "    #sc.pl.matrixplot(adata_inf, marker_genes_mglia, 'diet_celltype', dendrogram=False, swap_axes=True,\n",
    "    #                 colorbar_title='mean z-score', layer='scaled', vmin=-0.2, vmax=0.2, cmap='RdBu_r')\n",
    "    sc.pl.stacked_violin(adata_inf, \n",
    "                         marker_genes_mglia,\n",
    "                         'diet_celltype', \n",
    "                         dendrogram=False,\n",
    "                         swap_axes=True,\n",
    "                         colorbar_title='mean z-score',\n",
    "                         layer='scaled',\n",
    "                         vmin=-2,\n",
    "                         vmax=2,\n",
    "                         cmap='RdBu_r'\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Astorcytes only"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"Embedding\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Embedding and Clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "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_astro = np.asarray(adata_proc.obs_names)[\n",
    "        [x in ['astrocytes'] \n",
    "         for x in np.asarray(adata_proc.obs['celltypes'].values)]\n",
    "    ]\n",
    "    adata_astro=adata_raw[cell_ids_astro,:].copy()  # adata_raw\n",
    "    #dat = pd.DataFrame(adata_proc.X, index=adata_proc.obs.index, columns=adata_proc.var.index)\n",
    "    adata_astro.obs['n_genes'] = (adata_astro.X > 0).sum(1)\n",
    "    adata_astro.obs['n_counts'] = adata_astro.X.sum(1)\n",
    "    mt_gene_mask = [gene.startswith('mt-') for gene in adata_astro.var_names]\n",
    "    temp_mt_sum = adata_astro[:,mt_gene_mask].X.sum(1)\n",
    "    temp_mt_sum = np.squeeze(np.asarray(temp_mt_sum))\n",
    "    temp_n_counts = adata_astro.obs['n_counts']\n",
    "    adata_astro.obs['mt_frac'] = temp_mt_sum/adata_astro.obs['n_counts']\n",
    "    adata_astro.raw = adata_astro\n",
    "    sc.pp.normalize_per_cell(adata_astro)\n",
    "    sc.pp.log1p(adata_astro)\n",
    "    sc.pp.highly_variable_genes(adata_astro,n_top_genes=4000)\n",
    "    sc.pl.highly_variable_genes(adata_astro)\n",
    "    adata_astro.X = adata_astro.X.toarray()\n",
    "    \n",
    "    sc.pp.pca(adata_astro,\n",
    "              n_comps=50, \n",
    "              use_highly_variable=True,\n",
    "              random_state=0,\n",
    "              svd_solver='arpack'\n",
    "    )\n",
    "    sc.pp.neighbors(adata_astro,\n",
    "                    n_neighbors=100,\n",
    "                    knn=True,\n",
    "                    method='umap',\n",
    "                    n_pcs=50,\n",
    "                    random_state=0\n",
    "    )\n",
    "    sc.tl.umap(adata_astro)\n",
    "    if bool_recluster == True:\n",
    "        sc.tl.leiden(adata_astro, \n",
    "                     resolution=0.5\n",
    "        )\n",
    "        pd.DataFrame(adata_astro.obs).to_csv(path_or_buf=sc_settings_writedir+'obs_adata_astro.csv')\n",
    "    else:\n",
    "        obs=pd.read_csv(sc_settings_writedir+'obs_adata_astro.csv')\n",
    "        adata_astro.obs['leiden']=pd.Series(obs['leiden'].values,\n",
    "                                            dtype = 'category'\n",
    "        )\n",
    "    sc.write(sc_settings_writedir+'adata_astro.h5ad',adata_astro)\n",
    "else:\n",
    "    adata_astro=sc.read(sc_settings_writedir+'adata_astro.h5ad') \n",
    "sc.tl.paga(adata_astro)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#####################################################################################################################\n",
    "if bool_plot == True:\n",
    "    adata_a=adata_astro\n",
    "    adata_a.layers['scaled'] = sc.pp.scale(adata_a, copy=True).X\n",
    "    marker_pan_a1_a2 = {'PAN_reactive': ['Lcn2','Steap4','S1pr3','Timp1','Hsbp1','Cxcl10','Cd44','Osmr','Cp','Serpina3n',\n",
    "                                         'Aspg','Vim','Gfap'],\n",
    "                        'A1_specific': ['Serping1','H2-D1','Ggta1','Iigp1','Gbp2','Fbln5','Fkbp5','Psmb8',\n",
    "                                        'Srgn','Amigo2'], #'C3H2-T23',\n",
    "                        'A2_specific': ['Clcf1','Tgm1','Ptx3','S100a10','Sphk1','Cd109','Ptgs2','Emp1','Slc10a6','Tm4sf1',\n",
    "                                        'B3gnt5','Cd14','Stat3']}\n",
    "    sc.pl.matrixplot(adata_a,\n",
    "                     marker_pan_a1_a2,\n",
    "                     'diet', \n",
    "                     dendrogram=False, \n",
    "                     swap_axes=True, \n",
    "                     colorbar_title='mean z-score',\n",
    "                     layer='scaled',\n",
    "                     vmin=-0.2,\n",
    "                     vmax=0.2, \n",
    "                     cmap='RdBu_r'\n",
    "    )\n",
    "    sc.pl.heatmap(adata_a,\n",
    "                  marker_pan_a1_a2,\n",
    "                  swap_axes=True, \n",
    "                  show=False, \n",
    "                  use_raw=True,\n",
    "                  show_gene_labels=True, \n",
    "                  groupby='diet',\n",
    "                  dendrogram=False, \n",
    "                  layer='scaled', \n",
    "                  vmin=-2,\n",
    "                  vmax=2,\n",
    "                  cmap='RdBu_r'\n",
    "    )\n",
    "    sc.pl.matrixplot(adata=adata_a,\n",
    "                     var_names=marker_pan_a1_a2,\n",
    "                     groupby='diet',\n",
    "                     use_raw=False, \n",
    "                     log=False,  \n",
    "                     dendrogram=False,\n",
    "                     var_group_rotation=90,\n",
    "                     swap_axes=True, \n",
    "                     show=True\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Co-expression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp == True:\n",
    "    df = pd.DataFrame(adata_astro.raw.X.toarray(), \n",
    "                      columns = adata_astro.var.index\n",
    "    )\n",
    "    df.index=adata_astro.obs.index\n",
    "    df = df.assign(diet = adata_astro.obs['diet'])\n",
    "\n",
    "    df = df.assign(gfap_lepr = ((df['Gfap'] > 0) & (df['Lepr'] > 0) & (df['Aldh1l1'] == 0))*1)\n",
    "    df = df.assign(gfap_insr = ((df['Gfap'] > 0) & (df['Insr'] > 0) & (df['Aldh1l1'] == 0))*1)\n",
    "    df = df.assign(gfap_igf1r = ((df['Gfap'] > 0) & (df['Igf1r'] > 0) & (df['Aldh1l1'] == 0))*1)\n",
    "    df = df.assign(gfap_glp1r = ((df['Gfap'] > 0) & (df['Glp1r'] > 0) & (df['Aldh1l1'] == 0))*1)\n",
    "    df = df.assign(gfap_lpl = ((df['Gfap'] > 0) & (df['Lpl'] > 0) & (df['Aldh1l1'] == 0))*1)\n",
    "    df = df.assign(gfap_slc2a1 = ((df['Gfap'] > 0) & (df['Slc2a1'] > 0) & (df['Aldh1l1'] == 0))*1)\n",
    "    df = df.assign(gfap_ucp2 = ((df['Gfap'] > 0) & (df['Ucp2'] > 0) & (df['Aldh1l1'] == 0))*1)\n",
    "    df = df.assign(gfap_ghsr = ((df['Gfap'] > 0) & (df['Ghsr'] > 0) & (df['Aldh1l1'] == 0))*1)\n",
    "    df = df.assign(gfap_slc2a2 = ((df['Gfap'] > 0) & (df['Slc2a2'] > 0) & (df['Aldh1l1'] == 0))*1)\n",
    "\n",
    "    df = df.assign(aldh_lepr = ((df['Aldh1l1'] > 0) & (df['Lepr'] > 0) & (df['Gfap'] == 0))*1)\n",
    "    df = df.assign(aldh_insr = ((df['Aldh1l1'] > 0) & (df['Insr'] > 0) & (df['Gfap'] == 0))*1)\n",
    "    df = df.assign(aldh_igf1r = ((df['Aldh1l1'] > 0) & (df['Igf1r'] > 0) & (df['Gfap'] == 0))*1)\n",
    "    df = df.assign(aldh_glp1r = ((df['Aldh1l1'] > 0) & (df['Glp1r'] > 0) & (df['Gfap'] == 0))*1)\n",
    "    df = df.assign(aldh_lpl = ((df['Aldh1l1'] > 0) & (df['Lpl'] > 0) & (df['Gfap'] == 0))*1)\n",
    "    df = df.assign(aldh_slc2a1 = ((df['Aldh1l1'] > 0) & (df['Slc2a1'] > 0) & (df['Gfap'] == 0))*1)\n",
    "    df = df.assign(aldh_ucp2 = ((df['Aldh1l1'] > 0) & (df['Ucp2'] > 0) & (df['Gfap'] == 0))*1)\n",
    "    df = df.assign(aldh_ghsr = ((df['Aldh1l1'] > 0) & (df['Ghsr'] > 0) & (df['Gfap'] == 0))*1)\n",
    "    df = df.assign(aldh_slc2a2 = ((df['Aldh1l1'] > 0) & (df['Slc2a2'] > 0) & (df['Gfap'] == 0))*1)\n",
    "\n",
    "    df = df.assign(both_lepr = ((df['Aldh1l1'] > 0) & (df['Lepr'] > 0) & (df['Gfap'] > 0))*1)\n",
    "    df = df.assign(both_insr = ((df['Aldh1l1'] > 0) & (df['Insr'] > 0) & (df['Gfap'] > 0))*1)\n",
    "    df = df.assign(both_igf1r = ((df['Aldh1l1'] > 0) & (df['Igf1r'] > 0) & (df['Gfap'] > 0))*1)\n",
    "    df = df.assign(both_glp1r = ((df['Aldh1l1'] > 0) & (df['Glp1r'] > 0) & (df['Gfap'] > 0))*1)\n",
    "    df = df.assign(both_lpl = ((df['Aldh1l1'] > 0) & (df['Lpl'] > 0) & (df['Gfap'] > 0))*1)\n",
    "    df = df.assign(both_slc2a1 = ((df['Aldh1l1'] > 0) & (df['Slc2a1'] > 0) & (df['Gfap'] > 0))*1)\n",
    "\n",
    "    df = df.assign(both_ucp2 = ((df['Aldh1l1'] > 0) & (df['Ucp2'] > 0) & (df['Gfap'] > 0))*1)\n",
    "    df = df.assign(both_ghsr = ((df['Aldh1l1'] > 0) & (df['Ghsr'] > 0) & (df['Gfap'] > 0))*1)\n",
    "    df = df.assign(both_slc2a2 = ((df['Aldh1l1'] > 0) & (df['Slc2a2'] > 0) & (df['Gfap'] > 0))*1)\n",
    "    \n",
    "    df.iloc[:,31253:31281].to_csv(\"adata_astro_upset.obs.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RNA Velocitiy and clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "spmark=['Aldoc','Aqp4','S1pr1','Slc38a1','Npy','Pcsk1n','Ins2','Il1b','Ifit3b','Flt1',\n",
    "        'Ucp2','Vamp5','S100b','Gfap','Aldh1l1','Apoe','Ift43','Tpt1','Slc1a2','Slc3a2']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### chow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp == True:\n",
    "    astro0 = adata_astro[adata_astro.obs['diet']=='chow',]\n",
    "    sc.tl.leiden(astro0,\n",
    "                 resolution=0.5,\n",
    "                 key_added='cluster0'\n",
    "    )\n",
    "    astro0.obs['cluster0n'] = astro0.obs['cluster0'].replace({'0': 'b', '1': 'd', '2': 'a', '3': 'c'})\n",
    "    sc.write(sc_settings_writedir+'astro0.h5ad', astro0)\n",
    "else:\n",
    "    astro0 = sc.read(sc_settings_writedir+'astro0.h5ad')   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp == True:  \n",
    "    adata_loom0 = scv.read(dir_scv+'MUC26030/possorted_genome_bam_Z0I20.loom', \n",
    "                           cache=True\n",
    "    )\n",
    "    astro0_v = scv.utils.merge(astro0,\n",
    "                               adata_loom0\n",
    "    )\n",
    "    scv.pl.proportions(astro0_v)\n",
    "    scv.pp.filter_and_normalize(astro0_v,\n",
    "                                min_shared_counts=20,\n",
    "                                n_top_genes=2000\n",
    "    )\n",
    "    scv.pp.moments(astro0_v,\n",
    "                   n_pcs=30, \n",
    "                   n_neighbors=30\n",
    "    )\n",
    "    scv.tl.recover_dynamics(astro0_v)\n",
    "    scv.tl.velocity(astro0_v)\n",
    "    scv.tl.velocity_graph(astro0_v)\n",
    "    \n",
    "    sc.write(sc_settings_writedir+'astro0_v.h5ad', astro0_v)\n",
    "else:\n",
    "    astro0_v = sc.read(sc_settings_writedir+'astro0_v.h5ad') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:                              #'dimgray','gray','darkgray','lightgray'   'darkslategray','steelblue','lightsteelblue','bisque'\n",
    "    cf.plot_umap_marker(astro0,\n",
    "                        ['cluster0n'],\n",
    "                        palette=['firebrick','mediumblue','dodgerblue','aqua'],\n",
    "                        use_raw=False,\n",
    "                        size=60,\n",
    "                        frameon=False\n",
    "    )\n",
    "    scv.set_figure_params()\n",
    "    scv.pl.velocity_embedding_stream(astro0_v, \n",
    "                                     basis='umap',\n",
    "                                     color='cluster0n',\n",
    "                                     use_raw=False, \n",
    "                                     size=60, \n",
    "                                     figsize=(6.5,6.5), \n",
    "                                     legend_loc='right margin',\n",
    "                                     palette=['firebrick','mediumblue','dodgerblue','aqua']\n",
    "    ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True: \n",
    "    labels = 'b', 'd', 'a', 'c'\n",
    "    sizes = astro0.obs['cluster0n'].value_counts()\n",
    "\n",
    "    fig1, ax=plt.subplots()\n",
    "    ax1.pie(sizes,\n",
    "            labels=labels, \n",
    "            autopct='%1.1f%%', \n",
    "            startangle=90, \n",
    "            pctdistance=1.3, \n",
    "            labeldistance=None,\n",
    "            colors=['mediumblue','aqua','firebrick','dodgerblue'],\n",
    "            textprops={'fontsize': 20, 'color':\"black\"}\n",
    "    )\n",
    "    ax1.axis('equal')  # Equal aspect ratio ensures that pie is drawn as a circle.\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rank genes by likelihoods per cluster/regime.\n",
    "This ranks genes by their likelihood obtained from the dynamical model grouped by clusters specified in groupby."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp == True: \n",
    "    top_genes = astro0_v.var['fit_likelihood'].sort_values(ascending=False).index\n",
    "\n",
    "    scv.tl.rank_dynamical_genes(astro0_v, \n",
    "                                groupby='cluster0'\n",
    "    )\n",
    "    df = scv.get_df(astro0_v, \n",
    "                    'rank_dynamical_genes/names'\n",
    "    )\n",
    "    df.to_csv('chow_potential_drivers.csv')\n",
    "    print(df.shape)\n",
    "    df.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    scv.tl.latent_time(astro0_v)\n",
    "    top_genes = df['1']   #['Aldoc','Aqp4','S1pr1','Slc38a1']\n",
    "    scv.pl.heatmap(astro0_v,\n",
    "                   var_names=top_genes,\n",
    "                   sortby='latent_time', \n",
    "                   col_color='cluster0n',\n",
    "                   n_convolve=100,\n",
    "                   figsize=(10,15),\n",
    "                   font_scale=1,\n",
    "                   colorbar=True,\n",
    "                   yticklabels=True\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    scv.pl.velocity(astro0_v,\n",
    "                    var_names=['Pcsk1n'], \n",
    "                    figsize=(10,11),  #'Aldoc','Pcsk1n'\n",
    "                    colorbar=True, \n",
    "                    ncols=2,\n",
    "                    color='cluster0n'\n",
    "    )\n",
    "   # astro0_v_aldh = sc.read(sc_settings_writedir+'astro0_v_aldh.h5ad')\n",
    "   # scv.pl.velocity(astro0_v_aldh,\n",
    "   #                 var_names=['Gfap','Aldoc','Pcsk1n','Aqp4','S1pr1','Pcdh15','Pak3','Aldh1l1'],\n",
    "   #                 colorbar=True,\n",
    "   #                 ncols=1, \n",
    "   #                 color='cluster0n'\n",
    "   # )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### hfd_5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp == True:\n",
    "    astro5 = adata_astro[adata_astro.obs['diet'] == 'hfd_5',]\n",
    "    sc.tl.leiden(astro5,\n",
    "                 resolution=0.5, \n",
    "                 key_added='cluster5'\n",
    "    )\n",
    "    astro5.obs['cluster5n'] = astro5.obs['cluster5'].replace({'0': 'c', '1': 'a', '2': 'b'})\n",
    "    sc.write(sc_settings_writedir+'astro5.h5ad', \n",
    "             astro5\n",
    "    )\n",
    "else:\n",
    "    astro5=sc.read(sc_settings_writedir+'astro5.h5ad') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp == True: \n",
    "    adata_loom5 = scv.read(dir_scv+'MUC26031/possorted_genome_bam_VXMFJ.loom',\n",
    "                           cache=True\n",
    "    )\n",
    "    astro5_v = scv.utils.merge(astro5, \n",
    "                               adata_loom5\n",
    "    )\n",
    "    scv.pl.proportions(astro5_v)\n",
    "    scv.pp.filter_and_normalize(astro5_v,\n",
    "                                min_shared_counts=20,\n",
    "                                n_top_genes=2000\n",
    "    )\n",
    "    scv.pp.moments(astro5_v,\n",
    "                   n_pcs=30, \n",
    "                   n_neighbors=30\n",
    "    )\n",
    "    scv.tl.recover_dynamics(astro5_v)\n",
    "    scv.tl.velocity(astro5_v)\n",
    "    scv.tl.velocity_graph(astro5_v)\n",
    "\n",
    "    sc.write(sc_settings_writedir+'astro5_v.h5ad', \n",
    "             astro5_v\n",
    "    )\n",
    "else:\n",
    "    astro5_v=sc.read(sc_settings_writedir+'astro5_v.h5ad') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True: \n",
    "    cf.plot_umap_marker(astro5,\n",
    "                        ['cluster5n'],\n",
    "                        use_raw=False,\n",
    "                        size=60, \n",
    "                        frameon=False,\n",
    "                        palette=['darkorange', 'orangered', 'darkcyan']\n",
    "    )\n",
    "    scv.set_figure_params()\n",
    "    scv.pl.velocity_embedding_stream(astro5_v,\n",
    "                                     basis='umap',\n",
    "                                     color='cluster5n',\n",
    "                                     use_raw=False,\n",
    "                                     size=60, \n",
    "                                     figsize=(6.5,6.5), \n",
    "                                     legend_loc='right margin', \n",
    "                                     palette=['darkorange', 'orangered', 'darkcyan']\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True: \n",
    "    labels = 'c', 'a', 'b'\n",
    "    sizes = astro5.obs['cluster5n'].value_counts()\n",
    "\n",
    "    fig1, ax1 = plt.subplots()\n",
    "    ax1.pie(sizes, \n",
    "            labels=labels,\n",
    "            autopct='%1.1f%%',\n",
    "            startangle=90, \n",
    "            pctdistance=1.3,\n",
    "            labeldistance=None,\n",
    "            colors=['darkcyan', 'darkorange', 'orangered'], \n",
    "            textprops={'fontsize': 20, 'color':\"black\"}\n",
    "    )\n",
    "    ax1.axis('equal')  # Equal aspect ratio ensures that pie is drawn as a circle.\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp == True: \n",
    "    top_genes = astro5_v.var['fit_likelihood'].sort_values(ascending=False).index\n",
    "    scv.tl.rank_dynamical_genes(astro5_v,\n",
    "                                groupby='leiden',\n",
    "                                n_genes=400\n",
    "    )\n",
    "    df = scv.get_df(astro5_v, \n",
    "                    'rank_dynamical_genes/names'\n",
    "    )\n",
    "    df.to_csv('hfd5_potential_drivers.csv')\n",
    "    print(df.shape)\n",
    "    df.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    scv.tl.latent_time(astro5_v)\n",
    "    top_genes = df['1'] #['Npy','Pcsk1n','Ins2','Il1b','Ifit3b','Flt1','Ucp2','Vamp5','S100b','Gfap','Aldh1l1']\n",
    "    scv.pl.heatmap(astro5_v,\n",
    "                   var_names=top_genes,\n",
    "                   sortby='latent_time',\n",
    "                   col_color='cluster5n',\n",
    "                   n_convolve=100,\n",
    "                   col_cluster=False,\n",
    "                   figsize=(10,15),\n",
    "                   font_scale=1,\n",
    "                   yticklabels=True\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    scv.pl.velocity(astro5_v, \n",
    "                    var_names=['Ucp2'],\n",
    "                    figsize=(10,12), #'Gfap', 'Clu','Ucp2','Vamp5'\n",
    "                    colorbar=True, \n",
    "                    ncols=1, \n",
    "                    color='cluster5n'\n",
    "    )\n",
    "    #astro5_v_aldh = sc.read(sc_settings_writedir+'astro5_v_aldh.h5ad')\n",
    "    #scv.pl.velocity(astro5_v_aldh, var_names=['Aldh1l1'], colorbar=True, ncols=1, color='cluster5n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### hfd_15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp == True:\n",
    "    astro15 = adata_astro[adata_astro.obs['diet']=='hfd_15',]\n",
    "    sc.tl.leiden(astro15,\n",
    "                 resolution=0.5, \n",
    "                 key_added='cluster15'\n",
    "    )\n",
    "    astro15.obs['cluster15n']=astro15.obs['cluster15'].replace({'0': 'c', '1': 'b', '2': 'a', '3': 'd'})\n",
    "    sc.write(sc_settings_writedir+'astro15.h5ad', astro15)\n",
    "else:\n",
    "    astro15 = sc.read(sc_settings_writedir+'astro15.h5ad') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp == True: \n",
    "    adata_loom15 = scv.read(dir_scv+'MUC26032/possorted_genome_bam_2UK19.loom', \n",
    "                            cache=True\n",
    "    )\n",
    "    astro15_v = scv.utils.merge(astro15, \n",
    "                                adata_loom15\n",
    "    )\n",
    "    scv.pl.proportions(astro15_v)\n",
    "    scv.pp.filter_and_normalize(astro15_v, \n",
    "                                min_shared_counts=20,\n",
    "                                n_top_genes=2000\n",
    "    )\n",
    "    scv.pp.moments(astro15_v,\n",
    "                   n_pcs=30,\n",
    "                   n_neighbors=30\n",
    "    )\n",
    "    scv.tl.recover_dynamics(astro15_v)\n",
    "    scv.tl.velocity(astro15_v)\n",
    "    scv.tl.velocity_graph(astro15_v)\n",
    "    \n",
    "    sc.write(sc_settings_writedir+'astro15_v.h5ad', astro15_v)\n",
    "else:\n",
    "    astro15_v = sc.read(sc_settings_writedir+'astro15_v.h5ad') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    cf.plot_umap_marker(astro15, \n",
    "                        ['cluster15n'], \n",
    "                        use_raw=False, \n",
    "                        size=60,\n",
    "                        frameon=False, \n",
    "                        palette=['darkorange', 'orangered', 'darkcyan', 'lightseagreen']\n",
    "    )\n",
    "    scv.pl.velocity_embedding_stream(astro15_v,\n",
    "                                     basis='umap',\n",
    "                                     color='cluster15n',\n",
    "                                     use_raw=True,\n",
    "                                     save='check.png',\n",
    "                                     size=60, \n",
    "                                     legend_loc='right margin',\n",
    "                                     palette=['darkorange', 'orangered', 'darkcyan', 'lightseagreen']\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    labels = 'c', 'b', 'a', 'd'\n",
    "    sizes=astro15.obs['cluster15n'].value_counts()\n",
    "\n",
    "    fig1, ax1=plt.subplots()\n",
    "    ax1.pie(sizes, \n",
    "            labels=labels,\n",
    "            autopct='%1.1f%%',\n",
    "            startangle=90,\n",
    "            pctdistance=1.3,\n",
    "            labeldistance=None,\n",
    "            colors=['darkcyan', 'orangered', 'darkorange', 'lightseagreen'],\n",
    "            textprops={'fontsize': 20, 'color':\"black\"}\n",
    "    )\n",
    "    ax1.axis('equal')  # Equal aspect ratio ensures that pie is drawn as a circle.\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_recomp == True:\n",
    "    top_genes = astro15_v.var['fit_likelihood'].sort_values(ascending=False).index\n",
    "\n",
    "    scv.tl.rank_dynamical_genes(astro15_v, \n",
    "                                groupby='cluster15',\n",
    "                                n_genes=400\n",
    "    )\n",
    "    df = scv.get_df(astro15_v, \n",
    "                    'rank_dynamical_genes/names'\n",
    "    )\n",
    "    df.to_csv('hfd15_potential_drivers.csv')\n",
    "    print(df.shape)\n",
    "    df.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    scv.tl.latent_time(astro15_v)\n",
    "    top_genes=df['1'] # ['Apoe','Ift43','Cldn10','Tpt1','Slc1a2','Slc3a2']\n",
    "    scv.pl.heatmap(astro15_v, \n",
    "                   var_names=top_genes, \n",
    "                   sortby='latent_time',\n",
    "                   col_color='cluster15n',\n",
    "                   n_convolve=100,\n",
    "                   col_cluster=False,\n",
    "                   figsize=(10,15),\n",
    "                   font_scale=1,\n",
    "                   yticklabels=True\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    scv.pl.velocity(astro15_v,\n",
    "                    var_names=['Apoe'],\n",
    "                    figsize=(10,12),  # 'Ift43','Apoe'\n",
    "                    colorbar=True,\n",
    "                    ncols=1, \n",
    "                    color='cluster15n'\n",
    "    )\n",
    "    #astro15_v_aldh = sc.read(sc_settings_writedir+'astro15_v_aldh.h5ad')\n",
    "    #scv.pl.velocity(astro15_v_aldh, var_names=['Aldh1l1'], colorbar=True, ncols=1, color='cluster15n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    adata_astro.obs['cluster'] = list(pd.concat([astro0.obs['cluster0'],\n",
    "                                                 astro5.obs['cluster5'],\n",
    "                                                 astro15.obs['cluster15']]))\n",
    "    adata_astro.obs['diet_cluster'] = adata_astro.obs['diet'].astype(str) + '_' + adata_astro.obs['cluster'].astype(str)\n",
    "    adata_astro.obs['diet_location'] = adata_astro.obs['diet_cluster'].replace({'chow_2': 'chow_up',\n",
    "                                                                                'chow_0': 'chow_down',\n",
    "                                                                                'chow_1': 'chow_down',\n",
    "                                                                                'chow_3': 'chow_down',\n",
    "                                                                                'hfd_5_0': 'hfd_5_down',\n",
    "                                                                                'hfd_5_1': 'hfd_5_up',\n",
    "                                                                                'hfd_5_2': 'hfd_5_up',\n",
    "                                                                                'hfd_15_0': 'hfd_15_down',\n",
    "                                                                                'hfd_15_1': 'hfd_15_up',\n",
    "                                                                                'hfd_15_2': 'hfd_15_up',\n",
    "                                                                                'hfd_15_3': 'hfd_15_down'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### All diets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    adata_loom_all = adata_loom0.concatenate([adata_loom5,\n",
    "                                              adata_loom15]\n",
    "    )\n",
    "    astro_all_v = scv.utils.merge(adata_astro, \n",
    "                                  adata_loom_all\n",
    "    )\n",
    "    scv.pl.proportions(astro_all_v)\n",
    "    scv.pp.filter_and_normalize(astro_all_v,\n",
    "                                min_shared_counts=20,\n",
    "                                n_top_genes=2000\n",
    "    )\n",
    "    scv.pp.moments(astro_all_v,\n",
    "                   n_pcs=30,\n",
    "                   n_neighbors=30\n",
    "    )\n",
    "    scv.tl.recover_dynamics(astro_all_v)\n",
    "    scv.tl.velocity(astro_all_v)\n",
    "    scv.tl.velocity_graph(astro_all_v)\n",
    "    sc.write(sc_settings_writedir+'astro_all_v.h5ad', \n",
    "             astro_all_v\n",
    "    )\n",
    "else:\n",
    "    astro_all_v = sc.read(sc_settings_writedir+'astro_all_v.h5ad') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    cf.plot_umap_marker(adata_astro, \n",
    "                        ['leiden'],\n",
    "                        use_raw=False,\n",
    "                        size=60, \n",
    "                        frameon=False\n",
    "    )\n",
    "    scv.pl.velocity_embedding_stream(astro_all_v, \n",
    "                                     basis='umap',\n",
    "                                     color='leiden',\n",
    "                                     use_raw=True,\n",
    "                                     size=60, \n",
    "                                     legend_loc='right margin'\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    top_genes=astro_all_v.var['fit_likelihood'].sort_values(ascending=False).index\n",
    "\n",
    "    scv.tl.rank_dynamical_genes(astro_all_v,\n",
    "                                groupby='leiden',\n",
    "                                n_genes=400\n",
    "    )\n",
    "    df = scv.get_df(astro_all_v, \n",
    "                    'rank_dynamical_genes/names'\n",
    "    )\n",
    "    df.to_csv('astro_all_potential_drivers.csv')\n",
    "    print(df.shape)\n",
    "    df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differential gene expression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### diet effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:    \n",
    "    sc.tl.rank_genes_groups(adata_astro,\n",
    "                            'diet',\n",
    "                            groups=['hfd_15'], \n",
    "                            reference='chow',\n",
    "                            method='t-test_overestim_var'\n",
    "    )\n",
    "    result = adata_astro.uns['rank_genes_groups']\n",
    "    groups = result['names'].dtype.names\n",
    "    df = pd.DataFrame({group + '_' + key[:1]: result[key][group]\n",
    "                       for group in groups for key in ['names', 'logfoldchanges', 'pvals_adj']})\n",
    "\n",
    "    df.columns=['gene', 'log2fc', 'pval']\n",
    "    df.to_csv(path_or_buf = dir_tables+\"chow_vs_hfd_15.csv\", sep=\"\\t\")\n",
    "    df[df['gene'].isin(marker_genes_dict['Astrocytes'])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### diet effect of up and down clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.tl.rank_genes_groups(adata_astro, \n",
    "                            'diet_location',\n",
    "                            groups=['hfd_15_down'],\n",
    "                            reference='hfd_15_up',\n",
    "                            method='t-test_overestim_var'\n",
    "    ) \n",
    "    result = adata_astro.uns['rank_genes_groups']\n",
    "    groups = result['names'].dtype.names\n",
    "    df = pd.DataFrame({group + '_' + key[:1]: result[key][group]\n",
    "                      for group in groups for key in ['names', 'logfoldchanges', 'pvals_adj']})\n",
    "\n",
    "    df.columns=['gene', 'log2fc', 'pval']\n",
    "    df.to_csv(path_or_buf=dir_tables+\"hfd_15_up_vs_hfd_15_down.csv\", \n",
    "              sep=\"\\t\"\n",
    "    )\n",
    "    df[df['gene'].isin(marker_genes_dict['Astrocytes'])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### diet effect comparing clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    sc.tl.rank_genes_groups(adata_astro, \n",
    "                            'diet_cluster',\n",
    "                            groups=['hfd_15_2'],\n",
    "                            reference='chow_2',\n",
    "                            method='t-test_overestim_var'\n",
    "    ) \n",
    "    result = adata_astro.uns['rank_genes_groups']\n",
    "    groups = result['names'].dtype.names\n",
    "    df = pd.DataFrame({group + '_' + key[:1]: result[key][group]\n",
    "                      for group in groups for key in ['names', 'logfoldchanges', 'pvals_adj']})\n",
    "\n",
    "    df.columns=['gene', 'log2fc', 'pval']\n",
    "    df.to_csv(path_or_buf = dir_tables+\"chow_2_vs_hfd_15_2.csv\", sep=\"\\t\")\n",
    "    df[df['gene'].isin(marker_genes_dict['Astrocytes'])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Astorcyte markers per diet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### chow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    cf.plot_umap_marker(astro0, \n",
    "                        ['Gfap', 'Aldh1l1'],\n",
    "                        color_map=mymap,\n",
    "                        size=30, \n",
    "                        save=\"_Gfap_Aldh1l1_chow.png\",\n",
    "                        use_raw=False, \n",
    "                        vmax=3.4056416,\n",
    "                        frameon=False\n",
    "    )\n",
    "    cf.plot_umap_marker(astro0,\n",
    "                        ['Slc1a2', 'Slc1a3', 'Aqp4', 'Gja1', 'Gjb6', 'Atp1b2'], \n",
    "                        color_map=mymap,\n",
    "                        size=30,\n",
    "                        save=\"_astrocyte_markers_chow.png\", \n",
    "                        use_raw=False,\n",
    "                        vmax=4.1836486,\n",
    "                        frameon=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### hfd 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    cf.plot_umap_marker(astro5, \n",
    "                        ['Gfap', 'Aldh1l1'],\n",
    "                        color_map=mymap,\n",
    "                        size=30, \n",
    "                        save=\"_Gfap_Aldh1l1_hfd5.png\", \n",
    "                        use_raw=False, \n",
    "                        vmax=3.4056416,\n",
    "                        frameon=False\n",
    "    )\n",
    "    cf.plot_umap_marker(astro5, \n",
    "                        ['Slc1a2', 'Slc1a3', 'Aqp4', 'Gja1', 'Gjb6', 'Atp1b2'],\n",
    "                        color_map=mymap,\n",
    "                        size=30, \n",
    "                        save=\"_astrocyte_markers_hfd5.png\", \n",
    "                        use_raw=False,\n",
    "                        vmax=4.1836486, \n",
    "                        frameon=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### hfd 15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    cf.plot_umap_marker(astro15,\n",
    "                        ['Gfap', 'Aldh1l1'],\n",
    "                        color_map=mymap,\n",
    "                        size=30, \n",
    "                        save=\"_Gfap_Aldh1l1_hfd15.png\",\n",
    "                        use_raw=False, \n",
    "                        vmax=3.4056416,\n",
    "                        frameon=False\n",
    "    )\n",
    "    cf.plot_umap_marker(astro15, \n",
    "                        ['Slc1a2', 'Slc1a3', 'Aqp4', 'Gja1', 'Gjb6', 'Atp1b2'],\n",
    "                        color_map=mymap,\n",
    "                        size=30,\n",
    "                        save=\"_astrocyte_markers_hfd15.png\", \n",
    "                        use_raw=False, \n",
    "                        vmax=4.1836486, \n",
    "                        frameon=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Counting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    astro0.var['n_cells'] = np.squeeze(np.asarray((astro0.raw.X > 0).sum(0)))\n",
    "    astro5.var['n_cells'] = np.squeeze(np.asarray((astro5.raw.X > 0).sum(0)))\n",
    "    astro15.var['n_cells'] = np.squeeze(np.asarray((astro15.raw.X > 0).sum(0)))\n",
    "\n",
    "    p0_5 = (astro5.var[astro5.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']/astro0.var[astro0.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']-1)*100\n",
    "    p5_15 = (astro15.var[astro15.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']/astro5.var[astro5.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']-1)*100\n",
    "    \n",
    "    p0_5.to_csv('substraction_cells_hfd5-chow.csv')\n",
    "    p5_15.to_csv('substraction_cells_hfd15-hfd5.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    fig = plt.figure()\n",
    "    fig.set_size_inches(9, 5)\n",
    "    ax = fig.add_axes([0,0,1,1])\n",
    "    langs = p0_5.index\n",
    "    students = p0_5\n",
    "    ax.bar(langs,students)\n",
    "    plt.xlabel(\"astrocyte markers\")\n",
    "    plt.ylabel(\"percent of increase comparing to chow (%)\")\n",
    "    plt.title(\"Substraction HFD_5 - chow\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    fig = plt.figure()\n",
    "    fig.set_size_inches(9, 5)\n",
    "    ax = fig.add_axes([0,0,1,1])\n",
    "    langs = p5_15.index\n",
    "    students = p5_15\n",
    "    ax.bar(langs,students)\n",
    "    plt.xlabel(\"astrocyte markers\")\n",
    "    plt.ylabel(\"percentage of increase comparing to HFD_5 (%)\")\n",
    "    plt.title(\"Substraction HFD_15 - HFD_5\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### up cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    astro0up = adata_astro[adata_astro.obs['diet_location'] == 'chow_up',]\n",
    "    astro5up = adata_astro[adata_astro.obs['diet_location'] == 'hfd_5_up',]\n",
    "    astro15up = adata_astro[adata_astro.obs['diet_location'] == 'hfd_15_up',]\n",
    "    \n",
    "    astro0up.var['n_cells'] = (astro0up.X > 0).sum(0)\n",
    "    astro5up.var['n_cells'] = (astro5up.X > 0).sum(0)\n",
    "    astro15up.var['n_cells'] = (astro15up.X > 0).sum(0)\n",
    "\n",
    "    p0_5up = astro5up.var[astro5up.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']/astro0up.var[astro0up.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']\n",
    "    p0_15up = astro15up.var[astro15up.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']/astro0up.var[astro0up.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']\n",
    "    p5_15up = astro15up.var[astro15up.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']/astro5up.var[astro5up.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    fig = plt.figure()\n",
    "    fig.set_size_inches(9, 5)\n",
    "    ax = fig.add_axes([0,0,1,1])\n",
    "    langs = p0_5up.index\n",
    "    students = (p0_5up-1)*100\n",
    "    ax.bar(langs,students)\n",
    "    plt.xlabel(\"astrocyte markers\")\n",
    "    plt.ylabel(\"percent of increase comparing to chow (%)\")\n",
    "    plt.title(\"Substraction HFD_5-Chow - up\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    fig = plt.figure()\n",
    "    fig.set_size_inches(9, 5)\n",
    "    ax = fig.add_axes([0,0,1,1])\n",
    "    langs = p0_15up.index\n",
    "    students = (p0_15up-1)*100\n",
    "    ax.bar(langs,students)\n",
    "    plt.xlabel(\"astrocyte markers\")\n",
    "    plt.ylabel(\"percent of increase comparing to chow (%)\")\n",
    "    plt.title(\"Substraction HFD_15-Chow - up\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    fig = plt.figure()\n",
    "    fig.set_size_inches(9, 5)\n",
    "    ax = fig.add_axes([0,0,1,1])\n",
    "    langs = p5_15up.index\n",
    "    students = (p5_15up-1)*100\n",
    "    ax.bar(langs,students)\n",
    "    plt.xlabel(\"astrocyte markers\")\n",
    "    plt.ylabel(\"percentage of increase comparing to HFD_5 (%)\")\n",
    "    plt.title(\"Substraction HFD_15-HFD_5 - up\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### down clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    astro0down = adata_astro[adata_astro.obs['diet_location'] == 'chow_down',]\n",
    "    astro5down = adata_astro[adata_astro.obs['diet_location'] == 'hfd_5_down',]\n",
    "    astro15down = adata_astro[adata_astro.obs['diet_location'] == 'hfd_15_down',]\n",
    "    \n",
    "    astro0down.var['n_cells'] = (astro0down.X > 0).sum(0)\n",
    "    astro5down.var['n_cells'] = (astro5down.X > 0).sum(0)\n",
    "    astro15down.var['n_cells'] = (astro15down.X > 0).sum(0)\n",
    "\n",
    "    p0_5down = astro5down.var[astro5down.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']/astro0down.var[astro0down.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']\n",
    "    p0_15down = astro15down.var[astro15down.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']/astro0down.var[astro0down.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']\n",
    "    p5_15down = astro15down.var[astro15down.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']/astro5down.var[astro5down.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    fig = plt.figure()\n",
    "    fig.set_size_inches(9, 5)\n",
    "    ax = fig.add_axes([0,0,1,1])\n",
    "    langs = p0_5down.index\n",
    "    students = (p0_5down-1)*100\n",
    "    ax.bar(langs,students)\n",
    "    plt.xlabel(\"astrocyte markers\")\n",
    "    plt.ylabel(\"percent of increase comparing to chow (%)\")\n",
    "    plt.title(\"Substraction HFD_5-Chow - down\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    fig = plt.figure()\n",
    "    fig.set_size_inches(9, 5)\n",
    "    ax = fig.add_axes([0,0,1,1])\n",
    "    langs = p0_15down.index\n",
    "    students = (p0_15down-1)*100\n",
    "    ax.bar(langs,students)\n",
    "    plt.xlabel(\"astrocyte markers\")\n",
    "    plt.ylabel(\"percent of increase comparing to chow (%)\")\n",
    "    plt.title(\"Substraction HFD_15-Chow - down\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    fig = plt.figure()\n",
    "    fig.set_size_inches(9, 5)\n",
    "    ax = fig.add_axes([0,0,1,1])\n",
    "    langs = p5_15down.index\n",
    "    students = (p5_15down-1)*100\n",
    "    ax.bar(langs,students)\n",
    "    plt.xlabel(\"astrocyte markers\")\n",
    "    plt.ylabel(\"percentage of increase comparing to HFD_5 (%)\")\n",
    "    plt.title(\"Substraction HFD_15b-HFD_5 - down\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### up vs down clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    p0down_0up = astro0down.var[astro0down.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']/astro0up.var[astro0up.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']\n",
    "    p5down_5up = astro5down.var[astro0down.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']/astro5up.var[astro5up.var.index.isin(marker_genes_dict['Astrocytes'])]['n_cells']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    fig = plt.figure()\n",
    "    fig.set_size_inches(9, 5)\n",
    "    ax = fig.add_axes([0,0,1,1])\n",
    "    langs = p0down_0up.index\n",
    "    students = (p0down_0up-1)*100\n",
    "    ax.bar(langs,students)\n",
    "    plt.xlabel(\"astrocyte markers\")\n",
    "    plt.ylabel(\"percentage of increase comparing to chow up (%)\")\n",
    "    plt.title(\"Substraction chow_down-chow_up\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bool_plot == True:\n",
    "    fig = plt.figure()\n",
    "    fig.set_size_inches(9, 5)\n",
    "    ax = fig.add_axes([0,0,1,1])\n",
    "    langs = p5down_5up.index\n",
    "    students = (p5down_5up-1)*100\n",
    "    ax.bar(langs,students)\n",
    "    plt.xlabel(\"astrocyte markers\")\n",
    "    plt.ylabel(\"percentage of increase comparing to hfd5 up (%)\")\n",
    "    plt.title(\"Substraction hfd5_down-hfd5_up\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}