# Compute SLAKES score and SI-600 for results from 36 well format
# Claire Phillips
# Last updated July 2023

# About this script:
# This script requires the filepath to a folder named after the samples analyzed in one slaking experiment
# Inside that folder it expects:
#     -A sub-folder named "Images" that contains all the images imported into Image_J (and no extras)
#     -A spreadsheet named [SampleName]_SampleInfo.csv that identifies the sample in each region of interest (ROI).
#       It should contain columns named 'Well' (= ROI number), 'Sample' (= sample ID), and	'Rep' (= replicate number)
#     -A spreadsheet named [SampleName]_Results.csv that contains the ImageJ ROI Multi-Measure results.
#

# The function 'GetTimes' gets the time stamp for each image in 'Images', based on its file name, and writes  
# it to a spreadsheet called 'ImageInfo.csv'.

##### GetTimes ####
GetTimes <- function(FP){
  FP.Images <- paste(FP,"/Images/",sep="")
  FileNames <- list.files(FP.Images)
  if(grep("desktop.ini",FileNames)==1) FileNames <- FileNames[-1]
  LongFileNames<-paste(FP.Images,FileNames,sep="")
  
  CreatedTimes <- file.mtime(LongFileNames)
  class(CreatedTimes)
  
  Sec <- round(CreatedTimes-CreatedTimes[2],0)
  
  Summary <- as.data.frame(cbind(FileNames,Sec))
  Summary$Sec <- as.numeric(Summary$Sec)
  write.csv(Summary,paste(FP,"/ImageInfo.csv",sep=""),row.names=FALSE)
  return(Summary)
}

# The function 'CountImages' checks that there are the same number of images in 'Images' as rows in the ImageJ output.
# In some cases bad images were removed during the ImageJ process but not deleted from the 'Images' folder.
#### Count Images ####
CountImages <- function(FP, SampleName){
  # Get time stamp for each image in the Images subfolder
  Times<-GetTimes(FP)
  
  Results.FP<-paste(FP,"/",SampleName,"_Results.csv",sep="")
  Results.df<-read.csv(Results.FP)
  print(paste("There are", nrow(Times) ,"images, and" , nrow(Results.df), "were analyzed."))
}


# The function 'ComputeSI600' computes the slaking index following Flynn et al., 2020
# It creates a spreadsheet in FP called [SampleName]_SLAKES values.csv
# If Plot is TRUE it also plots the timeseries for each aggregate.

# N.B. I have not included much error handling.

#### ComputeSI600 ####
ComputeSI600 <- function(SampleName,StartWell, NumWell, FP){
  
  # Get the Image-J Multi Measure output
  Results.FP<-paste(FP,"/",SampleName,"_Results.csv",sep="")
  Results.df<-read.csv(Results.FP)
  #dim(Results.df);names(Results.df)
  
  #Rename first column as 'ROI'
  Results.df <- Results.df %>% rename(ROI = X)
  
  # Get the timestamp for each image
  Times<-GetTimes(FP)
  
  # Create a time series of areas for each ROI
  Areas <- dplyr::select(Results.df, contains('Area'))
  SLAKES.df <- cbind(Sec=Times$Sec,Areas)
  Reference <- slice_head(Areas,n=1)
  Slaking <- slice(Areas,-1)
  ChangeAreas <- sweep(Slaking,2,c(unlist(Reference)),'-')
  SI600Values <- sweep(ChangeAreas,2,c(unlist(Reference)),'/')
  SLAKES.df <- cbind(Sec=Times$Sec[-1],ChangeAreas)
  
  #Tidy data to plot with ggplot
  SLAKES.tidy <- pivot_longer(SLAKES.df,cols=-1,names_to="Cell",values_to="Area_Change")
  
  #Add column for well ID
  # Use this if wells are in 1A...6D format
  ROI <- as.numeric(str_split_fixed(SLAKES.tidy$Cell, "Area",2)[,2])
  if(max(ROI) < NumWell) print("Did the regions of interest start at well 1?")
  SLAKES.tidy$Well <- ROI + StartWell - 1
  
  #Get Sample Names and merge with well ID
  FP.SampleNames<-paste(FP,"/",SampleName,"_Samples.csv",sep="")
  SampleNames.df<-read.csv(FP.SampleNames)
  names(SampleNames.df)
  SLAKES.tidy<-merge(SLAKES.tidy,SampleNames.df[,c("Well","Sample","Rep")],by="Well")
  
  p <- ggplot(SLAKES.tidy,aes(x=Sec,y=Area_Change,color=as.factor(Rep)))+
    geom_point() +
    facet_wrap(~Sample) +
    xlab("Seconds") +
    ylab(expression(Slaking~index~group("(",mm^2~mm^'-2',")"))) +
    scale_color_discrete(name = "Replicate") +
    theme_classic(base_size = 16)
  print(p)
  #ggsave(paste(FigFP, "ExampleSlakingTimeseries.jpg",sep=""),width=5.4, height = 3.3, units="in")
  
  # Get the row closest to 600
  TimeDiff<-abs(SLAKES.df$Sec-600)
  ClosestTime<-SLAKES.df$Sec[which(TimeDiff==min(TimeDiff))][1]
  Report<-SLAKES.tidy %>% filter(Sec==ClosestTime) %>% 
    dplyr::select(Well, SI600 = Area_Change, Sample, Rep)
  Report$InitialArea <- unlist(Reference)
  
  # Output SI600 values and original area of each aggregate
  write.csv(Report, paste(FP,"/",SampleName,"_SLAKES values.csv",sep=""), row.names=FALSE)
  
  #N.B. Important to go into Report and flag wells that had problems, like hitting sidewall
}


# Load libraries ####
library(dplyr)
library(tidyr)
library(lubridate)
library(ggplot2)
library(stringr)

#Empty the object list environment, if desired
#rm(list=ls()) 

#### Enter user supplied info: plate name, filepath, the first ROI that actually contains an aggregate, and the numbe of ROIs. ####
SampleName<-"Kelly P16_P17_P20_Time4"
FP<-paste("G:/My Drive/Group Projects/Kelly Extracellular Polymeric Substances/SLAKES/",SampleName,sep="")
StartWell <- 1
NumWell <- 36

# Execute Functions with the user-supplied info, above
CountImages(FP, SampleName)
ComputeSI600(SampleName,StartWell,NumWell,FP)