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

#import data with columns "sample", "L" (ligand concentration), and "v" (initial reaction rate/initial reaction rate in absense of ligand)
setwd("/Users/andrewmorgenthaler/Google Drive/MCDB/Copley/AMorgenthaler_lab_notebook/Project_11_Adaptation2/kinetics/2018-08-27_CarABmutants_Kd_ornithine/")
df <- as.data.frame(read_excel("2018-08-27_CarABmutants_Kd_ornithine.xlsx",1))

#calculate Kd, a only for samples that fit model
samples <- unique(df$sample)
model.samples <- c("WT", "G369V", "K966E")
for (i in 1:length(model.samples)) {
  df.model <- df[df$sample == model.samples[i],]
  eval(parse(text = paste0("allosteric.model.", model.samples[i], "<- nls(v ~ (a*L + Kd)/(L + Kd), data = df.model, 
                           start = list(Kd = df.model[which.min(abs(df.model$v-max(df.model$v)/2)), 'L'], a = max(df.model$v)))")))
}  

#make combined table of all calculated Kd, a values
allost.const.comb <- data.frame()

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

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

#export Kd, a values and SE as excel file
write.xlsx(allost.const.comb.wide, file = "2018-08-27_CarABmutants_Kd_ornithine_modelfit.xlsx", row.names = FALSE)

#calculate allosteric ligand model predictive values for line of best fit
allost.calc.comb <- data.frame()

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

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

#plot allosteric ligand curves

#plot with different point shapes
# plot <- ggplot(data = df.mean, aes(x = L, y = v.mean, color = sample, shape = sample), log10 = "x") +
#           geom_point(size = 2) +
#           geom_line(data = allost.calc.comb, aes(x = L, y = v, color = sample), size = 1) +
#           xlab("log[ornithine] (mM)") +
#           ylab(substitute(paste(italic(v), "/", italic(v[0])))) +
#           scale_y_continuous(limits = c(0,12), expand = c(0, 0), breaks = pretty_breaks()) +
#           scale_x_log10(breaks = scales::trans_breaks("log10", function(x) 10^x, n = 8), 
#                         labels = scales::trans_format("log10", scales::math_format(.x))) +
#           scale_color_brewer("", labels = c("WT", "G369V", "L960P", "L964Q", "K966E", "∆12 bp (2906 nt)", "∆132 bp (2986 nt)", "∆12 bp (3108 nt)", "21 bp dup."),
#                              type = "qual", palette = "Set1") +
#           scale_shape_manual("", labels = c("WT", "G369V", "L960P", "L964Q", "K966E", "∆12 bp (2906 nt)", "∆132 bp (2986 nt)", "∆12 bp (3108 nt)", "21 bp dup."),
#                              values=c(1,15,19,19,17,19,19,19,6)) +
#           theme(text = element_text(family="Myriad Pro"),
#                 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")

#plot with same point shape for all samples
# plot <- ggplot(data = df.mean, aes(x = L, y = v.mean, color = sample), log10 = "x") +
#           geom_point(size = 2) +
#           geom_line(data = allost.calc.comb, aes(x = L, y = v, color = sample), size = 1) +
#           xlab("log[ornithine] (mM)") +
#           ylab(substitute(paste(italic(v), "/", italic(v[0])))) +
#           scale_y_continuous(limits = c(0,12), expand = c(0, 0), breaks = pretty_breaks()) +
#           scale_x_log10(breaks = scales::trans_breaks("log10", function(x) 10^x, n = 8), 
#                         labels = scales::trans_format("log10", scales::math_format(.x))) +
#           scale_color_brewer("", labels = c("WT", "G369V", "L960P", "L964Q", "K966E", "∆12 bp (2906 nt)", "∆132 bp (2986 nt)", "∆12 bp (3108 nt)", "21 bp dup."),
#                              type = "qual", palette = "Set1") +
#           theme(text = element_text(family="Myriad Pro"),
#                 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")

#plot with same point shape for all samples + black point outlines
plot <- ggplot(log10 = "x") +
          geom_point(data = df.mean, aes(x = L, y = v.mean, fill = sample), pch = 21, size = 2) +        
          geom_line(data = allost.calc.comb, aes(x = L, y = v, color = sample), size = 1) +        
          xlab("log[ornithine] (mM)") +
          ylab(substitute(paste(italic(v), "/", italic(v[0])))) +
          scale_y_continuous(limits = c(0,12), expand = c(0, 0), breaks = pretty_breaks()) +
          scale_x_log10(breaks = scales::trans_breaks("log10", function(x) 10^x, n = 8), 
                        labels = scales::trans_format("log10", scales::math_format(.x))) +
          scale_color_manual(values = c("#377eb8", "#ff7f00", "#e41a1c")) +
          scale_fill_brewer("", labels = c("WT", "G369V", "L960P", "L964Q", "K966E", "∆12 bp (at nt 2906)", "∆132 bp (at nt 2986)", "∆12 bp (at nt 3108)", "21 bp dup. (at nt 3130)"),
                            type = "qual", palette = "Set1") +
          guides(fill = guide_legend(override.aes = list(size=4))) +
          theme(text = element_text(family="Myriad Pro"),
                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.5, 'lines'),
                legend.justification = "top")

ggsave(plot, filename = "2018-08-27_CarABmutants_Kd_ornithine.tiff", width = 8.5, height = 5)
