---
title: "Image_analysis_Astrocytes"
output: html_document
---
```{r}
library(tidyverse)
library(ggforce)
library(ggsignif)
library(shiny)
library(ggpubr)
```


```{r}
# #setwd("P:/AxPD_Project/Results/ImageAnalysis/Yokogawa_results/GFAP_S100b_TUJ1/AllDataSet_Harmonized") #path to the files were the curate data is. 
# 
# #upload all datasets
# data9 <- readr::read_csv("AxPD9_Astrocytes_HarmonizedNaming.csv")
# data10 <- readr::read_csv("AxPD10_Astrocytes_HarmonizedNaming.csv")
# data11 <- readr::read_csv("AxPD11_Astrocytes_HarmonizedNaming.csv")
# data12 <- readr::read_csv("AxPD12_Astrocytes_HarmonizedNaming.csv")

```
```{r}
# data <- rbind(data9, data10, data11, data12) #merge datasets
# data %>% 
#   mutate(Cellline = case_when  (str_detect(AreaName, '375') ~ 'WT18075',
#                                str_detect(AreaName, '278') ~ 'WT580',
#                                str_detect(AreaName, '280') ~ 'WT581',
#                                str_detect(AreaName, '376') ~ 'Miro1R272Q',
#                                str_detect(AreaName, '377') ~ 'GC Miro1R272Q',
#                                TRUE ~ 'none'))%>%
#   
#   mutate(Condition = case_when (str_detect(AreaName, '375') ~ 'WT',
#                                str_detect(AreaName, '278') ~ 'WT',
#                                str_detect(AreaName, '280') ~ 'WT',
#                                str_detect(AreaName, '376') ~ 'PD',
#                                str_detect(AreaName, '377') ~ 'GC',
#                                 TRUE ~ 'none'))%>%
#   
#  mutate(Day = case_when (str_detect(AreaName, 'D20') ~ '20d',
#                          str_detect(AreaName, 'D30') ~ '30d',
#                          str_detect(AreaName, 'D60') ~ '60d',
#                          str_detect(AreaName, 'D90') ~ '90d',
#                          str_detect(AreaName, 'D120') ~ '120d',
#                          TRUE ~ 'none'))%>%
#  
#    mutate(Batch = case_when (str_detect(AreaName, 'B9') ~ 'AxPD9',
#                          str_detect(AreaName, 'B10') ~ 'AxPD10',
#                          str_detect(AreaName, 'B11') ~ 'AxPD11',
#                          str_detect(AreaName, 'B12') ~ 'AxPD12',
#                          TRUE ~ 'none'))%>%
# 
#   mutate(Replicate = case_when (str_detect(AreaName, 'R1') ~ 'R1',
#                                str_detect(AreaName, 'R2') ~ 'R2',
#                                TRUE ~ 'none'))-> Ctrldata
# 
# write.csv(Ctrldata, "C:\\Users\\claudia.saraiva\\Desktop\\Manuscript2\\Figure1\\Partial\\C\\staining_20to120d.csv", row.names=FALSE)

```
              
