# ==============================================================================
# AGS Summary Analysis
# Refactored: 2026-04-08

library(cowplot)
library(svglite)
library(ggalluvial)
library(tidyverse)
library(ggplot2)
library(ggrepel)
library(tidyr)
library(ggpubr)
library(ggdendro)
library(dendextend)   
library(pheatmap)
library(RColorBrewer)
library(viridis)
library(magrittr)
library(VennDiagram)
library(gridExtra)
library(grid)
library(plotrix)
library(eulerr)
library(stringr)
library(STRINGdb)
library(janitor)
library(emmeans)


AGS_Plotting_org <- read.csv('AGS_datasheet_ForPlotting.csv')


# Infer Sex from filename coding
add_sex <- function(df, col = Filename.Coding) {
  df %>% mutate(Sex = case_when(
    grepl('_sexM',  {{ col }}) | grepl('Male',   {{ col }}) ~ 'Male',
    grepl('_sexF',  {{ col }}) | grepl('Female', {{ col }}) ~ 'Female'
  ))
}

# Trim whitespace from Treatment and split into Genotype + Treat columns
split_treatment <- function(df) {
  df %>%
    mutate(Treatment = str_trim(Treatment)) %>%
    separate_wider_delim(cols = Treatment, delim = " ",
                         names = c("Genotype", "Treat"))
}

# Create an ordered Treatment factor from existing Genotype + Treat columns
# Uses local column names — safe inside mutate pipes
make_treatment_factor <- function(df) {
  df %>% mutate(
    Treatment = factor(
      paste0(Genotype, "_", Treat),
      levels = c("WT_SALINE", "WT_ESI", "KO_SALINE", "KO_ESI")
    )
  )
}

TREATMENT_LEVELS_SPACE  <- c("WT SALINE", "WT ESI", "KO SALINE", "KO ESI")
TREATMENT_LEVELS_UNDER  <- c("WT_SALINE", "WT_ESI", "KO_SALINE", "KO_ESI")

AGS_COLORS <- c(
  'KO_ESI'    = "#a31d47",
  'KO_SALINE' = "#ED2224",
  'WT_SALINE' = "grey0",
  'WT_ESI'    = "grey38"
)

BEHAVIOR_COLORS <- c(
  "Death"         = "#FFD0CD",
  "TonClon"       = "#E8B3AF",
  "Myoclonic"     = "#D19591",
  "Wildrun.multi" = "#BA7873",
  "Wildrun"       = "#A35B55",
  "NoSeize"       = "#616161"
)

ags_color_scale  <- scale_color_manual(values = AGS_COLORS)
ags_fill_scale   <- scale_fill_manual(values  = AGS_COLORS)
ags_shape_scale  <- scale_shape_manual(values = c('Female' = 1, 'Male' = 0))
ags_theme_clean  <- theme(panel.grid = element_blank(), legend.position = "none")
ags_theme_noaxis <- theme(
  panel.grid   = element_blank(), legend.position  = "none",
  axis.title.x = element_blank(), axis.title.y     = element_blank(),
  axis.text.x  = element_blank(), axis.text.y      = element_blank()
)
behavior_fill_scale <- scale_fill_manual(values = BEHAVIOR_COLORS)

# Shared stat layers for dot-plot style figures
stat_mean_tick   <- stat_summary(fun = mean, geom = "point", size = 10, shape = 95)
stat_se_errorbar <- stat_summary(fun.data = mean_se, geom = "errorbar",
                                 width = 0.1, size = 0.5)

# ==============================================================================
# SECTION 1 Raw outcome scores (RowSum)
# ==============================================================================

AGS_Plotting <- AGS_Plotting_org %>%
  mutate(
    RowSum  = rowSums(across(-c(Filename.Coding, ONSET..sec., Treatment))),
    Seized  = as.integer(RowSum > 0)
  ) %>%
  add_sex()

AGS_STATS <- AGS_Plotting %>%
  dplyr::select(Treatment, RowSum, Sex) %>%
  split_treatment()

# Stats
aov_outcomes <- aov(RowSum ~ Genotype * Treat, data = AGS_STATS)
summary(aov_outcomes)
emmeans(aov_outcomes, pairwise ~ Genotype * Treat)
AGS_STATS %>% filter(Genotype == "KO") %>% wilcox.test(RowSum ~ Treat, data = .)

