library(readxl)
library(ggplot2)
library(scales)
library(reshape2)
library(data.table)
library(xlsx)

#import data with columns "sample", "S" (substrate concentration), and "v" (initial reaction rate)
setwd("/Users/andrewmorgenthaler/Google Drive/MCDB/Copley/AMorgenthaler_lab_notebook/Project_11_Adaptation2/kinetics/2018-11-05_ProA_GSA/")
df <- as.data.frame(read_excel("2018-11-05_ProA_GSA.xlsx",1))

#calculate Km, Vmax
samples <- unique(df$sample)

for (i in 1:length(samples)) {
  df.model <- df[df$sample == samples[i],]
  eval(parse(text = paste0("mm.model.", samples[i], "<- nls(v ~ Vm * S/(Km+S), data = df.model, 
                  start = list(Km = df.model[which.min(abs(df.model$v-max(df.model$v)/2)), 'S'], Vm = max(df.model$v)))")))
}  

#make combined table of all calculated Km, Vmax values
MMconst.comb <- data.frame()

for (i in 1:length(samples)) {
  const <- as.data.frame(summary(eval(parse(text = paste0("mm.model.", samples[i]))))$coefficients)
  setDT(const, keep.rownames = TRUE)
  const$sample <- rep(samples[i], times = nrow(const))
  MMconst.comb <- rbind(MMconst.comb, const)
}

MMconst.comb.wide <- reshape(MMconst.comb, idvar = "sample", timevar = "rn", direction = "wide")

#export Km, Vmax values and SE as excel file
write.xlsx(MMconst.comb.wide, file = "2018-11-05_ProA_GSA_modelfit.xlsx", row.names = FALSE)

#calculate MM model predictive values for line of best fit
mm.calc.comb <- data.frame()

for (i in 1:length(samples)) {
  df.model <- df[df$sample == samples[i],]
  mm.calc <- data.frame(S = seq(0, max(df.model$S), length.out = 100))
  mm.calc$v <- predict(eval(parse(text = paste0("mm.model.", samples[i]))), newdata = mm.calc)
  mm.calc$sample <- rep(samples[i], times = 100)
  mm.calc.comb <- rbind(mm.calc.comb, mm.calc)
}

#summary df of observed data
df.mean <- aggregate(df$v, by = list(df$sample, df$S), FUN = mean)
colnames(df.mean) <- c('sample', 'S', 'v.mean')
df.mean$sample <- factor(df.mean$sample, levels = samples)

#plot MM curves
plot <- ggplot(data = df.mean[df.mean$sample=='WT',], aes(x = S, y = v.mean)) +
  geom_point(size = 2, color = '#f35e5a') +
  geom_line(data = mm.calc.comb[mm.calc.comb$sample=='WT',], aes(x = S, y = v), size = 1, color = '#f35e5a') +
  xlab("[GSA] (mM)") +
  ylab(substitute(paste(italic(v[0]), " (µM/s)"))) +
  scale_y_continuous(limits = c(0,1), expand = c(0, 0), breaks = pretty_breaks()) +
  scale_x_continuous(limits = c(0,0.4), expand = c(0, 0), breaks = pretty_breaks(n = 6)) +
  scale_color_discrete(labels = c("WT")) +
  ggtitle("WT") +
  theme(text = element_text(family="Myriad Pro"),
        plot.title = element_text(size=28, face="bold", hjust = 0.5),
        axis.title.y = element_text(size = 20),
        axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        panel.background = element_rect(fill = NA),
        panel.grid.major = element_line(colour = "grey95"),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, colour = "black"),
        legend.text = element_text(size = 16),
        legend.title = element_blank(),
        legend.key = element_blank(),
        legend.key.size = unit(1.75, 'lines'),
        legend.justification = "top")

ggsave(plot, filename = "2018-11-05_ProA_WT_GSA.tiff", width = 8, height = 5)

plot <- ggplot(data = df.mean[df.mean$sample=='E383A',], aes(x = S, y = v.mean)) +
  geom_point(size = 2, color = '#17b12b') +
  geom_line(data = mm.calc.comb[mm.calc.comb$sample=='E383A',], aes(x = S, y = v), size = 1, color = '#17b12b') +
  xlab("[GSA] (mM)") +
  ylab(substitute(paste(italic(v[0]), " (µM/s)"))) +
  scale_y_continuous(limits = c(0,0.03), expand = c(0, 0), breaks = pretty_breaks()) +
  scale_x_continuous(limits = c(0,0.4), expand = c(0, 0), breaks = pretty_breaks(n = 6)) +
  scale_color_discrete(labels = c("E383A")) +
  ggtitle("E383A") +
  theme(text = element_text(family="Myriad Pro"),
        plot.title = element_text(size=28, face="bold", hjust = 0.5),
        axis.title.y = element_text(size = 20),
        axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        panel.background = element_rect(fill = NA),
        panel.grid.major = element_line(colour = "grey95"),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, colour = "black"),
        legend.text = element_text(size = 16),
        legend.title = element_blank(),
        legend.key = element_blank(),
        legend.key.size = unit(1.75, 'lines'),
        legend.justification = "top")

ggsave(plot, filename = "2018-11-05_ProA_E383A_GSA.tiff", width = 8, height = 5)

plot <- ggplot(data = df.mean[df.mean$sample=='E383A_F372L',], aes(x = S, y = v.mean)) +
  geom_point(size = 2, color = '#5086ff') +
  geom_line(data = mm.calc.comb[mm.calc.comb$sample=='E383A_F372L',], aes(x = S, y = v), size = 1, color = '#5086ff') +
  xlab("[GSA] (mM)") +
  ylab(substitute(paste(italic(v[0]), " (µM/s)"))) +
  scale_y_continuous(limits = c(0,0.07), expand = c(0, 0), breaks = pretty_breaks(n = 7)) +
  scale_x_continuous(limits = c(0,0.4), expand = c(0, 0), breaks = pretty_breaks(n = 6)) +
  scale_color_discrete(labels = c("E383A F372L")) +
  ggtitle("E383A F372L") +
  theme(text = element_text(family="Myriad Pro"),
        plot.title = element_text(size=28, face="bold", hjust = 0.5),
        axis.title.y = element_text(size = 20),
        axis.title.x = element_text(size = 18),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        panel.background = element_rect(fill = NA),
        panel.grid.major = element_line(colour = "grey95"),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, colour = "black"),
        legend.text = element_text(size = 16),
        legend.title = element_blank(),
        legend.key = element_blank(),
        legend.key.size = unit(1.75, 'lines'),
        legend.justification = "top")

ggsave(plot, filename = "2018-11-05_ProA_E383A_F372L_GSA.tiff", width = 8, height = 5)