```{r}
data3c <- readr::read_csv("staining_20to90d_WT580&PD_final.csv") #Conditions without WT18075 & WT581

#subset data
#20 days
subset(data3c, Day == '20') -> data20d_WTPD
write.csv(data20d, "C:\\Users\\claudia.saraiva\\Desktop\\Manuscript2\\Figure1\\Partial\\A\\data20d.csv", row.names=FALSE)

#30 days
subset(data3c, Day == '30') -> data30d_WTPD
write.csv(data30d, "C:\\Users\\claudia.saraiva\\Desktop\\Manuscript2\\Figure1\\Partial\\A\\data30d.csv", row.names=FALSE)

#60 days
subset(data3c, Day == '60') -> data60d_WTPD
write.csv(data60d, "C:\\Users\\claudia.saraiva\\Desktop\\Manuscript2\\Figure1\\Partial\\A\\data60d.csv", row.names=FALSE)

#90 days
subset(data3c, Day == '90') -> data90d_WTPD
write.csv(data90d, "C:\\Users\\claudia.saraiva\\Desktop\\Manuscript2\\Figure1\\Partial\\A\\data90d.csv", row.names=FALSE)


#statistics:

#Normality test
#20days
shapiro.test(data20d_WTPD$GFAPoverlapS100betaByNuclei) # p-value = 0.02332 not normal
shapiro.test(data20d_WTPD$GFAPMeanInGFAP) # p-value = 0.0013 not normal
shapiro.test(data20d_WTPD$S100betaMeanInS100beta) #  p-value = 0.971 Normal

#30days
shapiro.test(data30d_WTPD$GFAPoverlapS100betaByNuclei) #p-value = 2.226e-05 not normal
shapiro.test(data30d_WTPD$GFAPMeanInGFAP) #  p-value = 0.6721 normal
shapiro.test(data30d_WTPD$S100betaMeanInS100beta) # p-value = 0.0006246 not Normal

#60days
shapiro.test(data60d_WTPD$GFAPoverlapS100betaByNuclei) # p-value = 1.397e-05 not normal
shapiro.test(data60d_WTPD$GFAPMeanInGFAP) #   p-value = 0.1489 normal
shapiro.test(data60d_WTPD$S100betaMeanInS100beta) #p-value = 0.01294 not Normal

#90days
shapiro.test(data90d_WTPD$GFAPoverlapS100betaByNuclei) #  p-value = 0.3782 normal
shapiro.test(data90d_WTPD$GFAPMeanInGFAP) # p-value = 0.8795 normal
shapiro.test(data90d_WTPD$S100betaMeanInS100beta) #p-value = 0.005229 not Normal

#Use t-test for all samples

#If shapiro test is significant, then data is not normal and Wilcoxon signed rank test is done


# #20days
# kruskal.test(GFAPoverlapS100betaByNuclei ~ Condition, data = data20d)
# stat.test.dunn.20d.1  <- dunn_test(data20d, GFAPoverlapS100betaByNuclei ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.20d.1
# 
# kruskal.test(GFAPMeanInGFAP ~ Condition, data = data20d)
# stat.test.dunn.20d.2  <- dunn_test(data20d, GFAPMeanInGFAP ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.20d.2
# 
# kruskal.test(S100betaMeanInS100beta ~ Condition, data = data20d)
# stat.test.dunn.20d.3  <- dunn_test(data20d, S100betaMeanInS100beta ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.20d.3
# 
# #30days
# kruskal.test(GFAPoverlapS100betaByNuclei ~ Condition, data = data30d)
# stat.test.dunn.30d.1  <- dunn_test(data30d, GFAPoverlapS100betaByNuclei ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.30d.1
# 
# kruskal.test(GFAPMeanInGFAP ~ Condition, data = data30d)
# stat.test.dunn.30d.2  <- dunn_test(data30d, GFAPMeanInGFAP ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.30d.2
# 
# kruskal.test(S100betaMeanInS100beta ~ Condition, data = data30d)
# stat.test.dunn.30d.3  <- dunn_test(data30d, S100betaMeanInS100beta ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.30d.3
# 
# #60days
# kruskal.test(GFAPoverlapS100betaByNuclei ~ Condition, data = data60d)
# stat.test.dunn.60d.1  <- dunn_test(data60d, GFAPoverlapS100betaByNuclei ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.60d.1
# 
# kruskal.test(GFAPMeanInGFAP ~ Condition, data = data60d)
# stat.test.dunn.60d.2  <- dunn_test(data60d, GFAPMeanInGFAP ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.60d.2
# 
# kruskal.test(S100betaMeanInS100beta ~ Condition, data = data60d)
# stat.test.dunn.60d.3  <- dunn_test(data60d, S100betaMeanInS100beta ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.60d.3
# 
# #90days
# kruskal.test(GFAPoverlapS100betaByNuclei ~ Condition, data = data90d)
# stat.test.dunn.90d.1  <- dunn_test(data90d, GFAPoverlapS100betaByNuclei ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.90d.1
# 
# kruskal.test(GFAPMeanInGFAP ~ Condition, data = data90d)
# stat.test.dunn.90d.2  <- dunn_test(data90d, GFAPMeanInGFAP ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.90d.2
# 
# kruskal.test(S100betaMeanInS100beta ~ Condition, data = data90d)
# stat.test.dunn.90d.3  <- dunn_test(data90d, S100betaMeanInS100beta ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.90d.3
# 
# #120days
# kruskal.test(GFAPoverlapS100betaByNuclei ~ Condition, data = data120d)
# stat.test.dunn.120d.1  <- dunn_test(data120d, GFAPoverlapS100betaByNuclei ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.120d.1
# 
# kruskal.test(GFAPMeanInGFAP ~ Condition, data = data120d)
# stat.test.dunn.120d.2  <- dunn_test(data120d, GFAPMeanInGFAP ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.120d.2
# 
# kruskal.test(S100betaMeanInS100beta ~ Condition, data = data120d)
# stat.test.dunn.120d.3  <- dunn_test(data120d, S100betaMeanInS100beta ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.120d.3
# 
# # merge statistics in 1 dataframe
# 
# stat.test.dunn.20d.1$Day <- c("20","20","20")
# stat.test.dunn.30d.1$Day <- c("30","30","30")
# stat.test.dunn.60d.1$Day <- c("60","60","60")
# stat.test.dunn.90d.1$Day <- c("90","90","90")
# #stat.test.dunn.120d.1$Day <- c("120","120","120")
# 
# stat.test.dunn.20d.2$Day <- c("20","20","20")
# stat.test.dunn.30d.2$Day <- c("30","30","30")
# stat.test.dunn.60d.2$Day <- c("60","60","60")
# stat.test.dunn.90d.2$Day <- c("90","90","90")
# #stat.test.dunn.120d.2$Day <- c("120","120","120")
# 
# stat.test.dunn.20d.3$Day <- c("20","20","20")
# stat.test.dunn.30d.3$Day <- c("30","30","30")
# stat.test.dunn.60d.3$Day <- c("60","60","60")
# stat.test.dunn.90d.3$Day <- c("90","90","90")
# #stat.test.dunn.120d.3$Day <- c("120","120","120")
# 
# stat.test.dunn.all.1 <- rbind.data.frame(stat.test.dunn.20d.1, stat.test.dunn.30d.1, stat.test.dunn.60d.1, stat.test.dunn.90d.1) #, stat.test.dunn.120d.1)
# 
# stat.test.dunn.all.2 <- rbind.data.frame(stat.test.dunn.20d.2, stat.test.dunn.30d.2, stat.test.dunn.60d.2, stat.test.dunn.90d.2)#, stat.test.dunn.120d.2)
# 
# stat.test.dunn.all.3 <- rbind.data.frame(stat.test.dunn.20d.3, stat.test.dunn.30d.3, stat.test.dunn.60d.3, stat.test.dunn.90d.3)#, stat.test.dunn.120d.3)

```

