---
title: "1_extract_tensometer_data"
author: "Matt George"
date: "2023-04-24"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Load packages; set ggplot theme
```{r message=TRUE, warning=TRUE}

## RUN as part of "PSMFC-mytilus-byssus-pilot.Rproj"


## clear workspace
rm(list=ls())

## Load Packages
load.lib<-c("readxl","tidyverse","openxlsx","EnvStats") # List of required packages
install.lib <- load.lib[!load.lib %in% installed.packages()] # Select missing packages
for(lib in install.lib) install.packages(lib,dependencies=TRUE) # Install missing packages + dependencies
sapply(load.lib,require,character=TRUE) # Load all packages.

## Set ggplot theme
my_theme <- theme(line              = element_line(size=1.2),
                  rect              = element_rect(size=1.2),
                  text              = element_text(size=14,color="black"),
                  panel.background  = element_blank(),
                  panel.grid.major  = element_blank(), 
                  panel.grid.minor  = element_blank(),
                  axis.text.x       = element_text(size=16,color="black"), #,angle=90),
                  axis.text.y       = element_text(size=16,color="black"),
                  axis.title.x      = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
                  axis.title.y      = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
                  axis.ticks.x      = element_line(color="black"),
                  axis.ticks.y      = element_line(color="black"),
                  # axis.line         = element_line(color = "black", size = 0.1),
                  panel.border      = element_rect(color = "black", fill=NA, size=1.2),
                  legend.key        = element_blank(), # removes background of legend bullets
                  legend.position   = "none"
                  )

```


## CONTROLS

```{r warning=FALSE}

# extract data from all .txt files in a folder located at thread_strength/tensometer_output/control/ and add it to a nested list object
# the data has three column headers: Time, displacement, force. 
# Add a third data column called "mussel" that is the first three characters of the file name.
getwd()
# Set the working directory to the folder containing the txt files
setwd("thread_strength/tensometer_output/control/")

# Create an empty list to hold the data
nested_list <- list()

# Get a vector of the names of all the txt files in the directory
txt_files <- list.files(pattern = "\\.txt$")

# Loop through each txt file, extract the data, and add it to the nested list
for (file_name in txt_files) {
    # Read in the data from the txt file
    data <- read.delim(file_name, 
                       header = FALSE,
                       sep = "\t", # how is the data separated?
                       dec = ".")
    
    # transpose
    data_transpose = t(data) 
    
    # Snip data, Fix headers
    data_snipped <- as.data.frame(data_transpose[,2:4]) # keep all rows, grab only column 2 through 4
    colnames(data_snipped) <- c('time', 'displacement', 'force')
    
    # Fix column data types
    data_snipped$time         <- as.numeric(data_snipped$time)
    data_snipped$displacement <- as.numeric(data_snipped$displacement)
    data_snipped$force        <- as.numeric(data_snipped$force)
    
    # Fix NA
    data_snipped$time[is.na(data_snipped$time)] <- 0
    data_snipped$displacement[is.na(data_snipped$displacement)] <- 0
    data_snipped$force[is.na(data_snipped$force)] <- 0
  
    # Extract mussel ID & thread test
    split_1 <- str_split(file_name, "_") # Split file name by _
    split_2 <- split_1[[1]][[2]] # Save thread test & .txt
    split_3 <- str_split(split_2, "\\.") # Split thread test and .txt
    data_snipped$mussel <- split_1[[1]][[1]] # Saving mussel ID
    data_snipped$thread <- split_3[[1]][[1]] # Saving thread number
  
    # Add the data to the nested list using the file name as the list name
    nested_list[[file_name]] <- data_snipped
}

```


```{r}
# Loop through the list object "nested_list"
# calculate the maximum value from the "force" column
# Save the maximum into a dataframe, along with the first values in the "mussel" column and "thread" column

# create a data frame with three columns: mussel, thread, and max_force
max_force_control <- data.frame(mussel = numeric(),
                        thread = numeric(),
                        max_force = numeric())

# loop through the nested list
for(i in seq_along(nested_list)){
  # create a data frame for each list in the nested list
  current_df <- nested_list[[i]]
  # add a row to the max_force data frame with the mussel, thread, and max_force
  max_force_control <- rbind(max_force_control, data.frame(mussel = current_df$mussel[1],
                                            thread = current_df$thread[1],
                                            max_force = max(current_df$force)))

}

```


