# Introduction ------------------------------------------------------------

# DISCLAIMER:
# This script was used to analyse .csv files output directly from the Zantiks
# MWP unit. Whilst it illustrates how R can be used to process data from the
# Zantiks unit, Zantiks Ltd. cannot guarantee this is how you want to run your
# experiments or analyse your data.

# The script below is for analysing adult Drosophila light-off startle data
# from the Zantiks Ltd MWP unit.

# To run this script you need the appropriate packages installed under 'Prerequisites'.

# Distance travelled is the main measure investigated. Activity is a measure
# being developed at the time of writing and relates to the pixel change overall
# in the video image.

# Prerequisites -----------------------------------------------------------

#install.packages("tidyverse")
library(tidyverse)
#install.packages("sjlabelled")
library(sjlabelled) # for data manipulation
#install.packages("ggpubr")
library(ggpubr) # for box plots
#install.packages("rstatix")
library(rstatix) # for statistics
#install.packages("car")
library(car) #for anova (aov) function
#install.packages("RColorBrewer")
library(RColorBrewer) # for colour palettes


# Data Loader -------------------------------------------------------------

# setting data directory
DATA_DIR="" #specify directory of zantiks file
setwd(DATA_DIR)
getwd()

#create list of files from directory
file_list <- list.files(pattern = "*.csv")

#create header from first file
df<-
  file_list[1] %>%
  read_csv(skip=3,col_names = TRUE, guess_max = 100) %>%
  head(0)

#create new list without demographic info
new_list<- c()

for (i in file_list){
  new_list[[i]]<-
    read_csv(i,skip=3, col_names = TRUE, guess_max = 100) %>%
    head(-1)
}

#append all files to df
for (i in new_list){
  df<-add_row(df,i)
}


# Data Formatting ---------------------------------------------------------

#convert variables to factors
df<-as_factor(df,BLOCK)
df<-as_factor(df,TYPE)
df<-as_factor(df,TIMESLOT)
df<-as_factor(df,UNIT)
df<-as_factor(df,FLASH_LENGTH)

#create file with total population columns
pop_data<-mutate(df,TOTAL_DIST = rowSums(select(df,!ends_with("MSD"), -RUNTIME,
                                                -UNIT, -TIMESLOT, -PLATE_ID, -TEMPERATURE, -TIME_BIN, -BLOCK,
                                                -TRIAL, -TYPE, -PRE_POST_COUNTER, -FLASH_LENGTH)))
pop_data<-mutate(pop_data,TOTAL_ACT = rowSums(select(df,ends_with("MSD"))))

#create file with well factor and distance dv
dfile_dist<-
  df %>% 
  select(!ends_with("MSD")) %>%
  gather(key = "WELL", value = "DISTANCE", -RUNTIME,
         -UNIT, -TIMESLOT, -PLATE_ID, -TEMPERATURE, -TIME_BIN, -BLOCK,
         -TRIAL, -TYPE, -PRE_POST_COUNTER, -FLASH_LENGTH) %>%
  convert_as_factor(WELL)

#create file with well factor and MSD activity dv only
dfile_act<-
  df %>% 
  select(ends_with("MSD")) %>%
  gather(key = "WELL", value = "ACTIVITY") %>%
  convert_as_factor(WELL)

#remove duplicate well variable before adding activity data
dfile_act<-select(dfile_act, -'WELL')

#add activity column to rest of data
df<-
  add_column(dfile_dist, dfile_act)

#create data file without acclimation data
no_acclimation<-
  df %>%
  filter(!(TYPE == "ACCLIMATION"))

#create data file with only startles and the pre-startle second
single_pre<-filter(df,TYPE=="STARTLE"|TYPE=="PRE" & PRE_POST_COUNTER==33)
pop_single_pre<-filter(pop_data,TYPE=="STARTLE"|TYPE=="PRE" & PRE_POST_COUNTER==33)

#create data file with only startles
startles_only<-filter(df, TYPE=="STARTLE")
pop_startles_only<-filter(pop_data,TYPE=="STARTLE")


# Graphs ------------------------------------------------------------------

### Distance data ###

#box plot - dist / light-off duration
dist_bxp <- ggboxplot(
  single_pre, x = "FLASH_LENGTH", y = "DISTANCE") +
  xlab("Light-off duration") +
  ylab("Distance travelled (pixels)")
dist_bxp