```{r}
#order the different parameters in a specific way
data3c$Condition <- factor (data3c$Condition, level=c("WT", "PD"))
data3c$Day <- factor (data3c$Day, level=c("20","30", "60", "90"))

data3c %>%
    ggplot(aes(x = factor(Condition), y=GFAPoverlapS100betaByNuclei) ) +
    geom_boxplot(aes(fill=Condition),width=0.5,outlier.shape=NA) +
      scale_fill_manual(values= c(alpha("skyblue4",0.7),"red4"),name = "Condition") +
    geom_jitter(, size = 1, colour= alpha("black",0.7), position=position_jitter(0), cex=2)+
  # #stat_compare_means()  + #default Kruskall-Wallis
  # stat_pvalue_manual(
  #   stat.test.dunn.all.1,
  #   #y.position = c(0.37, 0.4,0.5, 0.53, 0.33, 0.36, 0.33, 0.36)  , step.increase = 0,
  #   y.position = c(0.5, 0.53, 0.5, 0.53, 0.5, 0.53), step.increase = 0,
  #   label = "p.adj.signif", hide.ns =T, , tip.length = 0
  #   )+
      
  facet_grid(~fct_relevel(Day, "20", "30",  "60",  "90"), scales = "free") +
  labs(x     = "WT", "PD", "GC",
       y     = "GFAP&S100b / Nuclei (voxel ratio)",
       fill  = "Conditions",
       title = "") +
    # scale_y_continuous(limits=c(0, 0.7)) +
   ggpubr::stat_compare_means(comparisons = list(c("WT", "PD")), test='wilcox.test', label="p.signif", vjust=0.2)+
  
  
  theme_bw() +
  theme(
    axis.line = element_line(colour = 'black', size = 0.1) ,
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_text(size = 14),
    axis.text.y = element_text(size=12, color="black"),
    axis.ticks.y = element_line(),
    legend.position="right",
    legend.text = element_text(size = 12, face = "bold") ,
    legend.title = element_blank(),
    legend.key.size = unit(0.7, "cm"),
    legend.key.width = unit(0.6,"cm"),
    plot.title = element_text(size = 15, hjust=0.5, vjust= 1, face = "bold"),
    plot.subtitle = element_blank(),#element_text(size = 2, hjust=0.5)
    strip.text = element_text(size=12, vjust=0.5),
    strip.background = element_rect(fill="lightgray"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.spacing.y = unit(0.8, "lines"),
    strip.switch.pad.wrap=unit(20, "lines"),
    # panel.grid.major = element_blank(), 
    # #panel.grid.minor = element_blank(),
    # strip.background = element_rect(fill="lightgray"),
    # panel.border = element_rect(fill = NA, color = "black"),
    # panel.spacing.y = unit(0.8, "lines"),
    # strip.switch.pad.wrap=unit(20, "lines")
  ) -> p
 t<- cowplot::ggdraw(cowplot::add_sub(p, "Wilcox-test, *p<0.05",hjust=-1.3, size=10))
#+
#  theme(aspect.ratio = 1.25 / 0.75)
print(t)
ggsave("GFAP&S100bbyNuc_20t090days_WTPD.pdf", useDingbats=FALSE, plot=p, height = 4, width = 7)
ggsave("GFAP&S100bbyNuc_20t090days_WTPD.tiff", plot=p)

```

```{r}
#order the different parameters in a specific way
data3c$Condition <- factor (data3c$Condition, level=c("WT", "PD"))
data3c$Day <- factor (data3c$Day, level=c("20","30", "60", "90"))

data3c %>%
    ggplot(aes(x = factor(Condition), y=GFAPMeanInGFAP) ) +
    geom_boxplot(aes(fill=Condition),width=0.5,outlier.shape=NA) +
      scale_fill_manual(values= c(alpha("skyblue4",0.7),"red4"),name = "Condition") +
    geom_jitter(, size = 1, colour= alpha("black",0.7), position=position_jitter(0), cex=2)+
  # #stat_compare_means()  + #default Kruskall-Wallis
  # stat_pvalue_manual(
  #   stat.test.dunn.all.2,
  #   #y.position = c(0.37, 0.4,0.5, 0.53, 0.33, 0.36, 0.33, 0.36)  , step.increase = 0,
  #   y.position = c(7200, 7500, 7200, 7500), step.increase = 0,
  #   label = "p.adj.signif", hide.ns =T, , tip.length = 0
  #   )+
      
  facet_grid(~fct_relevel(Day, "20", "30",  "60",  "90"), scales = "free") +
  labs(x     = "WT", "PD", "GC",
       y     = "GFAP (MFI)",
       fill  = "Conditions",
       title = "") +
     #scale_y_continuous(limits=c(0, 8000)) +
  ggpubr::stat_compare_means(comparisons = list(c("WT", "PD")), test='wilcox.test', label="p.signif", vjust=0.2)+
  
  theme_bw() +
  theme(
    axis.line = element_line(colour = 'black', size = 0.1) ,
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_text(size = 14),
    axis.text.y = element_text(size=12, color="black"),
    axis.ticks.y = element_line(),
    legend.position="right",
    legend.text = element_text(size = 12, face = "bold") ,
    legend.title = element_blank(),
    legend.key.size = unit(0.7, "cm"),
    legend.key.width = unit(0.6,"cm"),
    plot.title = element_text(size = 15, hjust=0.5, vjust= 1, face = "bold"),
    plot.subtitle = element_blank(),#element_text(size = 2, hjust=0.5)
    strip.text = element_text(size=12, vjust=0.5),
    strip.background = element_rect(fill="lightgray"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.spacing.y = unit(0.8, "lines"),
    strip.switch.pad.wrap=unit(20, "lines"),
    # panel.grid.major = element_blank(), 
    # #panel.grid.minor = element_blank(),
    # strip.background = element_rect(fill="lightgray"),
    # panel.border = element_rect(fill = NA, color = "black"),
    # panel.spacing.y = unit(0.8, "lines"),
    # strip.switch.pad.wrap=unit(20, "lines")
  ) -> p
 t<- cowplot::ggdraw(cowplot::add_sub(p, "Kruskal-Wallis, **p<0.01, *p<0.05",hjust=-0.4, size=10))
#+
#  theme(aspect.ratio = 1.25 / 0.75)
print(t)
ggsave("GFAP-MFI_stats.pdf", useDingbats=FALSE, plot=p, height = 4, width = 7)
#ggsave("GFAP-MFI_20t0120days.tiff", plot=p)
```

