Cal_Clean.qmd

---
title: "Cal_Clean: A Script for Tidying Sable Indirect Calorimetry Data"
author: "Chelsea Faber"
format: 
   html:
     df-print: paged
     embed-resources: true
editor: visual
---

```{r}
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
```

## About

###### Inputs:

-   Date and run ID of calorimetry run
-   List of animal ID numbers to be excluded (NA if none)

###### Outputs:

-   .Rda file containing cleaned dataframe at original sampling rate (typically, 1-3 minutes), 1h-binned dataframe, and photoperiod-binned dataframes for further plotting and statistical analyses.

###### Requirements:

-   TimeSeries sheet from Macro-processed .exp data files (OneClickMacro \>= v2.52)
-   TimeSeries \<yyyy-mm-dd_run_m_timeseries.csv\>
-   Key \<yyyy-mm-dd_run_KEY.csv\> - must have same prefix as TimeSeries sheet, and be saved in the same directory

## Load required packages

```{r}
#| label: load-packages
#| warning: false
#| echo: false
#| include: false
library(plyr, include.only = 'mapvalues')
library(tidyverse)
library(lubridate)
```

## Define inputs to script for analysis

If everything is working properly, these variables are the only ones a user should have to change to clean a given file's calorimetry data. Eventually, this entire script will be converted to a function, such that the user would type these in-line when calling that function (e.g., `Cal_Clean(cohort = 'mon001',rundate = '2021-10-18',remove_animals = 274)`).

fpath \<- "C:/Users/kaspe/Dropbox (Barrow Neurological Institute)/Mirzadeh Lab Dropbox MAIN/Data/Calorimetry/macro_processed"

```{r}
# Run-specific user-defined parameters
filename        <- "2023-12-11_cal023_m_timeseries.csv"
fpath           <- "C:/Users/cfaber/Barrow Neurological Institute Dropbox/Chelsea Faber/Mirzadeh Lab Dropbox MAIN/Data/Calorimetry/macro_processed"


# list of strings if animals to be removed. Set remove_animals <- NA if none to be removed.
remove_animals  <- NA 

# Default parameters. For diet, temperature, or other interventions, set the value of intervention equal to TRUE.
NCDgcal         <- 3.35 #default conversion 
HFDgcal         <- 5.47 # default conversion
intervention    <- FALSE
segment.run     <- FALSE
cols2excl       <- c('Still_pct_M','Sleep_pct_M','Mass_g','AllMeters_M')
```

## Load data

To load your data, please ensure that it is located in the repository that the `here()` function points to. Otherwise, you will get an error.

A workaround, if you do not wish to load your data using the here() package, is to type the entire file path in the call to `read_csv()`.

