---
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")
```




###behavioural tests EPM (Fig1)
##Fig1 dataset description:
A. Experimental diagram
B. Open arms duration
C. Number of entry into closed arms
D. Latency time for first time entering open arms
E. Total distance moved
##FigS1 dataset description:
A. Closed arms duration
B. Number of entry into open arms

#Fig1B
```{r EPM_open arms duration}

df <- read_csv('assets/tables/EPM_open_arms_duration.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))

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)

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

Fig1A<-plot(dabest_hedges, palette = pal,  rawplot.ylabel = "duration in open arms (%)") 
ggsave('assets/figure/Fig1B.svg', plot = Fig1A, device = 'svg',dpi = 300)

Fig1A

dabest_hedges$result %>% mutate(p = p_tmp)
 
```

#Fig1C
```{r EPM_entry closed arms}
df <- read_csv('assets/tables/EPM_entry_closed_arms.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))

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)

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

Fig1C<-plot(dabest_hedges, palette = pal,  rawplot.ylabel = "number of entry") 
ggsave('assets/figure/EPM_entry_closed_arms.svg', plot = Fig1C, device = 'svg',dpi = 300)

Fig1C

dabest_hedges$result %>% mutate(p = p_tmp)
```

#Fig1D
```{r EPM_first entry latency}
df <- read_csv('assets/tables/EPM_first_time_entry.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))

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)

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

Fig1E<-plot(dabest_hedges, palette = pal,  rawplot.ylabel = "latency time (s)") 
ggsave('assets/figure/EPM_first_time_entry.svg', plot = Fig1E, device = 'svg',dpi = 300)

Fig1E

dabest_hedges$result %>% mutate(p = p_tmp)
```

#Fig1E
```{r EPM_total distance}
df <- read_csv('assets/tables/EPM_total_distance.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))

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)

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

Fig1F<-plot(dabest_hedges, palette = pal,  rawplot.ylabel = "total distance moved (cm)") 
ggsave('assets/figure/EPM total distance.svg', plot = Fig1F, device = 'svg',dpi = 300)

Fig1F

dabest_hedges$result %>% mutate(p = p_tmp)
```

#FigS1A
```{r EPM_closed arms duration}
df <- read_csv('assets/tables/EPM_closed_arms_duration.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))

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)

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

Fig1B<-plot(dabest_hedges, palette = pal,  rawplot.ylabel = "duration in closed arms (%)") 
ggsave('assets/figure/EPM_closed_arms_duration.svg', plot = Fig1B, device = 'svg',dpi = 300)

Fig1B

dabest_hedges$result %>% mutate(p = p_tmp)
```

#FigS1B
```{r EPM_entry open arms}
df <- read_csv('assets/tables/EPM_entry_open_arms.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))

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)

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

Fig1D<-plot(dabest_hedges, palette = pal,  rawplot.ylabel = "number of entry") 
ggsave('assets/figure/EPM_entry_open_arms.svg', plot = Fig1D, device = 'svg',dpi = 300)

Fig1D

dabest_hedges$result %>% mutate(p = p_tmp)
```

###behavioural tests Social interaction (Fig2)
##Fig2 dataset description:
A. Experimental diagram
B. Total no contact
C. Total mounting
D. Total aggressiveness
##FigS2 dataset description:
A. Total no contact at different time intervals
B. Total mounting at different time intervals
C. Total aggressiveness at different time intervals

#Fig1B
```{r SI_total no contact}
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))

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)

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

Fig1G<-plot(dabest_hedges, palette = pal,  rawplot.ylabel = "total no contact (%)") 
ggsave('assets/figure/SI total no contact.svg', plot = Fig1G, device = 'svg',dpi = 300)

Fig1G

dabest_hedges$result %>% mutate(p = p_tmp)
```

#Fig1C
```{r SI_total mounting}
df <- read_csv('assets/tables/SI_total_mounting.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))

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)

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

Fig1H<-plot(dabest_hedges, palette = pal,  rawplot.ylabel = "total mounting (%)") 
ggsave('assets/figure/SI total mounting.svg', plot = Fig1H, device = 'svg',dpi = 300)

Fig1H

dabest_hedges$result %>% mutate(p = p_tmp)
```

#Fig1D
```{r SI_total aggressiveness}
df <- read_csv('assets/tables/SI_total_aggressiveness.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))

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)

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