```{r}
#order the different parameters in a specific way
data3c$Condition <- factor (data3c$Condition, level=c("WT", "PD", "GC"))
data3c$Day <- factor (data3c$Day, level=c("20","30", "60", "90", "120"))

data3c %>%
    ggplot(aes(x = factor(Condition), y=S100betaMeanInS100beta) ) +
    geom_boxplot(aes(fill=Condition),width=0.5,outlier.shape=NA) +
      scale_fill_manual(values= c(alpha("skyblue4",0.7),"red4",alpha("coral",0.7)),name = "Condition") +
    geom_jitter(, size = 1, colour= alpha("black",0.7), position=position_jitter(0), cex=2)+
  #stat_compare_means()  + #default Kruskall-Wallis
  stat_pvalue_manual(
    stat.test.dunn.all.3,
    #y.position = c(0.37, 0.4,0.5, 0.53, 0.33, 0.36, 0.33, 0.36)  , step.increase = 0,
    y.position = c(7200, 7200, 7200, 7500, 7200), step.increase = 0,
    label = "p.adj.signif", hide.ns =T, , tip.length = 0
    )+
      
  facet_grid(~fct_relevel(Day, "20", "30",  "60",  "90", "120"), scales = "free") +
  labs(x     = "WT", "PD", "GC",
       y     = "S100b (MFI)",
       fill  = "Conditions",
       title = "") +
     scale_y_continuous(limits=c(0, 8000)) +
  
  
  theme_bw() +
  theme(
    axis.line = element_line(colour = 'black', size = 0.1) ,
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_text(size = 14),
    axis.text.y = element_text(size=12, color="black"),
    axis.ticks.y = element_line(),
    legend.position="right",
    legend.text = element_text(size = 12, face = "bold") ,
    legend.title = element_blank(),
    legend.key.size = unit(0.7, "cm"),
    legend.key.width = unit(0.6,"cm"),
    plot.title = element_text(size = 15, hjust=0.5, vjust= 1, face = "bold"),
    plot.subtitle = element_blank(),#element_text(size = 2, hjust=0.5)
    strip.text = element_text(size=12, vjust=0.5),
    strip.background = element_rect(fill="lightgray"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.spacing.y = unit(0.8, "lines"),
    strip.switch.pad.wrap=unit(20, "lines"),
    # panel.grid.major = element_blank(), 
    # #panel.grid.minor = element_blank(),
    # strip.background = element_rect(fill="lightgray"),
    # panel.border = element_rect(fill = NA, color = "black"),
    # panel.spacing.y = unit(0.8, "lines"),
    # strip.switch.pad.wrap=unit(20, "lines")
  ) -> p
 t<- cowplot::ggdraw(cowplot::add_sub(p, "Kruskal-Wallis, **p<0.01, *p<0.05",hjust=-0.4, size=10))
#+
#  theme(aspect.ratio = 1.25 / 0.75)
print(t)
ggsave("S100b-MFI_20t0120days.pdf", useDingbats=FALSE, plot=p, height = 4, width = 7)
ggsave("S100b-MFI_20t0120days.tiff", plot=p)
```



```{r}
data3c %>%
  ggplot(aes(x=Day, y=GFAPMeanInGFAP, group=Condition, color=Condition)) +
  stat_smooth(method='loess', formula ='y ~ x', size = 1)+
  scale_colour_manual(name="Condition", values= c("WT" = alpha("skyblue4",0.7), "PD" = "red4")) +
  # facet_wrap(~timepoint, scales="free", labeller = labeller(groupwrap = label_wrap_gen(10))) +
  # facet_grid(~fct_relevel(timepoint, "day15","day30",  "day70", "day90"), scales="free") +
  labs(x="Time (days)",
        y     = "GFAP (MFI)",
       fill  = "Condition",
       title = "") +
 scale_y_continuous(limits=c(0, 6500)) +
   theme_bw() +
     theme(
    axis.line = element_line(colour = 'black', size = 0.1) ,
    axis.title.x = element_blank(),
    axis.text.x = element_text(),
    axis.title.y = element_text(size = 14),
    axis.text.y = element_text(size=12, color="black"),
    axis.ticks.y = element_line(),
    legend.position="right",
    legend.text = element_text(size = 12, face = "bold") ,
    legend.title = element_blank(),
    legend.key.size = unit(0.7, "cm"),
    legend.key.width = unit(0.6,"cm"),
    plot.title = element_text(size = 15, hjust=0.5, vjust= 1, face = "bold"),
    plot.subtitle = element_blank(),#element_text(size = 2, hjust=0.5)
    strip.text = element_text(size=12, vjust=0.5),
    strip.background = element_rect(fill="lightgray"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.spacing.y = unit(0.8, "lines"),
    strip.switch.pad.wrap=unit(20, "lines"),
   ) -> p
  t<- cowplot::ggdraw(cowplot::add_sub(p, "Locally estimated scatterplot smoothing: y ~ Time \nWilcox-test, **p<0.01, *p<0.05",hjust=-0, size = 9))
  print(t)
ggsave("GFAP-MFI_timeline_WTPD.pdf", useDingbats=FALSE, plot=t, height = 4, width = 7)
ggsave("GFAP-MFI_timeline_WTPD.tiff", plot=p)

```

```{r}
data3c %>%
  ggplot(aes(x=Day, y=S100betaMeanInS100beta, group=Condition, color=Condition)) +
  stat_smooth(method='loess', formula ='y ~ x', size = 1)+
  scale_colour_manual(name="Condition", values= c("WT" = alpha("skyblue4",0.7), "PD" = "red4")) +
  # facet_wrap(~timepoint, scales="free", labeller = labeller(groupwrap = label_wrap_gen(10))) +
  # facet_grid(~fct_relevel(timepoint, "day15","day30",  "day70", "day90"), scales="free") +
  labs(x="Time (days)",
        y     = "S100b (MFI)",
       fill  = "Condition",
       title = "") +
  scale_y_continuous(limits=c(0, 4100)) +
   theme_bw() +
     theme(
    axis.line = element_line(colour = 'black', size = 0.1) ,
    axis.title.x = element_blank(),
    axis.text.x = element_text(),
    axis.title.y = element_text(size = 14),
    axis.text.y = element_text(size=12, color="black"),
    axis.ticks.y = element_line(),
    legend.position="right",
    legend.text = element_text(size = 12, face = "bold") ,
    legend.title = element_blank(),
    legend.key.size = unit(0.7, "cm"),
    legend.key.width = unit(0.6,"cm"),
    plot.title = element_text(size = 15, hjust=0.5, vjust= 1, face = "bold"),
    plot.subtitle = element_blank(),#element_text(size = 2, hjust=0.5)
    strip.text = element_text(size=12, vjust=0.5),
    strip.background = element_rect(fill="lightgray"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.spacing.y = unit(0.8, "lines"),
    strip.switch.pad.wrap=unit(20, "lines"),
   ) -> p
  t<- cowplot::ggdraw(cowplot::add_sub(p, "Locally estimated scatterplot smoothing: y ~ Time \nWilcox-test, ***p<0.001, *p<0.05",hjust=-0, size = 9))
  print(t)
 
 ggsave("S100B-MFI_timeline_WTPD.pdf", useDingbats=FALSE, plot=t, height = 4, width = 7)
ggsave("S100B-MFI_timeline_WTPD.tiff", plot=p)
```


