library(cowplot)
library(svglite)
library(tidyverse)
library(ggrepel)
library(ggpubr)
library(pheatmap)
library(RColorBrewer)
library(viridis)
library(janitor)
library(emmeans)

## ---------------------------------------------------------------------------
## WhiskPlotFacemap: longitudinal boxplot + diff plot for whisker stim data
## ---------------------------------------------------------------------------
WhiskPlotFacemap = function(dataframetorun, Variable, ylabel, xcategory, ggtitle, maincolor, Filename, ylimit){
  Mice_ID = c(unique(paste0(dataframetorun$Treatment, "_", dataframetorun$Mice_ID)), "ESI05", "Vehicle")
  dataframetorun$VariabletoPlot = Variable
  dataframetorun$Mice_ID_Treat = paste0(dataframetorun$Treatment, "_", dataframetorun$Mice_ID)
  colormap = data.frame(ID = Mice_ID[order(Mice_ID)], Color = rep(0, length(Mice_ID)))
  colormap$Color[grepl("ESI05",  colormap$ID)] = colorRampPalette(c("darksalmon", "indianred1"))(sum(grepl("ESI05",  colormap$ID)))
  colormap$Color[grepl("Vehicle", colormap$ID)] = colorRampPalette(c("grey50",    "grey30"))(sum(grepl("Vehicle", colormap$ID)))
  colormap$Color[colormap$ID == "ESI05"]   = maincolor[1]
  colormap$Color[colormap$ID == "Vehicle"] = maincolor[2]

  meandf = dataframetorun %>%
    group_by(Treatment_Day) %>%
    summarise(Median_Tot_Energy = median(VariabletoPlot)) %>%
    mutate(Treatment = stringr::str_extract(Treatment_Day, "ESI05|Vehicle")) %>%
    mutate(Day = stringr::str_extract(Treatment_Day, "Baseline|Acute|Chronic"))

  x_levels = c("Vehicle \n Baseline", "Vehicle \n Acute", "Vehicle \n Chronic",
                "ESI05 \n Baseline",   "ESI05 \n Acute",   "ESI05 \n Chronic")

  boxplt = ggplot(dataframetorun,
                  aes(x = factor(Treatment_Day, levels = x_levels),
                      y = VariabletoPlot,
                      color = factor(Treatment, levels = c("Vehicle", "ESI05")),
                      fill  = factor(Treatment, levels = c("Vehicle", "ESI05")))) +
    geom_boxplot(outlier.shape = NA, coef = 0, fatten = 0) +
    geom_point(aes(shape = factor(Sex)), fill = NA, size = 6) +
    geom_line(aes(group = factor(Mice_ID_Treat)), linetype = "twodash", linewidth = 0.75, alpha = 0.45) +
    geom_line(data = meandf,
              aes(x = factor(Treatment_Day, levels = x_levels), y = Median_Tot_Energy, group = Treatment),
              linewidth = 3, alpha = 1) +
    scale_color_manual(values = alpha(c(maincolor), 0.8)) +
    scale_fill_manual(values  = alpha(c(maincolor), 0.65)) +
    scale_shape_manual(values = c(2, 1)) +
    ylab(ylabel) + ggtitle(ggtitle) + theme_cowplot() +
    theme(axis.title.x = element_blank(),
          axis.text.x  = element_text(size = rel(1), angle = 0),
          axis.title.y = element_text(size = rel(1)),
          title        = element_text(size = rel(1)),
          legend.position = "none") +
    guides(color = guide_legend(title = ""), fill = "none", alpha = "none")

  vlnplt = ggplot(dataframetorun,
                  aes(x = factor(Treatment_Day, levels = x_levels),
                      y = VariabletoPlot,
                      color = factor(Treatment, levels = c("Vehicle", "ESI05")),
                      fill  = factor(Treatment, levels = c("Vehicle", "ESI05")))) +
    geom_point(aes(shape = factor(Sex)), fill = NA, size = 6) +
    geom_violin() +
    geom_line(aes(group = factor(Mice_ID_Treat)), linewidth = 1, linetype = "twodash", alpha = 0.35) +
    geom_line(data = meandf,
              aes(x = factor(Treatment_Day, levels = x_levels), y = Median_Tot_Energy, group = Treatment),
              linewidth = 3, alpha = 1) +
    scale_color_manual(values = alpha(c(maincolor), 0.8)) +
    scale_fill_manual(values  = alpha(c(maincolor), 0.65)) +
    scale_shape_manual(values = c(2, 1)) +
    ylab(ylabel) + ggtitle(ggtitle) + theme_cowplot() +
    theme(axis.title.x = element_blank(),
          axis.text.x  = element_text(size = rel(1), angle = 0),
          axis.title.y = element_text(size = rel(1)),
          title        = element_text(size = rel(1)),
          legend.position = "none") +
    guides(color = guide_legend(title = ""), fill = "none", alpha = "none")

  D0df = dataframetorun %>% subset(Day == "D0")
  D1df = dataframetorun %>% subset(Day == "D1")
  D5df = dataframetorun %>% subset(Day == "D5")
  D1df$Diff = D1df$VariabletoPlot - D0df$VariabletoPlot
  D5df$Diff = D5df$VariabletoPlot - D0df$VariabletoPlot
  DiffDf = rbind(D1df, D5df)

  DiffDf$Factortoplot = case_when(
    DiffDf$Treatment_Day == "ESI05 \n Acute"    ~ 'Acute-Baseline ESI05',
    DiffDf$Treatment_Day == "ESI05 \n Chronic"  ~ 'Chronic-Baseline ESI05',
    DiffDf$Treatment_Day == "Vehicle \n Acute"  ~ 'Acute-Baseline Vehicle',
    DiffDf$Treatment_Day == "Vehicle \n Chronic" ~ 'Chronic-Baseline Vehicle'
  )

  diff_levels = c('Acute-Baseline ESI05', 'Chronic-Baseline ESI05',
                  'Acute-Baseline Vehicle', 'Chronic-Baseline Vehicle')

  diffplt = ggplot(DiffDf,
                   aes(x = Diff,
                       y = factor(Factortoplot, levels = diff_levels),
                       color = factor(Treatment))) +
    geom_boxplot(aes(fill = factor(Treatment), alpha = 0.95), coef = 0, outlier.shape = NA, fatten = 0) +
    geom_point(size = 5, alpha = 0.5) +
    geom_vline(xintercept = 0, linewidth = 2, linetype = "dotted") +
    stat_summary(fun = "mean", geom = "point", shape = 20, size = 14, color = "grey0", fill = "grey0", alpha = 0.75) +
    scale_color_manual(values = c(ESI05 = maincolor[2], Vehicle = maincolor[1])) +
    scale_fill_manual(values  = c(ESI05 = maincolor[2], Vehicle = maincolor[1])) +
    xlab("Total Motion Energy - Total Motion Energy Baseline") +
    ggtitle(ggtitle) + theme_cowplot() +
    theme(axis.title.y = element_blank(),
          axis.text.x  = element_text(size = rel(1)),
          axis.title.x = element_text(size = rel(1)),
          title        = element_text(size = rel(1))) +
    guides(color = guide_legend(title = "Groups"), fill = "none", alpha = "none")

  diffplt_vln = ggplot(DiffDf,
                       aes(x = Diff,
                           y = factor(Factortoplot, levels = diff_levels),
                           color = factor(Treatment))) +
    geom_violin(aes(fill = factor(Treatment), alpha = 0.95)) +
    geom_point(size = 5, alpha = 0.5) +
    geom_vline(xintercept = 0, linewidth = 2, linetype = "dotted") +
    stat_summary(fun = "mean", geom = "point", shape = 20, size = 14, color = "grey0", fill = "grey0", alpha = 0.75) +
    scale_color_manual(values = c(ESI05 = maincolor[2], Vehicle = maincolor[1])) +
    scale_fill_manual(values  = c(ESI05 = maincolor[2], Vehicle = maincolor[1])) +
    ylab("Total Motion Energy - Total Motion Energy Baseline") +
    ggtitle(ggtitle) + theme_cowplot() +
    theme(axis.title.x = element_blank(),
          axis.text.y  = element_text(size = rel(1)),
          axis.title.y = element_text(size = rel(1)),
          title        = element_text(size = rel(1))) +
    guides(color = guide_legend(title = "Groups"), fill = "none", alpha = "none") +
    coord_flip()

  ggsave(boxplt,     filename = paste0(Filename, "StimPlots",        Sys.Date(), ".pdf"), width = 4, height = 5, units = "in")
  ggsave(boxplt,     filename = paste0(Filename, "StimPlots",        Sys.Date(), ".svg"), width = 4, height = 5, units = "in")
  ggsave(diffplt,    filename = paste0(Filename, "DiffPlots",        Sys.Date(), ".pdf"), width = 8, height = 4, units = "in")
  ggsave(diffplt,    filename = paste0(Filename, "DiffPlots",        Sys.Date(), ".svg"), width = 8, height = 4, units = "in")
  ggsave(diffplt_vln,filename = paste0(Filename, "Violin_DiffPlots", Sys.Date(), ".pdf"), width = 8, height = 5, units = "in")
  ggsave(diffplt_vln,filename = paste0(Filename, "Violin_DiffPlots", Sys.Date(), ".svg"), width = 8, height = 5, units = "in")

  boxplt + diffplt + plot_layout(width = c(3, 1))
  return(DiffDf)
}