#box plot - population dist / light-off duration
pop_dist_bxp <- ggboxplot(
  pop_single_pre, x = "FLASH_LENGTH", y = "TOTAL_DIST") +
  xlab("Light-off duration") +
  ylab("Population Distance travelled (pixels)")
pop_dist_bxp

### Activity data ###

#box plot - act / light-off duration
act_bxp <- ggboxplot(
  single_pre, x = "FLASH_LENGTH", y = "ACTIVITY") +
  xlab("Light-off duration") +
  ylab("Activity (MSD)")
act_bxp

#box plot - population act / light-off duration
act_bxp <- ggboxplot(
  pop_single_pre, x = "FLASH_LENGTH", y = "TOTAL_ACT") +
  xlab("Light-off duration") +
  ylab("Population Activity (MSD)")
act_bxp


# Neat Graphs -------------------------------------------------------------

### Distance data ###

#pairwise corrected (bonferroni) p values for boxplot
pwise_comp <- pop_startles_only %>%
  pairwise_t_test(TOTAL_DIST ~ FLASH_LENGTH) %>%
  adjust_pvalue(method = "bonferroni") %>%
  mutate(y.position = c(270,285,0,300,0,0,315,0,0,0))

# boxplot with sig markers
pop_dist_bxp <- ggboxplot(
  pop_single_pre, x = "FLASH_LENGTH", y = "TOTAL_DIST",
  xlab = "Light-off Duration (ms)", ylab = "Population Distance Travelled (pixels)") +
  stat_pvalue_manual(slice(pwise_comp,c(1,2,4,7)), label = "p.adj.signif") +
  stat_compare_means(method = "anova", label.y = 300)
pop_dist_bxp


# Distance Analysis -------------------------------------------------------

#means and sd
dist_means<-
  no_acclimation %>%
  group_by(FLASH_LENGTH) %>%
  get_summary_stats(DISTANCE, type = "mean_sd")
dist_means

#outliers
dist_outliers<-
  startles_only %>%
  group_by(FLASH_LENGTH) %>%
  identify_outliers(DISTANCE)

#normality test
dist_normality<-
  startles_only %>%
  group_by(FLASH_LENGTH) %>%
  shapiro_test(DISTANCE)

#normality plot
ggqqplot(startles_only, "DISTANCE", ggtheme = theme_bw()) +
  facet_wrap(startles_only$FLASH_LENGTH)

#anova
dist_anova <- aov(DISTANCE ~ FLASH_LENGTH, data = startles_only)
Anova(dist_anova, type = "II")

### Population data ###

#means and sd
pop_dist_means<-
  pop_single_pre %>%
  group_by(FLASH_LENGTH) %>%
  get_summary_stats(TOTAL_DIST, type = "mean_sd")
pop_dist_means

#outliers
pop_dist_outliers<-
  pop_startles_only %>%
  group_by(FLASH_LENGTH) %>%
  identify_outliers(TOTAL_DIST)

#normality test
pop_dist_normality<-
  pop_startles_only %>%
  group_by(FLASH_LENGTH) %>%
  shapiro_test(TOTAL_DIST)

#normality plot
ggqqplot(pop_startles_only, "TOTAL_DIST", ggtheme = theme_bw()) +
  facet_wrap(pop_startles_only$FLASH_LENGTH)

#anova
act_anova <- aov(ACTIVITY ~ FLASH_LENGTH, data = startles_only)
Anova(act_anova, type = "II")


# Activity Analysis -------------------------------------------------------

#means and sd
act_means<-
  no_acclimation %>%
  group_by(TYPE, UNIT, TIMESLOT) %>%
  get_summary_stats(ACTIVITY, type = "mean_sd")

#outliers
act_outliers<-
  no_acclimation %>%
  group_by(TYPE, UNIT, TIMESLOT) %>%
  identify_outliers(ACTIVITY)

#normality test
act_normality<-
  startles_only %>%
  group_by(FLASH_LENGTH) %>%
  shapiro_test(ACTIVITY)

#normality plot
ggqqplot(no_acclimation, "ACTIVITY", ggtheme = theme_bw()) +
  facet_grid(UNIT + TIMESLOT ~ TYPE , labeller = "label_both")

#anova
act_anova <- aov(ACTIVITY ~ TYPE * UNIT * TIMESLOT, data = no_acclimation)
Anova(act_anova, type = "III")