```{r warning=FALSE}

#order the different parameters in a specific way
data3c$Condition <- factor (data3c$Condition, level=c("WT", "PD"))
data3c$Day <- factor (data3c$Day, level=c("20","30", "60", "90"))

data3c %>%
    ggplot(aes(x = factor(Condition), y=S100betaMeanInS100beta) ) +
    geom_boxplot(aes(fill=Condition),width=0.5,outlier.shape=NA) +
      scale_fill_manual(values= c(alpha("skyblue4",0.7),"red4"),name = "Condition") +
    geom_jitter(, size = 1, colour= alpha("black",0.7), position=position_jitter(0), cex=2)+
  # #stat_compare_means()  + #default Kruskall-Wallis
  # stat_pvalue_manual(
  #   stat.test.dunn.all.2,
  #   #y.position = c(0.37, 0.4,0.5, 0.53, 0.33, 0.36, 0.33, 0.36)  , step.increase = 0,
  #   y.position = c(7200, 7500, 7200, 7500), step.increase = 0,
  #   label = "p.adj.signif", hide.ns =T, , tip.length = 0
  #   )+
      
  facet_grid(~fct_relevel(Day, "20", "30",  "60",  "90"), scales = "free") +
  labs(x     = "WT", "PD", "GC",
       y     = "S100B (MFI)",
       fill  = "Conditions",
       title = "") +
     #scale_y_continuous(limits=c(0, 8000)) +
  ggpubr::stat_compare_means(comparisons = list(c("WT", "PD")), test='wilcox.test', label="p.signif", vjust=0.2)+
  
  theme_bw() +
  theme(
    axis.line = element_line(colour = 'black', size = 0.1) ,
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_text(size = 14),
    axis.text.y = element_text(size=12, color="black"),
    axis.ticks.y = element_line(),
    legend.position="right",
    legend.text = element_text(size = 12, face = "bold") ,
    legend.title = element_blank(),
    legend.key.size = unit(0.7, "cm"),
    legend.key.width = unit(0.6,"cm"),
    plot.title = element_text(size = 15, hjust=0.5, vjust= 1, face = "bold"),
    plot.subtitle = element_blank(),#element_text(size = 2, hjust=0.5)
    strip.text = element_text(size=12, vjust=0.5),
    strip.background = element_rect(fill="lightgray"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.spacing.y = unit(0.8, "lines"),
    strip.switch.pad.wrap=unit(20, "lines"),
    # panel.grid.major = element_blank(), 
    # #panel.grid.minor = element_blank(),
    # strip.background = element_rect(fill="lightgray"),
    # panel.border = element_rect(fill = NA, color = "black"),
    # panel.spacing.y = unit(0.8, "lines"),
    # strip.switch.pad.wrap=unit(20, "lines")
  ) -> p
 t<- cowplot::ggdraw(cowplot::add_sub(p, "Kruskal-Wallis, **p<0.01, *p<0.05",hjust=-0.4, size=10))
#+
#  theme(aspect.ratio = 1.25 / 0.75)
print(t)
ggsave("S100B-MFI_stats.pdf", useDingbats=FALSE, plot=p, height = 4, width = 7)
#ggsave("GFAP-MFI_20t0120days.tiff", plot=p)
  
```
```{r}
#order the different parameters in a specific way
data4c <- rbind(data20d_WTPD, data30d_WTPD)
data4c$Condition <- factor (data4c$Condition, level=c("WT", "PD"))
data4c$Day <- factor (data4c$Day, level=c("20","30"))

data4c %>%
    ggplot(aes(x = factor(Condition), y=GFAPoverlapS100betaByNuclei) ) +
    geom_boxplot(aes(fill=Condition),width=0.5,outlier.shape=NA) +
      scale_fill_manual(values= c(alpha("skyblue4",0.7),"red4"),name = "Condition") +
    geom_jitter(, size = 1, colour= alpha("black",0.7), position=position_jitter(0), cex=2)+
  # #stat_compare_means()  + #default Kruskall-Wallis
  # stat_pvalue_manual(
  #   stat.test.dunn.all.1,
  #   #y.position = c(0.37, 0.4,0.5, 0.53, 0.33, 0.36, 0.33, 0.36)  , step.increase = 0,
  #   y.position = c(0.5, 0.53, 0.5, 0.53, 0.5, 0.53), step.increase = 0,
  #   label = "p.adj.signif", hide.ns =T, , tip.length = 0
  #   )+
      
  facet_grid(~fct_relevel(Day, "20", "30",  "60",  "90"), scales = "free") +
  labs(x     = "WT", "PD", "GC",
       y     = "GFAP&S100b / Nuclei (voxel ratio)",
       fill  = "Conditions",
       title = "") +
    # scale_y_continuous(limits=c(0, 0.7)) +
   ggpubr::stat_compare_means(comparisons = list(c("WT", "PD")), test='wilcox.test', label="p.signif", vjust=0.2)+
  
  
  theme_bw() +
  theme(
    axis.line = element_line(colour = 'black', size = 0.1) ,
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_text(size = 14),
    axis.text.y = element_text(size=12, color="black"),
    axis.ticks.y = element_line(),
    legend.position="right",
    legend.text = element_text(size = 12, face = "bold") ,
    legend.title = element_blank(),
    legend.key.size = unit(0.7, "cm"),
    legend.key.width = unit(0.6,"cm"),
    plot.title = element_text(size = 15, hjust=0.5, vjust= 1, face = "bold"),
    plot.subtitle = element_blank(),#element_text(size = 2, hjust=0.5)
    strip.text = element_text(size=12, vjust=0.5),
    strip.background = element_rect(fill="lightgray"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.spacing.y = unit(0.8, "lines"),
    strip.switch.pad.wrap=unit(20, "lines"),
    # panel.grid.major = element_blank(), 
    # #panel.grid.minor = element_blank(),
    # strip.background = element_rect(fill="lightgray"),
    # panel.border = element_rect(fill = NA, color = "black"),
    # panel.spacing.y = unit(0.8, "lines"),
    # strip.switch.pad.wrap=unit(20, "lines")
  ) -> p
 t<- cowplot::ggdraw(cowplot::add_sub(p, "Wilcox-test, *p<0.05",hjust=-1.3, size=10))
#+
#  theme(aspect.ratio = 1.25 / 0.75)
print(t)
ggsave("GFAP&S100bbyNuc_20&30days_WTPD.pdf", useDingbats=FALSE, plot=p, height = 4, width = 4.2)
ggsave("GFAP&S100bbyNuc_20&30days_WTPD.tiff", plot=p)
```