## ---------------------------------------------------------------------------
## Data loading
## ---------------------------------------------------------------------------

Grabbing_Sofia   = read.csv("Tactile_Defensiveness_Grabbing_Sofia_21Jan26.csv")
Grabbing_Index_KO = readxl::read_xlsx("IndexofMice_21Jan25.xlsx")

## ---------------------------------------------------------------------------
## Transpose and wrangle grabbing data
## ---------------------------------------------------------------------------
Grabbing_Sofia_Transposed = as.data.frame(t(Grabbing_Sofia))
Grabbing_Sofia_Transposed = Grabbing_Sofia_Transposed %>% row_to_names(row_number = 1)
Grabbing_Sofia_Transposed = Grabbing_Sofia_Transposed %>% mutate(ID = row.names(Grabbing_Sofia_Transposed))

Grabbing_Sofia_Transposed2 = Grabbing_Sofia_Transposed %>%
  mutate(Mice_ID = str_extract(ID, "^[A-Za-z]+\\d+")) %>%
  mutate(Day = case_when(
    grepl('Day0|D0',  ID) ~ "Day_0",
    grepl('Day6|_D5', ID) ~ "Day_6"
  )) %>%
  mutate(Genotype = case_when(
    grepl('ESI05',       ID) ~ "KO",
    grepl('Blue|Black',  ID) ~ "WT"
  ))