```{r}
#| warning: false
#| code-overflow: wrap
#| include: false

# Confirm file exists
if (!file.exists(paste0(fpath,"/",filename))) {
  stop("File not found! Please confirm file name and path are correct.")
} else {
  print(paste0("Found ", filename, ", extracting data."))
}

# This section assumes that files are named following yyyy-mm-dd_runID... .csv convention. If no underscores detected in filename, everything preceding .csv is taken as the name for all saved files.
fileparts       <- unlist(strsplit(filename,'[_.]+'))

if (length(fileparts) == 1) {
  savename <- unlist(strsplit(filename,'.csv'))
  code     <- paste(savename,"KEY.csv",sep = "_")
  
  if (!file.exists(paste0(fpath,"/",code))) {
      stop("Key file not found! Please confirm file name and path are correct.")
  }
  savename <- paste(fpath,"/r_cleaned/",savename,sep="")
} else {
  rundate  <- fileparts[1]
  runid    <- fileparts[2]
  savename <- paste(rundate,runid,"Clean.Rda",sep = "_")
  code     <- paste(rundate,runid,"KEY.csv",sep = "_")
  
  # Confirm key.csv file exists.
  if (!file.exists(paste0(fpath,"/",code))) {
      stop("Key file not found! Please confirm file name and path are correct.")
  }
  
  savename <- paste(fpath,"/r_cleaned/",savename,sep="")
}

# Check whether macro-processed file has already been run through script.
if (file.exists(savename)) {
  go <- menu(c("Re-run and overwrite existing Clean.Rda file","Re-run and save a new copy of Clean.Rda","Abort! Abort!"), title = "A .Rda file for this file already exists. How would you like to proceed?")
} else {
  go <- 0
}

if (go == 1) {
  print(paste0("Re-running Cal_Clean script and over-writing previous file at ", savename))
}else if (go == 2) {
  savename <- paste(rundate,cohort,"Clean_COPY.Rda",sep = "_")
  savename <- paste0(fpath,"/r_cleaned/",savename)
  print(paste0("A new Clean.Rda file will be generated at ", savename))
  print("Be sure to revise filename inputs to subsequent scripts to reflect this new naming convention.")
} else if (go == 3) {
      print("Aborting script, Clean.Rda file for this run already exists.")
      rm(list = ls())
}

# Load csv with run metadata
key             <- read_csv(paste(fpath,code,sep="/"),show_col_types=F) 

# Load cal.csv and merge metadata
df              <- read_csv(paste(fpath,filename,sep="/"),show_col_types = F) %>%
  left_join(key, by = "Animal") %>%                 #unblind by merging with decoding df 
  rename(Cage = Animal,
         Animal = ID_Code)%>%
  select(!(starts_with("Ped") 
           | all_of(cols2excl))) %>%
  mutate(across(c(Animal,Sex,Cage,Group,Cohort,Treatment), factor))
```

## Tidy data

Here, we remove any animals that we had previously specified, and perform some general data cleaning - most importantly, extraction of time of day information from the DateTime column.

```{r}
#| warning: false

# Remove dead/sick animals; if NA, don't run
# if (!anyNA(remove_animals)) {
#    df <- df[!df$Animal %in% remove_animals,]
# }

# Remove animals marked for removal from dataset within the _KEY.csv
if (all(key$Remove == "N") | all(is.na(key$Remove))) {
  print("No animals marked for removal. Keeping all animals in df")
} else if (any(key$Remove == "Y")) {
  
  print("Checking KEY for animals to remove from analysis")
  ct <- key %>% 
    count(Remove) %>% filter(Remove == "Y") %>% select(n) %>% pull()
  
  # Print how many and which animals are being removed
  print(paste0("Removing n = ", as.character(ct), " animals"))

  idx <- key %>%
    filter(Remove == "Y") %>% select(ID_Code) %>% pull()
  for (i in 1:length(idx)) {
  print(paste0("Removing subject ", idx[[i]]))
  }
  
  # Remove the specified animals
  df <- df %>%
  group_by(Animal) %>%
  filter(Remove == "N") %>%
  ungroup()
  print('Subject(s) removed successfully.')
  rm(list = c("ct","idx"))
}

# Parse date-times and sort by animal ID (ID_Code) and date-time
df <- df %>% 
  relocate(c(DateTime,Sex,Group,Treatment,Cohort,Animal)) %>%
  mutate(DateTime = round_date(ymd_hms(df$DateTime),unit = "minute"),
         month = as.numeric(month(DateTime)),
         day = as.numeric(day(DateTime)),
         hour = as.numeric(hour(DateTime)),
         minute = as.numeric(minute(DateTime))) %>%
  mutate(DietChangeDateTime = round_date(mdy_hm(df$DietChangeDateTime))) %>%
  mutate(LightChangeDateTime = round_date(mdy_hm(df$LightChangeDateTime)))

df <- df[order(df$DateTime),] %>%
  arrange(Animal)

# Remove any rows/columns with only NAs
df <- df[rowSums(is.na(df)) != ncol(df), ]
df <- df[, colSums(is.na(df)) != nrow(df)]
```

## Transform data

Here, we first filter the data so that all recordings begin and end at ZT12 (18:00). Next we compute a handful of new variables from the original DateTime, including Zeitgeber time, Photoperiod, and experimental day.