```{r}
data5c <- rbind(data60d_WTPD, data90d_WTPD)
#order the different parameters in a specific way
data5c$Condition <- factor (data5c$Condition, level=c("WT", "PD"))
data5c$Day <- factor (data5c$Day, level=c("60", "90"))

data5c %>%
    ggplot(aes(x = factor(Condition), y=GFAPoverlapS100betaByNuclei) ) +
    geom_boxplot(aes(fill=Condition),width=0.5,outlier.shape=NA) +
      scale_fill_manual(values= c(alpha("skyblue4",0.7),"red4"),name = "Condition") +
    geom_jitter(, size = 1, colour= alpha("black",0.7), position=position_jitter(0), cex=2)+
  # #stat_compare_means()  + #default Kruskall-Wallis
  # stat_pvalue_manual(
  #   stat.test.dunn.all.1,
  #   #y.position = c(0.37, 0.4,0.5, 0.53, 0.33, 0.36, 0.33, 0.36)  , step.increase = 0,
  #   y.position = c(0.5, 0.53, 0.5, 0.53, 0.5, 0.53), step.increase = 0,
  #   label = "p.adj.signif", hide.ns =T, , tip.length = 0
  #   )+
      
  facet_grid(~fct_relevel(Day, "20", "30",  "60",  "90"), scales = "free") +
  labs(x     = "WT", "PD", "GC",
       y     = "GFAP&S100b / Nuclei (voxel ratio)",
       fill  = "Conditions",
       title = "") +
    # scale_y_continuous(limits=c(0, 0.7)) +
   ggpubr::stat_compare_means(comparisons = list(c("WT", "PD")), test='wilcox.test', label="p.signif", vjust=0.2)+
  
  
  theme_bw() +
  theme(
    axis.line = element_line(colour = 'black', size = 0.1) ,
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_text(size = 14),
    axis.text.y = element_text(size=12, color="black"),
    axis.ticks.y = element_line(),
    legend.position="right",
    legend.text = element_text(size = 12, face = "bold") ,
    legend.title = element_blank(),
    legend.key.size = unit(0.7, "cm"),
    legend.key.width = unit(0.6,"cm"),
    plot.title = element_text(size = 15, hjust=0.5, vjust= 1, face = "bold"),
    plot.subtitle = element_blank(),#element_text(size = 2, hjust=0.5)
    strip.text = element_text(size=12, vjust=0.5),
    strip.background = element_rect(fill="lightgray"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.spacing.y = unit(0.8, "lines"),
    strip.switch.pad.wrap=unit(20, "lines"),
    # panel.grid.major = element_blank(), 
    # #panel.grid.minor = element_blank(),
    # strip.background = element_rect(fill="lightgray"),
    # panel.border = element_rect(fill = NA, color = "black"),
    # panel.spacing.y = unit(0.8, "lines"),
    # strip.switch.pad.wrap=unit(20, "lines")
  ) -> p
 t<- cowplot::ggdraw(cowplot::add_sub(p, "Wilcox-test, *p<0.05",hjust=-1.3, size=10))
#+
#  theme(aspect.ratio = 1.25 / 0.75)
print(t)
ggsave("GFAP&S100bbyNuc_60&90days_WTPD.pdf", useDingbats=FALSE, plot=p, height = 4, width = 4.2)
ggsave("GFAP&S100bbyNuc_60&90days_WTPD.tiff", plot=p)
```
```{r}
#order the different parameters in a specific way WT, PD, GC
data4c <- rbind(data20d_WTPD, data30d_WTPD)
data4c$Condition <- factor (data4c$Condition, level=c("WT", "PD", "GC"))
data4c$Day <- factor (data4c$Day, level=c("20","30"))

data4c %>%
    ggplot(aes(x = factor(Condition), y=GFAPoverlapS100betaByNuclei) ) +
    geom_boxplot(aes(fill=Condition),width=0.5,outlier.shape=NA) +
      scale_fill_manual(values= c(alpha("skyblue4",0.7),"red4",alpha("coral",0.7)),name = "Condition") +
    geom_jitter(, size = 1, colour= alpha("black",0.7), position=position_jitter(0), cex=2)+
 # stat_compare_means()  + #default Kruskall-Wallis
  stat_pvalue_manual(
  stat.test.dunn.all.1a,
   #y.position = c(0.37, 0.4,0.5, 0.53, 0.33, 0.36, 0.33, 0.36)  , step.increase = 0,
     y.position = c(0.0091, 0.0095), step.increase = 0,
    label = "p.adj.signif", hide.ns =T, , tip.length = 0
    )+
      
  facet_grid(~fct_relevel(Day, "20", "30"), scales = "free") +
  labs(x     = "WT", "PD", "GC",
       y     = "GFAP&S100b / Nuclei (voxel ratio)",
       fill  = "Conditions",
       title = "") +
    scale_y_continuous(limits=c(0, 0.01)) +
   #ggpubr::stat_compare_means(comparisons = list(c("WT", "PD")), test='wilcox.test', label="p.signif", vjust=0.2)+
  
  
  theme_bw() +
  theme(
    axis.line = element_line(colour = 'black', size = 0.1) ,
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_text(size = 14),
    axis.text.y = element_text(size=12, color="black"),
    axis.ticks.y = element_line(),
    legend.position="right",
    legend.text = element_text(size = 12, face = "bold") ,
    legend.title = element_blank(),
    legend.key.size = unit(0.7, "cm"),
    legend.key.width = unit(0.6,"cm"),
    plot.title = element_text(size = 15, hjust=0.5, vjust= 1, face = "bold"),
    plot.subtitle = element_blank(),#element_text(size = 2, hjust=0.5)
    strip.text = element_text(size=12, vjust=0.5),
    strip.background = element_rect(fill="lightgray"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.spacing.y = unit(0.8, "lines"),
    strip.switch.pad.wrap=unit(20, "lines"),
    # panel.grid.major = element_blank(), 
    # #panel.grid.minor = element_blank(),
    # strip.background = element_rect(fill="lightgray"),
    # panel.border = element_rect(fill = NA, color = "black"),
    # panel.spacing.y = unit(0.8, "lines"),
    # strip.switch.pad.wrap=unit(20, "lines")
  ) -> p
 t<- cowplot::ggdraw(cowplot::add_sub(p, "Wilcox-test, *p<0.05",hjust=-1.3, size=10))
#+
#  theme(aspect.ratio = 1.25 / 0.75)
print(t)
ggsave("GFAP&S100bbyNuc_20&30days_WTPDGC.pdf", useDingbats=FALSE, plot=p, height = 4, width = 4.2)
ggsave("GFAP&S100bbyNuc_20&30days_WTPDGC.tiff", plot=p)
```

