README.Rmd

---
title: "TPH2 KO statistics"
output: github_document
author: "Alex Meng, Joanes Grandjean"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval=FALSE,warning = FALSE,message=FALSE)
```


# Setup environement
Download DABEST
```{r setup env}
# you just need to run these once. 

install.packages('devtools')
install.pacakges("reshape2")
install.pacakges("wesanderson")
install.packages("svglite")
devtools::install_github("ACCLAB/dabestr")   
devtools::install_github("karthik/wesanderson")
```

# Load pacakges

```{r load env}
library(tidyverse)
library(glue)
library(dabestr)
library(wesanderson)  # see https://github.com/karthik/wesanderson
library(reshape2)
library(svglite)
pal <- wes_palette("Darjeeling1")
```


# Figure 1: Anxiety level in THP2. if necessary, show rest in table
A. experimental diagram
B. EPM time in open arm
C. EPM total distance

# Figure 2: Social behaviour. 
A. experimetnal diagram
B. Total no contact
C. Total no contact over time (as line plot)
D. Total mounting %
E. Total aggressive %

# Figure 3: SHow gene expression left, protein right, anat center? rest is in table. 
A. Brain punches location (color-code gene expression / protein)
B. IL/PL
C. centroid amygdala
D. ventral CA1
E. Dorsal raphe


# demo dabest in R

```{r demo dabest}

#load the table. for other cases, it might help to keep naming consistent between tables, see example. Also, avoid spaces or weird characters in table names. 
df <- read_csv('assets/tables/SI_total_no_contact.csv', col_types = cols()) %>% melt() %>% dplyr::rename(group = variable) %>% drop_na()

df %>% group_by(group) %>% summarise(mean=round(mean(value),2), sd=round(sd(value),2))

#Q-Q plot
qqnorm(df$value, pch = 1, frame = FALSE)
qqline(df$value, col = "steelblue", lwd = 2)

#shapiro test
shapiro.test(df$value)

# estimate the p-values using non-parametric test
p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2+/-'],var.equal=TRUE)$p.value, t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

# estimate hedges'g
dabest_hedges <- dabest(df, group, value, idx = c('Tph2+/+','Tph2+/-','Tph2-/-'), paired = FALSE)  %>% hedges_g() 
dabest_hedges


#make the plot and save to file
Fig1A<-plot(dabest_hedges, palette = pal,  rawplot.ylabel = "total no contact (%)") 
ggsave('assets/figure/SI total no contact.svg', plot = Fig1A, device = 'svg',dpi = 300)

#outputs the file to show in the doc. 
Fig1A

#outputs table with added p-values for good measure
dabest_hedges$result %>% mutate(p = p_tmp)
 
```

# Code for PCR and Oxytocin protein
```{r Code for two groups comparision}
#load table
df <- read_csv('assets/tables/OXY_PVN.csv', col_types = cols()) %>% melt() %>% dplyr::rename(group = variable) %>% drop_na()

#mean + SD

df %>% group_by(group) %>% summarise(mean=round(mean(value),2), sd=round(sd(value),2))

# estimate hedges'g

dabest_hedges <- dabest(df, group, value, idx = c('Tph2+/+','Tph2-/-'), paired = FALSE)  %>% hedges_g() 

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig2 <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "pg oxytocin/μg tissue protein/ml")
# estimate the p-values using non-parametric test
p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

#output table with p value
dabest_hedges$result %>% mutate(p = p_tmp)

ggsave('assets/figure/OXY_PVN.svg', plot = Fig2, device = 'svg',dpi = 300)
```


# what happens if dabest bugs? 

```{r demo dabest debug}

#load the table. for other cases, it might help to keep naming consistent between tables, see example. Also, avoid spaces or weird characters in table names. 
df <- read_csv('assets/tables/testname_measure.csv', col_types = cols()) %>% melt() %>% rename(group = variable) %>% drop_na()



dabest_hedges <- dabest(df, group, value, idx = c('Tph2+/+','Tph2+/-'), paired = FALSE)  %>% hedges_g() 

#figure will require post-pocessing to make do. 
plot(dabest_hedges, palette = pal)



dabest_hedges <- dabest(df, group, value, idx = c('Tph2+/+','Tph2-/-'), paired = FALSE)  %>% hedges_g() 

#figure will require post-pocessing to make do. 
plot(dabest_hedges, palette = pal)


 
```


Finally, you can consider https://rpkgs.datanovia.com/ggpubr/reference/ggarrange.html to make composite figures in R. Happy to help you get started.