# Plot
AGS_STATS_FIG <- AGS_STATS %>%
  make_treatment_factor() %>%
  ggplot(aes(x = Treatment, y = RowSum, color = Treatment, fill = Treatment)) +
  theme_cowplot() +
  stat_mean_tick + stat_se_errorbar +
  geom_jitter(aes(shape = factor(Sex)), width = 0.25, height = 0, size = 2) +
  ags_color_scale + ags_fill_scale + ags_shape_scale +
  ylab('Outcomes') + ags_theme_clean

ggsave("AGS_STATS_FIG_OUTCOMES.pdf", plot = AGS_STATS_FIG,
       width = 3, height = 3, dpi = 1000)

# ==============================================================================
# SECTION 2 — Osterweil incidence (Seized binary)
# ==============================================================================

AGS_STATS_OSTERWEIL <- AGS_Plotting_org %>%
  mutate(
    RowSum = rowSums(across(-c(Filename.Coding, ONSET..sec., Treatment))),
    Seized = as.integer(RowSum > 0)
  ) %>%
  add_sex() %>%
  dplyr::select(Treatment, Seized, Sex) %>%
  split_treatment()

Summarize_OSTERWEIL <- AGS_STATS_OSTERWEIL %>%
  group_by(Treat, Genotype) %>%
  summarize(
    Numofmice  = n(),
    Avg_Seized = mean(Seized),
    SE_Seized  = sd(Seized) / sqrt(n()),
    .groups    = "drop"
  )

write.csv(AGS_STATS_OSTERWEIL, file = 'AGS_STATS_OSTERWEIL_Claude.csv')

# Stats
aov_osterweil <- aov(Seized ~ Genotype * Treat, data = AGS_STATS_OSTERWEIL)
summary(aov_osterweil)
AGS_STATS_OSTERWEIL %>% filter(Genotype == "KO") %>% wilcox.test(Seized ~ Treat,    data = .)
AGS_STATS_OSTERWEIL %>% filter(Treat == "SALINE")  %>% wilcox.test(Seized ~ Genotype, data = .)
AGS_STATS_OSTERWEIL %>% filter(Genotype == "KO") %>% t.test(Seized ~ Treat,    data = .)
AGS_STATS_OSTERWEIL %>% filter(Treat == "SALINE")  %>% t.test(Seized ~ Genotype, data = .)

# Dot plot
AGS_BAR_FIG_2 <- AGS_STATS_OSTERWEIL %>%
  make_treatment_factor() %>%
  ggplot(aes(x = Treatment, y = Seized * 100, color = Treatment, fill = Treatment)) +
  theme_cowplot() +
  stat_mean_tick +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2, size = 0.5) +
  geom_jitter(aes(shape = factor(Sex)), width = 0.25, height = 0, size = 2, alpha = 0.5) +
  ags_color_scale + ags_fill_scale + ags_shape_scale +
  ylab('AGS Incidence (%)') + ags_theme_noaxis

ggsave("AGS_STATS_OSTERWEIL_DOTPLOT.pdf", plot = AGS_BAR_FIG_2,
       width = 1.5, height = 2, dpi = 1000)

# Bar plot
AGS_BAR_FIG <- AGS_STATS_OSTERWEIL %>%
  make_treatment_factor() %>%
  ggplot(aes(x = Treatment, y = Seized * 100, fill = Treatment)) +
  theme_cowplot() +
  stat_summary(fun = mean, geom = "bar", alpha = 0.5) +
  stat_se_errorbar +
  geom_jitter(aes(shape = factor(Sex), color = factor(Treatment)),
              width = 0.25, height = 0, size = 2, alpha = 0.5) +
  ags_fill_scale + ags_shape_scale + ags_color_scale +
  ylab('AGS Incidence') + ags_theme_noaxis

ggsave("AGS_STATS_OSTERWEIL_BAR.pdf", plot = AGS_BAR_FIG,
       width = 1.5, height = 2, dpi = 1000)

# ==============================================================================
# SECTION 3 Averaging method (unweighted RowAvg / RowSum)
# ==============================================================================