Grabbing_Sofia_Transposed2 = left_join(Grabbing_Sofia_Transposed2, Grabbing_Index_KO, by = 'Mice_ID')
Grabbing_Sofia_Transposed2$ROWID = paste0(Grabbing_Sofia_Transposed2$Genotype.x, "_", Grabbing_Sofia_Transposed2$Treatment)
Grabbing_Sofia_Transposed2 = Grabbing_Sofia_Transposed2 %>% arrange(desc(ROWID))

## Compute ROW_SUM from first 220 numeric-named columns
Grabbing_Sofia_Transposed2 = Grabbing_Sofia_Transposed2 %>%
  mutate(ROW_SUM = rowSums(dplyr::select(., dplyr::matches("^\\d"))[, 1:220], na.rm = TRUE))

## ---------------------------------------------------------------------------
## Heatmaps
## ---------------------------------------------------------------------------
day0_data = Grabbing_Sofia_Transposed2 %>% dplyr::filter(Day == 'Day_0')
Grabbing_Heatmap_Day0 = day0_data %>%
  dplyr::select(dplyr::matches("^\\d")) %>% dplyr::select(1:220) %>%
  pheatmap::pheatmap(cluster_cols = FALSE, cluster_rows = FALSE,
                     labels_row = day0_data$ROWID,
                     color = viridis(2), show_colnames = FALSE)
