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

## ---------------------------------------------------------------------------
## Data loading
## ---------------------------------------------------------------------------
Grabbing_Sofia    = read.csv("Tactile_Defensiveness_Grabbing_Sofia_21Jan25.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+")) %>%       # fixed: was $I, should be $ID
  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 time-bin columns
Grabbing_Sofia_Transposed2 = Grabbing_Sofia_Transposed2 %>%
  mutate(ROW_SUM = rowSums(dplyr::select(., dplyr::matches("^\\d"))[, 1:220], na.rm = TRUE))

## ---------------------------------------------------------------------------
## Heatmaps — filter first so labels_row matches the plotted rows
## ---------------------------------------------------------------------------
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"),
       plot = Grabbing_Heatmap_Day0, dpi = 1000, 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"),
       plot = Grabbing_Heatmap_Day6, dpi = 1000, width = 5, height = 3, units = "in")

## ---------------------------------------------------------------------------
## 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")) +  # fixed: keys now match Genotype.x
  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()))

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

## ---------------------------------------------------------------------------
## Day 0 vs Day 6 grabbing — KO mice only, longitudinal by Treatment
## ---------------------------------------------------------------------------
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(),                              # merged into single theme() call
        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()))