AGS_Plotting_MEAN <- AGS_Plotting_org %>%
  mutate(
    WR.1   = ifelse(X1.WR.1 == 1 | Multi.WR.2 == 1, 1, 0),   # FIX: or() -> |
    RowSum = rowSums(across(-c(Filename.Coding, X1.WR.1, Multi.WR.2, ONSET..sec., Treatment))),
    RowAvg = rowMeans(across(-c(Filename.Coding, X1.WR.1, Multi.WR.2, ONSET..sec., Treatment))),
    Seized = as.integer(RowSum > 0)
  ) %>%
  add_sex()

write.csv(AGS_Plotting_MEAN, "AGS_PLOTTING_MEAN.csv")

AGS_Plotting_UNWEIGHTED <- AGS_Plotting_MEAN %>%
  dplyr::select(Treatment, RowAvg, RowSum, Sex) %>%
  split_treatment()

# Stats
AGS_Plotting_UNWEIGHTED %>% filter(Genotype == "KO")    %>% wilcox.test(RowAvg ~ Treat,    data = .)
AGS_Plotting_UNWEIGHTED %>% filter(Genotype == "KO")    %>% wilcox.test(RowSum ~ Treat,    data = .)
AGS_Plotting_UNWEIGHTED %>% filter(Treat == "SALINE")   %>% wilcox.test(RowAvg ~ Genotype, data = .)
AGS_Plotting_UNWEIGHTED %>% filter(Treat == "SALINE")   %>% wilcox.test(RowSum ~ Genotype, data = .)

Summarize_UNWEIGHTED <- AGS_Plotting_UNWEIGHTED %>%
  group_by(Treat, Genotype) %>%
  summarize(Numofmice = n(), Avg_RowSum = mean(RowSum), Avg_RowAvg = mean(RowAvg),
            .groups = "drop")

aov_unweighted_mean <- aov(RowAvg ~ Genotype * Treat, data = AGS_Plotting_UNWEIGHTED)
summary(aov_unweighted_mean)
emmeans(aov_unweighted_mean, pairwise ~ Genotype * Treat)

# Plots use local Genotype/Treat columns, not AGS_Plotting_Long$ (FIX)
AGS_UNWEIGHTED_MEAN <- AGS_Plotting_UNWEIGHTED %>%
  make_treatment_factor() %>%
  ggplot(aes(x = Treatment, y = RowAvg, color = Treatment, fill = Treatment)) +
  theme_cowplot() +
  stat_mean_tick + stat_se_errorbar +
  geom_jitter(aes(shape = Sex), height = 0, alpha = 0.5, size = 2) +
  ags_color_scale + ags_fill_scale + ags_shape_scale +
  ylab('Unweighted AGS (Average score)') +
  theme(panel.grid = element_blank(), legend.position = "none",
        axis.title.x = element_blank())

ggsave("AGS_ROWAVG.pdf", plot = AGS_UNWEIGHTED_MEAN,
       width = 3, height = 3, dpi = 1000)

AGS_UNWEIGHTED_SUM <- AGS_Plotting_UNWEIGHTED %>%
  make_treatment_factor() %>%
  ggplot(aes(x = Treatment, y = RowSum, color = Treatment, fill = Treatment)) +
  theme_cowplot() +
  stat_mean_tick + stat_se_errorbar +
  geom_jitter(aes(shape = Sex), height = 0, alpha = 0.5, size = 2) +
  ags_color_scale + ags_fill_scale + ags_shape_scale +
  ylab('Unweighted AGS (Max score)') +
  theme(panel.grid = element_blank(), legend.position = "none",
        axis.title.x = element_blank())

ggsave("AGS_ROWSUM.pdf", plot = AGS_UNWEIGHTED_SUM,
       width = 3, height = 3, dpi = 1000)

# ==============================================================================
# SECTION 4 Racine scoring (CPC scheme: WR=1, Myo=2, TC=3, Death=4)
# ==============================================================================


AGS_Plotting_4_Racine_CPC <- AGS_Plotting_org %>%
  mutate(
    WR.1   = ifelse(X1.WR.1 == 1 | Multi.WR.2 == 1, 1, 0),   # FIX: or() -> |
    RowSum = rowSums(across(-c(Filename.Coding, ONSET..sec., Treatment))),
    NoSeize = as.integer(RowSum == 0)
  ) %>%
  dplyr::select(NoSeize, WR.1, Myoclonic.3, Tonic.Clonic.4, Death.5,
                Treatment, Filename.Coding) %>%
  mutate(Rac_Sex = case_when(
    grepl('_sexM', Filename.Coding) | grepl('Male',   Filename.Coding) ~ 'Male',
    grepl('_sexF', Filename.Coding) | grepl('Female', Filename.Coding) ~ 'Female'
  )) %>%
  dplyr::select(-Filename.Coding)