ggsave(paste0('SOFIA_GRABBING_HEATMAP_D0', Sys.Date(), ".pdf"), dpi = 1000, Grabbing_Heatmap_Day0, width = 5, height = 3, units = "in")

day6_data = Grabbing_Sofia_Transposed2 %>% dplyr::filter(Day == 'Day_6')
Grabbing_Heatmap_Day6 = day6_data %>%
  dplyr::select(dplyr::matches("^\\d")) %>% dplyr::select(1:220) %>%
  pheatmap::pheatmap(cluster_cols = FALSE, cluster_rows = FALSE,
                     labels_row = day6_data$ROWID,
                     color = viridis(2), show_colnames = FALSE)
ggsave(paste0('SOFIA_GRABBING_HEATMAP_D6', Sys.Date(), ".pdf"), dpi = 1000, Grabbing_Heatmap_Day6, width = 5, height = 3, units = "in")

## ---------------------------------------------------------------------------
## Build Table_Grabbing_Merged — mirrors CodeForPlotting_TACTILEDEF structure
## ---------------------------------------------------------------------------
## The grabbing CSV has 300 time-bin rows x 60 mouse-session columns.
## After transposing, each row is one mouse-session. We extract metadata from
## the column name (MouseID, Day, Sex, Treatment/Color) and merge with the
## index file, then compute Number_of_GrabsPerTime = ROW_SUM / 300 seconds.
## Day mapping: Day0/D0 -> "D0" (Baseline), Day6/D5 -> "D5" (Chronic),
## matching the D0/D1/D5 convention used in WhiskPlotFacemap.

## Extract Sex directly from the column ID string where it is encoded
## WT cohort:  MouseID_Color_Day_Sex_Freq_Cam  (e.g. ZD294_Blue_Day0_F_20Hz_Cam1)
## KO cohort:  MouseID_ESI05_Freq_Day_Rep_Date (no sex field — comes from index file)
Grabbing_Sofia_Transposed2 = Grabbing_Sofia_Transposed2 %>%
  mutate(Sex_fromID = str_extract(ID, "(?<=_)(F|M)(?=_)"))

## Where Sex is already in the joined index (Grabbing_Index_KO), prefer that;
## fall back to Sex_fromID for WT mice whose sex isn't in the index
Grabbing_Sofia_Transposed2 = Grabbing_Sofia_Transposed2 %>%
  mutate(Sex = coalesce(Sex, Sex_fromID))

## Remap Day to D0/D5 to match WhiskPlotFacemap internals (subset(Day=="D0") etc.)
## Note: there is no D1 (Acute) timepoint in this dataset — only Baseline and Chronic
Table_Grabbing_Merged = Grabbing_Sofia_Transposed2 %>%
  mutate(
    Day = case_when(
      Day == "Day_0" ~ "D0",
      Day == "Day_6" ~ "D5",
      TRUE           ~ NA_character_
    ),
    Number_of_GrabsPerTime = ROW_SUM / 300,           # grabs per second over 300s session
    LogValues              = log10(Number_of_GrabsPerTime + 1),
    Time = case_when(
      Day == "D0" ~ "Baseline",
      Day == "D5" ~ "Chronic",
      TRUE        ~ NA_character_
    ),
    Treatment_Day = paste0(Treatment, " \n ", Time)
  ) %>%
  filter(!is.na(Day), !is.na(Treatment), Treatment != "NA", Genotype.x != "NA")

## ---------------------------------------------------------------------------
## WhiskPlotFacemap calls on grabbing data — raw and log-transformed
## ---------------------------------------------------------------------------
## Note: xcategory labels cover Baseline/Acute/Chronic but only Baseline and
## Chronic timepoints exist here; Acute slots will be empty in the plots.