```{r}
data5c <- rbind(data60d_WTPD, data90d_WTPD)
#order the different parameters in a specific way
data5c$Condition <- factor (data5c$Condition, level=c("WT", "PD", "GC"))
data5c$Day <- factor (data5c$Day, level=c("60", "90"))

data5c %>%
   ggplot(aes(x = factor(Condition), y=GFAPoverlapS100betaByNuclei) ) +
    geom_boxplot(aes(fill=Condition),width=0.5,outlier.shape=NA) +
      scale_fill_manual(values= c(alpha("skyblue4",0.7),"red4",alpha("coral",0.7)),name = "Condition") +
    geom_jitter(, size = 1, colour= alpha("black",0.7), position=position_jitter(0), cex=2)+
 # stat_compare_means()  + #default Kruskall-Wallis
 stat_pvalue_manual(
 stat.test.dunn.all.1b,
  #y.position = c(0.37, 0.4,0.5, 0.53, 0.33, 0.36, 0.33, 0.36)  , step.increase = 0,
    y.position = c(0.525, 0.55), step.increase = 0,
   label = "p.adj.signif", hide.ns =T, , tip.length = 0
   )+
      
  facet_grid(~fct_relevel(Day, "60", "90"), scales = "free") +
  labs(x     = "WT", "PD", "GC",
       y     = "GFAP&S100b / Nuclei (voxel ratio)",
       fill  = "Conditions",
       title = "") +
    scale_y_continuous(limits=c(0, 0.6)) +
   #ggpubr::stat_compare_means(comparisons = list(c("WT", "PD")), test='wilcox.test', label="p.signif", vjust=0.2)+
  
  
  theme_bw() +
  theme(
    axis.line = element_line(colour = 'black', size = 0.1) ,
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_text(size = 14),
    axis.text.y = element_text(size=12, color="black"),
    axis.ticks.y = element_line(),
    legend.position="right",
    legend.text = element_text(size = 12, face = "bold") ,
    legend.title = element_blank(),
    legend.key.size = unit(0.7, "cm"),
    legend.key.width = unit(0.6,"cm"),
    plot.title = element_text(size = 15, hjust=0.5, vjust= 1, face = "bold"),
    plot.subtitle = element_blank(),#element_text(size = 2, hjust=0.5)
    strip.text = element_text(size=12, vjust=0.5),
    strip.background = element_rect(fill="lightgray"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.spacing.y = unit(0.8, "lines"),
    strip.switch.pad.wrap=unit(20, "lines"),
    # panel.grid.major = element_blank(), 
    # #panel.grid.minor = element_blank(),
    # strip.background = element_rect(fill="lightgray"),
    # panel.border = element_rect(fill = NA, color = "black"),
    # panel.spacing.y = unit(0.8, "lines"),
    # strip.switch.pad.wrap=unit(20, "lines")
  ) -> p
 t<- cowplot::ggdraw(cowplot::add_sub(p, " *p<0.05",hjust=-1.3, size=10))
#+
#  theme(aspect.ratio = 1.25 / 0.75)
print(t)
ggsave("GFAP&S100bbyNuc_60&90days_WTPDGC.pdf", useDingbats=FALSE, plot=p, height = 4, width = 4.2)
ggsave("GFAP&S100bbyNuc_60&90days_WTPDGC.tiff", plot=p)
```