AGS_Plotting_4_Racine_SUM <- AGS_Plotting_4_Racine_CPC %>%
  mutate(
    Rac_NoSeize    = NoSeize       * 0,
    Rac_WR_1       = WR.1          * 1,
    Rac_Myoclonic  = Myoclonic.3   * 2,
    Rac_TonicClonic= Tonic.Clonic.4* 3,
    Rac_Death      = Death.5       * 4,
    Rac_Treatment  = str_trim(Treatment),
    Rac_Sex        = Rac_Sex
  ) %>%
  dplyr::select(starts_with('Rac')) %>%
  mutate(
    RowMax  = pmax(Rac_NoSeize, Rac_WR_1, Rac_Myoclonic, Rac_TonicClonic,
                  Rac_Death, na.rm = TRUE),
    RowMean = rowMeans(across(-c(Rac_Treatment, Rac_Sex))),
    RowSum  = rowSums(across(-c(Rac_Treatment, Rac_Sex)))
  ) %>%
  separate_wider_delim(cols = Rac_Treatment, delim = " ",
                       names = c("Genotype", "Treat"))

AGS_Plotting_Long <- AGS_Plotting_4_Racine_SUM %>%
  dplyr::select(Treat, Genotype, RowMax, RowMean, RowSum, Rac_Sex)

# Stats
AGS_Plotting_Long %>% filter(Genotype == "KO")   %>% wilcox.test(RowSum  ~ Treat,    data = .)
AGS_Plotting_Long %>% filter(Genotype == "KO")   %>% wilcox.test(RowMax  ~ Treat,    data = .)
AGS_Plotting_Long %>% filter(Genotype == "KO")   %>% wilcox.test(RowMean ~ Treat,    data = .)
AGS_Plotting_Long %>% filter(Treat == "SALINE")  %>% wilcox.test(RowSum  ~ Genotype, data = .)
AGS_Plotting_Long %>% filter(Treat == "SALINE")  %>% wilcox.test(RowMax  ~ Genotype, data = .)

# Summary (renamed to avoid overwriting later objects)
AGS_SUMMARY_RACINE_CPC <- AGS_Plotting_Long %>%
  group_by(Treat, Genotype) %>%
  summarize(RowMax_Avg = mean(RowMax), RowMean_Avg = mean(RowMean),
            RowSum_Avg = mean(RowSum), .groups = "drop")

# Plots use local columns 
AGS_RACINE_SUM <- AGS_Plotting_Long %>%
  make_treatment_factor() %>%
  ggplot(aes(x = Treatment, y = RowSum, color = Treatment, fill = Treatment)) +
  theme_cowplot() +
  stat_mean_tick + stat_se_errorbar +
  geom_jitter(aes(shape = Rac_Sex), height = 0, alpha = 0.5, size = 2) +
  ags_color_scale + ags_fill_scale + ags_shape_scale +
  ylab('Total score (Racine Scores)') +
  theme(panel.grid = element_blank(), legend.position = "none",
        axis.title.x = element_blank())

ggsave("AGS_RACINE_SUM.pdf", plot = AGS_RACINE_SUM,
       width = 3, height = 3, dpi = 1000)

AGS_RACINE_MAX <- AGS_Plotting_Long %>%
  make_treatment_factor() %>%
  ggplot(aes(x = Treatment, y = RowMax, color = Treatment, fill = Treatment)) +
  theme_cowplot() +
  stat_mean_tick + stat_se_errorbar +
  geom_jitter(aes(shape = Rac_Sex), height = 0, alpha = 0.5, size = 2) +
  ags_color_scale + ags_fill_scale + ags_shape_scale +
  ylab('Max score (Racine Scores)') +
  theme(panel.grid = element_blank(), legend.position = "none",
        axis.title.x = element_blank())

ggsave("AGS_RACINE_MAX.pdf", plot = AGS_RACINE_MAX,
       width = 3, height = 3, dpi = 1000)