```{r}
# Extract original date-time information
og.start         <- head(df$DateTime, n=1)
og.end           <- last(df$DateTime)
og.dur           <- og.end - og.start

print(paste("Trimmming df to start and end at ZT18. Original recording length =", round(og.dur,2), "days."))

# Find index of first lights out (18:00) for each animal
df <- df %>% group_by(Animal)
idx.start        <- which(hour(df$DateTime) == 18)[1]-1
time.start       <- as.POSIXct(df$DateTime[idx.start])

# Trim everything before first lights out
df <- df %>%
  filter(DateTime > time.start)

# Find index of last lights out
idx.end          <- tail(which(hour(df$DateTime) == 17), n=1) + 1
time.end         <- as.POSIXct(df$DateTime[idx.end])

# Trim everything after last lights out
df <- df %>%
  filter(DateTime < time.end)

trim.dur         <- last(df$DateTime) - head(df$DateTime,n=1)
print(paste(round(og.dur - trim.dur,2), "days trimmed. Trimmed df length =", round(trim.dur,2), "days."))

# Compute ZT time, rename some variables
df <- df %>%
  mutate(ZT = plyr::mapvalues(hour, 
                              from = (0:23), 
                              to = c(18:23,0:17)),.before = Sex) %>% 
  mutate(Photoperiod = as.factor(
    plyr::mapvalues(ZT, 
                    from = c(0:23), 
                    to = c(rep("Light",12),rep("Dark",12)))),.before = ZT) %>%
  rename(VO2 = VO2_M,
         VCO2 = VCO2_M,
         EE = kcal_hr_M,
         RER = RER_M,
         FoodIn.cum = FoodInA_M,
         WaterIn.cum = WaterInA_M,
         AllMeters = AllMeters_R,
         BodyMass = BodyMass_Mnz,
         VH2O = VH2O_M)

# Compute experimental day (ExpDay) and elapsed time (Time) in hours
df <- df %>%
  mutate(ExpDay = as.integer(ceiling(difftime(DateTime,time.start, units = "days"))),
         Time = as.numeric(difftime(DateTime,time.start), units = "hours"),
         .after = DateTime) %>%
  ungroup()

# Calculate rolling Age, from the starting Age entered in the "KEY" metadata file.
if ("Age" %in% colnames(df)) {
  df <- df %>%
    group_by(ExpDay,Animal) %>% 
    mutate(Age = ExpDay-1 + first(Age)) %>% 
    ungroup()
} else {
  print("No Age values found in KEY If this is in error, please revised KEY sheet accordingly, save, and rerun script, otherwise, continue.")
}

# Recompute cumulative variables to begin at new start time
df <- df %>%
  group_by(Animal) %>% 
  mutate(FoodIn.cum = FoodIn.cum - first(FoodIn.cum),
         WaterIn.cum = WaterIn.cum - first(WaterIn.cum)) %>% 
  ungroup()
```

Next, we will compute some more parameters of interest, including binned food and water intake (from the original cumulative values), food intake in kcal (from grams), binned and cumulative energy expenditure (from the instantaneous estimate of hourly rate), cumulative distance traveled, and energy balance. The binned energy expenditure is questionable currently - please let me know if you find any issues.

