
#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 IMIDAZOLE FILES
df_p1 <- read_excel('Exp 20250603 - Quantitation of Imidazole with AMTA.xlsx','Plate 1 Transposed')
df_p2 <- read_excel('Exp 20250604 - Quantitation of Imidazole with AMTA Pt2.xlsx','Plate 1 Transposed')
df_p3 <- read_excel('Exp 20250616 - Quantitation of Imidazole with AMTA (Pt1 Repeat).xlsx','Plate 1 Transposed')

name_df_p1 <- read_excel('Exp 20250603 - Quantitation of Imidazole with AMTA.xlsx','Sample Names')
name_df_p2 <- read_excel('Exp 20250604 - Quantitation of Imidazole with AMTA Pt2.xlsx','Sample Names')
name_df_p3 <- read_excel('Exp 20250616 - Quantitation of Imidazole with AMTA (Pt1 Repeat).xlsx','Sample Names')

#OPEN 4-METHYLIMIDAZOLE FILES
df_p4 <- read_excel('Exp 20250820a - Quantitation of 4MEI with AMTA.xlsx','Plate 1 Transposed')
df_p5 <- read_excel('Exp 20250823 - Quantitation of 4MEI with AMTA and Selectivity.xlsx','Plate 1 Transposed')
df_p6 <- read_excel('Exp 20250824 - Quantitation of 4MEI with AMTA and Food Samples.xlsx','Plate 1 Transposed')

name_df_p4 <- read_excel('Exp 20250820a - Quantitation of 4MEI with AMTA.xlsx','Sample Names')
name_df_p5 <- read_excel('Exp 20250823 - Quantitation of 4MEI with AMTA and Selectivity.xlsx','Sample Names')
name_df_p6 <- read_excel('Exp 20250824 - Quantitation of 4MEI with AMTA and Food Samples.xlsx','Sample Names')

#RENAME FIRST COLUMN
data_list_Im <- list(df_p1, df_p2, df_p3)
data_list_4MEI <- list(df_p4, df_p5, df_p6)

rename_col1 <- function(df) {
  names(df)[1]<-"Well_ID"
  return(df)
}

data_list_Im <- lapply(data_list_Im, rename_col1)
data_list_4MEI <- lapply(data_list_4MEI, rename_col1)

#MERGE NAMES_DF WITH MAIN DATAFRAME
names_list_Im <- list(name_df_p1, name_df_p2, name_df_p3)
names_list_4MEI <- list(name_df_p4, name_df_p5, name_df_p6)

merge_names <- function(dfName, dfData, suffix, num){
  dfNew <- cbind(dfName, dfData[!names(dfData) %in% names(dfName)])
  dfNew <- cbind(Plate = num, dfNew)
  return(dfNew)
}

df1_list_Im <- Map(merge_names, names_list_Im, data_list_Im, MoreArgs = list(suffix = "Im"), num = seq_along(names_list_Im))
names(df1_list_Im) <- str_c("df1_p", seq_along(df1_list_Im))

df1_list_4MEI <- Map(merge_names, names_list_4MEI, data_list_4MEI, MoreArgs = list(suffix = "4MEI"), num = seq_along(names_list_4MEI)+3)
names(df1_list_4MEI) <- str_c("df1_p", seq_along(df1_list_4MEI)+3)

#MERGE DATAFRAMES IN THE LISTS
df1_Im <- bind_rows(df1_list_Im)
df1_4MEI <- bind_rows(df1_list_4MEI)

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

df2_4MEI <- melt(df1_4MEI, id=c("Well_ID", "SampleName", "Conc", "TechRep", "Plate"))
names(df2_4MEI)[names(df2_4MEI) == "variable"] <- "Wavelength"

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

df_waterControls_4MEI <- df2_4MEI %>%
  filter(SampleName == "Water Control") %>%
  group_by(Plate, Wavelength) %>%
  summarize_at(vars(value), list(MeanAbs = mean, SdAbs = sd))

#CALCULATE THE WATER (BLANK) CORRECTED ABSORBANCE VALUES
df3_Im <- inner_join(df2_Im, df_waterControls_Im, by=c("Plate", "Wavelength")) %>%
  mutate(BC_value = value - MeanAbs) %>%
  dplyr::select(Plate, Well_ID, SampleName, Conc, TechRep, Wavelength, value, BC_value) %>%
  filter(SampleName == "Im")