Fig1I<-plot(dabest_hedges, palette = pal,  rawplot.ylabel = "total aggressiveness (%)") 
ggsave('assets/figure/SI total aggressiveness.svg', plot = Fig1I, device = 'svg',dpi = 300)

Fig1I

dabest_hedges$result %>% mutate(p = p_tmp)
```

###molecular tests PCR for Oxytocin receptors (Fig3)
##Fig3 dataset description:
A: Diagram of brain punching position
B. Infralimbic cortex
C. Prelimbic cortex
D. Paraventricular nucleus
E. Dorsal raphe nucleus
F. Dorsal granular layer of dentate gyrus
G. Ventral CA1 region of hippocampus
H. Ventral CA3 region of hippocampus
I. Ventral granular layer of dentate gyrus

##FigS3 dataset description
A. Central amygdala
B. Dorsal CA1 region of hippocampus
C. Dorsal CA3 region of hippocampus

#Fig3B
```{r PCR_IF}

df <- read_csv('assets/tables/PCR_IF.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig2A <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "mRNA level (% vs Tph2+/+)")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#Fig3C
```{r PCR_PRL}
df <- read_csv('assets/tables/PCR_PRL.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig2B <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "mRNA level (% vs Tph2+/+)")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#Fig3D
```{r PCR_PVN}
df <- read_csv('assets/tables/PCR_PVN.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig2C <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "mRNA level (% vs Tph2+/+)")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#Fig3E
```{r PCR_DR}
df <- read_csv('assets/tables/PCR_DR.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig2D <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "mRNA level (% vs Tph2+/+)")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#Fig3F
```{r PCR_dGRDG}
df <- read_csv('assets/tables/PCR_dGRDG.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig2H <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "mRNA level (% vs Tph2+/+)")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#Fig3G
```{r PCR_vCA1}
df <- read_csv('assets/tables/PCR_vCA1.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig2I <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "mRNA level (% vs Tph2+/+)")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#Fig3H
```{r PCR_vCA3}
df <- read_csv('assets/tables/PCR_vCA3.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig2J <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "mRNA level (% vs Tph2+/+)")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#Fig3I
```{r PCR_vGRDG}
df <- read_csv('assets/tables/PCR_vGRDG.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig2K <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "mRNA level (% vs Tph2+/+)")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#FigS3A
```{r PCR_CA}
df <- read_csv('assets/tables/PCR_CA.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig2E <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "mRNA level (% vs Tph2+/+)")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#FigS3B
```{r PCR_dCA1}
df <- read_csv('assets/tables/PCR_dCA1.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig2F <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "mRNA level (% vs Tph2+/+)")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#FigS3C
```{r PCR_dCA3}
df <- read_csv('assets/tables/PCR_dCA3.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig2G <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "mRNA level (% vs Tph2+/+)")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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


###molecular tests ELISA for Oxytocin (Fig4)
##Fig3 dataset description:
A. Diagram of brain punching position
B. Medial frontal cortex
C. Paraventricular nucleus
D. Central amygdala
E. Hippocampus

##FigS4 dataset description:
A. Dorsal raphe nucleus

#Fig4B
```{r OXY_MFC}
df <- read_csv('assets/tables/OXY_MFC.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig3A <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "pg oxytocin/µg tissue protein/ml")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#Fig4C
```{r OXY_PVN}
df <- read_csv('assets/tables/OXY_PVN.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig3B <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "pg oxytocin/µg tissue protein/ml")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#Fig4D
```{r OXY_CA}
df <- read_csv('assets/tables/OXY_CA.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig3D <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "pg oxytocin/µg tissue protein/ml")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#Fig4E
```{r OXY_HIP}
df <- read_csv('assets/tables/OXY_HIP.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig3E <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "pg oxytocin/µg tissue protein/ml")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

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

#FigS4A
```{r OXY_DR}
df <- read_csv('assets/tables/OXY_DR.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))

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

dabest_hedges

plot(dabest_hedges, palette = pal)

Fig3C <- plot(dabest_hedges, palette = pal, rawplot.ylabel = "pg oxytocin/µg tissue protein/ml")

p_tmp<- c(t.test(df$value[df$group == 'Tph2+/+'], df$value[df$group == 'Tph2-/-'],var.equal=TRUE)$p.value)

dabest_hedges$result %>% mutate(p = p_tmp)

ggsave('assets/figure/OXY_DR.svg', plot = Fig3C, 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.