```{r}
# Compute time interval (in hours) for estimating binned energy expenditure
int <- df$Time[2] - df$Time[1]

# Compute daily average for body mass
df <- df %>%
  group_by(Animal, ExpDay) %>%
  mutate(meanBodyMass = mean(BodyMass),.after = BodyMass) %>%
  ungroup()

# Detect diet from key and establish conversion
if (sum(is.na(key$Diet1))==length(key$Diet1) & sum(is.na(key$Diet2)) ==length(key$Diet2)) {
  
  diet <- menu(c("NCD","HFD"), title = "No diet type specified in Key file. Which diet was provided?")
  
  df <- df %>%
    mutate(Diet = case_when(
      diet == 1 ~ "NCD",
      diet == 2 ~ "HFD"))
} else if (sum(is.na(key$Diet2)) == length(key$Diet2)) {

  df <- df %>%
    group_by(Animal) %>%
    rename(Diet = Diet1) %>%
    ungroup()
  
} else {
  
    df <- df %>% 
    group_by(Animal) %>%
    mutate(Diet = case_when(
      DateTime <= DietChangeDateTime ~ Diet1,
      DateTime > DietChangeDateTime ~ Diet2)) %>%
    ungroup()
  
}
  
# Detect light schedule from key 
if (sum(is.na(key$LightCycle1))==length(key$LightCycle1) & sum(is.na(key$LightCycle2)) ==length(key$LightCycle2)) {
  
  lt <- menu(c("LD","DD","LL"), title = "No light cycle  specified in Key file. Which light schedule was used?")
  
  df <- df %>%
    mutate(LightCycle = case_when(
      lt == 1 ~ "LD",
      lt == 2 ~ "DD",
      lt == 3 ~ "LL"))
} else if (sum(is.na(key$LightCycle2)) == length(key$LightCycle2)) {

  df <- df %>%
    group_by(Animal) %>%
    rename(LightCycle = LightCycle1) %>%
    ungroup()
  
} else {
  
    df <- df %>% 
    group_by(Animal) %>%
    mutate(LightCycle = case_when(
      DateTime <= LightChangeDateTime ~ LightCycle1,
      DateTime > LightChangeDateTime ~ LightCycle2)) %>%
    ungroup()
  
}

# Compute new columns
df <- df %>%
  group_by(Animal) %>%
  mutate(FoodIn.g = c(diff(FoodIn.cum),0),.before = FoodIn.cum) %>% # convert cumulative to binned food intake
  mutate(WaterIn.g = c(diff(WaterIn.cum),0),.before = WaterIn.cum) %>% # convert cumulative to binned water intake
  mutate(FoodIn.kcal = case_when(
    Diet == "NCD" ~ FoodIn.g * NCDgcal,
    Diet == "HFD" ~ FoodIn.g * HFDgcal),
    .before = FoodIn.g) %>% # convert g to kcal
  mutate(FoodIn.cum.kcal = cumsum(FoodIn.kcal),.after = FoodIn.cum) %>%
  mutate(EE.kcal.bin = EE * int, .before = EE) %>% # EE is kcal/hr, multiply by int (x hours between samples) to get kcal/bin
  mutate(EE.cum = cumsum(EE.kcal.bin), .before = EE.kcal.bin) %>%
  mutate(EBalance = FoodIn.kcal - EE.kcal.bin, .before = VO2) %>% # compute energy balance per bin
  mutate(EB.cum = cumsum(EBalance), .after = EBalance) %>%
  mutate(AllMeters.cum = cumsum(AllMeters),.before = AllMeters) %>% # compute cumulative distance traveled
  ungroup()

# Extract environmental sensor data to separate dataframe
df.env <- df %>% 
  select(DateTime,ExpDay,Time,Cage,starts_with("Enviro"),hour) %>%
  filter(Cage == 1 | Cage == 9) %>% 
  mutate(across(starts_with("Enviro"),as.numeric))

# Trim last row with mostly NAs, created during some transformations above
#df <- head(df,-1)
```

## Bin to hourly

Here, we compute hourly bins for each variable. Depending on variable, this is done by either summing all values within the hour (`cols2sum`: binned, non-cumulative measures), taking the mean of all values within the hour (`cols2avg`: rates), taking the maximum value within the hour (`cols4cum`: binned, cumulative measures), and taking the median value for assigning the new bin time (`cols2med`: dates and times only).

