
#Load libraries
setwd(dirname(rstudioapi::getSourceEditorContext()$path)) #path is where this code is saved - data must be in the same folder
library(readxl)
library(dplyr)
library(stringr)
library(reshape)
library(ggplot2)
library(RColorBrewer)
#library(scales)

#OPEN DATA FILES
df <- read_excel('Exp 20251201 - Quantitation of 4MEI with AMTA and Food Samples Pt4.xlsx','Plate 1 Transposed')

name_df <- read_excel('Exp 20251201 - Quantitation of 4MEI with AMTA and Food Samples Pt4.xlsx','Sample Names')

#RENAME FIRST COLUMN
names(df)[1]<-"Well_ID"

#MERGE NAMES_DF WITH MAIN DATAFRAME
df1 <- cbind(name_df, df[!names(df) %in% names(name_df)])

#WHERE THE INDIVIDUAL SPECTRA ARE IN COLUMNS, STACK THE DATA
df2 <- melt(df1, id=c("Well_ID", "SampleName", "Conc", "TechRep"))
names(df2)[names(df2) == "variable"] <- "Wavelength"

df2$value = as.numeric(df2$value)

#EXTRACT WATER CONTROLS AND COMPUTE THE AVERAGE OF TECHREPS
df_waterControls <- df2 %>%
  filter(SampleName == "Water Control") %>%
  group_by(Wavelength) %>%
  summarize_at(vars(value), list(MeanAbs = mean, SdAbs = sd))

#CALCULATE THE WATER (BLANK) CORRECTED ABSORBANCE VALUES
df3 <- inner_join(df2, df_waterControls, by=c("Wavelength")) %>%
  mutate(BC_value = value - MeanAbs) %>%
  dplyr::select(Well_ID, SampleName, Conc, TechRep, Wavelength, value, BC_value) %>%
  filter(SampleName != "Water Control")

#EXTRACT 0uM CONTROLS AND COMPUTE THE AVERAGE OF TECHREPS
df3_Blank <- df3 %>%
  filter(Conc == 0 & SampleName == "4MEI") %>%
  group_by(Wavelength) %>%
  summarize_at(vars(BC_value), list(BC_MeanAbs = mean, BC_SdAbs = sd))

#CALCULATE THE 0uM (REAGENT BLANK) CORRECTED ABSORBANCE VALUES
df4 <- inner_join(df3, df3_Blank, by=c("Wavelength")) %>%
  mutate(NC_value = BC_value - BC_MeanAbs) %>%
  dplyr::select(Well_ID, SampleName, Conc, TechRep, Wavelength, value, BC_value, NC_value)

#COMPUTE THE AVERAGE OF TECHREPS
df5 <- df4 %>%
  group_by(SampleName, Conc, Wavelength) %>%
  summarize_at(vars(NC_value), list(MeanAbs = mean, SdAbs = sd))

df5$Wavelength <- as.numeric(as.character(df5$Wavelength))

#ISOLATE THE DATA FOR THE LMAX PEAK
df5_peak <- df5 %>%
  filter(Wavelength == 506)

#SEPARATE THE FOOD SAMPLES FROM THE CALIBRATION SAMPLES
#df6_food <- df5_peak %>%
#  filter(!Conc %in% c(5,10,25,75,100))

df6_cal <- df5_peak %>%
  filter(SampleName == "4MEI")

#PERFORM LINEAR REGRESSION ON CALIBRATION STANDARDS
lm_4MEI <- lm(MeanAbs~Conc, data = df6_cal)
summary(lm_4MEI)
lm_int <- coef(lm_4MEI)[1]
lm_slope <- coef(lm_4MEI)[2]
MEI_r2 <- summary(lm_4MEI)$r.squared

#CALCUALTE CONCENTRATION OF 4MEI IN FOOD SAMPLES
df6_food_conc <- df4 %>%
  filter(Wavelength == 506 & !Conc %in% c(5,10,25,75,100)) %>%
  mutate(CalConc = (NC_value - lm_int)/lm_slope ) %>%
  group_by(SampleName, Conc, Wavelength) %>%
  summarize_at(vars(CalConc), list(MeanConc = mean, SdConc = sd))

df_rec <- df6_food_conc %>%
  inner_join(filter(df6_food_conc,Conc==0), filter(df6_food_conc, Conc==50), by="SampleName") %>%
  mutate(DiffConc=MeanConc.x-MeanConc.y) %>%
  mutate(Recovery=100*DiffConc/50) %>%
  filter(Recovery != 0 & SampleName != "4MEI") %>%
  select(SampleName, DiffConc, Recovery)

