R_visual_sc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "standing-gallery",
   "metadata": {},
   "source": [
    "### Load the libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "premier-integrity",
   "metadata": {},
   "outputs": [],
   "source": [
    "inst <- suppressMessages(lapply(c('SummarizedExperiment',\n",
    "                                  'pheatmap',\n",
    "                                  'ggplot2',\n",
    "                                  'EnhancedVolcano',\n",
    "                                  'VennDiagram',\n",
    "                                  'org.Mm.eg.db',\n",
    "                                  'clusterProfiler',\n",
    "                                  'readxl',\n",
    "                                  'tidyverse'),\n",
    "                                 library,\n",
    "                                character.only=TRUE)\n",
    ") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "balanced-hepatitis",
   "metadata": {},
   "source": [
    "### Global variable"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "opposed-electron",
   "metadata": {},
   "source": [
    "All embeddings and clusterings can be saved and loaded into this script. Be careful with over-writing cluster caches as soon as cell type annotation has started as cluster labels may be shuffled."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "clean-poker",
   "metadata": {},
   "source": [
    "Set whether to produce plots, set to False for test runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "willing-influence",
   "metadata": {},
   "outputs": [],
   "source": [
    "bool_plot = TRUE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "organized-geometry",
   "metadata": {},
   "source": [
    "### Set the working directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "activated-bennett",
   "metadata": {},
   "outputs": [],
   "source": [
    "setwd('~/Documents/consultation/Luiza/single_cell/results/figures')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "resistant-portrait",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 2 rows containing missing values (geom_vline).”\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=\"volcano_aldh_marker_chow.pdf\", width=10, height=10)\n",
    "    a0 <- read.csv('~/Documents/consultation/Luiza/single_cell/results/tables/aldh-gfap_both-chow.csv',\n",
    "                   header=TRUE, \n",
    "                   sep=\"\\t\"\n",
    "    )\n",
    "    a0 = a0[complete.cases(a0),]\n",
    "    a0 = a0[order(a0$pval),]\n",
    "    a0s = a0[a0$pval < 0.05,]\n",
    "    a0$log2fc <- ifelse(a0$log2fc > 4.5,\n",
    "                        4.5,\n",
    "                        a0$log2fc\n",
    "    )\n",
    "    a0$log2fc <- ifelse(a0$log2fc < -4.5,\n",
    "                        -4.5, \n",
    "                        a0$log2fc\n",
    "    )\n",
    "    EnhancedVolcano(a0,\n",
    "                    lab=a0$gene,\n",
    "                    pCutoff=0.05,\n",
    "                    FCcutoff=10,\n",
    "                    axisLabSize=30,\n",
    "                    title=NULL,\n",
    "                    subtitle=NULL,\n",
    "                    caption=NULL,\n",
    "                    x='log2fc',\n",
    "                    y='pval',\n",
    "                    ylim=c(0, 7.1),\n",
    "                    xlim=c(-4.5, 4.5),\n",
    "                    selectLab=a0$gene[1:25],\n",
    "                    labSize=8,\n",
    "                    col=c(\"grey30\", \"grey30\", \"red2\", \"red2\"),\n",
    "                    legendLabels=c('NS', expression(Log[2]~FC), 'p-value<0.05', expression(p-value~and~log[2]~FC)),\n",
    "                    legendLabSize=30\n",
    "    )\n",
    "dev.off()                    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "random-service",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 2 rows containing missing values (geom_vline).”\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=\"volcano_gfap_marker_chow.pdf\", width=10, height=10)\n",
    "    g0 <- read.csv('~/Documents/consultation/Luiza/single_cell/results/tables/gfap-aldh_both-chow.csv',\n",
    "                   header=TRUE,\n",
    "                   sep=\"\\t\"\n",
    "    )\n",
    "    g0 = g0[complete.cases(g0), ]\n",
    "    g0 = g0[order(g0$pval),]\n",
    "    g0s = g0[g0$pval < 0.05,]\n",
    "    g0$log2fc <- ifelse(g0$log2fc > 4.5,\n",
    "                        4.5,\n",
    "                        g0$log2fc\n",
    "    )\n",
    "    g0$log2fc <- ifelse(g0$log2fc < -4.5,\n",
    "                        -4.5,\n",
    "                        g0$log2fc\n",
    "    )\n",
    "    EnhancedVolcano(g0,\n",
    "                    lab=g0$gene,\n",
    "                    pCutoff=0.05,\n",
    "                    FCcutoff=10,\n",
    "                    axisLabSize=30,\n",
    "                    title=NULL,\n",
    "                    subtitle=NULL,\n",
    "                    caption=NULL,\n",
    "                    x='log2fc',\n",
    "                    y='pval',\n",
    "                    ylim=c(0, 7.1),\n",
    "                    xlim=c(-4.5, 4.5),\n",
    "                    selectLab=g0$gene[1:25],\n",
    "                    labSize=8,\n",
    "                    col=c(\"grey30\", \"grey30\", \"red2\", \"red2\"),\n",
    "                    legendLabels=c('NS', expression(Log[2]~FC), 'p-value<0.05', expression(p-value~and~log[2]~FC)),\n",
    "                    legendLabSize=30\n",
    "    )\n",
    "dev.off() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "pretty-environment",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 2 rows containing missing values (geom_vline).”\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=\"volcano_both_marker_chow.pdf\", width=10, height=10) \n",
    "    b0 <- read.csv('~/Documents/consultation/Luiza/single_cell/results/tables/both-aldh_gfap-chow.csv',\n",
    "                   header=TRUE,\n",
    "                   sep=\"\\t\"\n",
    "    )\n",
    "    b0 = b0[complete.cases(b0),]\n",
    "    b0 = b0[order(b0$pval),]\n",
    "    b0s = b0[b0$pval < 0.05,]\n",
    "    b0$log2fc <- ifelse(b0$log2fc > 4.5,\n",
    "                        4.5,\n",
    "                        b0$log2fc\n",
    "    )\n",
    "    b0$log2fc <- ifelse(b0$log2fc < -4.5,\n",
    "                        -4.5, \n",
    "                        b0$log2fc\n",
    "    )    \n",
    "    EnhancedVolcano(b0,\n",
    "                    lab=b0$gene,\n",
    "                    pCutoff=0.05,\n",
    "                    FCcutoff=10,\n",
    "                    axisLabSize=30,\n",
    "                    title=NULL,\n",
    "                    subtitle=NULL,\n",
    "                    caption=NULL,\n",
    "                    x='log2fc',\n",
    "                    y='pval',\n",
    "                    ylim=c(0, 7.1),\n",
    "                    xlim=c(-4.5, 4.5),\n",
    "                    selectLab=b0$gene[1:25],\n",
    "                    labSize=8,\n",
    "                    col=c(\"grey30\", \"grey30\", \"red2\", \"red2\"),\n",
    "                    legendLabels=c('NS', expression(Log[2]~FC), 'p-value<0.05', expression(p-value~and~log[2]~FC)),\n",
    "                    legendLabSize=30\n",
    "    )\n",
    "dev.off() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "dated-bulgarian",
   "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": [
    "pdf(file=\"venn_marker.pdf\", width=8, height=8)\n",
    "    vp <- venn.diagram(list(' '=a0s$gene,\n",
    "                            ' '=g0s$gene,\n",
    "                            ' '=b0s$gene),\n",
    "                       fill=c(\"green4\", \"darkmagenta\", \"blue\"),\n",
    "                       filename=NULL, \n",
    "                       cex=3.5,\n",
    "                       col=\"transparent\",\n",
    "                       cat.cex=0.01,\n",
    "                       margin=0.01,\n",
    "                       main=NULL,\n",
    "                       main.cex=3);\n",
    "    grid.draw(vp)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "dense-cylinder",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(setdiff(a0s$gene, union(g0s$gene, b0s$gene)), fromType = \"SYMBOL\", :\n",
      "“1.77% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(setdiff(g0s$gene, union(a0s$gene, b0s$gene)), fromType = \"SYMBOL\", :\n",
      "“3.34% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(setdiff(b0s$gene, union(a0s$gene, g0s$gene)), fromType = \"SYMBOL\", :\n",
      "“2.22% of input gene IDs are fail to map...”\n"
     ]
    }
   ],
   "source": [
    "sga <- bitr(setdiff(a0s$gene, union(g0s$gene, b0s$gene)), \n",
    "            fromType=\"SYMBOL\",\n",
    "            toType=\"ENTREZID\",\n",
    "            OrgDb=org.Mm.eg.db\n",
    ")\n",
    "ac_g <- enrichGO(sga[,2],\n",
    "                 'org.Mm.eg.db',\n",
    "                 ont=\"BP\", \n",
    "                 pvalueCutoff=0.1\n",
    ")\n",
    "sgg <- bitr(setdiff(g0s$gene, union(a0s$gene, b0s$gene)),\n",
    "            fromType=\"SYMBOL\",\n",
    "            toType=\"ENTREZID\",\n",
    "            OrgDb=org.Mm.eg.db\n",
    ")\n",
    "gc_g <- enrichGO(sgg[,2],\n",
    "                 'org.Mm.eg.db',\n",
    "                 ont=\"BP\",\n",
    "                 pvalueCutoff=0.1\n",
    ")\n",
    "sgb <- bitr(setdiff(b0s$gene, union(a0s$gene, g0s$gene)),\n",
    "            fromType=\"SYMBOL\",\n",
    "            toType=\"ENTREZID\",\n",
    "            OrgDb=org.Mm.eg.db\n",
    ")\n",
    "bc_g <- enrichGO(sgb[,2],\n",
    "                 'org.Mm.eg.db',\n",
    "                 ont=\"BP\",\n",
    "                 pvalueCutoff=0.1\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "divided-victoria",
   "metadata": {},
   "outputs": [],
   "source": [
    "go=c(\"regulation of TOR signaling\",\n",
    "     \"leptin-mediated signaling pathway\",\n",
    "     \"regulation of ERK1 and ERK2 cascade\",\n",
    "     \"glycogen catabolic process\",\n",
    "     \"polysaccharide catabolic process\",\n",
    "     \"carbohydrate homeostasis\",\n",
    "     \"regulation of lipid metabolic process\",\n",
    "     \"triglyceride metabolic process\",\n",
    "     \"intracellular lipid transport\",\n",
    "     \"reactive oxygen species metabolic process\",\n",
    "     \"response to oxidative stress\",\n",
    "     \"response to endoplasmic reticulum stress\",\n",
    "     \"reactive oxygen species biosynthetic process\",\n",
    "     \"mitochondrion organization\",\n",
    "     \"mitochondrial respiratory chain complex assembly\",\n",
    "     \"mitochondrial transmembrane transport\",\n",
    "     \"glial cell proliferation\",\n",
    "     \"interleukin-2 secretion\",\n",
    "     \"positive regulation of toll-like receptor signaling pathway\",\n",
    "     \"negative regulation of inflammatory response\",\n",
    "     \"positive regulation of cell-substrate adhesion\",\n",
    "     \"regulation of cell junction assembly\",\n",
    "     \"inositol phosphate metabolic process\",\n",
    "     \"calcium ion transmembrane transport\",\n",
    "     \"cytosolic calcium ion transport\",\n",
    "     \"axonogenesis\",\n",
    "     \"synapse organization\",\n",
    "     \"synapse assembly\",\n",
    "     \"regulation of glutamate receptor signaling pathway\",\n",
    "     \"positive regulation of synaptic transmission, glutamatergic\",\n",
    "     \"regulation of angiogenesis\",\n",
    "     \"regulation of vasculature development\",\n",
    "     \"response to hypoxia\",\n",
    "     \"blood vessel remodeling\",\n",
    "     \"sprouting angiogenesis\"\n",
    ") \n",
    "# \"negative regulation of cytokine production involved in inflammatory response\","
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "absent-simulation",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=15, repr.plot.height=15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "radical-failure",
   "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": [
    "pdf(file=\"enrich_path_marker.pdf\", width=15, height=15) \n",
    "    enr_a = ac_g@result[ac_g@result[,5] < 0.01, c(2, 5, 6, 8, 9)]\n",
    "    enr_a$Comparison = 'aldh'\n",
    "    enr_g = gc_g@result[gc_g@result[,5] < 0.01, c(2, 5, 6, 8, 9)]\n",
    "    enr_g$Comparison = 'gfap'\n",
    "    enr_b = bc_g@result[bc_g@result[,5] < 0.01, c(2, 5, 6, 8, 9)]\n",
    "    enr_b$Comparison = 'pair'\n",
    "\n",
    "    all = rbind(enr_a, enr_g, enr_b)\n",
    "    all = all[all$Description%in%go,]\n",
    "    #write.csv(all, file=\"GO_sc_marker_effect-chow.csv\")\n",
    "    \n",
    "    colnames(all) = c('Description', 'pvalue', 'p.adjust', 'geneID', 'count', 'Comparison')\n",
    "    ggplot(all, \n",
    "           aes(Comparison, y=factor(Description, level=go[length(go):1]), col=pvalue, size=count)) + \n",
    "           geom_point() + \n",
    "           scale_colour_gradient(low=\"red\", high=\"blue\") +\n",
    "           scale_size_continuous(range=c(3, 10)) +\n",
    "           xlab(\"\") + \n",
    "           ylab(\"\") + \n",
    "           theme(axis.text.x=element_blank()) +\n",
    "           theme(legend.position=\"none\") +\n",
    "           theme(text=element_text(size=35))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "intermediate-rolling",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(a0$gene, fromType = \"SYMBOL\", toType = \"ENTREZID\", OrgDb = org.Mm.eg.db):\n",
      "“3.59% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(g0$gene, fromType = \"SYMBOL\", toType = \"ENTREZID\", OrgDb = org.Mm.eg.db):\n",
      "“3.59% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(b0$gene, fromType = \"SYMBOL\", toType = \"ENTREZID\", OrgDb = org.Mm.eg.db):\n",
      "“3.59% of input gene IDs are fail to map...”\n"
     ]
    }
   ],
   "source": [
    "sg.a0.all <- bitr(a0$gene,\n",
    "                  fromType=\"SYMBOL\",\n",
    "                  toType=\"ENTREZID\",\n",
    "                  OrgDb=org.Mm.eg.db\n",
    ")\n",
    "sg.g0.all <- bitr(g0$gene,\n",
    "                  fromType=\"SYMBOL\",\n",
    "                  toType=\"ENTREZID\",\n",
    "                  OrgDb=org.Mm.eg.db\n",
    ")\n",
    "sg.b0.all <- bitr(b0$gene,\n",
    "                  fromType=\"SYMBOL\",\n",
    "                  toType=\"ENTREZID\",\n",
    "                  OrgDb=org.Mm.eg.db\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "bizarre-eating",
   "metadata": {},
   "outputs": [],
   "source": [
    "all = rbind(enr_a, enr_g, enr_b)\n",
    "all = all[all$Description%in%go,]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "worthy-chancellor",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=7, repr.plot.height=15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "compressed-academy",
   "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": [
    "pdf(file=\"enrichdot_aldh_marker.pdf\", width=7, height=15)   \n",
    "    df_total.a=data.frame()\n",
    "    for(i in 1:nrow(all)){   \n",
    "        p=as.data.frame(a0[which(a0$gene%in%\n",
    "                           sg.a0.all[sg.a0.all[,2]%in%unlist(strsplit(all$geneID[i],\"/\")), 1]), c(2, 3, 5)])\n",
    "        p$pathway=all$Description[i]\n",
    "        p$padjPath=all$p.adjust[i]\n",
    "        p$count=all$Count[i]\n",
    "        p$Comparison=all$Comparison[i]\n",
    "        df_total.a=rbind(df_total.a, p)\n",
    "    }\n",
    "    df_total.a$significance <- ifelse(df_total.a$pval < 0.05,\n",
    "                                      \"significant\",\n",
    "                                      \"not significant\"\n",
    "    )\n",
    "    ggplot(df_total.a, aes(log2fc, y=factor(pathway, level=go[length(go):1]))) + \n",
    "           geom_point(aes(colour=significance, size=-log10(pval)), alpha=0.5) + \n",
    "           scale_color_manual(values=c(\"gray85\", \"red\"))+ \n",
    "           scale_size_continuous(range=c(0, 7)) +\n",
    "           geom_vline(xintercept=0) +\n",
    "           theme(axis.title=element_text(size=20),\n",
    "                 title=element_text(size=20), \n",
    "                 text=element_text(size=35)) +\n",
    "           theme(legend.position=\"none\") + \n",
    "           theme(axis.text.y=element_blank()) +\n",
    "           ylab(\"\") +\n",
    "           xlab(\"\") +\n",
    "           xlim(-4.5, 4.5)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "disturbed-decision",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_a = df_total.a[,c(1:4,7)]\n",
    "df_a$pval <- ifelse(df_a$pval < 0.05,\n",
    "                    \"sig\",\n",
    "                    'no_sig'\n",
    ")\n",
    "df_a$log2fc <- ifelse(df_a$log2fc < 0,\n",
    "                      \"down\",\n",
    "                      'up'\n",
    ")\n",
    "df_a = df_a[order(df_a$pathway),]\n",
    "write.csv(df_a, 'marker_aldh1l1_chow_GO.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "resistant-winner",
   "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": [
    "pdf(file=\"enrichdot_gfap_marker.pdf\", width=7, height=15)   \n",
    "    df_total.g=data.frame()\n",
    "    for(i in 1:nrow(all)){   \n",
    "        p = as.data.frame(g0[which(g0$gene%in%\n",
    "                           sg.g0.all[sg.g0.all[,2]%in%unlist(strsplit(all$geneID[i],\"/\")), 1]), c(2, 3, 5)])\n",
    "        p$pathway = all$Description[i]\n",
    "        p$padjPath = all$p.adjust[i]\n",
    "        p$count = all$Count[i]\n",
    "        p$Comparison = all$Comparison[i]\n",
    "        df_total.g = rbind(df_total.g, p)\n",
    "    }\n",
    "    df_total.g$significance <- ifelse(df_total.g$pval < 0.05,\n",
    "                                      \"significant\",\n",
    "                                      \"not significant\"\n",
    "    )\n",
    "    ggplot(df_total.g, aes(log2fc, y=factor(pathway, level=go[length(go):1]))) + \n",
    "           geom_point(aes(colour=significance, size=-log10(pval)), alpha=0.5) + \n",
    "           scale_color_manual(values=c(\"gray85\", \"red\")) + \n",
    "           scale_size_continuous(range=c(0, 8)) +\n",
    "           geom_vline(xintercept=0) +\n",
    "           theme(axis.title=element_text(size=20),\n",
    "                 title=element_text(size=20), \n",
    "                 text=element_text(size=35)) +\n",
    "           theme(legend.position=\"none\") + \n",
    "           theme(axis.text.y=element_blank()) +\n",
    "           ylab(\"\") +\n",
    "           xlab(\"\") +\n",
    "           xlim(-4.5, 4.5)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "broke-behavior",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_g <- df_total.g[,c(1:4,7)]\n",
    "df_g$pval <- ifelse(df_g$pval < 0.05,\n",
    "                    \"sig\",\n",
    "                    'no_sig'\n",
    ")\n",
    "df_g$log2fc <- ifelse(df_g$log2fc < 0,\n",
    "                      \"down\",\n",
    "                      'up'\n",
    ")\n",
    "df_g = df_g[order(df_g$pathway),]\n",
    "write.csv(df_g, 'marker_gfap_chow_GO.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "noted-kazakhstan",
   "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": [
    "pdf(file=\"enrichdot_both_marker.pdf\", width=7, height=15)   \n",
    "    df_total.b=data.frame()\n",
    "    for(i in 1:nrow(all)){   \n",
    "        p = as.data.frame(b0[which(b0$gene%in%\n",
    "                          sg.b0.all[sg.b0.all[,2]%in%unlist(strsplit(all$geneID[i],\"/\")), 1]), c(2, 3, 5)])\n",
    "        p$pathway = all$Description[i]\n",
    "        p$padjPath = all$p.adjust[i]\n",
    "        p$count = all$Count[i]\n",
    "        p$Comparison = all$Comparison[i]\n",
    "        df_total.b = rbind(df_total.b, p)\n",
    "    }\n",
    "    df_total.b$significance <- ifelse(df_total.b$pval < 0.05,\n",
    "                                      \"significant\",\n",
    "                                      \"not significant\"\n",
    "    )\n",
    "    ggplot(df_total.b, aes(log2fc, y=factor(pathway, level=go[length(go):1]))) + \n",
    "           geom_point(aes(colour=significance, size=-log10(pval)), alpha=0.5) + \n",
    "           scale_color_manual(values=c(\"gray85\", \"red\")) + \n",
    "           scale_size_continuous(range=c(0, 8)) +\n",
    "           geom_vline(xintercept=0) +\n",
    "           theme(axis.title=element_text(size=20),\n",
    "                 title=element_text(size=20), \n",
    "                 text=element_text(size=35)) +\n",
    "           theme(legend.position=\"none\") + \n",
    "           theme(axis.text.y=element_blank()) +\n",
    "           ylab(\"\") +\n",
    "           xlab(\"\") + \n",
    "           xlim(-4.5, 4.5)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "stopped-disorder",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_b = df_total.b[,c(1:4,7)]\n",
    "df_b$pval <- ifelse(df_b$pval < 0.05, \n",
    "                    \"sig\",\n",
    "                    'no_sig'\n",
    ")\n",
    "df_b$log2fc <- ifelse(df_b$log2fc < 0,\n",
    "                      \"down\",\n",
    "                      'up'\n",
    ")\n",
    "df_b = df_b[order(df_b$pathway),]\n",
    "write.csv(df_b,'marker_both_chow_GO.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "moved-correlation",
   "metadata": {},
   "source": [
    "## Diet effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "reverse-transport",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=8, repr.plot.height=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "large-tablet",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 2 rows containing missing values (geom_vline).”\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=\"volcano_aldh_diet_chow-hfd5.pdf\", width=10, height=10)    \n",
    "    a0_5 <- read.csv('~/Documents/consultation/Luiza/single_cell/results/tables/diet_aldh_only_chow-hfd5.csv', \n",
    "                     header=TRUE, \n",
    "                     sep=\"\\t\"\n",
    "    )\n",
    "    a0_5 = a0_5[complete.cases(a0_5), ]\n",
    "    a0_5 = a0_5[order(a0_5$pval),]\n",
    "    a0_5s = a0_5[a0_5$pval < 0.05,]\n",
    "    a0_5$log2fc <- ifelse(a0_5$log2fc > 4,\n",
    "                          4, \n",
    "                          a0_5$log2fc\n",
    "    )\n",
    "    a0_5$log2fc <- ifelse(a0_5$log2fc < -4,\n",
    "                          -4,\n",
    "                          a0_5$log2fc\n",
    "    )\n",
    "    EnhancedVolcano(a0_5,\n",
    "                    lab=a0_5$gene,\n",
    "                    pCutoff=0.05,\n",
    "                    FCcutoff=10,\n",
    "                    axisLabSize=30,\n",
    "                    title=NULL,\n",
    "                    subtitle=NULL,\n",
    "                    caption=NULL,\n",
    "                    x='log2fc',\n",
    "                    y='pval',\n",
    "                    ylim=c(0, 9),\n",
    "                    xlim=c(-4, 4),\n",
    "                    selectLab=a0_5$gene[1:25],\n",
    "                    labSize=8,\n",
    "                    col=c(\"grey30\", \"grey30\", \"red2\", \"red2\"),\n",
    "                    legendLabels=c('NS', expression(Log[2]~FC), 'p-value<0.05', expression(p-value~and~log[2]~FC)),\n",
    "                    legendLabSize=30\n",
    "    )\n",
    "dev.off()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "inside-northern",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 2 rows containing missing values (geom_vline).”\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=\"volcano_aldh_diet_chow-hfd15.pdf\", width=10, height=10)       \n",
    "    a0_15 <- read.csv('~/Documents/consultation/Luiza/single_cell/results/tables/diet_aldh_only_chow-hfd15.csv',\n",
    "                      header=TRUE, \n",
    "                      sep=\"\\t\"\n",
    "    )\n",
    "    a0_15 = a0_15[complete.cases(a0_15), ]\n",
    "    a0_15 = a0_15[order(a0_15$pval),]\n",
    "    a0_15s = a0_15[a0_15$pval < 0.05,]\n",
    "    a0_15$log2fc <- ifelse(a0_15$log2fc > 4,\n",
    "                           4,\n",
    "                           a0_15$log2fc\n",
    "    )\n",
    "    a0_15$log2fc <- ifelse(a0_15$log2fc < -4,\n",
    "                           -4,\n",
    "                           a0_15$log2fc\n",
    "    )\n",
    "    EnhancedVolcano(a0_15,\n",
    "                    lab=a0_15$gene,\n",
    "                    pCutoff=0.05,\n",
    "                    FCcutoff=10,\n",
    "                    axisLabSize=30,\n",
    "                    title=NULL,\n",
    "                    subtitle=NULL,\n",
    "                    caption=NULL,\n",
    "                    x='log2fc',\n",
    "                    y='pval',\n",
    "                    #ylim=c(0, 8),\n",
    "                    xlim=c(-4, 4),\n",
    "                    selectLab=a0_15$gene[1:25],\n",
    "                    labSize=8,\n",
    "                    col=c(\"grey30\", \"grey30\", \"red2\", \"red2\"),\n",
    "                    legendLabels=c('NS', expression(Log[2]~FC), 'p-value<0.05', expression(p-value~and~log[2]~FC)),\n",
    "                    legendLabSize=30\n",
    "    )\n",
    "dev.off()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "parental-database",
   "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": [
    "pdf(file=\"venn_aldh_diet.pdf\", width=8, height=8)   \n",
    "    vp <- venn.diagram(list(' '=a0_5s$gene,\n",
    "                            ' '=a0_15s$gene),\n",
    "                       fill=c(\"limegreen\", \"green1\"),\n",
    "                       filename=NULL, \n",
    "                       cex=3.5,\n",
    "                       col=\"transparent\",\n",
    "                       cat.cex=0.01,\n",
    "                       margin=0.01,\n",
    "                       main=NULL,\n",
    "                       main.cex=3);\n",
    "    grid.draw(vp)   \n",
    "dev.off()   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "surgical-pointer",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 2 rows containing missing values (geom_vline).”\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=\"volcano_gfap_diet_chow-hfd5.pdf\", width=10, height=10) \n",
    "    g0_5 <- read.csv('~/Documents/consultation/Luiza/single_cell/results/tables/diet_gfap_only_chow-hfd5.csv', \n",
    "                     header=TRUE, \n",
    "                     sep=\"\\t\"\n",
    "    )\n",
    "    g0_5 = g0_5[complete.cases(g0_5), ]\n",
    "    g0_5 = g0_5[order(g0_5$pval),]\n",
    "    g0_5s = g0_5[g0_5$pval < 0.05,]\n",
    "    g0_5$log2fc <- ifelse(g0_5$log2fc > 4,\n",
    "                          4,\n",
    "                          g0_5$log2fc\n",
    "    )\n",
    "    g0_5$log2fc <- ifelse(g0_5$log2fc < -4,\n",
    "                          -4,\n",
    "                          g0_5$log2fc\n",
    "    )\n",
    "    EnhancedVolcano(g0_5,\n",
    "                    lab=g0_5$gene,\n",
    "                    pCutoff=0.05,\n",
    "                    FCcutoff=10,\n",
    "                    axisLabSize=30,\n",
    "                    title=NULL,\n",
    "                    subtitle=NULL,\n",
    "                    caption=NULL,\n",
    "                    x='log2fc',\n",
    "                    y='pval',\n",
    "                    ylim=c(0, 9),\n",
    "                    xlim=c(-4, 4),\n",
    "                    selectLab=g0_5$gene[1:25],\n",
    "                    labSize=8,\n",
    "                    col = c(\"grey30\", \"grey30\", \"red2\", \"red2\"),\n",
    "                    legendLabels=c('NS', expression(Log[2]~FC), 'p-value<0.05', expression(p-value~and~log[2]~FC)),\n",
    "                    legendLabSize=30\n",
    "    )\n",
    "dev.off()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "fewer-typing",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 2 rows containing missing values (geom_vline).”\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=\"volcano_gfap_diet_chow-hfd15.pdf\", width=10, height=10)  \n",
    "    g0_15 <- read.csv('~/Documents/consultation/Luiza/single_cell/results/tables/diet_gfap_only_chow-hfd15.csv',\n",
    "                      header=TRUE, \n",
    "                      sep=\"\\t\"\n",
    "    )\n",
    "    g0_15 = g0_15[complete.cases(g0_15), ]\n",
    "    g0_15 = g0_15[order(g0_15$pval),]\n",
    "    g0_15s = g0_15[g0_15$pval < 0.05,]\n",
    "    g0_15$log2fc <- ifelse(g0_15$log2fc > 4,\n",
    "                           4,\n",
    "                           g0_15$log2fc\n",
    "    )\n",
    "    g0_15$log2fc <- ifelse(g0_15$log2fc < -4,\n",
    "                           -4,\n",
    "                           g0_15$log2fc\n",
    "    )\n",
    "    EnhancedVolcano(g0_15,\n",
    "                    lab=g0_15$gene,\n",
    "                    pCutoff=0.05,\n",
    "                    FCcutoff=10,\n",
    "                    axisLabSize=30,\n",
    "                    title=NULL,\n",
    "                    subtitle=NULL,\n",
    "                    caption=NULL,\n",
    "                    x='log2fc',\n",
    "                    y='pval',\n",
    "                    #ylim=c(0, 8),\n",
    "                    xlim=c(-4, 4),\n",
    "                    selectLab=g0_15$gene[1:25],\n",
    "                    labSize=8,\n",
    "                    col = c(\"grey30\", \"grey30\", \"red2\", \"red2\"),\n",
    "                    legendLabels=c('NS', expression(Log[2]~FC), 'p-value<0.05', expression(p-value~and~log[2]~FC)),\n",
    "                    legendLabSize=30\n",
    "    )\n",
    "dev.off()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "loaded-archives",
   "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": [
    "pdf(file=\"venn_gfap_diet.pdf\", width=8, height=8) \n",
    "    vp <- venn.diagram(list(' '=g0_5s$gene,\n",
    "                            ' '=g0_15s$gene),\n",
    "                       fill=c(\"magenta\", \"violet\"),\n",
    "                       filename=NULL,\n",
    "                       cex=3.5,\n",
    "                       col=\"transparent\",\n",
    "                       cat.cex=0.01, \n",
    "                       margin=0.01, \n",
    "                       main=NULL, \n",
    "                       main.cex=3);\n",
    "    grid.draw(vp)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "floral-renaissance",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 2 rows containing missing values (geom_vline).”\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=\"volcano_both_diet_chow-hfd5.pdf\", width=10, height=10)     \n",
    "    b0_5 <- read.csv('~/Documents/consultation/Luiza/single_cell/results/tables/diet_both_chow-hfd5.csv',\n",
    "                     header=TRUE,\n",
    "                     sep=\"\\t\"\n",
    "    )\n",
    "    b0_5 = b0_5[complete.cases(b0_5), ]\n",
    "    b0_5 = b0_5[order(b0_5$pval),]\n",
    "    b0_5s = b0_5[b0_5$pval < 0.05,]\n",
    "    b0_5$log2fc <- ifelse(b0_5$log2fc > 4, \n",
    "                          4,\n",
    "                          b0_5$log2fc\n",
    "    )\n",
    "    b0_5$log2fc <- ifelse(b0_5$log2fc < -4,\n",
    "                          -4,\n",
    "                          b0_5$log2fc\n",
    "    )\n",
    "    EnhancedVolcano(b0_5,\n",
    "                    lab=b0_5$gene,\n",
    "                    pCutoff=0.05,\n",
    "                    FCcutoff=10,\n",
    "                    axisLabSize=30,\n",
    "                    title=NULL,\n",
    "                    subtitle=NULL,\n",
    "                    caption=NULL,\n",
    "                    x='log2fc',\n",
    "                    y='pval',\n",
    "                    ylim=c(0, 9),\n",
    "                    xlim=c(-4, 4),\n",
    "                    selectLab=b0_5$gene[1:25],\n",
    "                    labSize=8,\n",
    "                    col=c(\"grey30\", \"grey30\", \"red2\", \"red2\"),\n",
    "                    legendLabels=c('NS', expression(Log[2]~FC), 'p-value<0.05', expression(p-value~and~log[2]~FC)),\n",
    "                    legendLabSize=30)\n",
    "dev.off()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "damaged-distinction",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 2 rows containing missing values (geom_vline).”\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 = \"volcano_both_diet_chow-hfd15.pdf\", width = 10, height = 10)    \n",
    "    b0_15 <- read.csv('~/Documents/consultation/Luiza/single_cell/results/tables/diet_both_chow-hfd15.csv',\n",
    "                      header=TRUE, \n",
    "                      sep=\"\\t\"\n",
    "    )\n",
    "    b0_15 = b0_15[complete.cases(b0_15), ]\n",
    "    b0_15 = b0_15[order(b0_15$pval),]\n",
    "    b0_15s = b0_15[b0_15$pval<0.05,]\n",
    "    b0_15$log2fc <- ifelse(b0_15$log2fc>4, \n",
    "                           4,\n",
    "                           b0_15$log2fc\n",
    "    )\n",
    "    b0_15$log2fc <- ifelse(b0_15$log2fc< -4,\n",
    "                           -4,\n",
    "                           b0_15$log2fc\n",
    "    )\n",
    "    EnhancedVolcano(b0_15,\n",
    "                    lab=b0_15$gene,\n",
    "                    pCutoff=0.05,\n",
    "                    FCcutoff=10,\n",
    "                    axisLabSize=30,\n",
    "                    title=NULL,\n",
    "                    subtitle=NULL,\n",
    "                    caption=NULL,\n",
    "                    x='log2fc',\n",
    "                    y='pval',\n",
    "                    #ylim=c(0, 8),\n",
    "                    xlim=c(-4, 4),\n",
    "                    selectLab=b0_15$gene[1:25],\n",
    "                    labSize=8,\n",
    "                    col=c(\"grey30\", \"grey30\", \"red2\", \"red2\"),\n",
    "                    legendLabels=c('NS', expression(Log[2]~FC), 'p-value<0.05', expression(p-value~and~log[2]~FC)),\n",
    "                    legendLabSize=30\n",
    "    )\n",
    "dev.off()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "impressive-currency",
   "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": [
    "pdf(file=\"venn_both_diet.pdf\", width=8, height=8) \n",
    "    vp <- venn.diagram(list(' '=b0_5s$gene,\n",
    "                            ' '=b0_15s$gene),\n",
    "                       fill=c(\"dodgerblue\", \"turquoise3\"),\n",
    "                       filename=NULL, \n",
    "                       cex=3.5, \n",
    "                       col=\"transparent\",\n",
    "                       cat.cex=0.01, \n",
    "                       margin=0.01, \n",
    "                       main=NULL, \n",
    "                       main.cex=3);\n",
    "    grid.draw(vp)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "bibliographic-commercial",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(setdiff(a0_5s$gene, union(g0_5s$gene, b0_5s$gene)), fromType = \"SYMBOL\", :\n",
      "“2.06% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(setdiff(g0_5s$gene, union(a0_5s$gene, b0_5s$gene)), fromType = \"SYMBOL\", :\n",
      "“1.14% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(setdiff(b0_5s$gene, union(a0_5s$gene, g0_5s$gene)), fromType = \"SYMBOL\", :\n",
      "“2.04% of input gene IDs are fail to map...”\n"
     ]
    }
   ],
   "source": [
    "sga0_5 <- bitr(setdiff(a0_5s$gene, union(g0_5s$gene, b0_5s$gene)), \n",
    "               fromType=\"SYMBOL\",\n",
    "               toType=\"ENTREZID\",\n",
    "               OrgDb=org.Mm.eg.db\n",
    ")\n",
    "a0_5g <- enrichGO(sga0_5[,2],\n",
    "                  'org.Mm.eg.db',\n",
    "                  ont=\"BP\", \n",
    "                  pvalueCutoff=0.1\n",
    ")\n",
    "sgg0_5 <- bitr(setdiff(g0_5s$gene, union(a0_5s$gene, b0_5s$gene)),\n",
    "               fromType=\"SYMBOL\", \n",
    "               toType=\"ENTREZID\", \n",
    "               OrgDb=org.Mm.eg.db\n",
    ")\n",
    "g0_5g <- enrichGO(sgg0_5[,2],\n",
    "                  'org.Mm.eg.db',\n",
    "                  ont=\"BP\", \n",
    "                  pvalueCutoff=0.1\n",
    ")\n",
    "sgb0_5 <- bitr(setdiff(b0_5s$gene, union(a0_5s$gene, g0_5s$gene)),\n",
    "               fromType=\"SYMBOL\",\n",
    "               toType=\"ENTREZID\", \n",
    "               OrgDb=org.Mm.eg.db\n",
    ")\n",
    "b0_5g <- enrichGO(sgb0_5[,2],\n",
    "                  'org.Mm.eg.db',\n",
    "                  ont=\"BP\", \n",
    "                  pvalueCutoff=0.1\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "reported-bubble",
   "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": [
    "pdf(file=\"venn_all_diet_hfd5.pdf\", width=8, height=8)\n",
    "    vp <- venn.diagram(list(' '=a0_5s$gene,\n",
    "                            ' '=g0_5s$gene,\n",
    "                            ' '=b0_5s$gene),\n",
    "                       fill=c(\"limegreen\", \"magenta\", \"dodgerblue\"),\n",
    "                       filename=NULL, \n",
    "                       cex=3.5, \n",
    "                       col=\"transparent\",\n",
    "                       cat.cex = 0.01,\n",
    "                       margin = 0.01,\n",
    "                       main=NULL,\n",
    "                       main.cex=3);\n",
    "    grid.draw(vp)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "wound-azerbaijan",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(setdiff(a0_15s$gene, union(g0_15s$gene, b0_15s$gene)), fromType = \"SYMBOL\", :\n",
      "“1.98% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(setdiff(g0_15s$gene, union(a0_15s$gene, b0_15s$gene)), fromType = \"SYMBOL\", :\n",
      "“1.48% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(setdiff(b0_15s$gene, union(a0_15s$gene, g0_15s$gene)), fromType = \"SYMBOL\", :\n",
      "“1.92% of input gene IDs are fail to map...”\n"
     ]
    }
   ],
   "source": [
    "sga0_15 <- bitr(setdiff(a0_15s$gene,union(g0_15s$gene,b0_15s$gene)), \n",
    "                fromType=\"SYMBOL\", \n",
    "                toType=\"ENTREZID\",\n",
    "                OrgDb=org.Mm.eg.db\n",
    ")\n",
    "a0_15g <- enrichGO(sga0_15[,2], \n",
    "                   'org.Mm.eg.db',\n",
    "                   ont=\"BP\", \n",
    "                   pvalueCutoff=0.1\n",
    ")\n",
    "sgg0_15 <- bitr(setdiff(g0_15s$gene, union(a0_15s$gene, b0_15s$gene)), \n",
    "                fromType=\"SYMBOL\", \n",
    "                toType=\"ENTREZID\",\n",
    "                OrgDb=org.Mm.eg.db\n",
    ")\n",
    "g0_15g <- enrichGO(sgg0_15[,2], \n",
    "                   'org.Mm.eg.db', \n",
    "                   ont=\"BP\",\n",
    "                   pvalueCutoff=0.1\n",
    ")\n",
    "sgb0_15 <- bitr(setdiff(b0_15s$gene, union(a0_15s$gene, g0_15s$gene)), \n",
    "                fromType=\"SYMBOL\",\n",
    "                toType=\"ENTREZID\", \n",
    "                OrgDb=org.Mm.eg.db\n",
    ")\n",
    "b0_15g <- enrichGO(sgb0_15[,2], \n",
    "                   'org.Mm.eg.db',\n",
    "                   ont=\"BP\",\n",
    "                   pvalueCutoff=0.1\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "satisfied-front",
   "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": [
    "pdf(file=\"venn_all_diet_hfd15.pdf\", width=8, height=8)\n",
    "    vp <- venn.diagram(list(aldh_only=a0_15s$gene,\n",
    "                            gfap_only=g0_15s$gene,\n",
    "                            both=b0_15s$gene),\n",
    "                     fill=c(\"green1\", \"violet\", \"turquoise3\"),\n",
    "                     filename=NULL,\n",
    "                     cex=3.5,\n",
    "                     col=\"transparent\",\n",
    "                     cat.cex=0.01,\n",
    "                     margin=0.01,\n",
    "                     main=NULL,\n",
    "                     main.cex=3);\n",
    "    grid.draw(vp)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "patent-eagle",
   "metadata": {},
   "outputs": [],
   "source": [
    "go_diet = c(\"cellular response to peptide hormone stimulus\",\n",
    "            \"insulin receptor signaling pathway\",\n",
    "            \"regulation of TOR signaling\",\n",
    "            \"glucose metabolic process\",\n",
    "            \"glycogen metabolic process\",\n",
    "            \"carbohydrate catabolic process\",\"pentose-phosphate shunt\",\n",
    "            \"regulation of glucose import\",\n",
    "            \"fatty acid beta-oxidation\",\n",
    "            \"cellular ketone metabolic process\",\n",
    "            \"fatty acid metabolic process\",\n",
    "            \"lipid catabolic process\",\n",
    "            \"response to oxidative stress\",\n",
    "            \"reactive oxygen species metabolic process\",\n",
    "            \"regulation of stress-activated protein kinase signaling cascade\",\n",
    "            \"regulation of response to endoplasmic reticulum stress\",\n",
    "            \"mitochondrial ATP synthesis coupled electron transport\",\n",
    "            \"oxidative phosphorylation\",\n",
    "            \"mitochondrial respiratory chain complex assembly\",\n",
    "            \"NADH dehydrogenase complex assembly\",\n",
    "            \"smoothened signaling pathway\",\n",
    "            \"regulation of spindle assembly\",\n",
    "            \"mitotic sister chromatid segregation\",\n",
    "            \"cytokinesis\",\n",
    "            \"mitotic nuclear division\",\n",
    "            \"mitotic cell cycle phase transition\",\n",
    "            \"interleukin-17 production\",\n",
    "            \"transforming growth factor beta receptor signaling pathway\",\n",
    "            \"cell junction assembly\",\n",
    "            \"cell-substrate adhesion\",\n",
    "            \"adherens junction organization\",\n",
    "            \"COPII-coated vesicle budding\",\n",
    "            \"regulation of calcium ion transmembrane transporter activity\",\n",
    "            \"exocytosis\",\"phosphatidylinositol metabolic process\",\n",
    "            \"regulation of neurotransmitter levels\",\n",
    "            \"synapse organization\",\n",
    "            \"regulation of synapse structure or activity\",\n",
    "            \"axonogenesis\",\n",
    "            \"regulation of synaptic transmission, dopaminergic\",\n",
    "            \"gamma-aminobutyric acid secretion\",\n",
    "            \"synaptic transmission, glutamatergic\",\n",
    "            \"response to hypoxia\",\n",
    "            \"regulation of vasculature development\",\n",
    "            \"positive regulation of blood vessel endothelial cell migration\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "id": "developmental-cardiff",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=15, repr.plot.height=15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "id": "cheap-stock",
   "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": [
    "pdf(file=\"enrich_path_diet0_5.pdf\", width=15, height=15)   \n",
    "    cutoff = 0.01\n",
    "    enr_a0_5g = a0_5g@result[a0_5g@result[,5] < cutoff, c(2, 5, 6, 8, 9)]\n",
    "    enr_a0_5g$Comparison = 'aldh_chow-hfd5'\n",
    "    #enr_a0_15g = a0_15g@result[a0_15g@result[,5] < cutoff, c(2, 5, 6, 8, 9)]\n",
    "    #enr_a0_15g$Comparison = 'aldh_chow-hfd15'\n",
    "    \n",
    "    enr_g0_5g = g0_5g@result[g0_5g@result[,5] < cutoff, c(2, 5, 6, 8, 9)]\n",
    "    enr_g0_5g$Comparison = 'gfap_chow-hfd5'\n",
    "    #enr_g0_15g = g0_15g@result[g0_15g@result[,5] < cutoff, c(2, 5, 6, 8, 9)]\n",
    "    #enr_g0_15g$Comparison = 'gfap_chow-hfd15'\n",
    "    \n",
    "    enr_b0_5g = b0_5g@result[b0_5g@result[,5] < cutoff, c(2, 5, 6, 8, 9)]\n",
    "    enr_b0_5g$Comparison = 'pair_chow-hfd5'\n",
    "    #enr_b0_15g = b0_15g@result[b0_15g@result[,5] < cutoff, c(2, 5, 6, 8, 9)]\n",
    "    #enr_b0_15g$Comparison = 'pair_chow-hfd15'\n",
    "    \n",
    "    all = rbind(enr_a0_5g, enr_g0_5g, enr_b0_5g)# enr_a0_15g,enr_g0_15g,enr_b0_15g\n",
    "    all = all[all$Description%in%go_diet,]\n",
    "    #write.csv(all, file=\"GO_sc_diet_effect.csv\")\n",
    "    \n",
    "    colnames(all) = c('Description', 'pvalue', 'p.adjust', 'geneID', 'count', 'Comparison')\n",
    "    ggplot(all, \n",
    "           aes(Comparison, y=factor(Description, level=go_diet[length(go_diet):1]), col=pvalue, size=count)) + \n",
    "           geom_point()+ \n",
    "           scale_colour_gradient(low=\"red\", high=\"blue\") +\n",
    "           scale_size_continuous(range=c(3, 10)) +\n",
    "           xlab(\"\") + \n",
    "           ylab(\"\") + \n",
    "           #theme(axis.text.x=element_blank()) +\n",
    "           #theme(legend.position=\"none\") +\n",
    "           theme(text=element_text(size=15))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "welcome-catalyst",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(a0_5$gene, fromType = \"SYMBOL\", toType = \"ENTREZID\", OrgDb = org.Mm.eg.db):\n",
      "“3.36% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(a0_15$gene, fromType = \"SYMBOL\", toType = \"ENTREZID\", OrgDb = org.Mm.eg.db):\n",
      "“3.54% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(g0_5$gene, fromType = \"SYMBOL\", toType = \"ENTREZID\", OrgDb = org.Mm.eg.db):\n",
      "“3.17% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(g0_15$gene, fromType = \"SYMBOL\", toType = \"ENTREZID\", OrgDb = org.Mm.eg.db):\n",
      "“3.22% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(b0_5$gene, fromType = \"SYMBOL\", toType = \"ENTREZID\", OrgDb = org.Mm.eg.db):\n",
      "“3.42% of input gene IDs are fail to map...”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr(b0_15$gene, fromType = \"SYMBOL\", toType = \"ENTREZID\", OrgDb = org.Mm.eg.db):\n",
      "“3.22% of input gene IDs are fail to map...”\n"
     ]
    }
   ],
   "source": [
    "sg.a0_5 <- bitr(a0_5$gene,\n",
    "                fromType = \"SYMBOL\",\n",
    "                toType=\"ENTREZID\",\n",
    "                OrgDb=org.Mm.eg.db\n",
    ")\n",
    "sg.a0_15 <- bitr(a0_15$gene,\n",
    "                 fromType = \"SYMBOL\",\n",
    "                 toType=\"ENTREZID\",\n",
    "                 OrgDb=org.Mm.eg.db\n",
    ")\n",
    "sg.g0_5 <- bitr(g0_5$gene,\n",
    "                fromType=\"SYMBOL\",\n",
    "                toType=\"ENTREZID\",\n",
    "                OrgDb=org.Mm.eg.db)\n",
    "sg.g0_15 <- bitr(g0_15$gene,\n",
    "                 fromType=\"SYMBOL\",\n",
    "                 toType=\"ENTREZID\", \n",
    "                 OrgDb=org.Mm.eg.db\n",
    ")\n",
    "sg.b0_5 <- bitr(b0_5$gene,\n",
    "                fromType=\"SYMBOL\", \n",
    "                toType=\"ENTREZID\",\n",
    "                OrgDb=org.Mm.eg.db\n",
    ")\n",
    "sg.b0_15 <- bitr(b0_15$gene,\n",
    "                 fromType=\"SYMBOL\",\n",
    "                 toType=\"ENTREZID\",\n",
    "                 OrgDb=org.Mm.eg.db\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "id": "marine-classification",
   "metadata": {},
   "outputs": [],
   "source": [
    "all = rbind(enr_a0_5g, enr_g0_5g, enr_b0_5g)#,enr_a0_15g,enr_g0_15g,enr_b0_15g) #  \n",
    "all = all[all$Description%in%go_diet,]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "consistent-breeding",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=7, repr.plot.height=15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "minus-surgeon",
   "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": [
    "pdf(file=\"enrichdot_aldh0_5_diet.pdf\", width=7, height=15)   \n",
    "    df_total.a0_5 = data.frame()\n",
    "    for(i in 1:nrow(all)){   \n",
    "        p = as.data.frame(a0_5[which(a0_5$gene%in%\n",
    "                          sg.a0_5[sg.a0_5[,2]%in%unlist(strsplit(all$geneID[i],\"/\")),1]), c(2,3,5)])\n",
    "        p$pathway = all$Description[i]\n",
    "        p$padjPath = all$p.adjust[i]\n",
    "        p$count = all$Count[i]\n",
    "        df_total.a0_5 = rbind(df_total.a0_5, p)\n",
    "    }\n",
    "    df_total.a0_5$significance <- ifelse(df_total.a0_5$pval <= 0.05,\n",
    "                                         \"significant\",\n",
    "                                         \"not significant\"\n",
    "    )\n",
    "    ggplot(df_total.a0_5, aes(log2fc, y=factor(pathway, level=go_diet[length(go_diet):1]))) + \n",
    "           geom_point(aes(colour=significance, size=-log10(pval)), alpha=0.5) + \n",
    "           scale_color_manual(values=c(\"gray85\", \"red\"))+ \n",
    "           scale_size_continuous(range=c(0, 10)) +\n",
    "           geom_vline(xintercept=0) +\n",
    "           theme(axis.title=element_text(size=20),\n",
    "                 title=element_text(size=20), \n",
    "                 text=element_text(size=35)) +\n",
    "           theme(legend.position=\"none\") + \n",
    "           theme(axis.text.y=element_blank()) +\n",
    "           ylab(\"\") +\n",
    "           xlab(\"\") + \n",
    "           xlim(-4.5, 4.5)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "phantom-interest",
   "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": [
    "pdf(file=\"enrichdot_aldh0_15_diet.pdf\", width=4, height=15)   \n",
    "    df_total.a0_15 = data.frame()\n",
    "    for(i in 1:nrow(all)){   \n",
    "        p = as.data.frame(a0_15[which(a0_15$gene%in%\n",
    "                          sg.a0_15[sg.a0_15[,2]%in%unlist(strsplit(all$geneID[i],\"/\")),1]), c(3,5)])\n",
    "        p$pathway = all$Description[i]\n",
    "        p$padjPath = all$p.adjust[i]\n",
    "        p$count = all$Count[i]\n",
    "        df_total.a0_15 = rbind(df_total.a0_15, p)\n",
    "    }\n",
    "    df_total.a0_15$significance <- ifelse(df_total.a0_15$pval <= 0.05,\n",
    "                                          \"significant\",\n",
    "                                          \"not significant\"\n",
    "    )\n",
    "    ggplot(df_total.a0_15, aes(log2fc, y=factor(pathway, level=go_diet[length(go_diet):1]))) + \n",
    "           geom_point(aes(colour=significance, size=-log10(pval)), alpha=0.5) + \n",
    "           scale_color_manual(values=c(\"gray85\", \"red\")) + \n",
    "           scale_size_continuous(range=c(0, 8)) +\n",
    "           geom_vline(xintercept=0) +\n",
    "           theme(axis.title=element_text(size=20),\n",
    "                 title=element_text(size=20), \n",
    "                 text=element_text(size=35)) +\n",
    "           theme(legend.position=\"none\") + \n",
    "           theme(axis.text.y=element_blank()) +\n",
    "           ylab(\"\") + \n",
    "           xlab(\"\") +\n",
    "           xlim(-4.5, 4.5)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "forced-tribe",
   "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": [
    "pdf(file=\"enrichdot_gfap0_5_diet.pdf\", width=7, height=15)   \n",
    "    df_total.g0_5 = data.frame()\n",
    "    for(i in 1:nrow(all)){   \n",
    "        p = as.data.frame(g0_5[which(g0_5$gene%in%\n",
    "                          sg.g0_5[sg.g0_5[,2]%in%unlist(strsplit(all$geneID[i],\"/\")),1]), c(3,5)])\n",
    "        p$pathway = all$Description[i]\n",
    "        p$padjPath = all$p.adjust[i]\n",
    "        p$count = all$Count[i]\n",
    "        df_total.g0_5 = rbind(df_total.g0_5, p)\n",
    "    }\n",
    "    df_total.g0_5$significance <- ifelse(df_total.g0_5$pval <= 0.05,\n",
    "                                         \"significant\",\n",
    "                                         \"not significant\"\n",
    "    )\n",
    "    ggplot(df_total.g0_5, aes(log2fc, y=factor(pathway, level=go_diet[length(go_diet):1]))) + \n",
    "           geom_point(aes(colour=significance, size=-log10(pval)), alpha=0.5) + \n",
    "           scale_color_manual(values=c(\"gray85\", \"red\")) + \n",
    "           scale_size_continuous(range=c(0, 7)) +\n",
    "           geom_vline(xintercept=0) +\n",
    "           theme(axis.title=element_text(size=20),\n",
    "                 title=element_text(size=20), \n",
    "                 text=element_text(size=35)) +\n",
    "           theme(legend.position=\"none\") + \n",
    "           theme(axis.text.y=element_blank()) +\n",
    "           ylab(\"\") +\n",
    "           xlab(\"\") +\n",
    "           xlim(-4.5, 4.5)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "special-serve",
   "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": [
    "pdf(file=\"enrichdot_gfap0_15_diet.pdf\", width=4, height=15)   \n",
    "    df_total.g0_15 = data.frame()\n",
    "    for(i in 1:nrow(all)){   \n",
    "        p = as.data.frame(g0_15[which(g0_15$gene%in%\n",
    "                          sg.g0_15[sg.g0_15[,2]%in%unlist(strsplit(all$geneID[i],\"/\")), 1]), c(3,5)])\n",
    "        p$pathway = all$Description[i]\n",
    "        p$padjPath = all$p.adjust[i]\n",
    "        p$count = all$Count[i]\n",
    "        df_total.g0_15 = rbind(df_total.g0_15, p)\n",
    "    }\n",
    "    df_total.g0_15$significance <- ifelse(df_total.g0_15$pval <= 0.05, \n",
    "                                          \"significant\", \n",
    "                                          \"not significant\"\n",
    "    )\n",
    "    ggplot(df_total.g0_15, aes(log2fc, y=factor(pathway, level=go_diet[length(go_diet):1]))) + \n",
    "           geom_point(aes(colour=significance, size=-log10(pval)), alpha=0.5) + \n",
    "           scale_color_manual(values=c(\"gray85\", \"red\")) + \n",
    "           scale_size_continuous(range=c(0, 8)) +\n",
    "           geom_vline(xintercept=0) +\n",
    "           theme(axis.title=element_text(size=20),\n",
    "                 title=element_text(size=20), \n",
    "                 text=element_text(size=35)) +\n",
    "           #theme(legend.position=\"none\") + \n",
    "           theme(axis.text.y=element_blank()) +\n",
    "           ylab(\"\") +\n",
    "           xlab(\"\") +\n",
    "           xlim(-4.5, 4.5)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "id": "prerequisite-honey",
   "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": [
    "pdf(file=\"enrichdot_both0_5_diet.pdf\", width=7, height=15)   \n",
    "    df_total.b0_5 = data.frame()\n",
    "    for(i in 1:nrow(all)){   \n",
    "        p = as.data.frame(b0_5[which(b0_5$gene%in%\n",
    "                          sg.b0_5[sg.b0_5[,2]%in%unlist(strsplit(all$geneID[i],\"/\")), 1]), c(3,5)])\n",
    "        p$pathway = all$Description[i]\n",
    "        p$padjPath = all$p.adjust[i]\n",
    "        p$count = all$Count[i]\n",
    "        df_total.b0_5 = rbind(df_total.b0_5, p)\n",
    "    }\n",
    "    df_total.b0_5$significance <- ifelse(df_total.b0_5$pval <= 0.05, \n",
    "                                         \"significant\",\n",
    "                                         \"not significant\"\n",
    "    ) \n",
    "    ggplot(df_total.b0_5, aes(log2fc, y=factor(pathway, level=go_diet[length(go_diet):1]))) + \n",
    "           geom_point(aes(colour=significance, size=-log10(pval)), alpha=0.5) + \n",
    "           scale_color_manual(values=c(\"gray85\", \"red\"))+ \n",
    "           scale_size_continuous(range=c(0, 8)) +\n",
    "           geom_vline(xintercept=0) +\n",
    "           theme(axis.title=element_text(size=20),\n",
    "                 title=element_text(size=20), \n",
    "                 text=element_text(size=35)) +\n",
    "           theme(legend.position=\"none\") + \n",
    "           theme(axis.text.y=element_blank()) +\n",
    "           ylab(\"\") +\n",
    "           xlab(\"\") +\n",
    "           xlim(-4.5, 4.5)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "professional-steel",
   "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": [
    "pdf(file=\"enrichdot_both0_15_diet.pdf\", width=4, height=15)   \n",
    "    df_total.b0_15=data.frame()\n",
    "    for(i in 1:nrow(all)){   \n",
    "        p = as.data.frame(b0_15[which(b0_15$gene%in%\n",
    "                          sg.b0_15[sg.b0_15[,2]%in%unlist(strsplit(all$geneID[i],\"/\")), 1]), c(3,5)])\n",
    "        p$pathway = all$Description[i]\n",
    "        p$padjPath = all$p.adjust[i]\n",
    "        p$count = all$Count[i]\n",
    "        df_total.b0_15 = rbind(df_total.b0_15, p)\n",
    "    }\n",
    "    df_total.b0_15$significance <- ifelse(df_total.b0_15$pval <= 0.05,\n",
    "                                          \"significant\",\n",
    "                                          \"not significant\"\n",
    "    )\n",
    "    ggplot(df_total.b0_15, aes(log2fc, y=factor(pathway, level=go_diet[length(go_diet):1]))) + \n",
    "           geom_point(aes(colour=significance, size=-log10(pval)), alpha=0.5) + \n",
    "           scale_color_manual(values=c(\"gray85\", \"red\")) + \n",
    "           scale_size_continuous(range=c(0, 8)) +\n",
    "           geom_vline(xintercept=0) +\n",
    "           theme(axis.title=element_text(size=20),\n",
    "                 title=element_text(size=20), \n",
    "                 text=element_text(size=35)) +\n",
    "           theme(legend.position=\"none\") + \n",
    "           theme(axis.text.y=element_blank()) +\n",
    "           ylab(\"\") + \n",
    "           xlab(\"\") + \n",
    "           xlim(-4.5, 4.5)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "static-thanksgiving",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#if(bool_plot){  \n",
    "#    cutoff=0.05\n",
    "#    enr_a0_5k = a0_5k@result[a0_5k@result[,5] < cutoff, c(2,5,6,8,9)]\n",
    "#    enr_a0_5k$Comparison = 'aldh_only_chow-hfd5'\n",
    "#    enr_a0_15k = a0_15k@result[a0_15k@result[,5] < cutoff, c(2,5,6,8,9)]\n",
    "#    enr_a0_15k$Comparison = 'aldh_only_chow-hfd15'\n",
    "#    \n",
    "#    enr_g0_5k = g0_5k@result[g0_5k@result[,5] < cutoff, c(2,5,6,8,9)]\n",
    "#    enr_g0_5k$Comparison = 'gfap_only_chow-hfd5'\n",
    "#    enr_g0_15k = g0_15k@result[g0_15k@result[,5] < cutoff, c(2,5,6,8,9)]\n",
    "#    enr_g0_15k$Comparison = 'gfap_only_chow-hfd15'\n",
    "#    \n",
    "#    enr_b0_5k = b0_5k@result[b0_5k@result[,5] < cutoff, c(2,5,6,8,9)]\n",
    "#    enr_b0_5k$Comparison = 'both_chow-hfd5'\n",
    "#    enr_b0_15k = b0_15k@result[b0_15k@result[,5] < cutoff, c(2,5,6,8,9)]\n",
    "#    enr_b0_15k$Comparison = 'both_chow-hfd15'\n",
    "#    \n",
    "#    all = rbind(enr_a0_5k, enr_a0_15k, enr_g0_5k, enr_g0_15k, enr_b0_5k, enr_b0_15k)\n",
    "#    write.csv(all, file=\"KEGG_sc_diet_effect.csv\")\n",
    "#    options(repr.plot.width=15, repr.plot.height=30)\n",
    "#    ggplot(all, \n",
    "#           aes(Comparison, Description, col=pvalue, size=Count)) + \n",
    "#           geom_point()+ \n",
    "#           scale_colour_gradient(low=\"red\", high=\"blue\") +\n",
    "#           xlab(\"Tissue\") + \n",
    "#           ylab(\"KEGG enriched pathways\") + \n",
    "#           theme(text=element_text(size=25),axis.text.x=element_text(angle=45, hjust=1))\n",
    "#}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "white-registrar",
   "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": 5
}