# ==============================================================================
# NORT Texture (NORT) & Color (NovelObject) Analysis
# Open Field Behavior Genotype x Treatment
# ==============================================================================

library(tidyverse)
library(cowplot)
library(emmeans)

# ==============================================================================
# SECTION 0ata loading & shared constants
# ==============================================================================

# UPDATE this path to your local copy of the file
DATA_PATH <- "DataFrame_AveragedData_17Sep24.csv"

GT_LEVELS <- c('WT\nSaline', 'WT\nESI05', 'KO\nSaline', 'KO\nESI05')

GROUP_COLORS <- c(
  'WT\nSaline' = "grey0",
  'WT\nESI05'  = "grey38",
  'KO\nSaline' = "#ED2224",
  'KO\nESI05'  = "#a31d47"
)

color_scale <- scale_color_manual(values = GROUP_COLORS)
fill_scale  <- scale_fill_manual(values = GROUP_COLORS)
shape_scale <- scale_shape_manual(values = c(21, 22))  # Female = circle, Male = square

base_theme <- theme(
  legend.position = "none",
  axis.title.x    = element_blank(),
  axis.title.y    = element_blank()
)

# Shared stat layers for dot plots
stat_mean_tick   <- stat_summary(fun = mean, geom = "point", size = 10, shape = 95)
stat_sem_errbar  <- stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1, size = 0.5)
zero_hline       <- geom_hline(yintercept = 0, color = 'blue', linetype = 'dashed')

# Helper: pull a named vector of values from the data for t.test()
pull_values <- function(df, behavior, genotype, treatment, col) {
  df %>%
    filter(Behavior == behavior, Genotype == genotype, TreatmentDrug == treatment) %>%
    na.omit() %>%
    pull({{ col }})
}

# ==============================================================================
# SECTION 1 Load & prep data
# ==============================================================================

dat <- read.csv(DATA_PATH) %>%
  mutate(Genotype_Treatment = factor(Genotype_Treatment, levels = GT_LEVELS))

color_dat   <- dat %>% filter(Behavior == 'NovelObject')
texture_dat <- dat %>% filter(Behavior == 'NORT')

# ==============================================================================
# SECTION 2 NOVEL COLOR (NovelObject) Nose_NVOCCPY_INDX
# ==============================================================================

# --- 2a. Dot plot (mean ± SEM + jittered points) ---
ggplot(color_dat, aes(x = Genotype_Treatment, y = Nose_NVOCCPY_INDX, color = Genotype_Treatment)) +
  stat_mean_tick +
  stat_sem_errbar +
  geom_jitter(aes(shape = Sex), alpha = 0.5, size = 2, fill = NA, stroke = 0.6,
              width = 0.05, height = 0) +
  theme_cowplot() + color_scale + shape_scale +
  ylim(c(-0.91, 0.91)) +
  base_theme + theme(axis.text.x = element_blank()) +
  zero_hline
ggsave('Novel_Color_Day6_DotPlot.pdf', dpi = 1000, width = 1.5, height = 1.75)

# --- 2b. Dot plot with x-axis labels (mean ± SEM) ---
ggplot(color_dat, aes(x = Genotype_Treatment, y = Nose_NVOCCPY_INDX, color = Genotype_Treatment)) +
  stat_mean_tick +
  stat_sem_errbar +
  geom_point(aes(shape = Sex), size = 2, alpha = 0.5, fill = 'white', stroke = 0.6,
             position = position_jitter(width = 0.05, height = 0)) +
  theme_cowplot() + color_scale + shape_scale +
  ylim(c(-0.75, 0.75)) +
  base_theme +
  zero_hline
ggsave('Novel_COLOR_Day6_DotPlot_labeled.svg', dpi = 900, width = 1.75, height = 2)

# --- 2c. Boxplot ---
ggplot(color_dat, aes(x = Genotype_Treatment, y = Nose_NVOCCPY_INDX, color = Genotype_Treatment)) +
  geom_boxplot(aes(fill = Genotype_Treatment), alpha = 0.5, width = 0.6) +
  geom_point(aes(shape = Sex), size = 1.5, fill = NA, stroke = 0.5) +
  theme_cowplot() + color_scale + fill_scale + shape_scale +
  ylim(c(-0.75, 0.75)) +
  base_theme +
  zero_hline
