library(tidyverse)

if (!dir.exists("out_R/summary/stdev")) {
  dir.create("./out_R/summary/stdev")
}

if (!dir.exists("out_R/summary/std_err")) {
  dir.create("./out_R/summary/std_err")
}

if (!dir.exists("out_R/individual")) {
  dir.create("./out_R/individual")
}


## Read in datafile from skimage
d <- read_csv("./skimage_metrics.csv")

## Calculate relative values
rel <- d %>%
  group_by(Group, Fish, Genotype) %>%
  mutate(across(area:integrated_density, ~ .x/.x[Day == 0]),
         Group = fct_recode(Group, 
                            Control = "ctrl", 
                            Cisplatin = "cddp"),
         Group = fct_relevel(Group, "Control"),
         Fish = as_factor(Fish))

properties <- c("area", "filled_area", "major_axis_length", 
                "minor_axis_length", "mean_intensity", "integrated_density",
                "raw_integrated_density")
## Summary plots
for (prop in properties) {
  ggplot(filter(rel, Day < 30), aes_string("Day", prop, fill = "Group")) + 
    stat_summary(fun = mean, geom = "line") + 
    stat_summary(fun = mean, geom = "point", shape = 21) +
    stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.5) +
    geom_hline(aes(yintercept = 1), linetype = "dashed") + 
    facet_wrap(~ Genotype, nrow = 2)
  ggsave(paste0("./out_R/summary/std_err/", prop, ".png"), scale = 0.8)
}
  
## Individual plots
for (prop in properties) {
  ggplot(filter(rel, Genotype == "sw58"), aes_string("Day", prop)) + 
    geom_line(aes(color = Fish)) +
    geom_point(shape = 21, aes(fill = Fish)) +
    facet_wrap(~ Group) +
    coord_cartesian(ylim = c(0,2))
  ggsave(paste0("./out_R/individual/sw58_",prop,".png"), scale = 0.8)
}

for (prop in properties) {
  ggplot(filter(rel, Genotype == "sw70"), aes_string("Day", prop)) + 
    geom_line(aes(color = Fish)) +
    geom_point(shape = 21, aes(fill = Fish)) +
    facet_wrap(~ Group) +
    coord_cartesian(ylim = c(0,2.2))
  ggsave(paste0("./out_R/individual/sw70_",prop,".png"), scale = 0.8)
}

## Hypothesis testing
control_area <- rel %>%
  filter(Group == "Control", 
         Genotype == "sw58",
         Day == 28) %>%
  pull(area)
cddp_area <- rel %>%
  filter(Group == "Cisplatin", 
         Genotype == "sw58",
         Day == 28) %>%
  pull(area)
p_area_sw58 <- wilcox.test(control_area, cddp_area)$p.value

control_area <- rel %>%
  filter(Group == "Control", 
         Genotype == "sw70",
         Day == 28) %>%
  pull(area)
cddp_area <- rel %>%
  filter(Group == "Cisplatin", 
         Genotype == "sw70",
         Day == 28) %>%
  pull(area)
p_area_sw70 <- wilcox.test(control_area, cddp_area)$p.value

control_int_den <- rel %>%
  filter(Group == "Control", 
         Genotype == "sw58",
         Day == 28) %>%
  pull(integrated_density)
cddp_int_den <- rel %>%
  filter(Group == "Cisplatin", 
         Genotype == "sw58",
         Day == 28) %>%
  pull(integrated_density)
p_intDen_sw58 <- wilcox.test(control_int_den, cddp_int_den)$p.value

control_int_den <- rel %>%
  filter(Group == "Control", 
         Genotype == "sw70",
         Day == 28) %>%
  pull(integrated_density)
cddp_int_den <- rel %>%
  filter(Group == "Cisplatin", 
         Genotype == "sw70",
         Day == 28) %>%
  pull(integrated_density)
p_intDen_sw70 <- wilcox.test(control_int_den, cddp_int_den)$p.value
p_vals <- tibble(Metric = c("Area", "Integrated Density"),
                 sw58 = c(p_area_sw58, p_intDen_sw58),
                 sw70 = c(p_area_sw70, p_intDen_sw70))
write_csv(p_vals, "./out_R/p_values.csv")