AGS_RACINE_MEAN <- AGS_Plotting_Long %>%
  make_treatment_factor() %>%
  ggplot(aes(x = Treatment, y = RowMean, color = Treatment, fill = Treatment)) +
  theme_cowplot() +
  stat_mean_tick + stat_se_errorbar +
  geom_jitter(aes(shape = Rac_Sex), height = 0, alpha = 0.5, size = 2) +
  ags_color_scale + ags_fill_scale + ags_shape_scale +
  ylab('Mean score (Racine Scores)') +
  theme(panel.grid = element_blank(), legend.position = "none",
        axis.title.x = element_blank())

ggsave("AGS_RACINE_MEAN.pdf", plot = AGS_RACINE_MEAN,
       width = 3, height = 3, dpi = 1000)

# ==============================================================================
# SECTION 5 Racine scoring (standard: WR=1, MultiWR=2, Myo=3, TC=4, Death=5)
# ==============================================================================

AGS_Plotting_4_Racine <- AGS_Plotting_org %>%
  mutate(
    RowSum  = rowSums(across(-c(Filename.Coding, ONSET..sec., Treatment))),
    NoSeize = as.integer(RowSum == 0)
  )

AGS_Plotting_with_Racine <- AGS_Plotting_4_Racine %>%
  mutate(
    Rac_NoSeize     = NoSeize        * 0,
    Rac_WR_1        = X1.WR.1        * 1,
    Rac_Multi_WR    = Multi.WR.2     * 2,
    Rac_Myoclonic   = Myoclonic.3    * 3,
    Rac_TonicClonic = Tonic.Clonic.4 * 4,
    Rac_Death       = Death.5        * 5,
    RAC_Treatment   = str_trim(Treatment)
  ) %>%
  dplyr::select(starts_with('Rac')) %>%
  separate_wider_delim(cols = RAC_Treatment, delim = " ",
                       names = c("Genotype", "Treat"))

AGS_STATS_RACINE <- AGS_Plotting_with_Racine %>%
  mutate(RowSum = rowSums(across(-c(Genotype, Treat)))) %>%
  dplyr::select(Treat, Genotype, RowSum)

write.csv(AGS_STATS_RACINE, 'AGS_RACINE_STATS_TO_UPLOAD.csv')

# Stats
AGS_STATS_RACINE %>% filter(Genotype == "KO")  %>% wilcox.test(RowSum ~ Treat,    data = .)
AGS_STATS_RACINE %>% filter(Treat == "SALINE") %>% wilcox.test(RowSum ~ Genotype, data = .)
AGS_STATS_RACINE %>% filter(Genotype == "KO")  %>% t.test(RowSum ~ Treat,         data = .)

aov_racine <- aov(RowSum ~ Genotype * Treat, data = AGS_STATS_RACINE)
summary(aov_racine)
emmeans(aov_racine, pairwise ~ Genotype * Treat)

AGS_STATS_RACINE %>%
  make_treatment_factor() %>%
  ggplot(aes(x = Treatment, y = RowSum, color = Treatment, fill = Treatment)) +
  theme_cowplot() +
  stat_mean_tick + stat_se_errorbar +
  geom_jitter() +
  ags_color_scale + ags_fill_scale +
  ylab('Racine Scores')

# Modified Racine (WR scores zeroed out)
AGS_Plotting_with_Racine_Mod <- AGS_Plotting_4_Racine %>%
  mutate(
    Rac_NoSeize     = NoSeize        * 0,
    Rac_WR_1        = X1.WR.1        * 0,
    Rac_Multi_WR    = Multi.WR.2     * 0,
    Rac_Myoclonic   = Myoclonic.3    * 3,
    Rac_TonicClonic = Tonic.Clonic.4 * 4,
    Rac_Death       = Death.5        * 5,
    RAC_Treatment   = str_trim(Treatment)
  ) %>%
  dplyr::select(starts_with('Rac')) %>%
  separate_wider_delim(cols = RAC_Treatment, delim = " ",
                       names = c("Genotype", "Treat"))

AGS_STATS_RACINE_MOD <- AGS_Plotting_with_Racine_Mod %>%
  mutate(RowSum = rowSums(across(-c(Genotype, Treat)))) %>%
  dplyr::select(Treat, Genotype, RowSum)

AGS_STATS_RACINE_MOD %>% filter(Genotype == "KO") %>% wilcox.test(RowSum ~ Treat, data = .)