df3_4MEI <- inner_join(df2_4MEI, df_waterControls_4MEI, by=c("Plate", "Wavelength")) %>%
  mutate(BC_value = value - MeanAbs) %>%
  dplyr::select(Plate, Well_ID, SampleName, Conc, TechRep, Wavelength, value, BC_value) %>%
  filter(SampleName == "4MEI")

#COMPUTE THE AVERAGE OF TECHREPS
df4_Im <- df3_Im %>%
  group_by(SampleName, Conc, Wavelength) %>%
  summarize_at(vars(BC_value), list(MeanAbs = mean, SdAbs = sd))
df4_Im$Wavelength <- as.numeric(as.character(df4_Im$Wavelength))

df4_4MEI <- df3_4MEI %>%
  group_by(SampleName, Conc, Wavelength) %>%
  summarize_at(vars(BC_value), list(MeanAbs = mean, SdAbs = sd))
df4_4MEI$Wavelength <- as.numeric(as.character(df4_4MEI$Wavelength))

#PLOT DATA - Im
pal_Im <- brewer.pal(8, "Blues")
ggplot(df4_Im, aes(x=Wavelength, y=MeanAbs, color=as.factor(Conc))) + 
  geom_line(aes(color=as.factor(Conc))) +
  labs(x = "Wavelength (nm)",
       y = "Absorbance (a.u.)") +
  scale_y_continuous(breaks = seq(0, 1.7, by = 0.5), limits = c(-0.02,1.7)) +
  theme_classic() +
  theme(aspect.ratio=1,
        axis.title=element_text(size=18), 
        axis.text=element_text(color = 1, size = 14), 
        legend.title=element_text(color=1, size=14, face="bold", hjust = 0.5), 
        legend.justification=c(0,1),
        legend.position=c(0.7,0.97), 
        legend.text=element_text(size = 14),
        panel.border = element_rect(colour = "black", fill=NA)) +
  scale_color_manual(values = pal_Im[2:8], name = c("Imidazole"), labels = c("0 \U003BCM", "5 \U003BCM", "10 \U003BCM", "25 \U003BCM", "50 \U003BCM", "75 \U003BCM", "100 \U003BCM")) + 
  guides(color=guide_legend(byrow=TRUE, reverse = TRUE))

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

#PLOT DATA - 4MEI
pal_4MEI <- brewer.pal(8, "Reds")
ggplot(df4_4MEI, aes(x=Wavelength, y=MeanAbs, color=as.factor(Conc))) + 
  geom_line(aes(color=as.factor(Conc))) +
  labs(x = "Wavelength (nm)",
       y = "Absorbance (a.u.)") +
  scale_y_continuous(breaks = seq(0, 1.7, by = 0.5), limits = c(-0.02,1.7)) +
  theme_classic() +
  theme(aspect.ratio=1,
        axis.title=element_text(size=18), 
        axis.text=element_text(color = 1, size = 14), 
        legend.title=element_text(color=1, size=14, face="bold", hjust = 0.5), 
        legend.justification=c(0,1),
        legend.position=c(0.7,0.97), 
        legend.text=element_text(size = 14),
        panel.border = element_rect(colour = "black", fill=NA)) +
  scale_color_manual(values = pal_4MEI[2:8], name = c("4-MeI"), labels = c("0 \U003BCM", "5 \U003BCM", "10 \U003BCM", "25 \U003BCM", "50 \U003BCM", "75 \U003BCM", "100 \U003BCM")) + 
  guides(color=guide_legend(byrow=TRUE, reverse = TRUE))

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

#EXTRACT 0uM CONTROLS AND COMPUTE THE AVERAGE OF TECHREPS
df_Blank_Im <- df3_Im %>%
  filter(Conc == 0) %>%
  group_by(Plate, Wavelength) %>%
  summarize_at(vars(BC_value), list(BC_MeanAbs = mean, BC_SdAbs = sd))

df_Blank_4MEI <- df3_4MEI %>%
  filter(Conc == 0) %>%
  group_by(Plate, Wavelength) %>%
  summarize_at(vars(BC_value), list(BC_MeanAbs = mean, BC_SdAbs = sd))

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