Grabbing_DiffDf = WhiskPlotFacemap(
  Table_Grabbing_Merged,
  Table_Grabbing_Merged$Number_of_GrabsPerTime,
  "Number of Grabs/Second",
  c("Baseline", "Acute", "Chronic", "Baseline", "Acute", "Chronic"),
  "Number of Grabs/Second",
  c("grey10", "red2"),
  "GRAB_AS_MeanTotEnergyPerStim",
  c(0, 50)
)

Grabbing_Log_DiffDf = WhiskPlotFacemap(
  Table_Grabbing_Merged,
  Table_Grabbing_Merged$LogValues,
  "log10(Number of Grabs/Second + 1)",
  c("Baseline", "Acute", "Chronic", "Baseline", "Acute", "Chronic"),
  "Number of Grabs/Second (log10)",
  c("grey10", "red2"),
  "LN_GRAB_AS_MeanTotEnergyPerStim",
  c(0, 2)
)

## ---------------------------------------------------------------------------
## Statistics on grabbing data — matching CodeForPlotting_TACTILEDEF
## ---------------------------------------------------------------------------
levene_test(Number_of_GrabsPerTime ~ Time * Treatment, data = Table_Grabbing_Merged) # homogeneity of variance
levene_test(LogValues ~ Time * Treatment,              data = Table_Grabbing_Merged) # homogeneity of variance (log)

aov_Grab_LogValues = aov(LogValues ~ Treatment * Time, data = Table_Grabbing_Merged)
ggplot(Table_Grabbing_Merged, aes(x = LogValues)) + geom_density()   # check distribution
summary(aov_Grab_LogValues)
shapiro.test(x = residuals(object = aov_Grab_LogValues))              # normality of residuals
TukeyHSD(aov_Grab_LogValues)

## Chronic only (no Acute timepoint, but kept for structural consistency)
Table_Grabbing_Merged_Chronic = Table_Grabbing_Merged %>% filter(Day != "D1")
Grabbing_Chronic_DiffDf = WhiskPlotFacemap(
  Table_Grabbing_Merged_Chronic,
  Table_Grabbing_Merged_Chronic$LogValues,
  "log10(Number of Grabs/Second + 1)",
  c("Baseline", "Acute", "Chronic", "Baseline", "Acute", "Chronic"),
  "Number of Grabs/Second — Chronic only (log10)",
  c("grey10", "red2"),
  "LN_CHRONIC_GRAB_AS_MeanTotEnergyPerStim",
  c(0, 2)
)

levene_test(LogValues ~ Time * Treatment, data = Table_Grabbing_Merged_Chronic)
aov_Grab_Chronic = aov(LogValues ~ Treatment * Time, data = Table_Grabbing_Merged_Chronic)
summary(aov_Grab_Chronic)
shapiro.test(x = residuals(object = aov_Grab_Chronic))
TukeyHSD(aov_Grab_Chronic)

## ---------------------------------------------------------------------------
## Day 6 grabbing — Genotype x Treatment
## ---------------------------------------------------------------------------
Grabbing_Day6_plot = ggplot(
  Grabbing_Sofia_Transposed2 %>% filter(Day == 'Day_6'),
  aes(x = factor(ROWID, levels = c('WT_Vehicle', 'WT_ESI05', 'KO_Vehicle', 'KO_ESI05')),
      y = ROW_SUM, color = ROWID, fill = ROWID)) +
  stat_summary(fun = mean, geom = "point", size = 10, shape = 95) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1, linewidth = 0.5) +
  geom_jitter(aes(shape = factor(Sex)), fill = NA, size = 2, width = 0.1, height = 0.1) +
  theme_cowplot() +
  scale_shape_manual(values = c('F' = 1, 'M' = 0)) +
  scale_color_manual(values = c('WT_Vehicle' = "grey0", 'WT_ESI05' = "grey38", 'KO_Vehicle' = "#ED2224", 'KO_ESI05' = "#a31d47")) +
  scale_fill_manual(values  = c('WT_Vehicle' = "grey0", 'WT_ESI05' = "grey38", 'KO_Vehicle' = "#ED2224", 'KO_ESI05' = "#a31d47")) +
  theme(axis.title.x = element_blank(), legend.position = "none", plot.title = element_blank(),
        axis.title.y = element_blank(), axis.text.x = element_blank()) +
  ylab('Time grabbing (in seconds)') + xlab('')