```{r}
#/ warning: false

# DO NOT CHANGE WITHOUT EXPRESS PERMISSION FROM CHELSEA!!!!!!!!!!
cols2sum <- sort(c('FoodIn.g','FoodIn.kcal','WaterIn.g','EBalance','AllMeters','EE.kcal.bin'))
cols2avg <- sort(c('VO2','VCO2','VH2O','EE','RER','BodyMass','meanBodyMass','XBreak_R','YBreak_R'))
cols4cum <- sort(c('AllMeters.cum','FoodIn.cum','WaterIn.cum','EE.cum','EB.cum','FoodIn.cum.kcal'))
cols2med <- c('DateTime','Time','minute')

# Bin to hourly! 
df.hourly <- df %>%
  group_by(Animal,ExpDay,hour) %>%
  select(!starts_with("Enviro")) %>%
  mutate(
    across(all_of(cols2med),median)) %>% # assign middle of time bin to new bin
  mutate(across(
    all_of(cols2avg), mean)) %>% # rates get averaged
  mutate(across(
    all_of(cols2sum),sum)) %>% # intake, distances get summed
  mutate(across(
    all_of(cols4cum),last)) %>% # cumulative values just keep the maximum (total for the hour)
  ungroup() %>%
  distinct() %>% # squashes down to one observation per hour
  select(!c(hour,minute,month,day)) %>% 
  mutate(Animal = as.factor(Animal)) %>%
  group_by(Animal) %>%
  #slice(2:(n()-1)) %>% # trim incomplete hours at start and end
  mutate(DateTime = round_date(DateTime, "30 minutes")) %>% 
  mutate(Time = as.numeric(
    difftime(DateTime,first(DateTime)),units = "hours"),.after = DateTime) %>% # starts clock from 0 at start of recording
  ungroup()

# Removes duplicate row created at light transition when binning
if ("LightChangeDateTime" %in% colnames(df)) {
 df.hourly <- df.hourly %>% 
    group_by(Animal) %>%
    mutate(LightCycle = case_when(
      DateTime <= LightChangeDateTime ~ LightCycle1,
      DateTime > LightChangeDateTime ~ LightCycle2)) %>%
    ungroup() %>% 
    distinct() 
} 

df.env <- df.env %>%
  group_by(Cage,ExpDay,hour) %>%
  mutate(across(starts_with("Enviro"),mean)) %>%
  mutate(across(all_of(c("DateTime","Time")),median)) %>%
  ungroup() %>%
  distinct() %>%
  group_by(Cage) %>%
  mutate(DateTime = round_date(DateTime, "30 minutes")) %>% 
  mutate(Time = as.numeric(
    difftime(DateTime,first(DateTime)),units = "hours"),.after = DateTime) %>% # starts clock from 0 at start of recording
  ungroup()
```

## Photoperiod and daily means

Next, we compute the totals/averages of all variables within each experimental day.

```{r}
#/ warning: false
# Compute averages/totals for each day
df.dailytot <- df %>%
  group_by(Animal,ExpDay) %>% 
  summarize(
    across(all_of(cols2sum),sum),
    across(all_of(cols2avg),mean)) %>%
  mutate(Photoperiod = "Total") %>%
  ungroup() 

# Compute average within each photoperiod, append to total.avg.daily
df.daily <- df %>%
  group_by(Animal,ExpDay,Photoperiod) %>%
  summarize(
    across(all_of(cols2sum),sum),
    across(all_of(cols2avg),mean)) %>% 
  ungroup() %>%
  bind_rows(.,df.dailytot) %>% 
  left_join(key, by = c("Animal" = "ID_Code")) %>% 
  select(-Animal.y) %>%
  arrange(Animal,ExpDay,Photoperiod)
rm(df.dailytot)
```

Finally, we compute cumulative totals across the entire experiment.

```{r}
#/ warning: false
# Compute each animal's average/total per variable for entire experiment
df.avg.total <- df.daily %>%
  group_by(Animal, Photoperiod) %>%
  summarize(
    across(all_of(c(cols2sum,cols2avg)),mean)) %>%
  ungroup() %>%
  left_join(key, by = c("Animal" = "ID_Code")) %>%
  select(-c(Animal.y, Diet1,Diet2,DietChangeDateTime)) 
         
# Compute each animal's total per cumulative variable for entire experiment
df.cum.total <- df.daily %>%
  group_by(Animal, Photoperiod) %>%
  summarize(
    across(all_of(cols2sum),sum)) %>% 
  ungroup()  %>%
  left_join(key, by = c("Animal" = "ID_Code")) %>% 
  select(-c(Animal.y, Diet1,Diet2,DietChangeDateTime)) %>%
  rename(EE.cum = EE.kcal.bin)
```

## Save data

```{r}
save(df, df.hourly, df.env, df.daily, df.avg.total, df.cum.total, key, file = savename)
rm(list = ls())
```