library(lme4)
library(lmerTest)
library(ggplot2)
library(ggpubr)
library(tidyverse)
library(dplyr)
library(RColorBrewer)
library(multcomp)

#setwd("~/Dropbox/Glasgow/Projects/Reuniens CA1/CA1 pyramidal scRNAseq data Spruston")
#Load differential expression data from 2016 Neuron paper
#
#diff_exp<-read.delim(file="GSE67403_gene_exp.diff", header=TRUE, sep="\t")
#Load gene reads data from 2016 Neuron paper - data can be downloaded from https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE67403
#

Neuron2016_data<-read.delim(file="GSE67403_genes.read_group_tracking.txt", header=TRUE, sep="\t")
#Load gene reads dara from 2016 eLife paper
#eLife paper has exactly same number of reads for vCA1 and dCA1, so these are the same samples included in analysis 
eLife2016_data<-read.delim(file="GSE74985_genes.read_group_tracking.txt", header=TRUE, sep="\t")

#filter for Emx1
Emx1_Neuron_paper<-Neuron2016_data%>%
  filter(tracking_id=="ENSMUSG00000033726")%>%
  mutate(condition = case_when(
    condition == "dorsal" ~ "dCA1_whole",
    condition == "ventral" ~ "vCA1_whole",
    condition == "intermediate" ~ "iCA1_whole",
    condition == "deep" ~ "dCA1_deep",
    condition == "superficial" ~ "dCA1_superficial",
    condition == "proximal" ~ "dCA1_proximal",
    condition == "distal" ~ "dCA1_distal",
  ))
Emx1_data<-Emx1_Neuron_paper %>%
  na.omit

#Plot graph
my.pallet<-c(brewer.pal(7,"Greens")[6],brewer.pal(5,"Blues")[3],brewer.pal(6,"Set1")[6],brewer.pal(7,"Greens")[1:4])


ggbarplot(Emx1_data, x= "condition", y="FPKM", fill="condition",
          colour="condition", palette = my.pallet,
          add = c("mean_se","jitter"))+
  rremove("legend")+
  rotate_x_text(angle=45)+
  labs(x="Region",y="Emx1 FPKM")
 

#filter for Emx1 whole hippocampus only
Emx1_WholeHPC<-Neuron2016_data%>%
  filter(tracking_id=="ENSMUSG00000033726")%>%
  mutate(condition = case_when(
    condition == "dorsal" ~ "dCA1_whole",
    condition == "ventral" ~ "vCA1_whole",
    condition == "intermediate" ~ "iCA1_whole"
  ))
Emx1_WholeHPC<-Emx1_WholeHPC %>%
  na.omit

  Emx1_test<-aov(FPKM ~ condition, data=Emx1_WholeHPC)
  summary(Emx1_test)