Csppa_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“package ‘spatstat’ was built under R version 4.0.5”\n",
      "Warning message:\n",
      "“package ‘spatstat.geom’ was built under R version 4.0.5”\n",
      "Warning message:\n",
      "“package ‘spatstat.random’ was built under R version 4.0.5”\n",
      "Warning message:\n",
      "“package ‘spatstat.core’ was built under R version 4.0.5”\n",
      "Warning message:\n",
      "“package ‘nlme’ was built under R version 4.0.5”\n",
      "Warning message:\n",
      "“package ‘rpart’ was built under R version 4.0.5”\n",
      "Warning message:\n",
      "“package ‘spatstat.linnet’ was built under R version 4.0.5”\n",
      "Warning message:\n",
      "“package ‘mlr3learners’ was built under R version 4.0.5”\n",
      "Warning message:\n",
      "“package ‘raster’ was built under R version 4.0.5”\n",
      "Warning message:\n",
      "“package ‘tidyr’ was built under R version 4.0.5”\n",
      "Warning message:\n",
      "“package ‘readr’ was built under R version 4.0.5”\n",
      "Warning message:\n",
      "“package ‘dplyr’ was built under R version 4.0.5”\n",
      "Warning message:\n",
      "“package ‘MASS’ was built under R version 4.0.5”\n"
     ]
    }
   ],
   "source": [
    "inst <- suppressMessages(lapply(c(\"ggplot2\",\n",
    "                                  \"spatstat\",\n",
    "                                  \"randomForest\",\n",
    "                                  \"plotly\",\n",
    "                                  \"mlr3\",\n",
    "                                  \"mlr3learners\",\n",
    "                                  \"mlr3viz\",\n",
    "                                  \"vioplot\",\n",
    "                                  \"plot3D\",\n",
    "                                  \"raster\",\n",
    "                                  \"tidyverse\",\n",
    "                                  \"hrbrthemes\",\n",
    "                                  \"viridis\",\n",
    "                                  \"MASS\",\n",
    "                                  \"pscl\",\n",
    "                                  \"lme4\",\n",
    "                                  \"readxl\"\n",
    "                                ),\n",
    "                                library,\n",
    "                                character.only = TRUE)\n",
    ") \n",
    "\n",
    "# set the directory\n",
    "dir_out = '~/Documents/consultation/Luiza/SPPA_astrocytes/results'\n",
    "dir_tab = file.path(dir_out, 'tables')\n",
    "dir_fig = file.path(dir_out, 'figures')\n",
    "dir_dat = file.path(dir_out, 'anndata')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the location\n",
    "setwd(\"~/Documents/consultation/Luiza/SPPA_astrocytes/input/\")\n",
    "# load data and functions\n",
    "dat <- read.csv(\"spp_data_final.csv\", sep=\",\")# \"SPP_data_all.csv\"\n",
    "source(\"SPP_analysis_functions.R\")\n",
    "# define the window of the data\n",
    "window <- owin(poly = list(x = c(0, 500, 0), y = c(0, 0, 500)))\n",
    "num_size = 1.1\n",
    "num_size_sig = 1.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set the global variables"
   ]
  },
  {
   "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 = TRUE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set whether to produce plots, set to False for test runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "bool_plot = FALSE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the data objects of interest prepared for further analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    gfap = get.data(data=dat,\n",
    "                    marker='gfap_only',\n",
    "                    diet=c(0, 5, 15),\n",
    "                    mouse=1:4, \n",
    "                    arc=1:8,\n",
    "                    cal_arc=FALSE\n",
    "    ) \n",
    "    aldh = get.data(data=dat,\n",
    "                    marker='aldh_only',\n",
    "                    diet=c(0, 5, 15),\n",
    "                    mouse=1:4,\n",
    "                    arc=1:8, \n",
    "                    cal_arc=FALSE\n",
    "    )\n",
    "    both = get.data(data=dat, \n",
    "                    marker='both',\n",
    "                    diet=c(0, 5, 15), \n",
    "                    mouse=1:4, \n",
    "                    arc=1:8, \n",
    "                    cal_arc=FALSE\n",
    "    )\n",
    "    cfos = get.data(data=dat,\n",
    "                    marker='c-Fos',\n",
    "                    diet=c(0, 5, 15),\n",
    "                    mouse=1:5, \n",
    "                    arc=1:8, \n",
    "                    cal_arc=FALSE\n",
    "    )\n",
    "\n",
    "    gfap_sl_0 = get.data(data=dat,\n",
    "                         marker='gfap_only',\n",
    "                         diet=0,\n",
    "                         mouse=1:4, \n",
    "                         arc=1:8, \n",
    "                         cal_arc=TRUE\n",
    "    ) \n",
    "    gfap_sl_5 = get.data(data=dat,\n",
    "                         marker='gfap_only',\n",
    "                         diet=5, \n",
    "                         mouse=1:4,\n",
    "                         arc=1:8,\n",
    "                         cal_arc=TRUE\n",
    "    )\n",
    "    gfap_sl_15 = get.data(data=dat,\n",
    "                          marker='gfap_only',\n",
    "                          diet=15, \n",
    "                          mouse=1:4,\n",
    "                          arc=1:8,\n",
    "                          cal_arc=TRUE\n",
    "    ) \n",
    "\n",
    "    aldh_sl_0 = get.data(data=dat,\n",
    "                         marker='aldh_only', \n",
    "                         diet=0,\n",
    "                         mouse=1:4,\n",
    "                         arc=1:8,\n",
    "                         cal_arc=TRUE\n",
    "    ) \n",
    "    aldh_sl_5 = get.data(data=dat,\n",
    "                         marker='aldh_only',\n",
    "                         diet=5, \n",
    "                         mouse=1:4, \n",
    "                         arc=1:8,\n",
    "                         cal_arc=TRUE\n",
    "    )\n",
    "    aldh_sl_15 = get.data(data=dat,\n",
    "                          marker='aldh_only', \n",
    "                          diet=15,\n",
    "                          mouse=1:4, \n",
    "                          arc=1:8,\n",
    "                          cal_arc=TRUE\n",
    "    ) \n",
    "\n",
    "    both_sl_0 = get.data(data=dat, \n",
    "                         marker='both',\n",
    "                         diet=0, \n",
    "                         mouse=1:4,\n",
    "                         arc=1:8,\n",
    "                         cal_arc=TRUE\n",
    "    ) \n",
    "    both_sl_5 = get.data(data=dat,\n",
    "                         marker='both',\n",
    "                         diet=5, \n",
    "                         mouse=1:4,\n",
    "                         arc=1:8,\n",
    "                         cal_arc=TRUE\n",
    "    )\n",
    "    both_sl_15 = get.data(data=dat, \n",
    "                          marker='both',\n",
    "                          diet=15, \n",
    "                          mouse=1:4,\n",
    "                          arc=1:8,\n",
    "                          cal_arc=TRUE\n",
    "    ) \n",
    "\n",
    "    cfos_sl_0 = get.data(data=dat,\n",
    "                         marker='cfos', \n",
    "                         diet=0,\n",
    "                         mouse=1:5, \n",
    "                         arc=1:8, \n",
    "                         cal_arc=TRUE\n",
    "    ) \n",
    "    cfos_sl_5 = get.data(data=dat,\n",
    "                         marker='cfos',\n",
    "                         diet=5, \n",
    "                         mouse=1:5,\n",
    "                         arc=1:8,\n",
    "                         cal_arc=TRUE\n",
    "    )\n",
    "    cfos_sl_15 = get.data(data=dat,\n",
    "                          marker='cfos',\n",
    "                          diet=15, \n",
    "                          mouse=1:5,\n",
    "                          arc=1:8, \n",
    "                          cal_arc=TRUE\n",
    "    ) \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plots for initiall visuallization, based on the ggplot2 with four different options: \"hexagonal_heatmap\", \"scatter_plot\", \"density_contour\", \"polygon\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    gg.aspp(data=cfos[[1]],\n",
    "            type_plot=\"scatter_plot\",\n",
    "            title=\"cfos: chow diet\"\n",
    "    )\n",
    "    gg.aspp(data=cfos[[2]],\n",
    "            type_plot=\"scatter_plot\",\n",
    "            title=\"cfos: hfd_5 diet\"\n",
    "    )\n",
    "    gg.aspp(data=cfos[[3]],\n",
    "            type_plot=\"scatter_plot\",\n",
    "            title=\"cfos: hfd_15 diet\"\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"initial_visualization_possibilities.pdf\"), 8, 8)\n",
    "    gg.aspp(data=both[[1]],\n",
    "            type_plot=\"hexagonal_heatmap\",\n",
    "            title=\"Both: chow diet\"\n",
    "    )\n",
    "    gg.aspp(data=both[[1]],\n",
    "            type_plot=\"scatter_plot\",\n",
    "            title=\"Both: chow diet\"\n",
    "    )\n",
    "    gg.aspp(data=gfap[[3]],\n",
    "            type_plot=\"polygon\",\n",
    "            title=\"Gfap: hfd_15 diet\"\n",
    "    )\n",
    "    gg.aspp(data=gfap[[3]],\n",
    "            type_plot=\"density_contour\",\n",
    "            title=\"Gfap: hfd_15 diet\"\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test plots show kernel density structure, with contures, of each diet and marker combination, where on top of that we have number of cells in each square of 50 micromiters. The colours white and greean indicate (in the case of hfd5 and hfd15 days), whether particular square is non-significant or significant, respectively, comparing with chow diet. Also, coresponding overall comparison is performed using Freidman rank sum test and p-value is indicated in the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "type_plot=\"test_plot\"\n",
    "col=c(\"grey\", \"firebrick2\", \"yellow\")\n",
    "siz=6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    z_lim=c(0, 0.035)\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"test_plot_gfap_chow.pdf\"), siz, siz)\n",
    "    # Chow diet\n",
    "    dens.aspp(data=gfap[[1]], \n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window,\n",
    "              cex.axis=1.5,#0.55,\n",
    "              cex.lab=1.5,#0.6,\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              nr_slides=32\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"test_plot_gfap_hfd_5.pdf\"), siz, siz)\n",
    "    # HFD_5 diet\n",
    "    dens.aspp(data=gfap[[2]],\n",
    "              test_null=gfap_sl_0,\n",
    "              test_alt=gfap_sl_5,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              cex.axis=1.5,#0.55,\n",
    "              cex.lab=1.5,#0.6,\n",
    "              test_overall=FALSE,\n",
    "              col=col,\n",
    "              col.sig='seagreen1',\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig,\n",
    "              nr_slides=32\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"test_plot_gfap_hfd_15.pdf\"), siz, siz)\n",
    "    # HFD_15 diet\n",
    "    dens.aspp(data=gfap[[3]],\n",
    "              test_null=gfap_sl_0,\n",
    "              test_alt=gfap_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              cex.axis=1.5,#0.55,\n",
    "              cex.lab=1.5,#0.6,\n",
    "              test_overall=FALSE,\n",
    "              col=col,\n",
    "              col.sig='seagreen1',\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig,\n",
    "              nr_slides=32\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"test_plot_gfap_label.pdf\"), siz, siz)\n",
    "    # HFD_15 diet\n",
    "    dens.aspp(data=gfap[[3]],\n",
    "              test_null=gfap_sl_0,\n",
    "              test_alt=gfap_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              ribbon=TRUE,\n",
    "              cex.axis=0.9,#0.55,\n",
    "              cex.lab=.3,#0.6,\n",
    "              test_overall=FALSE,\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig,\n",
    "              nr_slides=32\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    z_lim=c(0, 0.06)\n",
    "    \n",
    "    pdf(file=file.path(dir_fig, \"test_plot_aldh_chow.pdf\"), siz, siz)\n",
    "    # Chow diet\n",
    "    dens.aspp(data=aldh[[1]], \n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window,\n",
    "              cex.axis=1.5,#0.55,\n",
    "              cex.lab=1.5,#0.6,\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              nr_slides=32\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file = file.path(dir_fig, \"test_plot_aldh_hfd_5.pdf\"), siz, siz)\n",
    "    # HFD_5 diet\n",
    "    dens.aspp(data=aldh[[2]], \n",
    "              test_null=aldh_sl_0,\n",
    "              test_alt=aldh_sl_5,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              cex.axis=1.5,#0.55,\n",
    "              cex.lab=1.5,#0.6,\n",
    "              test_overall=FALSE,\n",
    "              col=col,\n",
    "              col.sig='seagreen1',\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig,\n",
    "              nr_slides=32\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"test_plot_aldh_hfd_15.pdf\"), siz, siz)\n",
    "    # HFD_15 diet\n",
    "    dens.aspp(data=aldh[[3]], \n",
    "              test_null=aldh_sl_0,\n",
    "              test_alt=aldh_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              cex.axis=1.5,#0.55,\n",
    "              cex.lab=1.5,#0.6,\n",
    "              test_overall=FALSE,\n",
    "              col=col,\n",
    "              col.sig='seagreen1',\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig,\n",
    "              nr_slides=32\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"test_plot_aldh_label.pdf\"), siz, siz)\n",
    "    # HFD_15 diet\n",
    "    dens.aspp(data=aldh[[3]], \n",
    "              test_null=aldh_sl_0,\n",
    "              test_alt=aldh_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              ribbon=TRUE,\n",
    "              cex.axis=1,#0.55,\n",
    "              cex.lab=.3,#0.6,\n",
    "              test_overall=FALSE,\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig,\n",
    "              nr_slides=32\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    z_lim=c(0, 0.01)\n",
    "    \n",
    "    pdf(file=file.path(dir_fig, \"test_plot_both_chow.pdf\"), siz, siz)\n",
    "    # Chow diet\n",
    "    dens.aspp(data=both[[1]], \n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              cex.axis=1.5,#0.55,\n",
    "              cex.lab=1.5,#0.6,\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              nr_slides=32\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"test_plot_both_hfd_5.pdf\"), siz, siz)\n",
    "    # HFD_5 diet\n",
    "    dens.aspp(data=both[[2]], \n",
    "              test_null=both_sl_0,\n",
    "              test_alt=both_sl_5,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window,\n",
    "              cex.axis=1.5,#0.55,\n",
    "              cex.lab=1.5,#0.6,\n",
    "              test_overall=FALSE,\n",
    "              col=col,\n",
    "              col.sig='seagreen1',\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig,\n",
    "              nr_slides=32\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"test_plot_both_hfd_15.pdf\"), siz, siz)\n",
    "    # HFD_15 diet\n",
    "    dens.aspp(data=both[[3]], \n",
    "              test_null=both_sl_0,\n",
    "              test_alt=both_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              cex.axis=1.5,#0.55,\n",
    "              cex.lab=1.5,#0.6,\n",
    "              test_overall=FALSE,\n",
    "              col=col,\n",
    "              col.sig='seagreen1',\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig,\n",
    "              nr_slides=32\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"test_plot_both_label.pdf\"), siz, siz)\n",
    "    # HFD_15 diet\n",
    "    dens.aspp(data=both[[3]], \n",
    "              test_null=both_sl_0,\n",
    "              test_alt=both_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              ribbon=TRUE,\n",
    "              cex.axis=0.9,#0.55,\n",
    "              cex.lab=.3,#0.6,\n",
    "              test_overall=FALSE,\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig,\n",
    "              nr_slides=32\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subtraction plot kernel density"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plots indicate subtraction of kernel densities between all combinations of diets, always subtracting from longer to shorther hfd."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "type_plot=\"subtraction_dens_plot\"\n",
    "col=c(\"white\", \"grey\", \"firebrick2\") # blue  white\n",
    "siz=6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    z_lim=c(-0.015, 0.015)\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_dens_plot_gfap_5-0.pdf\"), siz, siz)\n",
    "    # HFD_5 - chow diet\n",
    "    dens.aspp(data=gfap[[1]],\n",
    "              data1=gfap[[2]],\n",
    "              test_null=gfap_sl_0,\n",
    "              test_alt=gfap_sl_5,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              #title=\"Gfap: hfd_5 vs. chow diet\",\n",
    "              col=col,\n",
    "              cex.axis=1.5,\n",
    "              cex.lab=1.5\n",
    "    )\n",
    "    dev.off()\n",
    "    \n",
    "    pdf(file=file.path(dir_fig, \"subtraction_dens_plot_gfap_15-0.pdf\"), siz, siz)\n",
    "    # HFD_15 - chow diet\n",
    "    dens.aspp(data=gfap[[1]],\n",
    "              data1=gfap[[3]],\n",
    "              test_null=gfap_sl_0,\n",
    "              test_alt=gfap_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              #title=\"Gfap: hfd_15 vs. chow diet\",\n",
    "              col=col,\n",
    "              cex.axis=1.5,\n",
    "              cex.lab=1.5\n",
    "    )\n",
    "    dev.off()\n",
    "    \n",
    "    pdf(file=file.path(dir_fig, \"subtraction_dens_plot_gfap_15-5.pdf\"), siz, siz)\n",
    "    dens.aspp(data=gfap[[2]],\n",
    "              data1=gfap[[3]],\n",
    "              test_null=gfap_sl_5,\n",
    "              test_alt=gfap_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              title=\"Gfap: hfd_15 vs. hfd_5\",\n",
    "              col=col,\n",
    "              cex.axis=1.5,\n",
    "              cex.lab=1.5\n",
    "    )\n",
    "    dev.off()\n",
    "    \n",
    "    pdf(file=file.path(dir_fig, \"subtraction_dens_plot_gfap_label.pdf\"), siz, siz)\n",
    "    # HFD_15 - chow diet\n",
    "    par(mar=c(5,6,4,4)+0.5, mgp=c(5,1,0))\n",
    "    dens.aspp(data=gfap[[1]],\n",
    "              data1=gfap[[3]],\n",
    "              test_null=gfap_sl_0,\n",
    "              test_alt=gfap_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              ribbon=TRUE,\n",
    "              cex.axis=1.5,\n",
    "              cex.lab=1.5,\n",
    "              #title=\"Gfap: hfd_15 vs. chow diet\",\n",
    "              col=col,\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    z_lim=c(-0.04,0.04)\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_dens_plot_aldh_5-0.pdf\"), siz, siz)\n",
    "    # HFD_5 - chow diet\n",
    "    dens.aspp(data=aldh[[1]],\n",
    "              data1=aldh[[2]],\n",
    "              test_null=aldh_sl_0,\n",
    "              test_alt=aldh_sl_5,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              #title=\"Aldh1l1: hfd_5 vs. chow diet\",\n",
    "              col=col,\n",
    "              cex.axis=1.5,\n",
    "              cex.lab=1.5\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_dens_plot_aldh_15-0.pdf\"), siz, siz)\n",
    "    # HFD_15 - chow diet\n",
    "    dens.aspp(data=aldh[[1]],\n",
    "              data1=aldh[[3]],\n",
    "              test_null=aldh_sl_0,\n",
    "              test_alt=aldh_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              #title=\"Aldh1l1: hfd_15 vs. chow diet\",\n",
    "              col=col,\n",
    "              cex.axis=1.5,\n",
    "              cex.lab=1.5\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_dens_plot_aldh_15-5.pdf\"), siz, siz)\n",
    "    # HFD_15 - chow diet\n",
    "    dens.aspp(data=aldh[[2]],\n",
    "              data1=aldh[[3]],\n",
    "              test_null=aldh_sl_5,\n",
    "              test_alt=aldh_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              title=\"Aldh1l1: hfd_15 vs. hfd_5\",\n",
    "              col=col\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_dens_plot_aldh_label.pdf\"), siz, siz)\n",
    "    par(mar=c(5,6,4,4)+0.5, mgp=c(5,1,0))\n",
    "    # HFD_15 - chow diet\n",
    "    dens.aspp(data=aldh[[1]],\n",
    "              data1=aldh[[3]],\n",
    "              test_null=aldh_sl_0,\n",
    "              test_alt=aldh_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              ribbon=TRUE,\n",
    "              cex.axis=1.5,\n",
    "              cex.lab=1.5,\n",
    "              title=\"Aldh1l1: hfd_15 vs. chow diet\",\n",
    "              col=col\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    z_lim=c(-0.00801,0.00801)\n",
    "    \n",
    "    pdf(file=file.path(dir_fig, \"subtraction_dens_plot_both_5-0.pdf\"), siz, siz)\n",
    "    # HFD_5 - chow diet\n",
    "    dens.aspp(data=both[[1]],\n",
    "              data1=both[[2]],\n",
    "              test_null=both_sl_0,\n",
    "              test_alt=both_sl_5,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              #title=\"Both: hfd_5 vs. chow diet\",\n",
    "              col=col,\n",
    "              cex.axis=1.5,\n",
    "              cex.lab=1.5\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_dens_plot_both_15-0.pdf\"), siz, siz)    \n",
    "    # HFD_15 - chow diet\n",
    "    dens.aspp(data=both[[1]],\n",
    "              data1=both[[3]],\n",
    "              test_null=both_sl_0,\n",
    "              test_alt=both_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              #title=\"Both: hfd_15 vs. chow diet\",\n",
    "              col=col,\n",
    "              cex.axis=1.5,\n",
    "              cex.lab=1.5\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_dens_plot_both_15-5.pdf\"), siz, siz)    \n",
    "    # HFD_15 - HFD_5\n",
    "    dens.aspp(data=both[[2]],\n",
    "              data1=both[[3]],\n",
    "              test_null=both_sl_5,\n",
    "              test_alt=both_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              title=\"Both: hfd_15 vs. HFD_5\",\n",
    "              col=col\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_dens_plot_both_label.pdf\"), siz, siz)    \n",
    "    # HFD_15 - chow diet\n",
    "    par(mar=c(5,6,4,4)+0.5, mgp=c(5,1,0))\n",
    "    dens.aspp(data=both[[1]],\n",
    "              data1=both[[3]],\n",
    "              test_null=both_sl_0,\n",
    "              test_alt=both_sl_15,\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              ribbon=TRUE,\n",
    "              cex.axis=1.5,\n",
    "              cex.lab=1.5,\n",
    "              title=\"Both: hfd_15 vs. chow diet\",\n",
    "              col=col\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3D kernel density plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive 3D plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=10, repr.plot.height=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    d3.aspp(data=gfap[[1]],\n",
    "            window=window,\n",
    "            type_plot=\"3D_plot_interactive\",\n",
    "            title=\"Gfap chow diet\",\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    d3.aspp(data=aldh[[1]],\n",
    "            window=window,\n",
    "            type_plot=\"3D_plot_interactive\",\n",
    "            title=\"Aldh1l1 chow diet\",\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    d3.aspp(data=both[[1]],\n",
    "            window=window,\n",
    "            type_plot=\"3D_plot_interactive\",\n",
    "            title=\"Both chow diet\",\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    d3.aspp(data=cfos[[1]],\n",
    "            window=window,\n",
    "            type_plot=\"3D_plot_interactive\",\n",
    "            title=\"Gfap chow diet\",\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    d3.aspp(data=cfos[[2]],\n",
    "            window=window,\n",
    "            type_plot=\"3D_plot_interactive\",\n",
    "            title=\"Gfap chow diet\",\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    d3.aspp(data=cfos[[3]],\n",
    "            window=window,\n",
    "            type_plot=\"3D_plot_interactive\",\n",
    "            title=\"Gfap chow diet\",\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standard 3D plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "type_plot=\"3D_plot\"\n",
    "theta=0#theta = 55\n",
    "phi=60#phi = 20\n",
    "siz=5#siz = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"3d_plot_gfap_chow.pdf\"), siz, siz)\n",
    "    d3.aspp(data=gfap[[1]],\n",
    "            type_plot=type_plot,\n",
    "            window=window,\n",
    "            title=\"Gfap chow diet\",\n",
    "            theta=theta,\n",
    "            phi=phi\n",
    "    )\n",
    "    dev.off()\n",
    "    pdf(file=file.path(dir_fig, \"3d_plot_gfap_hfd5.pdf\"), siz, siz)\n",
    "    d3.aspp(data=gfap[[2]],\n",
    "            type_plot=type_plot,\n",
    "            window=window,\n",
    "            title=\"Gfap hfd5 diet\",\n",
    "            theta=theta,\n",
    "            phi=phi\n",
    "    )\n",
    "    dev.off()\n",
    "    pdf(file=file.path(dir_fig, \"3d_plot_gfap_hfd15.pdf\"), siz, siz)\n",
    "    d3.aspp(data=gfap[[3]],\n",
    "            type_plot=type_plot,\n",
    "            window=window,\n",
    "            title=\"Gfap hfd15 diet\",\n",
    "            theta=theta,\n",
    "            phi=phi\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"3d_plot_aldh_chow.pdf\"), siz, siz)\n",
    "    d3.aspp(data=aldh[[1]],\n",
    "            type_plot=type_plot,\n",
    "            window=window,\n",
    "            title=\"Aldh1l1 chow diet\",\n",
    "            theta=theta,\n",
    "            phi=phi\n",
    "    )\n",
    "    dev.off()\n",
    "    pdf(file=file.path(dir_fig, \"3d_plot_aldh_hfd5.pdf\"), siz, siz)\n",
    "    d3.aspp(data=aldh[[2]],\n",
    "            type_plot=type_plot,\n",
    "            window=window,\n",
    "            title=\"Aldh1l1 hfd5 diet\",\n",
    "            theta=theta,\n",
    "            phi=phi\n",
    "    )\n",
    "    dev.off()\n",
    "    pdf(file=file.path(dir_fig, \"3d_plot_aldh_hfd15.pdf\"), siz, siz)\n",
    "    d3.aspp(data=aldh[[3]],\n",
    "            type_plot=type_plot,\n",
    "            window=window,\n",
    "            title=\"Aldh1l1 hfd15 diet\",\n",
    "            theta=theta,\n",
    "            phi=phi\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"3d_plot_both_chow.pdf\"), siz, siz)\n",
    "    d3.aspp(data=both[[1]],\n",
    "            type_plot=type_plot,\n",
    "            window=window,\n",
    "            title=\"Both chow diet\",\n",
    "            theta=theta,\n",
    "            phi=phi\n",
    "    )\n",
    "    dev.off()\n",
    "    pdf(file=file.path(dir_fig, \"3d_plot_both_hfd5.pdf\"), siz, siz)\n",
    "    d3.aspp(data=both[[2]],\n",
    "            type_plot=type_plot,\n",
    "            window=window,\n",
    "            title=\"Both hfd5 diet\",\n",
    "            theta=theta,\n",
    "            phi=phi\n",
    "    )\n",
    "    dev.off()\n",
    "    pdf(file=file.path(dir_fig, \"3d_plot_both_hfd15.pdf\"), siz, siz)\n",
    "    d3.aspp(data=both[[3]],\n",
    "            type_plot=type_plot,\n",
    "            window=window,\n",
    "            title=\"Both hfd15 diet\",\n",
    "            theta=theta,\n",
    "            phi=phi\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subtraction square counts plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we perform subtraction of the number of counts per each 50 micromiters square which is indicated as a color in the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "siz=4\n",
    "type_plot=\"subtraction_square_plot\"\n",
    "cex.axis=0.5\n",
    "cex.lab=0.55\n",
    "col=c(\"blue\", \"white\", \"red\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_square_plot_gfap_5-0.pdf\"), siz, siz)\n",
    "    # HFD_5 diet\n",
    "    dens.aspp(data=gfap[[2]],\n",
    "              data1=gfap[[1]],\n",
    "              type_plot=type_plot,\n",
    "              title=\"Gfap: hfd_5 vs. chow\",\n",
    "              window=window, \n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab,\n",
    "              col=col,\n",
    "              z_lim=c(-0.0082, 0.0082)\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_square_plot_gfap_15-0.pdf\"), siz, siz)\n",
    "    # HFD_15 diet\n",
    "    dens.aspp(data=gfap[[3]],\n",
    "              data1=gfap[[1]],\n",
    "              type_plot=type_plot,\n",
    "              title=\"Gfap: hfd_15 vs. chow\",\n",
    "              window=window, \n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab,\n",
    "              col=col,\n",
    "              z_lim=c(-0.0082, 0.0082)\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_square_plot_gfap_label.pdf\"), siz, siz)\n",
    "    # HFD_15 diet\n",
    "    dens.aspp(data=gfap[[3]],\n",
    "              data1=gfap[[1]],\n",
    "              type_plot=type_plot,\n",
    "              title=\"Gfap: hfd_15 vs. chow\",\n",
    "              window=window,\n",
    "              ribbon=TRUE,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab,\n",
    "              col=col,\n",
    "              z_lim=c(-0.0082, 0.0082)\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_square_plot_aldh_5-0.pdf\"), siz, siz)\n",
    "    # HFD_5 diet\n",
    "    dens.aspp(data=aldh[[2]],\n",
    "              data1=aldh[[1]],\n",
    "              type_plot=type_plot,\n",
    "              title=\"Aldh1l1: hfd_5 vs. chow\",\n",
    "              window=window, \n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab,\n",
    "              col=col,\n",
    "              z_lim=c(-0.019, 0.019)\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_square_plot_aldh_15-0.pdf\"), siz, siz)\n",
    "    # HFD_15 diet\n",
    "    dens.aspp(data=aldh[[3]],\n",
    "              data1=aldh[[1]],\n",
    "              type_plot=type_plot,\n",
    "              title=\"Aldh1l1: hfd_15 vs. chow\",\n",
    "              window=window, \n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab,\n",
    "              col=col,\n",
    "              z_lim=c(-0.019, 0.019)\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_square_plot_aldh_label.pdf\"), siz, siz)\n",
    "    # HFD_15 diet\n",
    "    dens.aspp(data=aldh[[3]],\n",
    "              data1=aldh[[1]],\n",
    "              type_plot=type_plot,\n",
    "              title=\"Aldh1l1: hfd_15 vs. chow\",\n",
    "              window=window,\n",
    "              ribbon=TRUE,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab,\n",
    "              col=col,\n",
    "              z_lim=c(-0.019, 0.019)\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_square_plot_both_5-0.pdf\"), siz, siz)\n",
    "    # HFD_5 diet\n",
    "    dens.aspp(data=both[[2]],\n",
    "              data1=both[[1]],\n",
    "              type_plot=type_plot,\n",
    "              title=\"Both: hfd_5 vs. chow\",\n",
    "              window=window, \n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab,\n",
    "              col=col,\n",
    "              z_lim=c(-0.007, 0.007)\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_square_plot_both_15-0.pdf\"), siz, siz)\n",
    "    # HFD_15 diet\n",
    "    dens.aspp(data=both[[3]],\n",
    "              data1=both[[1]],\n",
    "              type_plot=type_plot,\n",
    "              title=\"Both: hfd_15 vs. chow\",\n",
    "              window=window, \n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab,\n",
    "              col=col,\n",
    "              z_lim=c(-0.007, 0.007)\n",
    "    )\n",
    "    dev.off()\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"subtraction_square_plot_both_label.pdf\"), siz, siz)\n",
    "    # HFD_15 diet\n",
    "    dens.aspp(data=both[[3]],\n",
    "              data=both[[1]],\n",
    "              type_plot=type_plot,\n",
    "              title=\"Both: hfd_15 vs. chow\",\n",
    "              window=window,\n",
    "              ribbon=TRUE,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab,\n",
    "              col=col,\n",
    "              z_lim=c(-0.007, 0.007)\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing overall differences among the diets and markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "siz=11\n",
    "method=\"glmm\"\n",
    "cex.axis=0.9\n",
    "cex.lab=1\n",
    "col=c(\"grey\", \"firebrick2\", \"yellow\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aldh1l1 diet effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_aldh_0-5.pdf\"), siz, siz)\n",
    "    diff.aspp(null=aldh_sl_0, \n",
    "              alternative=aldh_sl_5, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"Aldh1l1: chow vs. hfd5\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab           \n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_aldh_0-15.pdf\"), siz, siz)\n",
    "    diff.aspp(null=aldh_sl_0, \n",
    "              alternative=aldh_sl_15, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"Aldh1l1: chow vs. hfd15\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_aldh_5-15.pdf\"), siz, siz)\n",
    "    diff.aspp(null=aldh_sl_5, \n",
    "              alternative=aldh_sl_15, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"Aldh1l1: hfd5 vs. hfd15\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gfap diet effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_gfap_0-5.pdf\"), siz, siz)\n",
    "    diff.aspp(null=gfap_sl_0, \n",
    "              alternative=gfap_sl_5, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"Gfap: chow vs. hfd5\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_gfap_0-15.pdf\"), siz, siz)\n",
    "    diff.aspp(null=gfap_sl_0, \n",
    "              alternative=gfap_sl_15, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"Gfap: chow vs. hfd15\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_gfap_5-15.pdf\"), siz, siz)\n",
    "    diff.aspp(null=gfap_sl_5, \n",
    "              alternative=gfap_sl_15, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"Gfap: hfd5 vs. hfd15\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Both diet effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_both_0-5.pdf\"), siz, siz)\n",
    "    diff.aspp(null=both_sl_0, \n",
    "              alternative=both_sl_5, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"Both: chow vs. hfd5\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_both_0-15.pdf\"), siz, siz)\n",
    "    diff.aspp(null=both_sl_0, \n",
    "              alternative=both_sl_15, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"Both: chow vs. hfd15\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_both_5-15.pdf\"), siz, siz)\n",
    "    diff.aspp(null=both_sl_5, \n",
    "              alternative=both_sl_15, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"Both: hfd5 vs. hfd15\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Chow - marker effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_chow_gfap-aldh.pdf\"), siz, siz)\n",
    "    diff.aspp(null=gfap_sl_0, \n",
    "              alternative=aldh_sl_0, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"chow: gfap vs. aldh\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_chow_both-gfap.pdf\"), siz, siz)\n",
    "    diff.aspp(null=both_sl_0, \n",
    "              alternative=gfap_sl_0, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"chow: both vs. gfap\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_chow_both-aldh.pdf\"), siz, siz)\n",
    "    diff.aspp(null=both_sl_0, \n",
    "              alternative=aldh_sl_0, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"chow: both vs. aldh\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HFD_5 - marker effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_hfd5_gfap-aldh.pdf\"), siz, siz)\n",
    "    diff.aspp(null=gfap_sl_5, \n",
    "              alternative=aldh_sl_5, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"hfd5: gfap vs. aldh\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_hfd5_both-gfap.pdf\"), siz, siz)\n",
    "    diff.aspp(null=both_sl_5, \n",
    "              alternative=gfap_sl_5, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"hfd5: both vs. gfap\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_hfd5_both-aldh.pdf\"), siz, siz)\n",
    "    diff.aspp(null=both_sl_5, \n",
    "              alternative=aldh_sl_5, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"hfd5: both vs. aldh\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HFD_15 - marker effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_hfd15_gfap-aldh.pdf\"), siz, siz)\n",
    "    diff.aspp(null=gfap_sl_15, \n",
    "              alternative=aldh_sl_15, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"hfd15: gfap vs. aldh\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_hfd15_both-gfap.pdf\"), siz, siz)\n",
    "    diff.aspp(null=both_sl_15, \n",
    "              alternative=gfap_sl_15, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"hfd15: both vs. gfap\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"freidman_test_plot_hfd15_both-aldh.pdf\"), siz, siz)\n",
    "    diff.aspp(null=both_sl_15, \n",
    "              alternative=aldh_sl_15, \n",
    "              window=window,\n",
    "              method=method,\n",
    "              title=\"hfd15: both vs. aldh\",\n",
    "              col=col,\n",
    "              cex.axis=cex.axis,\n",
    "              cex.lab=cex.lab\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classificatioin plots - random forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "ename": "ERROR",
     "evalue": "Error in rbind(gfap[[1]], aldh[[1]], both[[1]]): объект 'gfap' не найден\n",
     "output_type": "error",
     "traceback": [
      "Error in rbind(gfap[[1]], aldh[[1]], both[[1]]): объект 'gfap' не найден\nTraceback:\n",
      "1. rbind(gfap[[1]], aldh[[1]], both[[1]])"
     ]
    }
   ],
   "source": [
    "all_0 = rbind(gfap[[1]],aldh[[1]],both[[1]]) #,cfos[[1]]\n",
    "all_5 = rbind(gfap[[2]],aldh[[2]],both[[2]]) # ,cfos[[2]]\n",
    "all_15 = rbind(gfap[[3]],aldh[[3]],both[[3]]) # ,cfos[[3]]\n",
    "siz=30\n",
    "siz_lab=80"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    p = test.saptcorr(dataX=all_0$X,\n",
    "                      dataY=all_0$Y, \n",
    "                      dataG=all_0$Gene, \n",
    "                      method=\"Moran.I\"\n",
    "    )\n",
    "    pdf(file=file.path(dir_fig, \"classification_chow_plot.pdf\"), siz, siz)\n",
    "    theme_update(text=element_text(size=siz_lab))\n",
    "    rf0 = randF.class(dataX=all_0$X,\n",
    "                dataY=all_0$Y, \n",
    "                dataG=all_0$Gene) + \n",
    "                xlab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "                ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) + \n",
    "                scale_fill_manual(values=c(\"limegreen\", \"dodgerblue\", \"magenta\")) + #\"green4\",\"blue\",\"darkmagenta\"\n",
    "                xlim(0,520) +\n",
    "                ylim(0,525) +\n",
    "                #theme(axis.text = element_text(size=25), axis.title = element_text(size=30)) +\n",
    "                geom_text(x=400, \n",
    "                          y=400,\n",
    "                          label=paste(\"Moran.I = \", as.character(round(p$observed, 3))),\n",
    "                          size=20)\n",
    "    rf0\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"classification_chow_pie.pdf\"), siz, siz)\n",
    "    grid=rf0$data\n",
    "    rem=numeric()\n",
    "    for(i in 2:100){\n",
    "      rem=c(rem, (i*100-(i-2)):(i*100))\n",
    "    }\n",
    "    Chow=as.numeric(table(grid[!rem,5]))\n",
    "    pct=as.numeric(round(table(grid[!rem,5])/sum(table(grid[!rem,5]))*100, 2))\n",
    "    pct=paste(pct,\"%\", sep=\"\")\n",
    "    pie(Chow,\n",
    "        labels=pct,\n",
    "        col=c(\"limegreen\", \"dodgerblue\", \"magenta\"),\n",
    "        cex=13,\n",
    "        radius=0.5\n",
    "    ) # \"green4\",\"blue\",\"darkmagenta\"\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    p = test.saptcorr(dataX=all_5$X, \n",
    "                      dataY=all_5$Y, \n",
    "                      dataG=all_5$Gene, \n",
    "                      method=\"Moran.I\"\n",
    "    )\n",
    "\n",
    "    pdf(file=file.path(dir_fig, \"classification_hfd5_plot.pdf\"), siz, siz)\n",
    "    theme_update(text=element_text(size=siz_lab))\n",
    "    rf5 = randF.class(dataX=all_5$X,\n",
    "                      dataY=all_5$Y, \n",
    "                      dataG=all_5$Gene) + \n",
    "                      xlab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "                      ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) + \n",
    "                      scale_fill_manual(values=c(\"limegreen\", \"dodgerblue\", \"magenta\")) +\n",
    "                      xlim(0,520) +\n",
    "                      ylim(0,525) +\n",
    "                      geom_text(x=400, \n",
    "                                y=400,\n",
    "                                label=paste(\"Moran.I = \", as.character(round(p$observed, 3))),\n",
    "                                size=20)\n",
    "    rf5\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"classification_hfd5_pie.pdf\"), siz, siz)\n",
    "    grid=rf5$data\n",
    "    rem=numeric()\n",
    "    for(i in 2:100){\n",
    "      rem=c(rem, (i*100-(i-2)):(i*100))\n",
    "    }\n",
    "    HFD_5=as.numeric(table(grid[!rem,5]))\n",
    "    pct=as.numeric(round(table(grid[!rem,5])/sum(table(grid[!rem,5]))*100, 2))\n",
    "    pct=paste(pct,\"%\", sep=\"\")\n",
    "    pie(HFD_5,\n",
    "        labels=pct,\n",
    "        col=c(\"limegreen\", \"dodgerblue\", \"magenta\"),\n",
    "        cex=13,\n",
    "        radius=0.5\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    p = test.saptcorr(dataX=all_15$X,\n",
    "                      dataY=all_15$Y,\n",
    "                      dataG=all_15$Gene,\n",
    "                      method=\"Moran.I\"\n",
    "    )\n",
    "    pdf(file=file.path(dir_fig, \"classification_hfd15_plot.pdf\"), siz, siz)\n",
    "    theme_update(text = element_text(size=siz_lab))\n",
    "    rf15 = randF.class(dataX=all_15$X,\n",
    "                       dataY=all_15$Y, \n",
    "                       dataG=all_15$Gene) +\n",
    "                       xlab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "                       ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) + \n",
    "                       scale_fill_manual(values=c(\"limegreen\", \"dodgerblue\", \"magenta\")) + # \"green1\",\"turquoise3\",\"violet\"\n",
    "                       xlim(0, 520) +\n",
    "                       ylim(0, 525) + \n",
    "                       geom_text(x=400, \n",
    "                                 y=400, \n",
    "                                 label=paste(\"Moran.I = \", as.character(round(p$observed, 3))), \n",
    "                                 size=20)\n",
    "    rf15\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    pdf(file=file.path(dir_fig, \"classification_hfd15_pie.pdf\"), siz, siz)\n",
    "    grid =rf15$data\n",
    "    rem=numeric()\n",
    "    for(i in 2:100){\n",
    "      rem=c(rem, (i*100-(i-2)):(i*100))\n",
    "    }\n",
    "    HFD_15=as.numeric(table(grid[!rem,5]))\n",
    "    pct=as.numeric(round(table(grid[!rem,5])/sum(table(grid[!rem,5]))*100,2))\n",
    "    pct=paste(pct,\"%\", sep=\"\")\n",
    "    pie(HFD_15, \n",
    "        labels=pct,\n",
    "        col=c(\"limegreen\",\"dodgerblue\",\"magenta\"),\n",
    "        cex=13,\n",
    "        radius=0.5\n",
    "    ) # \"green1\",\"turquoise3\",\"violet\"\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    dat = data.frame(value=c(Chow, HFD_5, HFD_15),\n",
    "                      diet=c(rep(\"Chow\",3), rep(\"HFD_5\",3), rep(\"HFD_15\",3)), \n",
    "                      marker=rep(c(\"Aldh1l1_only\", \"Both\", \"Gfap_only\"),3))\n",
    "\n",
    "    ggplot(dat, aes(fill=marker, y=value, x=diet)) + \n",
    "           geom_bar(position=\"fill\", stat=\"identity\") +\n",
    "           scale_fill_manual(values=c(\"green2\", \"blue\", \"magenta\")\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classification using k-NN networks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "ename": "ERROR",
     "evalue": "Error in is(train_labels, \"character\"): объект 'all_0' не найден\n",
     "output_type": "error",
     "traceback": [
      "Error in is(train_labels, \"character\"): объект 'all_0' не найден\nTraceback:\n",
      "1. get.knn(train_labels = all_0$Gene, train_data = cbind(all_0$X, \n .     all_0$Y), to_classify_data = cbind(all_0$X, all_0$Y), k = 5, \n .     standardize = FALSE)",
      "2. is(train_labels, \"character\")"
     ]
    }
   ],
   "source": [
    "# classify the cell using k-NN\n",
    "to_check_result = get.knn(train_labels=all_0$Gene,\n",
    "                          train_data=cbind(all_0$X, all_0$Y),\n",
    "                          to_classify_data=cbind(all_0$X, all_0$Y),\n",
    "                          k=5, \n",
    "                          standardize=FALSE\n",
    ")\n",
    "table(actual=all_0$Gene, \n",
    "      predicted=to_check_result$prediction\n",
    ")\n",
    "\n",
    "correct = (all_0$Gene == to_check_result$prediction) \n",
    "print(paste( 100 * round(mean(correct), 4),\n",
    "             \"% correctly classified\"\n",
    "))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    t = test.saptcorr(dataX=all_0$X,\n",
    "                      dataY=all_0$Y, \n",
    "                      dataG=all_0$Gene, \n",
    "                      method=\"Moran.I\"\n",
    "    )\n",
    "    pdf(file=file.path(dir_fig, \"knn_classification_plot_chow.pdf\"), 10, 10)\n",
    "    theme_update(text=element_text(size=20))\n",
    "    plot.2d.knn(train_labels=all_0$Gene, \n",
    "                train_data=all_0[,6:7],\n",
    "                k=5,\n",
    "                X1=\"X\",\n",
    "                X2=\"Y\") + \n",
    "                annotate(\"text\", \n",
    "                         x=400,\n",
    "                         y=400, \n",
    "                         label=paste(\"Moran.I = \", as.character(round(t$observed, 3))),\n",
    "                         size=6) + \n",
    "                xlab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "                ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\")))\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    t = test.saptcorr(dataX=all_5$X,\n",
    "                      dataY=all_5$Y, \n",
    "                      dataG=all_5$Gene, \n",
    "                      method=\"Moran.I\"\n",
    "    )\n",
    "    pdf(file=file.path(dir_fig, \"knn_classification_plot_hfd5.pdf\"), 10, 10)\n",
    "    theme_update(text=element_text(size=20))\n",
    "    plot.2d.knn(train_labels=all_5$Gene, \n",
    "                train_data=all_5[,6:7],\n",
    "                k=5,\n",
    "                X1=\"X\",\n",
    "                X2=\"Y\") + \n",
    "                annotate(\"text\", \n",
    "                         x=400,\n",
    "                         y=400, \n",
    "                         label=paste(\"Moran.I = \", as.character(round(t$observed, 3))),\n",
    "                         size=6) + \n",
    "                xlab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "                ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\")))\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    t = test.saptcorr(dataX=all_15$X,\n",
    "                      dataY=all_15$Y, \n",
    "                      dataG=all_15$Gene, \n",
    "                      method=\"Moran.I\"\n",
    "    )\n",
    "    pdf(file=file.path(dir_fig, \"knn_classification_plot_hfd15.pdf\"), 10, 10)\n",
    "    theme_update(text=element_text(size=20))\n",
    "    plot.2d.knn(train_labels=all_15$Gene, \n",
    "                train_data=all_15[,6:7],\n",
    "                k=5,\n",
    "                X1=\"X\",\n",
    "                X2=\"Y\") + \n",
    "                annotate(\"text\", \n",
    "                         x=400,\n",
    "                         y=400, \n",
    "                         label=paste(\"Moran.I = \", as.character(round(t$observed, 3))),\n",
    "                         size=6) + \n",
    "                xlab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "                ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\")))\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Violin plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){\n",
    "    out0 = neighbor.ratio(dataX=all_0$X,\n",
    "                           dataY=all_0$Y,\n",
    "                           dataG=all_0$Gene,\n",
    "                           diam=100\n",
    "    )\n",
    "    out5 = neighbor.ratio(dataX=all_5$X,\n",
    "                           dataY=all_5$Y,\n",
    "                           dataG=all_5$Gene,\n",
    "                           diam=100\n",
    "    )\n",
    "    out15 = neighbor.ratio(dataX=all_15$X,\n",
    "                           dataY=all_15$Y,\n",
    "                           dataG=all_15$Gene,\n",
    "                           diam=100\n",
    "    )\n",
    "    pdf(file = file.path(dir_fig, \"violin_plot_neighbourhood_percentage_per_marker.pdf\"), 10, 10)\n",
    "    vioplot(out0[,1], out0[,2], out0[,3],\n",
    "            out5[,1], out5[,2], out5[,3], \n",
    "            out15[,1], out15[,2], out15[,3],\n",
    "            names=c(\"gfap_only\", \"aldh_only\", \"both\",\n",
    "                    \"gfap_only\", \"aldh_only\", \"both\",\n",
    "                    \"gfap_only\", \"aldh_only\", \"both\"),\n",
    "            col = c(rep(\"darkorange\",3), rep(\"darkorange3\",3), rep(\"darkorange4\",3))\n",
    "    )\n",
    "    title(\"Violin plot of the neighbourhood ratios\")\n",
    "\n",
    "    legend(8,0.85, \n",
    "           legend=c(\"Chow\", \"HFHSD_5\", \"HFHSD15\"), \n",
    "           col=c(\"darkorange\", \"darkorange3\", \"darkorange4\"),  \n",
    "           pch=c(19,19,19), \n",
    "           bty=\"n\", \n",
    "           pt.cex=2, \n",
    "           cex=1.2, \n",
    "           text.col=\"black\", \n",
    "           horiz=FALSE, \n",
    "           inset=c(0.1, 0.1)\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "ename": "ERROR",
     "evalue": "Error in rbind(gfap[[1]], aldh[[1]], both[[1]], cfos[[1]]): объект 'gfap' не найден\n",
     "output_type": "error",
     "traceback": [
      "Error in rbind(gfap[[1]], aldh[[1]], both[[1]], cfos[[1]]): объект 'gfap' не найден\nTraceback:\n",
      "1. rbind(gfap[[1]], aldh[[1]], both[[1]], cfos[[1]])"
     ]
    }
   ],
   "source": [
    "all_0 = rbind(gfap[[1]],aldh[[1]],both[[1]],cfos[[1]]) \n",
    "all_5 = rbind(gfap[[2]],aldh[[2]],both[[2]],cfos[[2]]) \n",
    "all_15 = rbind(gfap[[3]],aldh[[3]],both[[3]],cfos[[3]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## c-Fos -  chow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){  \n",
    "    astro=all_0[all_0$Gene != \"c-Fos\",]\n",
    "    neuro=all_0[all_0$Gene == \"c-Fos\",]\n",
    "\n",
    "    neb=character()\n",
    "    for(i in 1:nrow(neuro)){\n",
    "        h=cbind(astro$Gene, pointDistance(neuro[i,c(6,7)], astro[,c(6,7)], lonlat=F))\n",
    "        neb=cbind(neb,names(sort(table(h[order(h[,2]),][1:k])/k, decreasing=T))[1])\n",
    "    }\n",
    "    neuro$col=as.character(ifelse(neb==\"aldh_only\", \"green\", \"magenta\"))\n",
    "    plot(neuro$X,\n",
    "         neuro$Y,\n",
    "         col=neuro$col,\n",
    "         pch=19,\n",
    "         xlim=c(0,500),\n",
    "         ylim=c(0,500),\n",
    "         xlab=expression(paste(\"Distance \", \"(\", mu, \"m)\")),\n",
    "         ylab=expression(paste(\"Distance \", \"(\", mu, \"m)\")),\n",
    "         cex=1.3, \n",
    "         cex.axis=1.6,\n",
    "         cex.lab=1.4\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){  \n",
    "    options(repr.plot.width=8, repr.plot.height=8)\n",
    "    pdf(file=file.path(dir_fig, \"cfos_neighbour_chow.pdf\"), 8, 8)\n",
    "    ggplot(neuro, aes(x=X, y=Y, color=col, size=10)) +\n",
    "           geom_point() + \n",
    "           theme(axis.text=element_text(size=30), \n",
    "                 axis.title=element_text(size=30)) +\n",
    "           scale_color_manual(values = c(\"green4\",\"darkmagenta\")) + \n",
    "           xlim(0,500) +\n",
    "           ylim(0,500) + \n",
    "           xlab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) + \n",
    "           ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "           theme(legend.position=\"none\")\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){  \n",
    "    pdf(file=file.path(dir_fig, \"pie_cfos_neighbour_chow.pdf\"), 8, 8)\n",
    "    pctc=as.numeric(round(table(neuro$col)/sum(table(neuro$col))*100,2))\n",
    "    pctc=paste(pctc,\"%\", sep=\"\")\n",
    "    chow=as.numeric(table(neuro$col))\n",
    "    pie(chow,\n",
    "        labels=pctc, \n",
    "        col=c(\"green4\", \"darkmagenta\"), \n",
    "        cex=4,\n",
    "        radius=0.5\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){  \n",
    "    neb=data.frame()\n",
    "    for(i in 1:nrow(neuro)){\n",
    "        h=data.frame(astro$Gene, pointDistance(neuro[i,c(6,7)], astro[,c(6,7)], lonlat=FALSE))\n",
    "        hk=h[order(h[,2]),]\n",
    "        hk=hk[1:k,]\n",
    "        neb=rbind(neb,hk)\n",
    "    }\n",
    "    colnames(neb)=c(\"gene\",\"distance\")\n",
    "    neb[neb[,1] == 'both',1]='pair'\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){  \n",
    "    options(repr.plot.width=10, repr.plot.height=10)\n",
    "    pdf(file=file.path(dir_fig, \"boxplot_cfos_chow.pdf\"), 8, 8)\n",
    "    ggplot(neb, aes(x=gene, y=distance, fill=gene)) + \n",
    "           geom_jitter(color=\"black\", size=0.4, alpha=0.9) +\n",
    "           geom_boxplot() + \n",
    "           theme(axis.text=element_text(size=30),\n",
    "                 axis.text.x=element_blank(),\n",
    "                 axis.title=element_text(size=30)) +\n",
    "           scale_fill_manual(values=c(\"green4\",\"darkmagenta\",\"blue\")) +\n",
    "           theme(legend.position=\"none\") +\n",
    "           xlab(\"\") + \n",
    "           ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "           ylim(0,18)\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## c-Fos -  HFD_5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){  \n",
    "    astro=all_5[all_5$Gene != \"c-Fos\",]\n",
    "    neuro=all_5[all_5$Gene == \"c-Fos\",]\n",
    "\n",
    "    neb=character()\n",
    "    for(i in 1:nrow(neuro)){\n",
    "        h=cbind(astro$Gene, pointDistance(neuro[i,c(6,7)], astro[,c(6,7)], lonlat=FALSE))\n",
    "        neb=cbind(neb,names(sort(table(h[order(h[,2]),][1:k])/k, decreasing=T))[1])\n",
    "    }\n",
    "    neuro$col=as.character(ifelse(neb==\"aldh_only\", \"green\", \"magenta\"))\n",
    "    plot(neuro$X, \n",
    "         neuro$Y,\n",
    "         col=neuro$col,\n",
    "         pch=19,\n",
    "         main=\"Cfos: HFD5 diet\",\n",
    "         xlim=c(0,500),\n",
    "         ylim=c(0,500),\n",
    "         xlab=\"Distance\",\n",
    "         ylab=\"Distance\"\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){ \n",
    "    options(repr.plot.width=8, repr.plot.height=8)\n",
    "    pdf(file=file.path(dir_fig, \"cfos_neighbour_hfd5.pdf\"), 8, 8)\n",
    "    ggplot(neuro, aes(x=X, y=Y, color=col, size=10)) +\n",
    "           geom_point() + \n",
    "           theme(axis.text=element_text(size=30), \n",
    "                 axis.title=element_text(size=30)) +\n",
    "           scale_color_manual(values = c(\"limegreen\", \"magenta\")) + \n",
    "           xlim(0,500) +\n",
    "           ylim(0,500) + \n",
    "           xlab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) + \n",
    "           ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) + \n",
    "           theme(legend.position=\"none\")\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){ \n",
    "    pdf(file=file.path(dir_fig, \"pie_cfos_neighbour_hfd5.pdf\"), 8, 8)\n",
    "    pctc=as.numeric(round(table(neuro$col)/sum(table(neuro$col))*100,2))\n",
    "    pctc=paste(pctc,\"%\",sep=\"\")\n",
    "    chow=as.numeric(table(neuro$col))\n",
    "    pie(chow, \n",
    "        labels=pctc,\n",
    "        col=c(\"limegreen\",\"magenta\"),\n",
    "        cex=4,\n",
    "        radius=0.5\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){ \n",
    "    neb=data.frame()\n",
    "    for(i in 1:nrow(neuro)){\n",
    "        h=data.frame(astro$Gene, pointDistance(neuro[i,c(6,7)], astro[,c(6,7)], lonlat=FALSE))\n",
    "        hk=h[order(h[,2]),]\n",
    "        hk=hk[1:k,]\n",
    "        neb=rbind(neb,hk)\n",
    "    }\n",
    "    colnames(neb)=c(\"gene\", \"distance\")\n",
    "    neb[neb[,1] == 'both',1]='pair'\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){ \n",
    "    options(repr.plot.width=10, repr.plot.height=10)\n",
    "    pdf(file=file.path(dir_fig, \"boxplot_cfos_hfd5.pdf\"), 8, 8)\n",
    "    ggplot(neb, aes(x=gene, y=distance, fill=gene)) + \n",
    "           geom_jitter(color=\"black\", size=0.4, alpha=0.9) +\n",
    "           geom_boxplot() + \n",
    "           theme(axis.text = element_text(size=30),\n",
    "                 axis.text.x=element_blank(),\n",
    "                 axis.title = element_text(size=30)) +\n",
    "           scale_fill_manual(values = c(\"limegreen\",\"magenta\",\"dodgerblue\")) +\n",
    "           theme(legend.position = \"none\") + \n",
    "           xlab(\"\") +\n",
    "           ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "           ylim(0,18)\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## c-Fos -  HFD_15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){ \n",
    "    astro=all_15[all_15$Gene != \"c-Fos\",]\n",
    "    neuro=all_15[all_15$Gene == \"c-Fos\",]\n",
    "\n",
    "    neb=character()\n",
    "    for(i in 1:nrow(neuro)){\n",
    "        h=cbind(astro$Gene, pointDistance(neuro[i,c(6,7)], astro[,c(6,7)], lonlat=FALSE))\n",
    "        neb=cbind(neb,names(sort(table(h[order(h[,2]),][1:k])/k, decreasing=T))[1])\n",
    "    }\n",
    "    neuro$col=as.character(ifelse(neb==\"aldh_only\", \"green\", \"magenta\"))\n",
    "    plot(neuro$X, \n",
    "         neuro$Y, \n",
    "         col=neuro$col,\n",
    "         pch=19,\n",
    "         main=\"Cfos: HFD15 diet\", \n",
    "         xlim=c(0,500),\n",
    "         ylim=c(0,500),\n",
    "         xlab=\"Distance\",\n",
    "         ylab=\"Distance\"\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){ \n",
    "    options(repr.plot.width=8, repr.plot.height=8)\n",
    "    pdf(file=file.path(dir_fig, \"cfos_neighbour_hfd15.pdf\"), 8, 8)\n",
    "    ggplot(neuro, aes(x=X, y=Y, color=col, size=10)) +\n",
    "           geom_point() + \n",
    "           theme(axis.text=element_text(size=30), \n",
    "                 axis.title=element_text(size=30)) +\n",
    "           scale_color_manual(values=c(\"green1\",\"violet\")) + \n",
    "           xlim(0,500) +\n",
    "           ylim(0,500) + \n",
    "           xlab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) + \n",
    "           ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "    theme(legend.position=\"none\")\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){ \n",
    "    pdf(file=file.path(dir_fig, \"pie_cfos_neighbour_hfd15.pdf\"), 8, 8)\n",
    "    pctc=as.numeric(round(table(neuro$col)/sum(table(neuro$col))*100,2))\n",
    "    pctc=paste(pctc,\"%\",sep=\"\")\n",
    "    chow=as.numeric(table(neuro$col))\n",
    "    pie(chow,\n",
    "        labels=pctc,\n",
    "        col=c(\"green1\",\"violet\"),\n",
    "        cex=4,\n",
    "        radius=0.5\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){ \n",
    "    neb=data.frame()\n",
    "    for(i in 1:nrow(neuro)){\n",
    "        h=data.frame(astro$Gene, pointDistance(neuro[i,c(6,7)], astro[,c(6,7)], lonlat=FALSE))\n",
    "        hk=h[order(h[,2]),]\n",
    "        hk=hk[1:k,]\n",
    "        neb=rbind(neb,hk)\n",
    "    }\n",
    "    colnames(neb)=c(\"gene\",\"distance\")\n",
    "    neb[neb[,1] == 'both',1]='pair'\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){ \n",
    "    options(repr.plot.width=10, repr.plot.height=10)\n",
    "    pdf(file=file.path(dir_fig, \"boxplot_cfos_hfd15.pdf\"), 8, 8)\n",
    "    ggplot(neb, aes(x=gene, y=distance, fill=gene)) + \n",
    "           geom_jitter(color=\"black\", size=0.4, alpha=0.9) +\n",
    "           geom_boxplot() + \n",
    "           theme(axis.text=element_text(size=30),\n",
    "                 axis.text.x=element_blank(),\n",
    "                 axis.title = element_text(size=30)) +\n",
    "           scale_fill_manual(values = c(\"green1\",\"violet\",\"turquoise3\")) +\n",
    "           theme(legend.position=\"none\") +\n",
    "           xlab(\"\") +\n",
    "           ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "           ylim(0,18)\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_plot){     \n",
    "    res0 <- neighbor.ratio(dataX=all_0$X,\n",
    "                           dataY=all_0$Y,\n",
    "                           dataG=all_0$Gene,\n",
    "                           diam=60,\n",
    "                           per_marker=FALSE\n",
    "    )\n",
    "    res5 <- neighbor.ratio(dataX=all_5$X,\n",
    "                           dataY=all_5$Y,\n",
    "                           dataG=all_5$Gene,\n",
    "                           diam=60,\n",
    "                           per_marker=FALSE\n",
    "    )\n",
    "    res15 <- neighbor.ratio(dataX=all_15$X,\n",
    "                            dataY=all_15$Y,\n",
    "                            dataG=all_15$Gene,\n",
    "                            diam=60,\n",
    "                            per_marker=FALSE\n",
    "    )\n",
    "    pdf(file=file.path(dir_fig, \"violin_plot_neighbourhood_percentage_overall.pdf\"), 10, 10)\n",
    "    vioplot(res0, \n",
    "            res5,\n",
    "            res15, \n",
    "            names=c(\"Chow data\", \"5 days HFHS data\", \"15 days HFHS data\"),\n",
    "            col=\"deepskyblue\")\n",
    "    title(\"Violin plot of the neighbourhood ratios\")\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis of data for Figure 3 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "bool_recomp=FALSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "###############################################################################\n",
    "#        ARC\n",
    "###############################################################################\n",
    "if(bool_recomp){\n",
    "    dat = read_excel('~/Documents/consultation/Luiza/manuscript/revision_1/number of cells for statistics.xlsx',\n",
    "                     sheet=1\n",
    "    )\n",
    "\n",
    "    dat$Diet = as.factor(dat$Diet)\n",
    "    dat$Marker = as.factor(dat$Marker)\n",
    "    dat$Mouse = as.factor(dat$Mouse)\n",
    "    dat$ARC = as.factor(dat$ARC)\n",
    "    dat$Value = as.integer(dat$Value)\n",
    "\n",
    "    dat1 = dat[(dat$Diet == \"SC diet\" & dat$Marker == \"Aldh1L1_only\")|\n",
    "               (dat$Diet == \"5d HFHS diet\" & dat$Marker == \"Aldh1L1_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse), \n",
    "                        data=dat1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]   #odTest(glmFitNB)\n",
    "\n",
    "    dat1 = dat[(dat$Diet == \"SC diet\" & dat$Marker == \"Aldh1L1_only\")|\n",
    "               (dat$Diet == \"15d HFHS diet\" & dat$Marker == \"Aldh1L1_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=dat1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "\n",
    "    dat1 = dat[(dat$Diet == \"5d HFHS diet\" & dat$Marker == \"Aldh1L1_only\")|\n",
    "               (dat$Diet == \"15d HFHS diet\" & dat$Marker == \"Aldh1L1_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=dat1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "\n",
    "    ##################\n",
    "    dat1 = dat[(dat$Diet == \"SC diet\" & dat$Marker == \"Gfap_only\")|\n",
    "               (dat$Diet == \"5d HFHS diet\" & dat$Marker == \"Gfap_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=dat1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "\n",
    "    dat1 = dat[(dat$Diet == \"SC diet\" & dat$Marker == \"Gfap_only\")|\n",
    "               (dat$Diet == \"15d HFHS diet\" & dat$Marker == \"Gfap_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse), \n",
    "                        data=dat1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "    ##################\n",
    "    dat1 = dat[(dat$Diet == \"SC diet\" & dat$Marker == \"Double_positive\")|\n",
    "               (dat$Diet == \"5d HFHS diet\" & dat$Marker == \"Double_positive\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=dat1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "\n",
    "    dat1 = dat[(dat$Diet == \"SC diet\" & dat$Marker == \"Double_positive\")|\n",
    "               (dat$Diet == \"15d HFHS diet\" & dat$Marker == \"Double_positive\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=dat1\n",
    "    )   \n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "###############################################################################\n",
    "#           RNA scope\n",
    "###############################################################################\n",
    "if(bool_recomp){\n",
    "    scope = read_excel('~/Documents/consultation/Luiza/manuscript/revision_1/number of cells for statistics.xlsx',\n",
    "                       sheet=4\n",
    "    )\n",
    "\n",
    "    scope$Diet = as.factor(scope$Diet)\n",
    "    scope$Mouse = as.factor(scope$Mouse)\n",
    "    scope$ARC = as.factor(scope$ARC)\n",
    "\n",
    "    scope1 = scope[(scope$Diet == \"SC diet\")|(scope$Diet == \"5d HFHS diet\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=scope1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]   \n",
    "\n",
    "    scope1 = scope[(scope$Diet == \"SC diet\" )|(scope$Diet == \"15d HFHS diet\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse), \n",
    "                        data=scope1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]  \n",
    "\n",
    "    scope1 = scope[(scope$Diet == \"5d HFHS diet\")|(scope$Diet == \"15d HFHS diet\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=scope1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]  \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "###############################################################################\n",
    "#           CTX\n",
    "###############################################################################\n",
    "if(bool_recomp){\n",
    "    ctx = read_excel('~/Documents/consultation/Luiza/manuscript/revision_1/number of cells for statistics.xlsx',\n",
    "                     sheet=2\n",
    "    )\n",
    "\n",
    "    ctx$Diet = as.factor(ctx$Diet)\n",
    "    ctx$Marker = as.factor(ctx$Marker)\n",
    "    ctx$Mouse = as.factor(ctx$Mouse)\n",
    "    ctx$CTX = as.factor(ctx$CTX)\n",
    "    ctx$Value = as.integer(ctx$Value)\n",
    "\n",
    "    ctx1 = ctx[(ctx$Diet == \"SC diet\" & ctx$Marker == \"Aldh1L1_only\")|\n",
    "               (ctx$Diet == \"5d HFHS diet\" & ctx$Marker == \"Aldh1L1_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=ctx1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]   #odTest(glmFitNB)  + CTX\n",
    "\n",
    "    ctx1 = ctx[(ctx$Diet == \"SC diet\" & ctx$Marker == \"Aldh1L1_only\")|\n",
    "               (ctx$Diet == \"15d HFHS diet\" & ctx$Marker == \"Aldh1L1_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse), \n",
    "                        data=ctx1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "\n",
    "    ctx1 = ctx[(ctx$Diet == \"5d HFHS diet\" & ctx$Marker == \"Aldh1L1_only\")|\n",
    "               (ctx$Diet == \"15d HFHS diet\" & ctx$Marker == \"Aldh1L1_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse), \n",
    "                        data=ctx1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "    ##################\n",
    "    ctx1 = ctx[(ctx$Diet == \"SC diet\" & ctx$Marker == \"Gfap_only\")|\n",
    "               (ctx$Diet == \"5d HFHS diet\" & ctx$Marker == \"Gfap_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=ctx1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]   #odTest(glmFitNB)\n",
    "\n",
    "    ctx1 = ctx[(ctx$Diet == \"SC diet\" & ctx$Marker == \"Gfap_only\")|\n",
    "               (ctx$Diet == \"15d HFHS diet\" & ctx$Marker == \"Gfap_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse), \n",
    "                        data=ctx1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "\n",
    "    ##################\n",
    "    ctx1 = ctx[(ctx$Diet == \"SC diet\" & ctx$Marker == \"Double_positive\")|\n",
    "               (ctx$Diet == \"5d HFHS diet\" & ctx$Marker == \"Gfap_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=ctx1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]  #odTest(glmFitNB)\n",
    "\n",
    "    ctx1 = ctx[(ctx$Diet == \"SC diet\" & ctx$Marker == \"Double_positive\")|\n",
    "               (ctx$Diet == \"15d HFHS diet\" & ctx$Marker == \"Gfap_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=ctx1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "###############################################################################\n",
    "#           Hippo\n",
    "###############################################################################\n",
    "if(bool_recomp){\n",
    "    hippo = read_excel('~/Documents/consultation/Luiza/manuscript/revision_1/number of cells for statistics.xlsx',\n",
    "                       sheet=3\n",
    "    )\n",
    "\n",
    "    hippo$Diet = as.factor(hippo$Diet)\n",
    "    hippo$Marker = as.factor(hippo$Marker)\n",
    "    hippo$Mouse = as.factor(hippo$Mouse)\n",
    "    hippo$HIPPO = as.factor(hippo$HIPPO)\n",
    "\n",
    "    hippo1 = hippo[(hippo$Diet == \"SC diet\" & hippo$Marker == \"Aldh1L1_only\")|\n",
    "                   (hippo$Diet == \"5d HFHS diet\" & hippo$Marker == \"Aldh1L1_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=hippo1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]   #odTest(glmFitNB)\n",
    "\n",
    "    hippo1 = hippo[(hippo$Diet == \"SC diet\" & hippo$Marker == \"Aldh1L1_only\")|\n",
    "                   (hippo$Diet == \"15d HFHS diet\" & hippo$Marker == \"Aldh1L1_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=hippo1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "\n",
    "    hippo1 = hippo[(hippo$Diet == \"5d HFHS diet\" & hippo$Marker == \"Aldh1L1_only\")|\n",
    "                   (hippo$Diet == \"15d HFHS diet\" & hippo$Marker == \"Aldh1L1_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse), \n",
    "                        data=hippo1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "    ##################\n",
    "    hippo1 = hippo[(hippo$Diet == \"SC diet\" & hippo$Marker == \"Gfap_only\")|\n",
    "                   (hippo$Diet == \"5d HFHS diet\" & hippo$Marker == \"Gfap_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=hippo1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]   #odTest(glmFitNB)\n",
    "\n",
    "    hippo1 = hippo[(hippo$Diet == \"SC diet\" & hippo$Marker == \"Gfap_only\")|\n",
    "                   (hippo$Diet == \"15d HFHS diet\" & hippo$Marker == \"Gfap_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse),\n",
    "                        data=hippo1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "    ##################\n",
    "    hippo1 = hippo[(hippo$Diet == \"SC diet\" & hippo$Marker == \"Double_positive\")|\n",
    "                   (hippo$Diet == \"5d HFHS diet\" & hippo$Marker == \"Gfap_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse), \n",
    "                        data=hippo1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]   #odTest(glmFitNB)\n",
    "\n",
    "    hippo1 = hippo[(hippo$Diet == \"SC diet\" & hippo$Marker == \"Double_positive\")|\n",
    "                   (hippo$Diet == \"15d HFHS diet\" & hippo$Marker == \"Gfap_only\"),]\n",
    "    glmFitNB = glmer.nb(Value ~ Diet + (1|Mouse), \n",
    "                        data=hippo1\n",
    "    )\n",
    "    coef(summary(glmFitNB))[2,4]\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R 4.0.3",
   "language": "R",
   "name": "ir32"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}