notebooks/asppa_code.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "inst <- suppressMessages(lapply(c(\"ggplot2\",\n",
    "                                  \"spatstat\",\n",
    "                                  \"randomForest\",\n",
    "                                  \"plotly\",\n",
    "                                  \"mlr3\",\n",
    "                                  \"mlr3learners\",\n",
    "                                  \"mlr3viz\",\n",
    "                                  \"vioplot\"\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_all.csv\",\n",
    "                sep=\";\"\n",
    ")\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.2\n",
    "num_size_sig = 1.5"
   ]
  },
  {
   "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 = TRUE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the data objects of interest prepared for further analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gfap positive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aldh1l1 positive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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, 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",
    ") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Double positive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    ") "
   ]
  },
  {
   "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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 396 rows containing missing values (geom_raster).”\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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()"
   ]
  },
  {
   "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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "type_plot = \"test_plot\"\n",
    "col = c(\"grey\", \"firebrick2\", \"yellow\")\n",
    "siz = 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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",
    "              title=\"Gfap: chow diet\",\n",
    "              col=col,\n",
    "              num_size=num_size,\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",
    "              title=\"Gfap: hfd_5 vs. chow diet\",\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig\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",
    "              title=\"Gfap: hfd_15 vs. chow diet\",\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig\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.55,\n",
    "              cex.lab=0.6,\n",
    "              title=\"Gfap: hfd_15 vs. chow diet\",\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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",
    "              title=\"Aldh1l1: chow diet\",\n",
    "              col=col,\n",
    "              num_size=num_size,\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",
    "              title=\"Aldh1l1: hfd_5 vs. chow diet\",\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig\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",
    "              title=\"Aldh1l1: hfd_15 vs. chow diet\",\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig\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",
    "              title=\"Aldh1l1: hfd_15 vs. chow diet\",\n",
    "              ribbon=TRUE,\n",
    "              cex.axis=0.55,\n",
    "              cex.lab=0.6,\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Some expected counts are small; chi^2 approximation may be inaccurate”\n",
      "Warning message:\n",
      "“Some expected counts are small; chi^2 approximation may be inaccurate”\n",
      "Warning message:\n",
      "“Some expected counts are small; chi^2 approximation may be inaccurate”\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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",
    "              title=\"Both: chow diet\",\n",
    "              col=col,\n",
    "              num_size=num_size\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",
    "              title=\"Both: hfd_5 vs. chow diet\",\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig\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",
    "              title=\"Both: hfd_15 vs. chow diet\",\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig\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",
    "              title=\"Both: hfd_15 vs. chow diet\",\n",
    "              ribbon=TRUE,\n",
    "              cex.axis=0.55,\n",
    "              cex.lab=0.6,\n",
    "              col=col,\n",
    "              num_size=num_size,\n",
    "              num_size_sig=num_size_sig\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": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "type_plot = \"subtraction_dens_plot\"\n",
    "col = c(\"blue\", \"grey\", \"firebrick2\")\n",
    "siz = 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              title=\"Gfap: hfd_5 vs. chow diet\",\n",
    "              col=col,\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",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              title=\"Gfap: hfd_15 vs. chow diet\",\n",
    "              col=col,\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",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              title=\"Gfap: hfd_15 vs. hfd_5\",\n",
    "              col=col,\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",
    "    dens.aspp(data=gfap[[1]],\n",
    "              data1=gfap[[3]],\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              ribbon=TRUE,\n",
    "              cex.axis=0.6,\n",
    "              cex.lab=0.7,\n",
    "              title=\"Gfap: hfd_15 vs. chow diet\",\n",
    "              col=col,\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              title=\"Aldh1l1: hfd_5 vs. chow diet\",\n",
    "              col=col\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",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              title=\"Aldh1l1: hfd_15 vs. chow diet\",\n",
    "              col=col\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[[1]],\n",
    "              data1=aldh[[3]],\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",
    "    # HFD_15 - chow diet\n",
    "    dens.aspp(data=aldh[[1]],\n",
    "              data1=aldh[[3]],\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              ribbon=TRUE,\n",
    "              cex.axis=0.6,\n",
    "              cex.lab=0.7,\n",
    "              title=\"Aldh1l1: hfd_15 vs. chow diet\",\n",
    "              col=col\n",
    "    )\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              title=\"Both: hfd_5 vs. chow diet\",\n",
    "              col=col\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",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              title=\"Both: hfd_15 vs. chow diet\",\n",
    "              col=col\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[[1]],\n",
    "              data1=both[[3]],\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",
    "    dens.aspp(data=both[[1]],\n",
    "              data1=both[[3]],\n",
    "              z_lim=z_lim,\n",
    "              type_plot=type_plot,\n",
    "              window=window, \n",
    "              ribbon=TRUE,\n",
    "              cex.axis=0.6,\n",
    "              cex.lab=0.7,\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": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "type_plot = \"3D_plot\"\n",
    "theta = 55\n",
    "phi = 20\n",
    "siz = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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()"
   ]
  },
  {
   "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": 19,
   "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": 20,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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",
    "       data1=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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing differences among the diets and markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "siz = 11\n",
    "method = \"friedman\"\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": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Some expected counts are small; chi^2 approximation may be inaccurate”\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Some expected counts are small; chi^2 approximation may be inaccurate”\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Some expected counts are small; chi^2 approximation may be inaccurate”\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_0 = rbind(gfap[[1]],aldh[[1]],both[[1]])\n",
    "all_5 = rbind(gfap[[2]],aldh[[2]],both[[2]])\n",
    "all_15 = rbind(gfap[[3]],aldh[[3]],both[[3]])\n",
    "siz = 30\n",
    "siz_lab = 40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO  [16:02:06.977] Applying learner 'classif.rpart' on task 'spirals' (iter 1/1) \n",
      "INFO  [16:02:07.479] Applying learner 'classif.ranger' on task 'spirals' (iter 1/1) \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 199 rows containing missing values (geom_raster).”\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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",
    "randF.class(dataX=all_0$X, dataY=all_0$Y, dataG=all_0$Gene) + \n",
    "            xlab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "            ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) + \n",
    "            scale_fill_manual(values=c(\"green1\", \"yellow\", \"magenta\")) +\n",
    "            xlim(0,520) +\n",
    "            ylim(0,500) +\n",
    "            geom_text(x=400, \n",
    "                      y=400,\n",
    "                      label=paste(\"Moran.I=\", as.character(round(p$observed, 3))),\n",
    "                      size=20)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO  [16:02:11.654] Applying learner 'classif.rpart' on task 'spirals' (iter 1/1) \n",
      "INFO  [16:02:11.972] Applying learner 'classif.ranger' on task 'spirals' (iter 1/1) \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 199 rows containing missing values (geom_raster).”\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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",
    "randF.class(dataX=all_5$X, dataY=all_5$Y, dataG=all_5$Gene) + \n",
    "            xlab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) +\n",
    "            ylab(expression(paste(\"Distance \", \"(\", mu, \"m)\"))) + \n",
    "            scale_fill_manual(values=c(\"green1\", \"yellow\", \"magenta\")) +\n",
    "            xlim(0,520) +\n",
    "            ylim(0,500) +\n",
    "            geom_text(x=400, \n",
    "                      y=400,\n",
    "                      label=paste(\"Moran.I=\", as.character(round(p$observed, 3))),\n",
    "                      size=20)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO  [16:02:15.556] Applying learner 'classif.rpart' on task 'spirals' (iter 1/1) \n",
      "INFO  [16:02:15.814] Applying learner 'classif.ranger' on task 'spirals' (iter 1/1) \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 199 rows containing missing values (geom_raster).”\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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",
    "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(\"green1\", \"yellow\", \"magenta\")) +\n",
    "            xlim(0, 520) +\n",
    "            ylim(0, 500) + \n",
    "            geom_text(x=400, \n",
    "                      y=400, \n",
    "                      label=paste(\"Moran.I=\", as.character(round(p$observed, 3))), \n",
    "                      size=20)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classification using k-NN networks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "           predicted\n",
       "actual      aldh_only both gfap_only\n",
       "  aldh_only       938   22       221\n",
       "  both             89   27        50\n",
       "  gfap_only       250   16       548"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"70.01 % correctly classified\"\n"
     ]
    }
   ],
   "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": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{pdf:} 2"
      ],
      "text/markdown": [
       "**pdf:** 2"
      ],
      "text/plain": [
       "pdf \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Violin plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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=F , \n",
    "           inset=c(0.1, 0.1))\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "    )\n",
    "    title(\"Violin plot of the neighbourhood ratios\")\n",
    "    dev.off()\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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
}