df5_4MEI <- inner_join(df3_4MEI, df_Blank_4MEI, by=c("Plate", "Wavelength")) %>%
  mutate(NC_value = BC_value - BC_MeanAbs) %>%
  dplyr::select(Plate, Well_ID, SampleName, Conc, TechRep, Wavelength, value, BC_value, NC_value)

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

df6_4MEI <- df5_4MEI %>%
  group_by(Plate, SampleName, Conc, Wavelength) %>%
  summarize_at(vars(NC_value), list(MeanAbs = mean, SdAbs = sd))

#COMPUTE THE AVERAGE OF EXPERIMENTAL REPEATS
df7_Im <- df6_Im %>%
  group_by(SampleName, Conc, Wavelength) %>%
  summarize_at(vars(MeanAbs), list(MeanAbs2 = mean, SdAbs2 = sd))
df7_Im$Wavelength <- as.numeric(as.character(df7_Im$Wavelength))

df7_4MEI <- df6_4MEI %>%
  group_by(SampleName, Conc, Wavelength) %>%
  summarize_at(vars(MeanAbs), list(MeanAbs2 = mean, SdAbs2 = sd))
df7_4MEI$Wavelength <- as.numeric(as.character(df7_4MEI$Wavelength))

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

df5_4MEI$Wavelength <- as.numeric(as.character(df5_4MEI$Wavelength))
lmax_4MEI <- df5_4MEI %>%
  filter(Conc > 0) %>%
  group_by(SampleName, Conc, TechRep) %>% 
  slice(which.max(NC_value)) %>%
  group_by(SampleName, Conc) %>%
  summarize_at(vars(Wavelength), list(MeanLmax = mean, SdLmax = sd))

#PLOT DATA - Im
ggplot(df7_Im, aes(x=Wavelength, y=MeanAbs2, color=as.factor(Conc))) + 
  geom_line(aes(color=as.factor(Conc))) +
  labs(x = "Wavelength (nm)",
       y = "Δ Absorbance (a.u.)") +
  scale_y_continuous(breaks = seq(0, 0.7, by = 0.2), limits = c(-0.02,0.7)) +
  theme_classic() +
  theme(aspect.ratio=1,
        axis.title=element_text(size=18), 
        axis.text=element_text(color = 1, size = 14), 
        legend.title=element_text(color=1, size=14, face="bold", hjust = 0.5), 
        legend.justification=c(0,1),
        legend.position=c(0.7,0.97), 
        legend.text=element_text(size = 14),
        panel.border = element_rect(colour = "black", fill=NA)) +
  scale_color_manual(values = pal_Im[2:8], name = c("Imidazole"), labels = c("0 \U003BCM", "5 \U003BCM", "10 \U003BCM", "25 \U003BCM", "50 \U003BCM", "75 \U003BCM", "100 \U003BCM")) + 
  #annotate("text", x = 470, y = 0.75, label = expression(lambda["max"]~"494 nm")) +
  guides(color=guide_legend(byrow=TRUE, reverse = TRUE))

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

#PLOT DATA - 4MEI
ggplot(df7_4MEI, aes(x=Wavelength, y=MeanAbs2, color=as.factor(Conc))) + 
  geom_line(aes(color=as.factor(Conc))) +
  labs(x = "Wavelength (nm)",
       y = "Δ Absorbance (a.u.)") +
  scale_y_continuous(breaks = seq(0, 0.6, by = 0.2), limits = c(-0.02,0.6)) +
  theme_classic() +
  theme(aspect.ratio=1,
        axis.title=element_text(size=18), 
        axis.text=element_text(color = 1, size = 14), 
        legend.title=element_text(color=1, size=14, face="bold", hjust = 0.5), 
        legend.justification=c(0,1),
        legend.position=c(0.7,0.97), 
        legend.text=element_text(size = 14),
        panel.border = element_rect(colour = "black", fill=NA)) +
  scale_color_manual(values = pal_4MEI[2:8], name = c("4-MeI"), labels = c("0 \U003BCM", "5 \U003BCM", "10 \U003BCM", "25 \U003BCM", "50 \U003BCM", "75 \U003BCM", "100 \U003BCM")) + 
  guides(color=guide_legend(byrow=TRUE, reverse = TRUE))

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

