
#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_p1 <- read_excel('Exp 20250823 - Quantitation of 4MEI with AMTA and Selectivity.xlsx','Plate 1 Transposed')
df_p2 <- read_excel('Exp 20250924 - Quantitation of 4MEI with AMTA and Food Samples Pt2.xlsx','Plate 1 Transposed')

name_df_p1 <- read_excel('Exp 20250823 - Quantitation of 4MEI with AMTA and Selectivity.xlsx','Sample Names')
name_df_p2 <- read_excel('Exp 20250924 - Quantitation of 4MEI with AMTA and Food Samples Pt2.xlsx','Sample Names')

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

#MERGE NAMES_DF WITH MAIN DATAFRAME
df1_p1 <- cbind(name_df_p1, df_p1[!names(df_p1) %in% names(name_df_p1)])
df1_p2 <- cbind(name_df_p2, df_p2[!names(df_p2) %in% names(name_df_p2)])

#ADD A COLUMN FOR DATE OF EXPERIMENT
df1_p1_sub <- cbind(Plate=rep("1", nrow(df1_p1)), df1_p1)

df1_p2_sub <- cbind(Plate=rep("2", nrow(df1_p2)), df1_p2)

#COMBINE DATASETS
df2 <- rbind(df1_p1_sub, df1_p2_sub)

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

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

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

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

#ISOLATE THE SELECTIVITY SAMPLES
df4_p1 <- df4 %>%
  filter(Plate == 1 & Conc %in% c(0,50))

df4_p2 <- df4 %>%
  filter(Plate == 2 & SampleName %in% c("4MEIControl", "Histamine", "NegControl")) %>%
  mutate(SampleName = ifelse(as.character(SampleName) == "4MEIControl", "4MEI", 
                             ifelse(as.character(SampleName) == "NegControl", "4MEI", as.character(SampleName)
                                    )))

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

df4_Blank_p2 <- df4_p2 %>%
  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_p1 <- inner_join(df4_p1, df4_Blank_p1, 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_p2 <- inner_join(df4_p2, df4_Blank_p2, 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 <- rbind(df5_p1, df5_p2)

#COMPUTE THE AVERAGE OF TECHREPS
df6 <- df5 %>%
  group_by(SampleName, Conc, Wavelength) %>%
  summarize_at(vars(NC_value), list(MeanAbs = mean, SdAbs = sd)) %>%
  filter(Conc >0)
df6$Wavelength <- as.numeric(as.character(df6$Wavelength))

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

df6$SampleName <- factor(df6$SampleName, levels = c("Im", "4MEI", "2-MEI", "His", "Histamine", "Fur", "5-HMF", "2-FP"))

#PLOT DATA 
ggplot(df6, 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.36, by = 0.05), limits = c(-0.02,0.36)) +
  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.97), 
        legend.text=element_text(size = 14),
        panel.border = element_rect(colour = "black", fill=NA)) +
  scale_color_brewer(palette = "Dark2",labels = c("Imidazole", "4-MeI", "2-MeI", "L-Histidine", "Histamine", "Furfural", "5-HMF", "2-FP")) +
  guides(color=guide_legend(byrow=TRUE))

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

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

df6_peak$SampleName <- factor(df6_peak$SampleName, levels = c("Im", "4MEI", "2-MEI", "His", "Histamine", "Fur", "5-HMF", "2-FP"))

#PLOT DATA 
pal <- brewer.pal(6, "Reds")
ggplot(df6_peak, aes(x=SampleName, y=MeanAbs)) + 
  geom_col(fill=pal[6], width = 0.75, color="black") +
  geom_errorbar(aes(ymin=(MeanAbs - SdAbs), ymax=(MeanAbs + SdAbs)), color = "black", width = 0.5) + 
  labs(y = "Δ Absorbance 506 nm (a.u.)") +
  scale_y_continuous(breaks = seq(0, 0.36, by = 0.05), limits = c(-0.01,0.36)) +
  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(angle=45, vjust=1, hjust=1, size = 16),
        legend.title=element_blank(), 
        legend.justification=c(0,1),
        legend.position=c(0.6,0.97), 
        legend.text=element_text(size = 14),
        panel.border = element_rect(colour = "black", fill=NA)) +
  scale_x_discrete(labels = c("Imidazole", "4-MeI", "2-MeI", "L-Histidine", "Histamine", "Furfural", "5-HMF", "2-FP")) +
  guides(color=guide_legend(byrow=TRUE))

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