```{r}
# path: C:\Users\claudia.saraiva\Desktop\Manuscript2\Figure1\Partial\C\allconditions

library(readr)
data3c <- read_csv("allconditions/staining_20to9d_WT580-PD-GC.csv")
View(data3c)

#data3c <- readr::read_csv("staining_20to90d_WT580-PD-GC.csv") #Conditions without WT18075 & WT581

#subset data
#20 days
subset(data3c, Day == '20') -> data20d_WTPD
write.csv(data20d_WTPD, "C:\\Users\\claudia.saraiva\\Desktop\\Manuscript2\\Figure1\\Partial\\C\\allconditions\\data20d.csv", row.names=FALSE)

#30 days
subset(data3c, Day == '30') -> data30d_WTPD
write.csv(data30d_WTPD, "C:\\Users\\claudia.saraiva\\Desktop\\Manuscript2\\Figure1\\Partial\\C\\allconditions\\data30d.csv", row.names=FALSE)

#60 days
subset(data3c, Day == '60') -> data60d_WTPD
write.csv(data60d_WTPD, "C:\\Users\\claudia.saraiva\\Desktop\\Manuscript2\\Figure1\\Partial\\C\\allconditions\\data60d.csv", row.names=FALSE)

#90 days
subset(data3c, Day == '90') -> data90d_WTPD
write.csv(data90d_WTPD, "C:\\Users\\claudia.saraiva\\Desktop\\Manuscript2\\Figure1\\Partial\\C\\allconditions\\data90d.csv", row.names=FALSE)


#statistics:

#Normality test
#20days
shapiro.test(data20d_WTPD$GFAPoverlapS100betaByNuclei) # p-value = 2.505e-07 not normal
shapiro.test(data20d_WTPD$GFAPMeanInGFAP) # p-value = 0.0125 not normal
shapiro.test(data20d_WTPD$S100betaMeanInS100beta) #  p-value = 0.4223 not Normal

#30days
shapiro.test(data30d_WTPD$GFAPoverlapS100betaByNuclei) #p-value = 4.108e-07 not normal
shapiro.test(data30d_WTPD$GFAPMeanInGFAP) #  p-value = 0.1543 normal
shapiro.test(data30d_WTPD$S100betaMeanInS100beta) # p-value = 2.633e-05 not Normal

#60days
shapiro.test(data60d_WTPD$GFAPoverlapS100betaByNuclei) # p-value = 0.001013 not normal
shapiro.test(data60d_WTPD$GFAPMeanInGFAP) #   p-value = 0.1489 normal
shapiro.test(data60d_WTPD$S100betaMeanInS100beta) #p-value = 0.01294 not Normal

#90days
shapiro.test(data90d_WTPD$GFAPoverlapS100betaByNuclei) #  p-value = 0.6388 normal
shapiro.test(data90d_WTPD$GFAPMeanInGFAP) # p-value = 0.83 normal
shapiro.test(data90d_WTPD$S100betaMeanInS100beta) #p-value = 0.004304 not Normal

#Use t-test for all samples

#If shapiro test is significant, then data is not normal and Wilcoxon signed rank test is done


#20days
kruskal.test(GFAPoverlapS100betaByNuclei ~ Condition, data = data20d_WTPD)
stat.test.dunn.20d.1  <- dunn_test(data20d_WTPD, GFAPoverlapS100betaByNuclei ~ Condition, p.adjust.method = "BH")
stat.test.dunn.20d.1
# 
# kruskal.test(GFAPMeanInGFAP ~ Condition, data = data20d)
# stat.test.dunn.20d.2  <- dunn_test(data20d, GFAPMeanInGFAP ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.20d.2
# 
# kruskal.test(S100betaMeanInS100beta ~ Condition, data = data20d)
# stat.test.dunn.20d.3  <- dunn_test(data20d, S100betaMeanInS100beta ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.20d.3
# 
#30days
kruskal.test(GFAPoverlapS100betaByNuclei ~ Condition, data = data30d_WTPD)
stat.test.dunn.30d.1  <- dunn_test(data30d_WTPD, GFAPoverlapS100betaByNuclei ~ Condition, p.adjust.method = "BH")
stat.test.dunn.30d.1
# 
# kruskal.test(GFAPMeanInGFAP ~ Condition, data = data30d)
# stat.test.dunn.30d.2  <- dunn_test(data30d, GFAPMeanInGFAP ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.30d.2
# 
# kruskal.test(S100betaMeanInS100beta ~ Condition, data = data30d)
# stat.test.dunn.30d.3  <- dunn_test(data30d, S100betaMeanInS100beta ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.30d.3
# 
#60days
kruskal.test(GFAPoverlapS100betaByNuclei ~ Condition, data = data60d_WTPD)
stat.test.dunn.60d.1  <- dunn_test(data60d_WTPD, GFAPoverlapS100betaByNuclei ~ Condition, p.adjust.method = "BH")
stat.test.dunn.60d.1
# 
# kruskal.test(GFAPMeanInGFAP ~ Condition, data = data60d)
# stat.test.dunn.60d.2  <- dunn_test(data60d, GFAPMeanInGFAP ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.60d.2
# 
# kruskal.test(S100betaMeanInS100beta ~ Condition, data = data60d)
# stat.test.dunn.60d.3  <- dunn_test(data60d, S100betaMeanInS100beta ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.60d.3
# 
#90days
kruskal.test(GFAPoverlapS100betaByNuclei ~ Condition, data = data90d_WTPD)
stat.test.dunn.90d.1  <- dunn_test(data90d_WTPD, GFAPoverlapS100betaByNuclei ~ Condition, p.adjust.method = "BH")
stat.test.dunn.90d.1
# 
# kruskal.test(GFAPMeanInGFAP ~ Condition, data = data90d)
# stat.test.dunn.90d.2  <- dunn_test(data90d, GFAPMeanInGFAP ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.90d.2
# 
# kruskal.test(S100betaMeanInS100beta ~ Condition, data = data90d)
# stat.test.dunn.90d.3  <- dunn_test(data90d, S100betaMeanInS100beta ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.90d.3
# 
# #120days
# kruskal.test(GFAPoverlapS100betaByNuclei ~ Condition, data = data120d)
# stat.test.dunn.120d.1  <- dunn_test(data120d, GFAPoverlapS100betaByNuclei ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.120d.1
# 
# kruskal.test(GFAPMeanInGFAP ~ Condition, data = data120d)
# stat.test.dunn.120d.2  <- dunn_test(data120d, GFAPMeanInGFAP ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.120d.2
# 
# kruskal.test(S100betaMeanInS100beta ~ Condition, data = data120d)
# stat.test.dunn.120d.3  <- dunn_test(data120d, S100betaMeanInS100beta ~ Condition, p.adjust.method = "BH")
# stat.test.dunn.120d.3
# 
# # merge statistics in 1 dataframe
# 
stat.test.dunn.20d.1$Day <- c("20","20","20")
stat.test.dunn.30d.1$Day <- c("30","30","30")
stat.test.dunn.60d.1$Day <- c("60","60","60")
stat.test.dunn.90d.1$Day <- c("90","90","90")
#stat.test.dunn.120d.1$Day <- c("120","120","120")
# 
# stat.test.dunn.20d.2$Day <- c("20","20","20")
# stat.test.dunn.30d.2$Day <- c("30","30","30")
# stat.test.dunn.60d.2$Day <- c("60","60","60")
# stat.test.dunn.90d.2$Day <- c("90","90","90")
# #stat.test.dunn.120d.2$Day <- c("120","120","120")
# 
# stat.test.dunn.20d.3$Day <- c("20","20","20")
# stat.test.dunn.30d.3$Day <- c("30","30","30")
# stat.test.dunn.60d.3$Day <- c("60","60","60")
# stat.test.dunn.90d.3$Day <- c("90","90","90")
# #stat.test.dunn.120d.3$Day <- c("120","120","120")
# 
stat.test.dunn.all.1a <- rbind.data.frame(stat.test.dunn.20d.1, stat.test.dunn.30d.1) #, stat.test.dunn.120d.1)
stat.test.dunn.all.1b <- rbind.data.frame(stat.test.dunn.60d.1, stat.test.dunn.90d.1) #, stat.test.dunn.120d.1)
# 
# stat.test.dunn.all.2 <- rbind.data.frame(stat.test.dunn.20d.2, stat.test.dunn.30d.2, stat.test.dunn.60d.2, stat.test.dunn.90d.2)#, stat.test.dunn.120d.2)
# 
# stat.test.dunn.all.3 <- rbind.data.frame(stat.test.dunn.20d.3, stat.test.dunn.30d.3, stat.test.dunn.60d.3, stat.test.dunn.90d.3)#, stat.test.dunn.120d.3)
```