AGS_STATS_RACINE_MOD %>%
  make_treatment_factor() %>%
  ggplot(aes(x = Treatment, y = RowSum, color = Treatment, fill = Treatment)) +
  theme_cowplot() +
  stat_mean_tick + stat_se_errorbar +
  geom_jitter() +
  ags_color_scale + ags_fill_scale +
  ylab('Racine Scores (Modified)')

# ==============================================================================
# SECTION 6 All-behavior stacked / proportion plots (normalized by N mice)
# ==============================================================================

# Renamed AGS_SUMMARY -> AGS_SUMMARY_BEHAVIOR_RAW (FIX: avoid silent overwrite)
AGS_SUMMARY_BEHAVIOR_RAW <- AGS_Plotting_4_Racine %>%
  mutate(Treatment = str_trim(Treatment)) %>%
  group_by(Treatment) %>%
  summarize(
    Num_mice      = n(),
    Wildrun       = sum(X1.WR.1),
    Wildrun.multi = sum(Multi.WR.2),
    Myoclonic     = sum(Myoclonic.3),
    TonClon       = sum(Tonic.Clonic.4),
    Death         = sum(Death.5),
    NoSeize       = sum(NoSeize),
    .groups       = "drop"
  )

AGS_Proportion_Raw <- AGS_SUMMARY_BEHAVIOR_RAW %>%
  mutate(across(-c(Treatment, Num_mice), ~ . / Num_mice))

AGS_Proportion_Long_Raw <- AGS_Proportion_Raw %>%
  pivot_longer(cols = -c(Treatment, Num_mice),
               values_to = 'Proportion', names_to = 'Behaviors')

# Stats
AGS_Proportion_Long_Raw %>%
  filter(Behaviors == "NoSeize") %>%   # adjust to the behavior of interest
  mutate(Treatment = str_trim(Treatment)) %>%
  separate_wider_delim(cols = Treatment, delim = " ", names = c("Genotype", "Treat")) %>%
  filter(Genotype == "KO")  # note: kept from original; check variable name

# Dodge bar plot
BEHAVIOR_LEVELS <- c("Death", "TonClon", "Myoclonic", "Wildrun.multi", "Wildrun", "NoSeize")

prep_behavior_long <- function(df) {
  df %>%
    filter(!Behaviors %in% c('Sumofbehavior', 'Num_mice')) %>%
    mutate(
      Treatment = factor(Treatment, levels = TREATMENT_LEVELS_SPACE),
      Behaviors = factor(Behaviors, levels = BEHAVIOR_LEVELS)
    )
}

AGS_Proportion_Normalized_2_mice_DODGE <- ggplot(
  prep_behavior_long(AGS_Proportion_Long_Raw),
  aes(x = Treatment, y = Proportion, fill = Behaviors)
) +
  geom_col(position = "dodge") +
  theme_cowplot() +
  theme(panel.grid = element_blank(), legend.position = "none",
        axis.title.x = element_blank(), axis.title.y = element_blank()) +
  behavior_fill_scale

ggsave("AGS_Proportion_Normalized_2_mice_DODGE.pdf",
       plot = AGS_Proportion_Normalized_2_mice_DODGE,
       width = 3, height = 3, dpi = 1000)

AGS_Proportion_Normalized_2_mice_STACK <- ggplot(
  prep_behavior_long(AGS_Proportion_Long_Raw),
  aes(x = Treatment, y = Proportion, fill = Behaviors)
) +
  geom_col(position = "stack") +    # FIX: was "stacked"
  theme_cowplot() +
  theme(panel.grid = element_blank(), legend.position = "none",
        axis.title.x = element_blank(), axis.title.y = element_blank()) +
  behavior_fill_scale

ggsave("AGS_Proportion_Stacked_2_mice_DODGE.pdf",
       plot = AGS_Proportion_Normalized_2_mice_STACK,
       width = 3, height = 3, dpi = 1000)

# Not-normalised stacked (counts per mouse)
AGS_SUMMARY_Long_Raw <- AGS_SUMMARY_BEHAVIOR_RAW %>%
  pivot_longer(cols = -Treatment, values_to = 'Proportion', names_to = 'Behaviors')