```{r}

require(stats)

# Loop through the list object "nested_list"
# extract the "time" and "force" column
# remove rows that contain zeros before first non-zero row. Apply to time column as well.
# remove rows that contain zeros after last non-zero row. Apply to time column as well.
# fit loess function to "time" and "force" after trimmed. Time is x, force is y.
# calculate area under the curve and save in dataframe, along with the first values in the "mussel" column and "thread" column from "nested_list"
# use a for loop to iterate through the nested list object "nested_list"

# Set the working directory to the folder containing the txt files
setwd("thread_strength/QC_plots/control/")

# create a data frame with three columns: mussel, thread, and max_force
integral_control <- data.frame(mussel = numeric(),
                        thread = numeric(),
                        auc = numeric())

for(i in seq_along(nested_list)){
  # extract time and force column from the current nested list element using $ notation
  current_df <- nested_list[[i]][c("time", "force")]
  
  # fit loess function to "time" and "force" using loess() and save as an object
  current_loess <- loess(force ~ time, data = current_df)
  
  # plot the loess function with the corresponding time and force data points
  jpeg(paste(nested_list[[i]]$mussel[1],"_",nested_list[[i]]$thread[1],".jpg",sep=""))
  plot(current_df$time, current_df$force, main= "Force over Time Plot with Loess Function", 
       xlab = "Time", ylab = "Force", col = "blue")
  lines(current_loess, col = "red") # add the loess function to the plot in red
  dev.off()
  
  # calculate area under the curve using approx() to approximate the integral and save to a variable
  auc <- sum(diff(current_loess$x) * (approx(current_loess$x, current_loess$y, n = length(current_loess$x)))$y[-1])
  
  # save the first values in the "mussel" and "thread" columns to variables
  mussel_first_val <- nested_list[[i]]$mussel[1]
  thread_first_val <- nested_list[[i]]$thread[1]
  
  # create a dataframe to store the results, and save the variables to the corresponding columns in the dataframe
  integral_control <- rbind(integral_control, data.frame(mussel = nested_list[[i]]$mussel[1],
                                           thread = nested_list[[i]]$thread[1],
                                           auc = auc))
}


```

```{r}
library(ggplot2)
hist(max_force_control$max_force)
max_force_control$ID <- paste(max_force_control$mussel, max_force_control$thread)
  ggplot(data=max_force_control, aes(x=1, y=max_force, label = ID))+
  geom_boxplot(width=.5)+
  # jittered text with geom_text
  geom_text(check_overlap = TRUE,
            position=position_jitter(width=0.15))+
  theme(legend.position="none")
```
```{r}
library(EnvStats)
test <- rosnerTest(max_force_control$max_force,
  k = 5
)
test
```

```{r}
write.xlsx(max_force_control, file="thread_strength/summarized_data/max_force_control.xlsx", colNames = TRUE)
write.xlsx(integral_control, file="thread_strength/summarized_data/integral_control.xlsx", colNames = TRUE)
print("done")
```




## TREATMENT

```{r warning=FALSE}

# extract data from all .txt files in a folder located at thread_strength/tensometer_output/control/ and add it to a nested list object
# the data has three column headers: Time, displacement, force. 
# Add a third data column called "mussel" that is the first three characters of the file name.

# Set the working directory to the folder containing the txt files
setwd("thread_strength/tensometer_output/treatment/")

# Create an empty list to hold the data
nested_list <- list()

# Get a vector of the names of all the txt files in the directory
txt_files <- list.files()

# Loop through each txt file, extract the data, and add it to the nested list
for (file_name in txt_files) {
    # Read in the data from the txt file
    data <- read.delim(file_name, 
                       header = FALSE,
                       sep = "\t", # how is the data separated?
                       dec = ".")
    
    # transpose
    data_transpose = t(data) 
    
    # Snip data, Fix headers
    data_snipped <- as.data.frame(data_transpose[,2:4]) # keep all rows, grab only column 2 through 4
    colnames(data_snipped) <- c('time', 'displacement', 'force')
    
    # Fix column data types
    data_snipped$time         <- as.numeric(data_snipped$time)
    data_snipped$displacement <- as.numeric(data_snipped$displacement)
    data_snipped$force        <- as.numeric(data_snipped$force)
    
    # Fix NA
    data_snipped$time[is.na(data_snipped$time)] <- 0
    data_snipped$displacement[is.na(data_snipped$displacement)] <- 0
    data_snipped$force[is.na(data_snipped$force)] <- 0
  
    # Add a new column called "mussel" that contains the first FOUR characters of the file name
    data_snipped$mussel <- substr(file_name, start = 1, stop = 4)
    
    # Add new column that adds thread run
    data_snipped$thread <- substr(file_name, start = 6, stop = 7)
  
    # Add the data to the nested list using the file name as the list name
    nested_list[[file_name]] <- data_snipped
}

```