ggsave(paste0('SOFIA_GRABBING_D6', Sys.Date(), ".pdf"), plot = Grabbing_Day6_plot, dpi = 1000, width = 1.5, height = 1.75, units = "in")

aov_Grab_D6 = aov(ROW_SUM ~ Genotype.x * Treatment, data = Grabbing_Sofia_Transposed2 %>% filter(Day == 'Day_6'))
summary(aov_Grab_D6)
emmeans(aov_Grab_D6, pairwise ~ Genotype.x * Treatment)
Grabbing_Sofia_Transposed2 %>% filter(Day == 'Day_6') %>%
  group_by(Genotype.x, Treatment) %>%
  summarize(mean = mean(ROW_SUM), SEM = sd(ROW_SUM) / sqrt(n()))

## ---------------------------------------------------------------------------
## Day 0 grabbing — Genotype x Treatment
## ---------------------------------------------------------------------------
Grabbing_Day0_plot = ggplot(
  Grabbing_Sofia_Transposed2 %>% filter(Day == 'Day_0', ROWID != 'WT_NA', ROWID != 'NA_NA'),
  aes(x = factor(ROWID, levels = c('WT_Vehicle', 'WT_ESI05', 'KO_Vehicle', 'KO_ESI05')),
      y = ROW_SUM, color = ROWID, fill = ROWID)) +
  stat_summary(fun = mean, geom = "point", size = 10, shape = 95) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1, linewidth = 0.5) +
  geom_jitter(aes(shape = factor(Sex)), fill = NA, size = 2, width = 0.1, height = 0.1, alpha = 0.5) +
  theme_cowplot() +
  scale_shape_manual(values = c('F' = 1, 'M' = 0)) +
  scale_color_manual(values = c('WT_Vehicle' = "grey0", 'WT_ESI05' = "grey38", 'KO_Vehicle' = "#ED2224", 'KO_ESI05' = "#a31d47")) +
  scale_fill_manual(values  = c('WT_Vehicle' = "grey0", 'WT_ESI05' = "grey38", 'KO_Vehicle' = "#ED2224", 'KO_ESI05' = "#a31d47")) +
  theme(axis.title.x = element_blank(), legend.position = "none", plot.title = element_blank(),
        axis.title.y = element_blank(), axis.text.x = element_blank()) +
  ylab('Time grabbing (in seconds)') + xlab('')

ggsave(paste0('SOFIA_GRABBING_ALLPLOTs_Day0', Sys.Date(), ".pdf"), plot = Grabbing_Day0_plot, dpi = 1000, width = 1.5, height = 1.75, units = "in")

aov_Grab_D0_GenoTreat = aov(ROW_SUM ~ Genotype.x * Treatment,
  data = Grabbing_Sofia_Transposed2 %>% filter(Day == 'Day_0', ROWID != 'WT_NA', ROWID != 'NA_NA'))
summary(aov_Grab_D0_GenoTreat)
emmeans(aov_Grab_D0_GenoTreat, pairwise ~ Genotype.x * Treatment)

## ---------------------------------------------------------------------------
## Day 0 grabbing — Genotype only
## ---------------------------------------------------------------------------
Grabbing_Day0_Genotype = ggplot(
  Grabbing_Sofia_Transposed2 %>% filter(Day == 'Day_0', ROWID != 'WT_NA', ROWID != 'NA_NA'),
  aes(x = factor(Genotype.x, levels = c('WT', 'KO')),
      y = ROW_SUM, color = Genotype.x, fill = Genotype.x)) +
  stat_summary(fun = mean, geom = "point", size = 10, shape = 95) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1, linewidth = 0.5) +
  geom_jitter(aes(shape = factor(Sex)), fill = NA, size = 2, width = 0.1, height = 0.1, alpha = 0.5) +
  theme_cowplot() +
  scale_shape_manual(values = c('F' = 1, 'M' = 0)) +
  scale_color_manual(values = c('WT' = "grey0", 'KO' = "#ED2224")) +
  scale_fill_manual(values  = c('WT' = "grey0", 'KO' = "#ED2224")) +
  theme(axis.title.x = element_blank(), legend.position = "none", plot.title = element_blank(),
        axis.title.y = element_blank(), axis.text.x = element_blank()) +
  ylab('Time grabbing (in seconds)') + xlab('')