MeanRecovery <- mean(df_rec$Recovery)
SdRecovery <- sd(df_rec$Recovery)


df6_food_conc$SampleName <- factor(df6_food_conc$SampleName, levels = c("4MEI", "E150A", "E150D", "Cola 1", "WpTe", "WpTe+E150D", "WpAu", "WpAu+E150D", "PreD","PreD+E150D"))

#PLOT DATA 
pal <- brewer.pal(6, "Reds")
ggplot(df6_food_conc, aes(x=SampleName, y=MeanConc, group=as.factor(Conc))) + 
  geom_col(aes(fill=as.factor(Conc)), color="black", width = 0.7, position = position_dodge(width = 0.7)) +
  geom_errorbar(aes(ymin=(MeanConc - SdConc), ymax=(MeanConc + SdConc)), color = "black", width = 0.45, position = position_dodge(width = 0.7)) + 
  labs(y = "[4-MeI] (\U003BCM)") +
  #scale_y_continuous(breaks = seq(0, 200, by = 50), limits = c(-0.01,200)) +
  theme_classic() +
  theme(aspect.ratio=1,
        axis.title.y=element_text(size=18), 
        axis.title.x=element_blank(),
        axis.text=element_text(color = 1, size = 14),
        axis.text.x=element_text(size = 12),
        legend.title=element_blank(), 
        legend.justification=c(0,1),
        legend.position=c(0.05,0.95), 
        legend.text=element_text(size = 14),
        panel.border = element_rect(colour = "black", fill=NA)) +
  scale_fill_manual(values = c(pal[2],pal[6]), labels = c("Added 0 \U003BCM 4-MeI", "Added 50 \U003BCM 4-MeI")) +
  #scale_x_discrete(labels = c("Water\nControl", "E150A\n(0.2%)", "E150D\n(0.2%)", "Cola", "Whisky", "Whisky\n+ E150D")) +
  guides(color=guide_legend(byrow=TRUE))

#PRINT/SAVE IMAGES, CHANGE FILE NAMES AS APPROPRIATE
ggsave("Food_bar_BGCorr.pdf", height=13, width=13, units="cm", device = cairo_pdf)

#COMPUTE THE MAXIMUM ABSORBANCE WAVELENGTH
df4$Wavelength <- as.numeric(as.character(df4$Wavelength))
lmax <- df4 %>%
  filter(Conc == 0 & Wavelength > 490) %>%
  group_by(SampleName, Conc, TechRep) %>% 
  slice(which.max(NC_value)) %>%
  group_by(SampleName) %>%
  summarize_at(vars(Wavelength), list(MeanLmax = mean, SdLmax = sd))

df5_food <- df5 %>%
  filter(Conc == 0 & SampleName !="4MEI")

df5_control <- df5 %>%
  filter(Conc == 50 & SampleName == "4MEI")

df5_food <- rbind(df5_food, df5_control)

df5_food$SampleName <- factor(df5_food$SampleName, levels = c("4MEI", "E150A", "E150D", "Cola 1", "PreD+E150D", "PreD", "WpAu+E150D", "WpAu", "WpTe+E150D", "WpTe"))

#PLOT DATA
pal2 <- brewer.pal(8, "Dark2")
ggplot(df5_food, aes(x=Wavelength, y=MeanAbs, color=as.factor(SampleName))) + 
  geom_line(aes(color=as.factor(SampleName), linetype=as.factor(SampleName))) +
  labs(x = "Wavelength (nm)",
       y = "Δ Absorbance (a.u.)") +
  scale_y_continuous(breaks = seq(0, 0.3, by = 0.05), limits = c(-0.02,0.3)) +
  theme_classic() +
  theme(aspect.ratio=1,
        axis.title=element_text(size=18), 
        axis.text=element_text(color = 1, size = 14), 
        legend.title=element_blank(), 
        legend.justification=c(0,1),
        legend.position=c(0.5,0.97), 
        legend.text=element_text(size = 12),
        panel.border = element_rect(colour = "black", fill=NA)) +
  scale_color_manual(values=c(pal2[8], pal2[6], pal2[5],pal2[4],pal2[3],pal2[3],pal2[2],pal2[2],pal2[1],pal2[1]),labels = c("4-MeI Control", "E150A (0.2%)", "E150D (0.2%)", "Cola", "Bourbon + E150D", "Bourbon", "Scotch B + E150D", "Scotch B", "Scotch A + E150D", "Scotch A")) +
  scale_linetype_manual(values = c("solid", "solid", "solid", "solid", "dashed", "solid", "dashed", "solid", "dashed", "solid"), labels = c("4-MeI Control", "E150A (0.2%)", "E150D (0.2%)", "Cola", "Bourbon + E150D", "Bourbon", "Scotch B + E150D", "Scotch B", "Scotch A + E150D", "Scotch A")) +
  guides(color=guide_legend(byrow=TRUE, reverse=TRUE), linetype=guide_legend(byrow=TRUE, reverse=TRUE)) 