ggsave('Novel_Color_Day6_Boxplot.pdf', dpi = 900, width = 1.5, height = 2)

# --- 2d. Descriptive stats ---
color_dat %>%
  group_by(Genotype_Treatment) %>%
  summarize(N = n(),
            Mean = mean(Nose_NVOCCPY_INDX),
            SEM  = sd(Nose_NVOCCPY_INDX) / sqrt(n()))

# --- 2e. 2-way ANOVA + post-hoc ---
anova_color <- aov(Nose_NVOCCPY_INDX ~ Genotype * TreatmentDrug, data = color_dat)
summary(anova_color)
TukeyHSD(anova_color)
emmeans(anova_color, pairwise ~ Genotype * TreatmentDrug)

# --- 2f. Pairwise t-tests (NovelObject only) ---
# KO Saline vs WT Saline
t.test(
  pull_values(dat, 'NovelObject', 'KO', 'Saline',  Nose_NVOCCPY_INDX),
  pull_values(dat, 'NovelObject', 'WT', 'Saline',  Nose_NVOCCPY_INDX)
)

# KO Saline vs KO ESI05
t.test(
  pull_values(dat, 'NovelObject', 'KO', 'Saline',  Nose_NVOCCPY_INDX),
  pull_values(dat, 'NovelObject', 'KO', 'ESI05',   Nose_NVOCCPY_INDX)
)

# ==============================================================================
# SECTION 3  NOVEL TEXTURE (NORT) Nose_NROCCPY_INDX
# ==============================================================================

# --- 3a. Dot plot (mean ± SEM + jittered points) ---
ggplot(texture_dat, aes(x = Genotype_Treatment, y = Nose_NROCCPY_INDX, color = Genotype_Treatment)) +
  stat_mean_tick +
  stat_sem_errbar +
  geom_point(aes(shape = Sex), size = 2, alpha = 0.5, fill = 'white', stroke = 0.6,
             position = position_jitter(width = 0.05, height = 0)) +
  theme_cowplot() + color_scale + shape_scale +
  ylim(c(-0.15, 0.8)) +
  base_theme + theme(axis.text.x = element_blank()) +
  zero_hline
ggsave('Novel_Texture_Day6_DotPlot.svg', dpi = 900, width = 1.75, height = 2)

# --- 3b. Boxplot ---
ggplot(texture_dat, aes(x = Genotype_Treatment, y = Nose_NROCCPY_INDX, color = Genotype_Treatment)) +
  geom_boxplot(aes(fill = Genotype_Treatment), alpha = 0.5, width = 0.6) +
  geom_point(aes(shape = Sex), size = 1.5, fill = NA, stroke = 0.5) +
  theme_cowplot() + color_scale + fill_scale + shape_scale +
  ylim(c(-0.2, 0.6)) +
  base_theme +
  zero_hline
ggsave('Novel_Texture_Day6_Boxplot.pdf', dpi = 900, width = 2, height = 2.5)

# --- 3c. Descriptive stats ---
texture_dat %>%
  group_by(Genotype_Treatment) %>%
  summarize(N = n(),
            Mean = mean(Nose_NROCCPY_INDX),
            SEM  = sd(Nose_NROCCPY_INDX) / sqrt(n()))

# --- 3d. 2-way ANOVA + post-hoc ---
anova_texture <- aov(Nose_NROCCPY_INDX ~ Genotype * TreatmentDrug, data = texture_dat)
summary(anova_texture)
TukeyHSD(anova_texture)
emmeans(anova_texture, pairwise ~ Genotype * TreatmentDrug)

# --- 3e. Pairwise t-tests (NORT only) ---
# KO Saline vs WT Saline
t.test(
  pull_values(dat, 'NORT', 'KO', 'Saline', Nose_NROCCPY_INDX),
  pull_values(dat, 'NORT', 'WT', 'Saline', Nose_NROCCPY_INDX)
)

# KO Saline vs KO ESI05
t.test(
  pull_values(dat, 'NORT', 'KO', 'Saline', Nose_NROCCPY_INDX),
  pull_values(dat, 'NORT', 'KO', 'ESI05',  Nose_NROCCPY_INDX)
)