#SUBSET TO EXAMPLE CURVES OF Im AND 4MEI
df7_100uM_Im <- df7_Im %>%
  filter(Conc == 100)

df7_100uM_4MEI <- df7_4MEI %>%
  filter(Conc == 100)

df7_100uM <- rbind(df7_100uM_Im, df7_100uM_4MEI)

#PLOT DATA - 100 uM Im AND 4MEI
ImCol <- brewer.pal(7, "Blues")
MEICol <- brewer.pal(7, "Reds")
ggplot(df7_100uM, aes(x=Wavelength, y=MeanAbs2, color=as.factor(SampleName))) +
  geom_vline(xintercept = 494, color=ImCol[3], linetype=2) +
  geom_vline(xintercept = 506, color=MEICol[3], linetype=2) +
  geom_line(aes(color=as.factor(SampleName))) +
  labs(x = "Wavelength (nm)",
       y = "Δ Absorbance (a.u.)") +
  scale_y_continuous(breaks = seq(0, 0.7, by = 0.2), limits = c(-0.02,0.7)) +
  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.6,0.95), 
        legend.text=element_text(size = 14),
        panel.border = element_rect(colour = "black", fill=NA)) +
  scale_color_manual(values = c(MEICol[5], ImCol[5]), labels = c("4-MeI", "Imidazole")) +
  guides(color=guide_legend(byrow=TRUE, reverse =TRUE))
  
#PRINT/SAVE IMAGES, CHANGE FILE NAMES AS APPROPRIATE
ggsave("4MEI_Im_100uM_Spectra_BGCorr.pdf", height=13, width=13, units="cm", device = cairo_pdf)

#ISOLATE THE DATA FOR THE LMAX PEAK
df7_Im_peak <- df7_Im %>%
  filter(Wavelength == 494) %>%
  mutate(RSD2 = 100*SdAbs2/MeanAbs2)

df7_4MEI_peak <- df7_4MEI %>%
  filter(Wavelength == 506) %>%
  mutate(RSD2 = 100*SdAbs2/MeanAbs2)

#PERFORM LINEAR REGRESSION ON CALIBRATION STANDARDS
lm_Im <- lm(MeanAbs2~Conc, data = df7_Im_peak)
summary(lm_Im)
Im_r2 <- summary(lm_Im)$r.squared

df6_Im_peak_blank <- df6_Im %>%
  filter(Conc==0 & Wavelength == 494) %>%
  mutate(lod = 3.3*SdAbs/lm_Im$coefficients[2]) %>%
  mutate(loq = 10*SdAbs/lm_Im$coefficients[2])

mean(df6_Im_peak_blank$lod)
sd(df6_Im_peak_blank$lod)

lod_Im <- 0.0111284*3.3/lm_Im$coefficients[2]
loq_Im <- 0.0111284*10/lm_Im$coefficients[2]

lm_4MEI <- lm(MeanAbs2~Conc, data = df7_4MEI_peak)
summary(lm_4MEI)
MEI_r2 <- summary(lm_4MEI)$r.squared

df6_4MEI_peak_blank <- df6_4MEI %>%
  filter(Conc==0 & Wavelength == 506) %>%
  mutate(lod = 3.3*SdAbs/lm_4MEI$coefficients[2]) %>%
  mutate(loq = 10*SdAbs/lm_4MEI$coefficients[2])

mean(df6_4MEI_peak_blank$lod)
sd(df6_4MEI_peak_blank$lod)

lod_4MEI <- 0.0094417*3.3/lm_4MEI$coefficients[2]
loq_4MEI <- 0.0094417*10/lm_4MEI$coefficients[2]