#PRINT/SAVE IMAGES, CHANGE FILE NAMES AS APPROPRIATE
ggsave("Food_Spectra_BGCorr2.pdf", height=13, width=13, units="cm", device = cairo_pdf)


# REMOVE THE E150D SPIKED FOODS
df7_food <- df5_food %>%
  filter(!SampleName %in% c("PreD+E150D", "WpAu+E150D", "WpTe+E150D"))
df7_food$SampleName <- factor(df7_food$SampleName, levels = c("4MEI", "E150A", "E150D", "Cola 1", "PreD",  "WpAu",  "WpTe"))

#PLOT DATA
pal2 <- brewer.pal(8, "Dark2")
ggplot(df7_food, aes(x=Wavelength, y=MeanAbs, color=as.factor(SampleName))) + 
  geom_line(aes(color=as.factor(SampleName))) +
  labs(x = "Wavelength (nm)",
       y = "Δ Absorbance (a.u.)") +
  scale_y_continuous(breaks = seq(0, 0.3, by = 0.05), limits = c(-0.02,0.3)) +
  theme_classic() +
  theme(aspect.ratio=1,
        axis.title=element_text(size=18), 
        axis.text=element_text(color = 1, size = 14), 
        legend.title=element_blank(), 
        legend.justification=c(0,1),
        legend.position=c(0.5,0.97), 
        legend.text=element_text(size = 12),
        panel.border = element_rect(colour = "black", fill=NA)) +
  scale_color_manual(values=c(pal2[8], pal2[6], pal2[5],pal2[4],pal2[3],pal2[2],pal2[1]),labels = c("4-MeI Control", "E150A (0.2%)", "E150D (0.2%)", "Cola", "Bourbon", "Scotch B", "Scotch A")) +
  guides(color=guide_legend(byrow=TRUE, reverse=TRUE)) 

#PRINT/SAVE IMAGES, CHANGE FILE NAMES AS APPROPRIATE
ggsave("Food_Spectra_BGCorr3.pdf", height=13, width=13, units="cm", device = cairo_pdf)

ggplot(df5_food, aes(x=Wavelength, y=MeanAbs, color=as.factor(SampleName))) + 
  geom_line(aes(color=as.factor(SampleName), linetype=as.factor(SampleName))) +
  #geom_smooth(aes(color=as.factor(SampleName)), stat = "smooth", method="gam") +
  labs(x = "Wavelength (nm)",
       y = "Δ Absorbance (a.u.)") +
  scale_y_continuous(breaks = seq(0, 0.25, by = 0.05), limits = c(0,0.25)) +
  theme_classic() +
  theme(aspect.ratio=1,
        axis.title=element_text(size=18), 
        axis.text=element_text(color = 1, size = 14), 
        legend.title=element_blank(), 
        legend.justification=c(0,1),
        legend.position=c(0.5,0.97), 
        legend.text=element_text(size = 14),
        panel.border = element_rect(colour = "black", fill=NA)) +
  scale_color_manual(values=c(pal2[8], pal2[6], pal2[5],pal2[4],pal2[3],pal2[3],pal2[2],pal2[2],pal2[1],pal2[1]),labels = c("4-MeI Control", "E150A (0.2%)", "E150D (0.2%)", "Cola", "Bourbon + E150D", "Bourbon", "Scotch B + E150D", "Scotch B", "Scotch A + E150D", "Scotch A")) +
  scale_linetype_manual(values = c("solid", "solid", "solid", "solid", "dashed", "solid", "dashed", "solid", "dashed", "solid"), labels = c("4-MeI Control", "E150A (0.2%)", "E150D (0.2%)", "Cola", "Bourbon + E150D", "Bourbon", "Scotch B + E150D", "Scotch B", "Scotch A + E150D", "Scotch A")) +
  guides(color=guide_legend(byrow=TRUE, reverse=TRUE), linetype=guide_legend(byrow=TRUE, reverse=TRUE)) +
  geom_vline(xintercept=510, color="magenta", linetype="dotted", size=1)+
  #geom_vline(xintercept=485, color="green3", linetype="dotted", size=1)+
  #geom_vline(xintercept=485, color="magenta", linetype="dotted", size=1)+
  geom_vline(xintercept=470, color="purple3", linetype="dotted", size=1)