```{r}
# Loop through the list object "nested_list"
# calculate the maximum value from the "force" column
# Save the maximum into a dataframe, along with the first values in the "mussel" column and "thread" column

# create a data frame with three columns: mussel, thread, and max_force
max_force_treatment <- data.frame(mussel = numeric(),
                        thread = numeric(),
                        max_force = numeric())

# loop through the nested list
for(i in seq_along(nested_list)){
  # create a data frame for each list in the nested list
  current_df <- nested_list[[i]]
  # add a row to the max_force data frame with the mussel, thread, and max_force
  max_force_treatment <- rbind(max_force_treatment, data.frame(mussel = current_df$mussel[1],
                                            thread = current_df$thread[1],
                                            max_force = max(current_df$force)))

}

```


```{r}

require(stats)

# Loop through the list object "nested_list"
# extract the "time" and "force" column
# remove rows that contain zeros before first non-zero row. Apply to time column as well.
# remove rows that contain zeros after last non-zero row. Apply to time column as well.
# fit loess function to "time" and "force" after trimmed. Time is x, force is y.
# calculate area under the curve and save in dataframe, along with the first values in the "mussel" column and "thread" column from "nested_list"
# use a for loop to iterate through the nested list object "nested_list"

# Set the working directory to the folder containing the txt files
setwd("thread_strength/QC_plots/treatment/")

# create a data frame with three columns: mussel, thread, and max_force
integral <- data.frame(mussel = numeric(),
                        thread = numeric(),
                        auc = numeric())

# create a data frame with three columns: mussel, thread, and max_force
integral_treatment <- data.frame(mussel = numeric(),
                        thread = numeric(),
                        auc = numeric())

for(i in seq_along(nested_list)){
  # extract time and force column from the current nested list element using $ notation
  current_df <- nested_list[[i]][c("time", "force")]
  
  # fit loess function to "time" and "force" using loess() and save as an object
  current_loess <- loess(force ~ time, data = current_df)
  
  # plot the loess function with the corresponding time and force data points
  jpeg(paste(nested_list[[i]]$mussel[1],"_",nested_list[[i]]$thread[1],".jpg",sep=""))
  plot(current_df$time, current_df$force, main= "Force over Time Plot with Loess Function", 
       xlab = "Time", ylab = "Force", col = "blue")
  lines(current_loess, col = "red") # add the loess function to the plot in red
  dev.off()
  
  # calculate area under the curve using approx() to approximate the integral and save to a variable
  auc <- sum(diff(current_loess$x) * (approx(current_loess$x, current_loess$y, n = length(current_loess$x)))$y[-1])
  
  # save the first values in the "mussel" and "thread" columns to variables
  mussel_first_val <- nested_list[[i]]$mussel[1]
  thread_first_val <- nested_list[[i]]$thread[1]
  
  # create a dataframe to store the results, and save the variables to the corresponding columns in the dataframe
  integral_treatment <- rbind(integral_treatment, data.frame(mussel = nested_list[[i]]$mussel[1],
                                           thread = nested_list[[i]]$thread[1],
                                           auc = auc))
}


```

```{r}
library(ggplot2)
hist(max_force_treatment$max_force)
max_force_treatment$ID <- paste(max_force_treatment$mussel, max_force_treatment$thread)
  ggplot(data=max_force_treatment, aes(x=1, y=max_force, label = ID))+
  geom_boxplot(width=.5)+
  # jittered text with geom_text
  geom_text(check_overlap = TRUE,
            position=position_jitter(width=0.15))+
  theme(legend.position="none")
```

```{r}
library(EnvStats)
test <- rosnerTest(max_force_treatment$max_force,
  k = 5
)
test
```

```{r}
write.xlsx(max_force_treatment, file="thread_strength/summarized_data/max_force_treatment.xlsx", colNames = TRUE)
write.xlsx(integral_treatment, file="thread_strength/summarized_data/integral_treatment.xlsx", colNames = TRUE)
print("done")
```