#PLOT DATA - Im
ggplot(df7_Im_peak, aes(x=Conc, y=MeanAbs2)) +
  geom_abline(intercept = coef(lm_Im)[1], slope = coef(lm_Im)[2], color = ImCol[5], linetype=2, size=1) +
  geom_point(color= ImCol[6], size=3) +
  geom_errorbar(aes(ymin=(MeanAbs2 - SdAbs2), ymax=(MeanAbs2 + SdAbs2)), color = ImCol[6], width = 2) +
  labs(x = "[Imidazole] (\U003BCM)",
       y = "Δ Absorbance 494 nm (a.u.)") +
  scale_y_continuous(breaks = seq(0, 0.8, by = 0.2), limits = c(-0.02,0.8)) +
  theme_classic() +
  theme(aspect.ratio=1,
        axis.title=element_text(size=18), 
        axis.text=element_text(color = 1, size = 14), 
        legend.title=element_text(color=1, size=14, face="bold", hjust = 0.5), 
        legend.justification=c(0,1),
        legend.position=c(0.7,0.97), 
        legend.text=element_text(size = 14),
        panel.border = element_rect(color = "black", fill=NA)) +
  annotate("text", color=ImCol[6], label = paste0("R^2 == ", round(Im_r2, 4)), parse = TRUE, x=15,y=0.6, size=5)


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

#PLOT DATA - 4MEI
ggplot(df7_4MEI_peak, aes(x=Conc, y=MeanAbs2)) +
  geom_abline(intercept = coef(lm_4MEI)[1], slope = coef(lm_4MEI)[2], color = MEICol[5], linetype=2, size=1) +
  geom_point(color= MEICol[6], size=3) +
  geom_errorbar(aes(ymin=(MeanAbs2 - SdAbs2), ymax=(MeanAbs2 + SdAbs2)), color = MEICol[6], width = 2) +
  labs(x = "[4-MeI] (\U003BCM)",
       y = "Δ Absorbance 506 nm (a.u.)") +
  scale_y_continuous(breaks = seq(0, 0.7, by = 0.2), limits = c(-0.02,0.7)) +
  theme_classic() +
  theme(aspect.ratio=1,
        axis.title=element_text(size=18), 
        axis.text=element_text(color = 1, size = 14), 
        legend.title=element_text(color=1, size=14, face="bold", hjust = 0.5), 
        legend.justification=c(0,1),
        legend.position=c(0.7,0.97), 
        legend.text=element_text(size = 14),
        panel.border = element_rect(color = "black", fill=NA)) +
  annotate("text", color=MEICol[6], label = paste0("R^2 == ", round(MEI_r2, 4)), parse = TRUE, x=15,y=0.5, size=5)


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

#PLOT DATA - BOTH Im AND 4MEI
df7_peak <- rbind(df7_Im_peak, df7_4MEI_peak)

ggplot(df7_peak, aes(x=Conc, y=MeanAbs2, group=as.factor(SampleName))) +
  geom_abline(intercept = coef(lm_Im)[1], slope = coef(lm_Im)[2], color = ImCol[5], linetype=2, size=1) +
  geom_abline(intercept = coef(lm_4MEI)[1], slope = coef(lm_4MEI)[2], color = MEICol[5], linetype=2, size=1) +
  geom_point(aes(color=SampleName, shape=SampleName), size=3) +
  geom_errorbar(aes(ymin=(MeanAbs2 - SdAbs2), ymax=(MeanAbs2 + SdAbs2), color=SampleName), width = 2) +
  labs(x = "Concentration (\U003BCM)",
       y = "Δ Absorbance (a.u.)") +
  scale_y_continuous(breaks = seq(0, 0.8, by = 0.2), limits = c(-0.02,0.8)) +
  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.05,0.95), 
        legend.text=element_text(size = 14),
        panel.border = element_rect(colour = "black", fill=NA)) +
  scale_color_manual(values = c(MEICol[5], ImCol[5]), labels = c("4-MeI (506 nm)", "Imidazole (494 nm)")) +
  scale_shape_manual(values = c(17,16), labels = c("4-MeI (506 nm)", "Imidazole (494 nm)")) +
  guides(color=guide_legend(byrow=TRUE, reverse =TRUE), 
         shape=guide_legend(byrow=TRUE, reverse =TRUE)) +
  annotate("text", color=MEICol[5], label = paste0("R^2 == ", round(MEI_r2, 4)), parse = TRUE, x=75,y=0.71, size=5) +
  annotate("text", color=ImCol[5], label = paste0("R^2 == ", round(Im_r2, 4)), parse = TRUE, x=75,y=0.76, size=5)

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