transcriptomics_analysis_v1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of the astrocytes RNA-seq data from HPT, CTX, and Hippo  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of contents:\n",
    "\n",
    "* <a href=#Load>Load packages and set global variables</a>\n",
    "* <a href=#Dataloading>Loading data, preprocessing and make DEseq object</a>\n",
    "* <a href=#Visualisataion>Visualisation of the data</a>\n",
    "* <a href=#DGE>Differential gene expression - diet effect</a>\n",
    "* <a href=#Celltype>Heatmap marker genes</a>\n",
    "* <a href=#HeatmapDGE>Heatmap DGE</a>\n",
    "* <a href=#GOenrichment>GO pathway enrichment</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load packages and set global variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "inst <- suppressMessages(lapply(c('DESeq2',\n",
    "                                  'vidger',\n",
    "                                  'pheatmap',\n",
    "                                  'ggplot2',\n",
    "                                  'EnhancedVolcano',\n",
    "                                  'VennDiagram',\n",
    "                                  'org.Mm.eg.db',\n",
    "                                  'clusterProfiler',\n",
    "                                  'UpSetR',\n",
    "                                  'readxl',\n",
    "                                  'tidyverse',\n",
    "                                  'KEGGREST'),\n",
    "                                library,\n",
    "                                character.only=TRUE)\n",
    ") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "prot_input_file = '~/Documents/consultation/Luiza/omics_data/data/astrocyteprotoeme12102017-1.xlsx'\n",
    "options(repr.plot.width=8, repr.plot.height=8)\n",
    "my_palette = colorRampPalette(c(\"blue\", \"white\", \"red\"))(n=255)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set whether anndata objects are recomputed or loaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "bool_recomp = FALSE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set whether to produce plots, set to False for test runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "bool_plot = FALSE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading data, preprocessing and make DEseq object "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "dat <- read.csv(file=\"/Users/viktorian.miok/Documents/data/RNA-seq_Data-NEW/new_counts1.csv\", \n",
    "                header=TRUE,\n",
    "                sep=\",\"\n",
    ")\n",
    "rownames(dat) = dat[,1]\n",
    "dat = dat[,-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_recomp){  \n",
    "    # make a meta data\n",
    "    group = c(rep(\"HPT_Chow\",6),rep(\"HPT_HFD\",7),rep(\"Hippo_Chow\",6),\n",
    "              rep(\"Hippo_HFD\",7),rep(\"CTX_Chow\",6),rep(\"CTX_HFD\",7)\n",
    "    )\n",
    "    lab = colnames(dat)\n",
    "    loc = c(rep(\"Hypothalamus\",13), rep(\"Hippocampus\",13), rep(\"Cortex\",13))\n",
    "    diet = c(rep(\"chow\",6), rep(\"hfd\",7), rep(\"chow\",6),\n",
    "             rep(\"hfd\",7), rep(\"chow\",6), rep(\"hfd\",7)\n",
    "    )\n",
    "    id = colnames(dat)\n",
    "    metaData = data.frame(id, group, loc,diet)\n",
    "\n",
    "    # remove probleatic samples\n",
    "    dat = dat[,-c(10,13,19,34)] #dat[,-c(10,13,19,34,36)]    19,,36  #10,13,19,32,39\n",
    "    metaData = metaData[-c(10,13,19,34),] #metaData[-c(10,13,19,34,36),]    19,36\n",
    "\n",
    "    # convert from ensembl to symbol gene id\n",
    "    rownames(dat) = gsub(\"\\\\..*\",\"\",rownames(dat))\n",
    "    sig.gene <- bitr(gsub(\"\\\\..*\",\"\",rownames(dat)),\n",
    "                     fromType=\"ENSEMBL\",\n",
    "                     toType=\"SYMBOL\", \n",
    "                     OrgDb=org.Mm.eg.db\n",
    "    )\n",
    "    symb = sig.gene$SYMBOL[ match(rownames(dat), sig.gene$ENSEMBL) ]\n",
    "    dat = dat[!is.na(symb),]\n",
    "    symb = symb[!is.na(symb)]\n",
    "    rownames(dat) = make.names(symb, TRUE)\n",
    "\n",
    "    # filtering\n",
    "    keep = rowSums(dat) > 10\n",
    "    dat = dat[keep,]\n",
    "\n",
    "    keep = rowSums(dat >= 10) >= 5\n",
    "    dat = dat[keep,]\n",
    "\n",
    "    #make a deseq2 object\n",
    "    dds <- DESeqDataSetFromMatrix(countData=dat, \n",
    "                                  colData=metaData, \n",
    "                                  design=~ group\n",
    "    )\n",
    "    dds <- DESeq(dds)\n",
    "    save(dds, file=\"dds.RData\")\n",
    "} else {\n",
    "    load(\"dds.RData\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualisation of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# variance stabilizing transformation data\n",
    "vsd <- assay(vst(dds,\n",
    "                 blind=TRUE)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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=\"SM/PCA_hpt.pdf\", width=10, height=10)\n",
    "    pcaData <- plotPCA(vst(dds[,1:11]),\n",
    "                       intgroup=c(\"group\"),\n",
    "                       returnData=TRUE\n",
    "    )\n",
    "    gr1 = dds@colData$group[1:11]\n",
    "    percentVar = round(100 * attr(pcaData, \"percentVar\"))\n",
    "\n",
    "    ggplot(pcaData, aes(x=PC1, y=PC2, color=gr1)) +\n",
    "           geom_point(aes(color=gr1), size=5) + \n",
    "           geom_point(shape=1, color=\"black\", size=5) + \n",
    "           xlab(paste0(\"PC1: \",percentVar[1],\"% variance\")) + \n",
    "           ylab(paste0(\"PC2: \",percentVar[2],\"% variance\")) +\n",
    "           theme(axis.text=element_text(size=20), \n",
    "                 axis.title=element_text(size=20), \n",
    "                 legend.title=element_text(siz=20),\n",
    "                 legend.text=element_text(size=20),\n",
    "                 plot.title=element_text(size=20)) + \n",
    "           scale_color_manual(values=c(\"chocolate\", \"tan1\")) +\n",
    "           ylim(-25,15) +\n",
    "           xlim(-20,20)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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=\"SM/PCA_hippo.pdf\", width=10, height=10)\n",
    "    pcaData <- plotPCA(vst(dds[,12:23]),\n",
    "                       intgroup=c(\"group\"),\n",
    "                       returnData=TRUE\n",
    "    )\n",
    "    gr1 = dds@colData$group[12:23]\n",
    "    percentVar = round(100 * attr(pcaData, \"percentVar\"))\n",
    "\n",
    "    p <- ggplot(pcaData, aes(x=PC1, y=PC2, color=gr1)) +\n",
    "                geom_point(aes(color=gr1), size=5) + \n",
    "                geom_point(shape=1, color=\"black\", size=5) + \n",
    "                xlab(paste0(\"PC1: \",percentVar[1], \"% variance\")) + \n",
    "                ylab(paste0(\"PC2: \",percentVar[2], \"% variance\")) +\n",
    "                theme(axis.text=element_text(size=20), \n",
    "                      axis.title=element_text(size=20), \n",
    "                      legend.title=element_text(siz=20),\n",
    "                      legend.text=element_text(size=20),\n",
    "                      plot.title=element_text(size=20)) + \n",
    "                scale_color_manual(values=c(\"mediumseagreen\", \"lightgreen\"))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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=\"SM/PCA_ctx.pdf\", width=10, height=10)\n",
    "    pcaData <- plotPCA(vst(dds[,24:35]),\n",
    "                       intgroup=c(\"group\"),\n",
    "                       returnData=TRUE\n",
    "    )\n",
    "    gr1 = dds@colData$group[24:35]\n",
    "    percentVar = round(100 * attr(pcaData, \"percentVar\"))\n",
    "\n",
    "    ggplot(pcaData, aes(x=PC1, y=PC2, color=gr1)) +\n",
    "           geom_point(aes(color=gr1), size=5) + \n",
    "           geom_point(shape=1, color=\"black\", size=5) + \n",
    "           xlab(paste0(\"PC1: \",percentVar[1], \"% variance\")) + \n",
    "           ylab(paste0(\"PC2: \",percentVar[2], \"% variance\")) +\n",
    "           theme(axis.text=element_text(size=20), \n",
    "                 axis.title=element_text(size=20), \n",
    "                 legend.title=element_text(siz=20),\n",
    "                 legend.text=element_text(size=20),\n",
    "                 plot.title=element_text(size=20)) + \n",
    "           scale_color_manual(values=c(\"lightblue4\", \"lightblue\"))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if(bool_plot){ \n",
    "    df_pca <- prcomp(t(vsd))\n",
    "\n",
    "    df_out = as.data.frame(df_pca$x)\n",
    "    df_out$tissue = as.factor(dds@colData$group) \n",
    "\n",
    "    ggplot(df_out, aes(x=PC1, y=PC2, color=tissue)) + \n",
    "           geom_point(aes(color=tissue), size=5) + \n",
    "           geom_point(shape=1, color=\"black\", size=5) + \n",
    "           xlab(paste0(\"PC1: \",round(df_pca$sdev[1]), \"% variance\")) + \n",
    "           ylab(paste0(\"PC2: \",round(df_pca$sdev[2]), \"% variance\")) +\n",
    "           theme(axis.text=element_text(size=20), \n",
    "                 axis.title=element_text(size=20), \n",
    "                 legend.title=element_text(siz=20),\n",
    "                 legend.text=element_text(size=20),\n",
    "                 plot.title=element_text(size=20)) + \n",
    "           scale_color_manual(values=c(\"lightblue\", \"lightblue4\", \"tan1\",\n",
    "                                       \"chocolate\", \"lightgreen\", \"mediumseagreen\"))\n",
    "}"
   ]
  },
  {
   "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": [
    "pdf(file=\"PCA.pdf\", width=10, height=10)\n",
    "    pcaData <- plotPCA(vst(dds),\n",
    "                       intgroup=c(\"group\"),\n",
    "                       returnData=TRUE\n",
    "    )\n",
    "    percentVar = round(100 * attr(pcaData, \"percentVar\"))\n",
    "\n",
    "    ggplot(pcaData, aes(x=PC1, y=PC2, color=dds@colData$group)) +\n",
    "           geom_point(aes(color=dds@colData$group), size=10) + \n",
    "           geom_point(shape=1, color=\"black\", size=10) + \n",
    "           xlab(paste0(\"PC1: \",percentVar[1],\"% variance\")) + \n",
    "           ylab(paste0(\"PC2: \",percentVar[2],\"% variance\")) +\n",
    "           theme(axis.text=element_text(size=25), \n",
    "                 axis.title=element_text(size=30)) +\n",
    "           scale_color_manual(values=c(\"aquamarine4\", \"aquamarine2\", \"darkorange2\",\n",
    "                                       \"orange\", \"blueviolet\", \"mediumpurple1\")) +\n",
    "           theme(legend.position=\"none\")\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>pdf:</strong> 3"
      ],
      "text/latex": [
       "\\textbf{pdf:} 3"
      ],
      "text/markdown": [
       "**pdf:** 3"
      ],
      "text/plain": [
       "pdf \n",
       "  3 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Heatmap of all samples and genes\n",
    "pdf(file=\"SM/heatmap_alldat.pdf\", width=10, height=10)\n",
    "    agsh <- pheatmap(vsd,\n",
    "                     cluster_rows=T, \n",
    "                     scale=\"row\", \n",
    "                     show_rownames=FALSE, \n",
    "                     color=my_palette,\n",
    "                     cluster_cols=F, \n",
    "                     annotation_col=as.data.frame(colData(dds)[,c(\"diet\", \"loc\")]),\n",
    "                     fontsize=10,\n",
    "                     breaks=seq(-3, 3, length.out=255)\n",
    "    )\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differential gene expression - diet effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_recomp){ \n",
    "    hpt <- results(dds, \n",
    "                   contrast=c(\"group\", \"HPT_HFD\", \"HPT_Chow\")\n",
    "    )\n",
    "    hippo <- results(dds, \n",
    "                     contrast=c(\"group\", \"Hippo_HFD\", \"Hippo_Chow\")\n",
    "    )\n",
    "    ctx <- results(dds,\n",
    "                   contrast=c(\"group\", \"CTX_HFD\", \"CTX_Chow\")\n",
    "    )\n",
    "\n",
    "    dge = list(hpt=as.data.frame(hpt[complete.cases(hpt), ]),\n",
    "               hippo=as.data.frame(hippo[complete.cases(hippo), ]),\n",
    "               ctx=as.data.frame(ctx[complete.cases(ctx), ])\n",
    "    )\n",
    "    save(dge, file=\"dge.RData\")\n",
    "} else {\n",
    "    load(\"dge.RData\")\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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_hpt.pdf\", width=10, height=10)\n",
    "    EnhancedVolcano(dge[['hpt']],\n",
    "                    lab=rownames(dge[['hpt']]),\n",
    "                    pCutoff=0.05,\n",
    "                    FCcutoff=100,\n",
    "                    axisLabSize=30,\n",
    "                    x='log2FoldChange',\n",
    "                    y='pvalue',\n",
    "                    title=NULL,\n",
    "                    subtitle=NULL,\n",
    "                    caption=NULL,\n",
    "                    ylim=c(0, 8),\n",
    "                    xlim=c(-8, 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": 16,
   "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_hippo.pdf\", width=10, height=10)\n",
    "    EnhancedVolcano(dge[['hippo']],\n",
    "                    lab=rownames(dge[['hippo']]),\n",
    "                    pCutoff=0.05,\n",
    "                    FCcutoff=100,\n",
    "                    axisLabSize=30,\n",
    "                    x='log2FoldChange',\n",
    "                    y='pvalue',\n",
    "                    title=NULL,\n",
    "                    subtitle=NULL,\n",
    "                    caption=NULL,\n",
    "                    ylim=c(0, 8),\n",
    "                    xlim=c(-8, 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": 17,
   "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_ctx.pdf\", width=10, height=10)\n",
    "    EnhancedVolcano(dge[['ctx']],\n",
    "                    lab = rownames(dge[['ctx']]),\n",
    "                    pCutoff=0.05,\n",
    "                    FCcutoff=100,\n",
    "                    axisLabSize=30,\n",
    "                    x='log2FoldChange',\n",
    "                    y='pvalue',\n",
    "                    title=NULL,\n",
    "                    subtitle=NULL,\n",
    "                    caption=NULL,\n",
    "                    ylim=c(0, 8),\n",
    "                    xlim=c(-8, 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": 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"
    }
   ],
   "source": [
    "pdf(file=\"SM/Nr_sig_0.05.pdf\", width=10, height=10)\n",
    "    df = data.frame(Tissue=c(\"CTX\", \"HPT\", \"Hippo\"),\n",
    "                    Nr_significant_genes=c(sum(dge[['ctx']][,5] < 0.05, na.rm=T),\n",
    "                                           sum(dge[['hpt']][,5] < 0.05, na.rm=T), \n",
    "                                           sum(dge[['hippo']][,5] < 0.05, na.rm=T))\n",
    "    )\n",
    "    ggplot(df, aes(x=Tissue, y=Nr_significant_genes, fill=Tissue))  +\n",
    "           geom_bar(stat=\"identity\") + \n",
    "            theme(axis.text=element_text(size=20), \n",
    "                  axis.title=element_text(size=20), \n",
    "                  legend.title=element_text(siz=20),\n",
    "                  legend.text=element_text(size=20),\n",
    "                  plot.title=element_text(size=20)) +\n",
    "           theme(legend.position=\"none\") +\n",
    "           ggtitle(\"Number significant genes (p<0.05)\") +\n",
    "           scale_fill_manual(values=c(\"lightblue\", \"tan1\", \"lightgreen\"))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "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_tran.pdf\", width=8, height=8)\n",
    "    vp <- venn.diagram(list(' '=rownames(dge[['hpt']][dge[['hpt']][,5] < 0.05,]), \n",
    "                            ' '=rownames(dge[['hippo']][dge[['hippo']][,5] < 0.05,]),\n",
    "                            ' '=rownames(dge[['ctx']][dge[['ctx']][,5] < 0.05,])),\n",
    "                       fill=c(\"orange\", \"mediumpurple1\", \"aquamarine2\"),\n",
    "                       filename=NULL,\n",
    "                       cex=3.5, \n",
    "                       col=\"transparent\",\n",
    "                       cat.cex=0.01,\n",
    "                       margin=0.01,\n",
    "                       main.cex=3,\n",
    "                       main=NULL);\n",
    "    grid.draw(vp)  \n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heatmap marker genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_recomp){ \n",
    "    vsd1 = as.data.frame(vsd[,c(24:35, 1:23)])\n",
    "    datfin = data.frame(CTX=apply(vsd1 %>% select(starts_with(\"ctx_\")), 1, mean),\n",
    "                        HPT=apply(vsd1 %>% select(starts_with(\"hpt_\")), 1, mean),\n",
    "                        Hippo=apply(vsd1 %>% select(starts_with(\"hippo_\")), 1, mean))\n",
    "\n",
    "    astro = datfin[rownames(datfin)%in%c('Slc1a2', 'Slc1a3', 'Aqp4', 'Gfap', 'Aldh1l1', 'Gja1', 'Gjb6', 'Atp1b2'),]\n",
    "    neuro = datfin[rownames(datfin)%in%c('Syp', 'Tubb3', 'Snap25', 'Syt1'),]\n",
    "    mglia = datfin[rownames(datfin)%in%c('Itgam', 'Aif1'),]\n",
    "    oligo = datfin[rownames(datfin)%in%c('Mog', 'Mag'),]\n",
    "    endo = datfin[rownames(datfin)%in%c('Slco1c1'),]\n",
    "    mural = datfin[rownames(datfin)%in%c('Mustn1', 'Pdgfrb', 'Des'),]\n",
    "    tcyes = datfin[rownames(datfin)%in%c('Crym'),]\n",
    "    npc = datfin[rownames(datfin)%in%c('Nes', 'Sox2', 'Prom1'),]\n",
    "    vlmc = datfin[rownames(datfin)%in%c('Col1a1'),]\n",
    "\n",
    "    astro = astro[order(row.names(astro)),]\n",
    "    neuro = neuro[order(row.names(neuro)),]\n",
    "    mglia = mglia[order(row.names(mglia)),]\n",
    "    oligo = oligo[order(row.names(oligo)),]\n",
    "    mural = mural[order(row.names(mural)),]\n",
    "    tcyes = tcyes[order(row.names(tcyes)),]\n",
    "    mat = rbind(astro, neuro, mglia, oligo, endo, mural) # npc,vlmc,tcyes\n",
    "    \n",
    "    write.csv(mat,'celltype_heatmap.csv', row.names=T)\n",
    "} else{\n",
    "    mat=read.csv('celltype_heatmap.csv', row.names=1)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Omics = c(rep(\"transcriptomics\", 3))\n",
    "#Acsa2 = rep(\"Acsa2_pos\", 3)\n",
    "#Tissue = c(\"cortex\", \"hypothalamus\", \"hippocampus\")\n",
    "#dcol = as.data.frame(cbind(Tissue,Acsa2,Omics))\n",
    "#rownames(dcol) = colnames(mat)\n",
    "mat = as.matrix(mat)\n",
    "pheatmap(mat, \n",
    "         scale=\"column\",\n",
    "         show_rownames=F, \n",
    "         show_colnames=F,\n",
    "         color=my_palette, \n",
    "         cluster_cols=F,\n",
    "         cluster_rows=F, \n",
    "         gaps_col=3,\n",
    "         gaps_row=c(8,12,14,16,17,20),#(11,15,17,19,20,23), \n",
    "         fontsize=18,\n",
    "         filename=\"heatmap_ACSA2.pdf\",\n",
    "         legend=F,\n",
    "         width=3.3, \n",
    "         height=8\n",
    ") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heatmap DGE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### HPT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "vsd_hpt = vsd[rownames(vsd)%in%rownames(dge[['hpt']][dge[['hpt']][,5] < 0.05,]),] \n",
    "pheatmap(vsd_hpt[,1:11], \n",
    "         scale=\"row\",\n",
    "         show_rownames=FALSE,\n",
    "         show_colnames=FALSE,\n",
    "         color=my_palette, \n",
    "         cluster_cols=F,\n",
    "         legend=FALSE,\n",
    "         fontsize=15,\n",
    "         filename=\"heatmap_hpt.pdf\",\n",
    "         width=5,\n",
    "         height=8,\n",
    "         breaks=seq(-2, 2, length.out=255)\n",
    ") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hippo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "vsd_hip = vsd[rownames(vsd)%in%rownames(dge[['hippo']][dge[['hippo']][,5] < 0.05,]),] \n",
    "pheatmap(vsd_hip[,12:23], \n",
    "         scale=\"row\", \n",
    "         show_rownames=FALSE,\n",
    "         show_colnames=FALSE,\n",
    "         color=my_palette,\n",
    "         cluster_cols=F, \n",
    "         legend=FALSE,\n",
    "         fontsize=15, \n",
    "         filename=\"heatmap_hippo.pdf\",\n",
    "         width=5, \n",
    "         height=8,\n",
    "         breaks=seq(-2, 2, length.out=255)\n",
    ") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### CTX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "vsd_ctx = vsd[rownames(vsd)%in%rownames(dge[['ctx']][dge[['ctx']][,5] < 0.05,]),] \n",
    "pheatmap(vsd_ctx[,24:35],\n",
    "         scale=\"row\",\n",
    "         show_rownames=FALSE,\n",
    "         show_colnames=FALSE,\n",
    "         color=my_palette, \n",
    "         cluster_cols=F,\n",
    "         legend=F,\n",
    "         fontsize=15,\n",
    "         filename=\"heatmap_ctx.pdf\",\n",
    "         width=5, \n",
    "         height=8,\n",
    "         breaks=seq(-2, 2, length.out=255)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GO - enrichment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(bool_recomp){ \n",
    "    sg.ctx <- bitr(rownames(vsd_ctx),\n",
    "                   fromType=\"SYMBOL\", \n",
    "                   toType=\"ENTREZID\", \n",
    "                   OrgDb=org.Mm.eg.db\n",
    "    )\n",
    "    go.ctx <- enrichGO(sg.ctx[,2], \n",
    "                       'org.Mm.eg.db', \n",
    "                       ont=\"BP\",\n",
    "                       pvalueCutoff=0.1\n",
    "    )\n",
    "\n",
    "    sg.hpt <- bitr(rownames(vsd_hpt),\n",
    "                   fromType=\"SYMBOL\",\n",
    "                   toType=\"ENTREZID\",\n",
    "                   OrgDb=org.Mm.eg.db\n",
    "    )\n",
    "    go.hpt <- enrichGO(sg.hpt[,2], \n",
    "                       'org.Mm.eg.db',\n",
    "                       ont=\"BP\",\n",
    "                       pvalueCutoff=0.1\n",
    "    )\n",
    "\n",
    "    sg.hip <- bitr(rownames(vsd_hip),\n",
    "                   fromType=\"SYMBOL\",\n",
    "                   toType=\"ENTREZID\",\n",
    "                   OrgDb=org.Mm.eg.db\n",
    "    )\n",
    "    go.hip <- enrichGO(sg.hip[,2],\n",
    "                       'org.Mm.eg.db',\n",
    "                       ont=\"BP\",\n",
    "                       pvalueCutoff=0.1\n",
    "    )\n",
    "\n",
    "    go = c(\"insulin receptor signaling pathway\",\n",
    "           \"response to insulin\",\n",
    "           \"TOR signaling\",\n",
    "           \"pyruvate metabolic process\",\n",
    "           \"intracellular lipid transport\",\n",
    "           \"fatty acid metabolic process\",\n",
    "           \"fatty acid oxidation\",\n",
    "           \"response to reactive oxygen species\",\n",
    "           \"response to endoplasmic reticulum stress\",\n",
    "           \"regulation of mitotic cell cycle\",\n",
    "           \"transforming growth factor beta receptor signaling pathway\",\n",
    "           \"I-kappaB kinase/NF-kappaB signaling\",\n",
    "           \"regulation of inflammatory response\",\n",
    "           \"regulation of cell junction assembly\",\n",
    "           \"cell-matrix adhesion\",\n",
    "           \"regulation of cytosolic calcium ion concentration\",\n",
    "           \"vesicle budding from membrane\",\n",
    "           \"axon ensheathment\",\n",
    "           \"maintenance of synapse structure\",\n",
    "           \"regulation of synapse structure or activity\",\n",
    "           \"synaptic transmission, GABAergic\",\n",
    "           \"AMPA glutamate receptor clustering\",\n",
    "           \"positive regulation of blood vessel endothelial cell migration\",\n",
    "           \"sprouting angiogenesis\",\n",
    "           \"regulation of angiogenesis\"\n",
    "    )\n",
    "    \n",
    "    enr_ctx.go = go.ctx@result[go.ctx@result[,6] < 0.05, c(2,5,6,8,9)]\n",
    "    enr_ctx.go$Tissue = 'CTX'\n",
    "    enr_hpt.go = go.hpt@result[go.hpt@result[,6] < 0.05, c(2,5,6,8,9)]\n",
    "    enr_hpt.go$Tissue = 'HPT'\n",
    "    enr_hip.go = go.hip@result[go.hip@result[,6] < 0.05, c(2,5,6,8,9)]\n",
    "    enr_hip.go$Tissue = 'Hippo'\n",
    "    \n",
    "    all = rbind(enr_ctx.go, enr_hpt.go, enr_hip.go)\n",
    "    #write.csv(all, file=\"GO_transcriptome.csv\")\n",
    "    all = all[all$Description%in%go,]\n",
    "    options(repr.plot.width=20, repr.plot.height=15)\n",
    "    colnames(all) = c('Description', 'pvalue', 'p.adjust', 'geneID', 'count', 'Tissue')\n",
    "\n",
    "    all=rbind(enr_ctx.go, enr_hpt.go, enr_hip.go)\n",
    "    all1=all[all$Description%in%go,]\n",
    "\n",
    "    sig.gene.hpt.all <- bitr(rownames(dge[['hpt']]),\n",
    "                             fromType=\"SYMBOL\",\n",
    "                             toType=\"ENTREZID\", \n",
    "                             OrgDb=org.Mm.eg.db\n",
    "    )\n",
    "    sig.gene.ctx.all <- bitr(rownames(dge[['ctx']]), \n",
    "                             fromType=\"SYMBOL\",\n",
    "                             toType=\"ENTREZID\",\n",
    "                             OrgDb=org.Mm.eg.db\n",
    "    )\n",
    "    sig.gene.hip.all <- bitr(rownames(dge[['hippo']]),\n",
    "                             fromType=\"SYMBOL\",\n",
    "                             toType=\"ENTREZID\",\n",
    "                             OrgDb=org.Mm.eg.db\n",
    "    )\n",
    "    \n",
    "    # CTX\n",
    "    df_total.ctx = data.frame()\n",
    "    for(i in 1:nrow(all1)){   \n",
    "        p = as.data.frame(dge[['ctx']][which(rownames(dge[['ctx']])%in%\n",
    "                          sig.gene.ctx.all[sig.gene.ctx.all[,2]%in%unlist(strsplit(all1$geneID[i],\"/\")),1]), c(2,5)])\n",
    "        p$pathway = all1$Description[i]\n",
    "        p$padjPath = all1$p.adjust[i]\n",
    "        p$count = all1$Count[i]\n",
    "        df_total.ctx = rbind(df_total.ctx, p)    \n",
    "    }\n",
    "    df_total.ctx$significance <- ifelse(df_total.ctx$pvalue < 0.05, \n",
    "                                        \"significant\", \n",
    "                                        \"not significant\"\n",
    "    )\n",
    "    # HPT\n",
    "    df_total.hpt = data.frame()\n",
    "    for(i in 1:nrow(all1)){\n",
    "        p = as.data.frame(dge[['hpt']][which(rownames(dge[['hpt']])%in%\n",
    "                          sig.gene.hpt.all[sig.gene.hpt.all[,2]%in%unlist(strsplit(all1$geneID[i],\"/\")),1]),c(2,5)])\n",
    "        p$pathway = all1$Description[i]\n",
    "        p$padjPath = all1$p.adjust[i]\n",
    "        p$count = all1$Count[i]\n",
    "        df_total.hpt = rbind(df_total.hpt, p)    \n",
    "    } \n",
    "    df_total.hpt$col <- ifelse(df_total.hpt$pvalue < 0.05, \n",
    "                               \"red\", \n",
    "                               \"grey\"\n",
    "    )\n",
    "    \n",
    "    # Hippo\n",
    "    df_total.hip = data.frame()\n",
    "    for(i in 1:nrow(all1)){\n",
    "        p = as.data.frame(dge[['hippo']][which(rownames(dge[['hippo']])%in%\n",
    "                          sig.gene.hip.all[sig.gene.hip.all[,2]%in%unlist(strsplit(all1$geneID[i],\"/\")),1]), c(2,5)])\n",
    "        p$pathway = all1$Description[i]\n",
    "        p$padjPath = all1$p.adjust[i]\n",
    "        p$count = all1$Count[i]\n",
    "        df_total.hip = rbind(df_total.hip, p)    \n",
    "    }    \n",
    "    df_total.hip$col <- ifelse(df_total.hip$pvalue < 0.05,\n",
    "                               \"red\",\n",
    "                               \"grey\"\n",
    "    )\n",
    "\n",
    "    gsea = list(all=all,\n",
    "                ctx=df_total.ctx,\n",
    "                hpt=df_total.hpt, \n",
    "                hippo=df_total.hip,\n",
    "                go=go\n",
    "    )\n",
    "    save(gsea, file=\"gsea.RData\")\n",
    "} else {\n",
    "    load(\"gsea.RData\")\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "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.pdf\", width=15, height=15)\n",
    "    ggplot(gsea[['all']],   \n",
    "           aes(Tissue, y=factor(Description, level=gsea[['go']][length(gsea[['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()) + theme(legend.position=\"none\") +\n",
    "           theme(text=element_text(size=40))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### CTX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 12 rows containing missing values (geom_point).”\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=\"enrichdot_ctx.pdf\", width=7, height=15)\n",
    "    ggplot(gsea[['ctx']], aes(log2FoldChange, y=factor(pathway, level=gsea[['go']][length(gsea[['go']]):1]))) + \n",
    "           geom_point(aes(colour=significance, size=-log10(pvalue)), alpha=0.5) + \n",
    "           scale_size_continuous(range=c(0,8)) +\n",
    "           scale_color_manual(values=c(\"gray85\", \"red\")) +\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)) + #guides(color = guide_legend(override.aes = list(size=5) ) ) +\n",
    "           theme(legend.position=\"none\") + \n",
    "           theme(axis.text.y=element_blank()) +\n",
    "           ylab(\"\") + \n",
    "           xlab(\"\") + \n",
    "           xlim(-2.68, 3.82)\n",
    "dev.off() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### HPT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 5 rows containing missing values (geom_point).”\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=\"enrichdot_hpt.pdf\", width=7, height=15)\n",
    "\n",
    "    ggplot(gsea[['hpt']], aes(log2FoldChange, y=factor(pathway, level=gsea[['go']][length(gsea[['go']]):1]))) + \n",
    "           geom_point(aes(colour=col, size=-log10(pvalue)), alpha=0.5) + \n",
    "           scale_size_continuous(range=c(0, 8)) +\n",
    "           scale_color_manual(values=c(\"gray85\", \"red\")) + \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(-2.68,3.82)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hippo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Removed 15 rows containing missing values (geom_point).”\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=\"enrichdot_hippo.pdf\", width=7, height=15)\n",
    "\n",
    "    ggplot(gsea[['hippo']], aes(log2FoldChange, y=factor(pathway, level=gsea[['go']][length(gsea[['go']]):1]))) + \n",
    "           geom_point(aes(colour=col, size=-log10(pvalue)), alpha=0.5) + \n",
    "           scale_color_manual(values=c(\"gray85\", \"red\")) + \n",
    "           scale_size_continuous(range=c(0, 6)) + \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(-2.68, 3.82)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "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": 2
}