AGS_Proportion_Not_Normalized_BYMICE <- ggplot(
  AGS_SUMMARY_Long_Raw %>%
    filter(Behaviors != 'Num_mice') %>%
    mutate(
      Treatment = factor(Treatment, levels = TREATMENT_LEVELS_SPACE),
      Behaviors = factor(Behaviors, levels = BEHAVIOR_LEVELS)
    ),
  aes(x = Treatment, y = Proportion, fill = Behaviors)
) +
  geom_col(position = "stack", color = 'white', linewidth = 0.5) +
  theme_cowplot() +
  theme(panel.grid = element_blank(), legend.position = "none",
        axis.title.x = element_blank(), axis.title.y = element_blank()) +
  behavior_fill_scale

ggsave("AGS_Proportion_Not_Normalized_BYMICE.pdf",
       plot = AGS_Proportion_Not_Normalized_BYMICE,
       width = 3, height = 3, dpi = 1000)

# ==============================================================================
# SECTION 7 Normalized behavior (proportion of total behavior events)
# ==============================================================================

# Renamed AGS_SUMMARY -> AGS_SUMMARY_BEHAVIOR_NORM (FIX: avoid silent overwrite)
AGS_SUMMARY_BEHAVIOR_NORM <- AGS_Plotting_4_Racine %>%
  group_by(Treatment) %>%
  summarize(
    Wildrun       = sum(X1.WR.1),
    Wildrun.multi = sum(Multi.WR.2),
    Myoclonic     = sum(Myoclonic.3),
    TonClon       = sum(Tonic.Clonic.4),
    Death         = sum(Death.5),
    NoSeize       = sum(NoSeize),
    .groups       = "drop"
  ) %>%
  mutate(Sumofbehavior = rowSums(across(-Treatment)))

AGS_Proportion_Norm <- AGS_SUMMARY_BEHAVIOR_NORM %>%
  mutate(across(-c(Treatment, Sumofbehavior), ~ . / Sumofbehavior)) %>%
  mutate(No_Seizure = 1 - rowSums(across(-c(Treatment, Sumofbehavior))))

AGS_Proportion_Long_Norm <- AGS_Proportion_Norm %>%
  pivot_longer(cols = -c(Treatment, Sumofbehavior),
               values_to = 'Proportion', names_to = 'Behaviors')

# Quick exploratory stacked bar
ggplot(AGS_Proportion_Long_Norm,
       aes(x = Treatment, y = Proportion, fill = Behaviors)) +
  geom_col(position = "stack", color = 'white') +
  theme_cowplot() +
  behavior_fill_scale

# Donut charts (normalized)
DONUTS_NORMALIZED_BEHAVIOR <- ggplot(
  AGS_Proportion_Long_Norm %>%
    mutate(
      Treatment = factor(Treatment, levels = TREATMENT_LEVELS_SPACE),
      Behaviors = factor(Behaviors, levels = BEHAVIOR_LEVELS)
    ),
  aes(x = 1, y = Proportion, fill = Behaviors)
) +
  geom_col(width = 0.5, color = 'white') +
  coord_polar(theta = "y") +
  xlim(0.5, 1.5) +
  facet_wrap(~Treatment, scales = "free") +
  theme_cowplot() +
  theme(
    axis.text = element_blank(), axis.title = element_blank(),
    panel.grid = element_blank(), strip.background = element_blank(),
    axis.ticks = element_blank(), axis.line = element_blank(),
    legend.position = 'none'
  ) +
  behavior_fill_scale

ggsave("AGS_DONUTS_NORMALIZED_BEHAVIOR.pdf", plot = DONUTS_NORMALIZED_BEHAVIOR,
       width = 3, height = 3, dpi = 1000)

# Wide-hole donut variant
ggplot(
  AGS_Proportion_Long_Norm %>%
    mutate(
      Treatment = factor(Treatment, levels = TREATMENT_LEVELS_SPACE),
      Behaviors = factor(Behaviors, levels = BEHAVIOR_LEVELS)
    ),
  aes(x = 2, y = Proportion, fill = Behaviors)
) +
  geom_col(position = "stack", color = 'white', linewidth = 0.5) +
  coord_polar(theta = "y") +
  facet_wrap(~Treatment) +
  xlim(0.5, 2.5) +
  theme_cowplot() +
  behavior_fill_scale +
  theme(
    axis.title = element_blank(), axis.text  = element_blank(),
    axis.ticks = element_blank(), axis.line  = element_blank(),
    strip.background = element_blank()
  )