ggsave(paste0('Grabbing_Day0_Genotype', Sys.Date(), ".pdf"), plot = Grabbing_Day0_Genotype, dpi = 1000, width = 1.5, height = 1.75, units = "in")

aov_Grab_D0_Geno = aov(ROW_SUM ~ Genotype.x,
  data = Grabbing_Sofia_Transposed2 %>% filter(Day == 'Day_0', ROWID != 'WT_NA', ROWID != 'NA_NA'))
summary(aov_Grab_D0_Geno)
emmeans(aov_Grab_D0_Geno, pairwise ~ Genotype.x)
Grabbing_Sofia_Transposed2 %>%
  filter(Day == 'Day_0', ROWID != 'WT_NA', ROWID != 'NA_NA') %>%
  group_by(Genotype.x) %>%
  summarize(mean = mean(ROW_SUM), SEM = sd(ROW_SUM) / sqrt(n()))

## ---------------------------------------------------------------------------
## Day 0 vs Day 6 grabbing — KO mice, by Treatment (longitudinal)
## ---------------------------------------------------------------------------
Grabbing_Day0vsDay6_KO = ggplot(
  Grabbing_Sofia_Transposed2 %>%
    mutate(GROUP = paste(Treatment, '\n', Day)) %>%
    filter(Genotype.x == 'KO'),
  aes(x = factor(GROUP, levels = c("Vehicle \n Day_0", "Vehicle \n Day_6", "ESI05 \n Day_0", "ESI05 \n Day_6")),
      y = ROW_SUM, color = GROUP, fill = GROUP)) +
  stat_summary(fun = mean, geom = "point", size = 10, shape = 95) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1, linewidth = 0.5) +
  geom_point(aes(shape = factor(Sex)), fill = NA, size = 2, alpha = 0.5) +
  geom_line(aes(group = factor(Mice_ID)), linewidth = 0.5, alpha = 0.45) +
  theme_cowplot() +
  scale_shape_manual(values = c('F' = 1, 'M' = 0)) +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size = rel(1), angle = 0),
        axis.title.y = element_text(size = rel(1)), title = element_text(size = rel(1)),
        legend.position = "none") +
  scale_color_manual(values = c('ESI05 \n Day_0'   = "#a31d47", 'ESI05 \n Day_6'   = "#a31d47",
                                 'Vehicle \n Day_0' = "grey38",  'Vehicle \n Day_6' = "grey0")) +
  scale_fill_manual(values  = c('ESI05 \n Day_0'   = "#a31d47", 'ESI05 \n Day_6'   = "#a31d47",
                                 'Vehicle \n Day_0' = "grey38",  'Vehicle \n Day_6' = "grey0")) +
  ylab('Time grabbing (in seconds)') + xlab('')

ggsave(paste0('Grabbing_Day0vsDay6_KO', Sys.Date(), ".pdf"), plot = Grabbing_Day0vsDay6_KO, dpi = 1000, width = 1.5, height = 1.75, units = "in")

aov_Grab_D0_D6_KO = aov(log10(ROW_SUM + 1) ~ Treatment * Day,
  data = Grabbing_Sofia_Transposed2 %>% filter(Genotype.x == 'KO'))
summary(aov_Grab_D0_D6_KO)
emmeans(aov_Grab_D0_D6_KO, pairwise ~ Treatment * Day)
Grabbing_Sofia_Transposed2 %>% filter(Genotype.x == 'KO') %>%
  group_by(Treatment, Day) %>%
  summarize(mean = mean(ROW_SUM), SEM = sd(ROW_SUM) / sqrt(n()))
