library(data.table)
library(ggplot2)
library(reshape2)
library(stringr)
library(gridExtra)
library(tcpl)
source("L:/Lab/NCCT_Lab/Simmons_Lab/R_Git_Repos/tcpl/LoadTcpl.R")
source("L:/Lab/NCCT_Lab/Simmons_Lab/R_Git_Repos/tcpl_lite/tcplFit.R")
source("L:/Lab/NCCT_Lab/Simmons_Lab/R_Git_Repos/tcpl_lite/tcpl_Fit_Lite_Sample_Data.R")
source("L:/Lab/NCCT_Lab/Simmons_Lab/R_Git_Repos/tcpl_lite/tcplObjCnst.R")
source("L:/Lab/NCCT_Lab/Simmons_Lab/R_Git_Repos/tcpl_lite/tcplObjHill.R")
source("L:/Lab/NCCT_Lab/Simmons_Lab/R_Git_Repos/tcpl_lite/tcplObjGnls.R")

match.function <- function(x, a, b){
  return(x[match(a[1:length(a)],b)])
}

hill_curve <- function(hill_tp, hill_ga, hill_gw, lconc){
  return(hill_tp/(1+10^((hill_ga - lconc)*hill_gw)))
}
##########
#Sample data is a data frame with 5 chemicals. Change to data table, then loop over chemical name through tcplFit
#tcplFit needs something called bmad, but it doesn't impact the results, 
#but can limit the number of curves attempted to be fit, so set smaller than you think it really is
bmad <- .001
##########




####### Tier 1 Single-Concentration Testing ########


### read-in data table
dt0 <- read.csv(file = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2/sc_seahorse.lvl0.merged.data.csv", stringsAsFactors = FALSE)
dt0 <- as.data.table(dt0)
dt0$time <- factor(dt0$time, levels = c("T1", "T2", "T3", "T4", "T5", "T6", "T7", "T8", "T9", "T10", "T11", "T12"))
  
### correct chem stock concnetrations and logc from ChemTracker info
samples <- tcplQuery("SELECT * FROM sample")
dt0$stck.correct <- ifelse(is.na(samples$stkc[match(dt0$spid[1:length(dt0$spid)], samples$spid)]), dt0$stock.conc.mM, samples$stkc[match(dt0$spid[1:length(dt0$spid)], samples$spid)])
dt0$conc.correction <- ifelse(is.na(dt0$stck.correct), dt0$Assay.Conc.uM, dt0$Assay.Conc.uM * dt0$stck.correct / dt0$stock.conc.mM)
dt0$logc <- log10(dt0$conc.correction)
dt0$trim.conc <- signif(dt0$conc.correction, digits =3)

#remove blank wells, convert rval to %norm, give DMSO a logc name, and define a dt w/o DMSO
dt <- dt0[spid != "blank", ]
dt$cval <- dt$rval * 100
dt$trim.conc[dt$spid == "DMSO"] <- "DMSO"

dt2 <- dt[spid != "DMSO", ]

### match chmn to spid
dt.chem <- tcplLoadChem(field = "spid", val = unique(dt0$spid), include.spid = TRUE)
dt2$chmn <- ifelse(is.na(dt.chem$chnm[match(dt2$spid[1:length(dt2$spid)], dt.chem$spid)]), as.character(dt2$spid), dt.chem$chnm[match(dt2$spid[1:length(dt2$spid)], dt.chem$spid)])

#NCCT_MITO_basal_respiration_state_OCR_up and down
dt.basal <- NULL
dt.basal <- dt0[time == "T6",]
dt.basal[, bval := median(rval[wllt == "n"]), by = apid]
dt.basal$fc.up <- dt.basal$rval / dt.basal$bval - 1
dt.basal$fc.down <- -1 * dt.basal$fc.up

dt.basal.med <- dt.basal[, .(med.up = median(fc.up), med.down = median(fc.down), mad.down = mad(fc.down)), by = spid]

T6_coeff <- 0.2 #20% activity threshold
basal.up.actives <- dt.basal.med[med.up > T6_coeff& spid != "DNP", ]
basal.down.actives <- dt.basal.med[med.down > T6_coeff & spid != "Fenpyroximate", ]


#NCCT_MITO_maximal_respiration_state_OCR_up and down
dt.max <- NULL
dt.max <- dt0[time == "T9",]
dt.max[, bval := median(rval[wllt == "n"]), by = apid]
dt.max$fc.up <- dt.max$rval / dt.max$bval - 1
dt.max$fc.down <- -1 * dt.max$fc.up

dt.max.med <- dt.max[, .(med.up = median(fc.up), med.down = median(fc.down), mad.down = mad(fc.down)), by = spid]

T9_coeff <- 0.2 #20% activity threshold
max.up.actives <- dt.max.med[med.up > T9_coeff & spid != "DNP", ]
max.down.actives <- dt.max.med[med.down > T9_coeff & spid != "Fenpyroximate", ]


#NCCT_MITO_inhibited_respiration_state_OCR_up and down
dt.inhib <- NULL
dt.inhib <- dt0[time == "T12",]
dt.inhib[, bval := median(rval[wllt == "n"]), by = apid]
dt.inhib$fc.up <- dt.inhib$rval / dt.inhib$bval - 1
dt.inhib$fc.down <- -1 * dt.inhib$fc.up

dt.inhib.med <- dt.inhib[, .(med.up = median(fc.up), med.down = median(fc.down)), by = spid]

T12_coeff <- 0.3 #30% activity threshold
inhib.up.actives <- dt.inhib.med[med.up > T12_coeff & spid != "DNP", ]
inhib.down.actives <- dt.inhib.med[med.down > T12_coeff & spid != "Fenpyroximate", ]

#combine unique, active spids into single list and save

hit_spids <- unique(c(as.character(basal.up.actives$spid), as.character(basal.down.actives$spid), as.character(max.up.actives$spid), 
                      as.character(max.down.actives$spid), as.character(inhib.up.actives$spid), as.character(inhib.down.actives$spid)))
write.csv(hit_spids, file = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2/SC-hit-spids2.csv", row.names = FALSE)

### supplemental table of tested spids- by chmn with test conc and activity
test_spids <- unique(dt0$spid[!dt0$spid %in% c("Fenpyroximate", "DNP", "blank", "DMSO")])
inactive_spids <- setdiff(test_spids, hit_spids)

SC.table <- data.table(spid = test_spids, chnm = match.function(dt2$chmn, test_spids, dt2$spid), 
                       conc = match.function(dt2$trim.conc, test_spids, dt2$spid), 
                       activity = ifelse(test_spids %in% hit_spids, "active", "inactive"))

SC.summ.table <- SC.table[, .(test_conc = mean(as.numeric(conc)), activity = unique(activity)), by  = chnm]


write.csv(SC.summ.table, file = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2/SC.summ.table.csv", row.names = FALSE)

### subset active spids and controls for plotting (sort by chmn)
SC.active.spids <- hit_spids
dt.hits <- data.table(spid = SC.active.spids$x)
dt.hits$chmn <- dt2$chmn[match(dt.hits$spid[1:length(dt.hits$spid)], dt2$spid)]
setkey(dt.hits, chmn)
dt.controls <- data.table(spid = c("Fenpyroximate", "DNP"), chmn = c("Fenpyroximate", "DNP"))
dt.plot <- rbind(dt.controls, dt.hits)

### backcalculate coff to rval and then %OCR for seahorse plots
basal.upper <- median(dt.basal$rval[dt.basal$spid == "DMSO"]) *1.2 *100
basal.lower <- median(dt.basal$rval[dt.basal$spid == "DMSO"]) *0.8 *100
max.upper <- median(dt.max$rval[dt.max$spid == "DMSO"]) *1.2 *100
max.lower <- median(dt.max$rval[dt.max$spid == "DMSO"]) *0.8 *100
inhib.upper <- median(dt.inhib$rval[dt.inhib$spid == "DMSO"]) *1.3 *100
inhib.lower <- median(dt.inhib$rval[dt.inhib$spid == "DMSO"]) *0.7 *100

### Seahorse-style plots by test chemical/control plotted with DMSO

plot_spids <- unique(dt0$spid[!dt0$spid %in% c("blank", "DMSO")])

for (i in 1:length(plot_spids)) {
  
  dt3 <- dt[spid == "DMSO" | spid == plot_spids[i], ]
  concs <- sort(as.numeric(unique(dt3[spid != "DMSO", trim.conc])))
  concs.order <- c("DMSO", concs)
  dt3$trim.conc <- factor(dt3$trim.conc, levels = concs.order)
  dt4 <- dt3[ , .(median(rval)*100, mad(rval)*100), by = .(trim.conc, time)]
  
  title <- ifelse(plot_spids[i] %in% dt2$spid, unique(dt2[spid == plot_spids[i], chmn]), unique(dt3$spid)[2])
  title.font <- ifelse(nchar(title) < 50, 18, 14)
  
  summplot <- ggplot(dt4, aes(x = time, y = V1, group = trim.conc, color = trim.conc)) + geom_line(aes(group = trim.conc)) +
    geom_errorbar(aes(ymin = dt4$V1 - dt4$V2, ymax = dt4$V1 + dt4$V2), width = 0.05) +
    geom_point(data = dt4, aes(x = time, y = V1, color = trim.conc), size = 2.5) +
    scale_color_manual(values = c("#000000", "#0000ff", "#ff0000", "#ff8800", "#aaff00", "#00dd00", "#00aaff", "#1500ff"), 
                       labels = concs.order, name = "[Compound] uM") +
    scale_x_discrete(breaks = unique(dt4$time), labels= as.factor(seq(from = 6, to = 72, by = 6))) +
    coord_cartesian(ylim=c(0,275)) +
    geom_segment(aes(x = 3.2, y = 0, xend = 3.2, yend = 266), color = "#555555", linetype = 5) +
    geom_text(aes(x = 3.2, y = 275, label = "Test Compound"), size = 4, color = "#444444", show.legend = FALSE) +
    geom_segment(aes(x = 6.2, y = 0, xend = 6.2, yend = 266), color = "#555555", linetype = 5) +
    geom_text(aes(x = 6.2, y = 275, label = "FCCP"), size = 4, color = "#444444", show.legend = FALSE) +
    geom_segment(aes(x = 9.2, y = 0, xend = 9.2, yend = 266), color = "#555555", linetype = 5) +
    geom_text(aes(x = 9.2, y = 275, label = "ROT/AA"), size = 4, color = "#444444", show.legend = FALSE) +
    geom_text(aes(x = 1.7, y = 245, label = "Pre-Injection"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_text(aes(x = 4.7, y = 245, label = "Basal"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_text(aes(x = 7.7, y = 245, label = "Maximal"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_text(aes(x = 10.7, y = 245, label = "Inhibited"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_segment(aes(x = 0.8, y = 245, xend = 0.4, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 2.6, y = 245, xend = 3.1, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 4.2, y = 245, xend = 3.3, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 5.2, y = 245, xend = 6.1, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 7.0, y = 245, xend = 6.3, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 8.4, y = 245, xend = 9.1, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 10.0, y = 245, xend = 9.3, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 11.4, y = 245, xend = 12.2, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 5.7, y = basal.upper, xend = 6.3, yend = basal.upper), color = "#ff0000", linetype = 1) +
    geom_segment(aes(x = 5.7, y = basal.lower, xend = 6.3, yend = basal.lower), color = "#ff0000", linetype = 1) +
    geom_segment(aes(x = 8.7, y = max.upper, xend = 9.3, yend = max.upper), color = "#ff0000", linetype = 1) +
    geom_segment(aes(x = 8.7, y = max.lower, xend = 9.3, yend = max.lower), color = "#ff0000", linetype = 1) +
    geom_segment(aes(x = 11.7, y = inhib.upper, xend = 12.3, yend = inhib.upper), color = "#ff0000", linetype = 1) +
    geom_segment(aes(x = 11.7, y = inhib.lower, xend = 12.3, yend = inhib.lower), color = "#ff0000", linetype = 1) +
    theme_bw() +
    ggtitle(title) +
    theme(panel.grid.major = element_line(colour = "grey"), panel.grid.major.x = element_blank()) +
    theme(panel.grid.minor = element_line(colour = "lightgrey")) +
    ylab("% Initial Respiration") +
    xlab("Time (min)") +
    theme (axis.title=element_text(size= 16,face="bold"),
           plot.title=element_text(size= title.font,face="bold", hjust = 0.5),
           legend.text=element_text(size=12), 
           legend.title=element_text(color="#000000"), 
           axis.text = element_text(size = 12, color = "#000000")) 

  file.name <- paste(plot_spids[i], ".png", sep = "")
  
  ggsave(summplot, filename = file.name, path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2/Seahorse_plots", 
         width = 7.5, height = 6.5, units = "in", device = "png", dpi = 300)
}

### duplicate chems plot for Figure 2
unique.spids <- unique(dt2$spid[grep(dt2$spid, pattern = "EPAPLT")])
spid.chmn <- ifelse(is.na(dt.chem$chnm[match(unique.spids[1:length(unique.spids)], dt.chem$spid)]), as.character(unique.spids), dt.chem$chnm[match(unique.spids[1:length(unique.spids)], dt.chem$spid)])
tested.chmn <- data.table(spid = unique.spids, chmn = spid.chmn)

chmn.occurance <- tested.chmn[, .(RowCount = .N), by = chmn]
duplicate.chmn <- chmn.occurance[RowCount > 1, chmn]
duplicate.chmn <- sort(duplicate.chmn)

tested.duplicates <- tested.chmn[chmn %in% duplicate.chmn, ]
tested.duplicates$basal.med.fc.down <- dt.basal.med$med.down[match(tested.duplicates$spid[1:length(tested.duplicates$spid)], dt.basal.med$spid)]
tested.duplicates$basal.mad.fc.down <- dt.basal.med$mad.down[match(tested.duplicates$spid[1:length(tested.duplicates$spid)], dt.basal.med$spid)]
tested.duplicates$max.med.fc.down <- dt.max.med$med.down[match(tested.duplicates$spid[1:length(tested.duplicates$spid)], dt.max.med$spid)]
tested.duplicates$max.mad.fc.down <- dt.max.med$mad.down[match(tested.duplicates$spid[1:length(tested.duplicates$spid)], dt.max.med$spid)]

tested.duplicates$chmn[tested.duplicates$chmn == "Perfluorooctanesulfonic acid"] <- "PFOS"

#Supp Figure 1A
SFigure1A.plot <- NULL
SFigure1A.plot <- ggplot(tested.duplicates, aes(x = chmn, y = basal.med.fc.down, group = spid)) + 
                    geom_point(aes(x = chmn, y = basal.med.fc.down, group = spid), size = 4, shape = 16, position = position_dodge(width = 0.5)) + 
                    geom_errorbar(aes(ymin = tested.duplicates$basal.med.fc.down - tested.duplicates$basal.mad.fc.down, 
                                      ymax = tested.duplicates$basal.med.fc.down + tested.duplicates$basal.mad.fc.down), width = 0.4,
                                      position = position_dodge(width = 0.5)) +
                    geom_hline(yintercept = 0.2, color = "#ff0000", linetype = 5) +
                    theme_bw() +  
                    ggtitle("Basal Respiration") +
                    ylab("Fold Decrease") +
                    xlab("") +
                    theme(axis.title=element_text(size=16,face="bold", color = "#000000"), plot.title=element_text(size=title.font,face="bold", hjust = 0.5)) +
                    theme(text = element_text(size=16), axis.text.x = element_text(hjust= 1, vjust = 1, angle = 45, color = "#000000"), 
                          axis.text.y = element_text(color = "#000000"))

ggsave(SFigure1A.plot, filename = "basal.duplicates.png", path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2", 
       width = 7.5, height = 6.5, units = "in", device = "png", dpi = 300)

#Supp Figure 1B
SFigure1B.plot <- NULL
SFigure1B.plot <- ggplot(tested.duplicates, aes(x = chmn, y = max.med.fc.down, group = spid)) + 
                    geom_point(aes(x = chmn, y = max.med.fc.down, group = spid), size = 4, shape = 16, position = position_dodge(width = 0.5)) + 
                    geom_errorbar(aes(ymin = tested.duplicates$max.med.fc.down - tested.duplicates$max.mad.fc.down, 
                                      ymax = tested.duplicates$max.med.fc.down + tested.duplicates$max.mad.fc.down), width = 0.4,
                                      position = position_dodge(width = 0.5)) +
                    geom_hline(yintercept = 0.2, color = "#ff0000", linetype = 5) +
                    theme_bw() + 
                    ggtitle("Maximal Respiration") +
                    ylab("Fold Decrease") +
                    xlab("") +
                    theme (axis.title=element_text(size=16,face="bold", color = "#000000"), plot.title=element_text(size=title.font,face="bold", hjust = 0.5)) +
                    theme(text = element_text(size=16), axis.text.x = element_text(hjust= 1, vjust = 1, angle = 45, color = "#000000"), 
                          axis.text.y = element_text(color = "#000000"))

ggsave(SFigure1B.plot, filename = "max.duplicates.png", path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2", 
       width = 7.5, height = 6.5, units = "in", device = "png", dpi = 300)

### Activity range plots for Figure 3
basal.activity <- dt.basal.med[!dt.basal.med$spid %in% c("blank", "DMSO", "DNP", "Fenpyroximate"), ]
basal.activity$hitc <- ifelse(basal.activity$med.up > 0.2 | basal.activity$med.down > 0.2, "active", "inactive")
basal.activity$order <- order(basal.activity$med.up)
max.activity <- dt.max.med[!dt.max.med$spid %in% c("blank", "DMSO", "DNP", "Fenpyroximate"), ]
max.activity$hitc <- ifelse(max.activity$med.up > 0.2 | max.activity$med.down > 0.2, "active", "inactive")
max.activity$order <- order(max.activity$med.up)
inhib.activity <- dt.inhib.med[!dt.inhib.med$spid %in% c("blank", "DMSO", "DNP", "Fenpyroximate"), ]
inhib.activity$hitc <- ifelse(inhib.activity$med.up > 0.3 | inhib.activity$med.down > 0.3, "active", "inactive")
inhib.activity$order <- order(inhib.activity$med.up)

#Figure 3A
Figure3A.plot <- NULL
Figure3A.plot <- ggplot(basal.activity, aes(y = order(basal.activity$order), x = med.up, color = hitc)) + 
  geom_point() + 
  scale_y_continuous(breaks = c(0, 200, 400, 600, 800, 1000), limits = c(1,1050)) + 
  scale_x_continuous(breaks = c(seq(-0.8,0.8,0.2)), limits = c(-0.8,0.8)) + 
  scale_color_manual(values = c("#000000", "#999999"), guide = FALSE) +
  geom_vline(xintercept = -0.2, linetype = 5, color = "red", size = 0.4) + 
  geom_vline(xintercept = 0.2, linetype = 5, color = "red", size = 0.4) + 
  geom_vline(xintercept = 0, linetype = 1, color = "blue", size = 0.4) +
  ggtitle("Basal Respiration") +
  xlab("Fold Change") + ylab("Chemical Sample Rank") + 
  theme_bw() +
  theme(axis.title=element_text(size=16,face="bold", color = "#000000"), 
        plot.title=element_text(size=title.font,face="bold", hjust = 0.5),
        text = element_text(size=16), 
        axis.text.x = element_text(color = "#000000"), 
        axis.text.y = element_text(color = "#000000"),
        axis.line = element_line(color = "#000000", size = 0.5)) 

ggsave(Figure3A.plot, filename = "basal.ranged.png", path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2", 
       width = 4.5, height = 7.5, units = "in", device = "png", dpi = 300)

#Figure 3B
Figure3B.plot <- NULL
Figure3B.plot <- ggplot(max.activity, aes(y = order(max.activity$order), x = med.up, color = hitc)) + 
  geom_point() + 
  scale_y_continuous(breaks = c(0, 200, 400, 600, 800, 1000), limits = c(1,1050)) + 
  scale_x_continuous(breaks = c(seq(-0.9,0.3,0.2)), limits = c(-0.9,0.3)) + 
  scale_color_manual(values = c("#000000", "#999999"), guide = FALSE) +
  geom_vline(xintercept = -0.2, linetype = 5, color = "red", size = 0.4) + 
  geom_vline(xintercept = 0.2, linetype = 5, color = "red", size = 0.4) + 
  geom_vline(xintercept = 0, linetype = 1, color = "blue", size = 0.4) +
  ggtitle("Maximal Respiration") +
  xlab("Fold Change") + ylab("") + 
  theme_bw() +
  theme(axis.title=element_text(size=16,face="bold", color = "#000000"), 
        plot.title=element_text(size=title.font,face="bold", hjust = 0.5),
        text = element_text(size=16), 
        axis.text.x = element_text(color = "#000000"), 
        axis.text.y = element_text(color = "#000000"),
        axis.line = element_line(color = "#000000", size = 0.5))

ggsave(Figure3B.plot, filename = "max.ranged.png", path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2", 
       width = 4.5, height = 7.5, units = "in", device = "png", dpi = 300)

#Figure 3C
Figure3C.plot <- NULL
Figure3C.plot <- ggplot(inhib.activity, aes(y = order(inhib.activity$order), x = med.up, color = hitc)) + 
  geom_point() + 
  scale_y_continuous(breaks = c(0, 200, 400, 600, 800, 1000), limits = c(1,1050)) + 
  scale_x_continuous(breaks = c(seq(-0.7, 0.9, 0.2)), limits = c(-0.7,0.9)) + 
  scale_color_manual(values = c("#000000", "#999999"), guide = FALSE) +
  geom_vline(xintercept = -0.3, linetype = 5, color = "red", size = 0.4) + 
  geom_vline(xintercept = 0.3, linetype = 5, color = "red", size = 0.4) + 
  geom_vline(xintercept = 0, linetype = 1, color = "blue", size = 0.4) +
  ggtitle("Inhibited Respiration") +
  xlab("Fold Change") + ylab("") + 
  theme_bw() +
  theme(axis.title=element_text(size=16,face="bold", color = "#000000"), 
        plot.title=element_text(size=title.font,face="bold", hjust = 0.5),
        text = element_text(size=16), 
        axis.text.x = element_text(color = "#000000"), 
        axis.text.y = element_text(color = "#000000"),
        axis.line = element_line(color = "#000000", size = 0.5))

ggsave(Figure3C.plot, filename = "inhib.ranged.png", path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2", 
       width = 4.5, height = 7.5, units = "in", device = "png", dpi = 300)

### Tier 1 (single-conc) QC ###

#Global Plot for Fenpyroximate (all apid together)
dt.control <- dt[wllt %in% c("n", "o", "m"), ] #subset control data
dt.control[, bval := median(rval[time == "T6" & spid == "DMSO"]), by = apid]
dt.control$resp <- (1 - dt.control$rval / dt.control$bval) * 100

dt.plot <- dt.control[spid != "DMSO" & time == "T6", ]
dt.plot$logc <- as.numeric(dt.plot$logc)

for (i in 1:length(unique(dt.plot$spid))){
  
  chemtitle <- unique(dt.plot$spid)[i]
  
  sub.dt <- dt.plot[spid == unique(dt.plot$spid)[i],]
  
  dat_steve <- sub.dt
  setkey(dat_steve, spid)
  dat_steve <- dat_steve[,spid := as.factor(spid)]
  
  #Get hill parameters for each spid
  dat_hill <- data.table(spid = unique(dat_steve[,spid]))
  
  for(chem in dat_hill[,spid]){
    pipefit <- tcplFit_Lite(dat_steve[spid == chem,logc], dat_steve[spid == chem,resp], bmad)
    dat_hill[spid == chem, hill_tp := pipefit$hill_tp]
    dat_hill[spid == chem, hill_ga := pipefit$hill_ga]
    dat_hill[spid == chem, hill_gw := pipefit$hill_gw]
  }
  
  spid.ac50 <- paste("AC50 = ", signif(10^dat_hill$hill_ga, digits = 3), "uM", sep = "")
  
  ymax <- max(0.5, max(sub.dt$resp)) + 0.1
  ymin <- min(-0.1, min(sub.dt$resp))
  xval <- min(dat_steve$logc) + 0.1*(max(dat_steve$logc) - min(dat_steve$logc))
  yval <- ymax - 0.1*(ymax-ymin)
  
  # The palette with black:
  cbbPalette <- c("#000000", "#ff0000")
  
  curves_plot <- ggplot(dat_steve, aes(x=logc, y=resp, color = spid))
  
  #plot the fitted curves, color by spid. This code may need to be modified if you have more colors than in cbbPalette
  i <- 1L
  for(chem in dat_hill[,spid]){
    chemical_fit_hill_tp <- dat_hill[spid == chem, hill_tp]
    chemical_fit_hill_ga <- dat_hill[spid == chem, hill_ga]
    chemical_fit_hill_gw <- dat_hill[spid == chem, hill_gw]
    curves_plot <- curves_plot +
      stat_function(fun = hill_curve, args=list(hill_tp = chemical_fit_hill_tp, 
                                                hill_ga = chemical_fit_hill_ga, 
                                                hill_gw = chemical_fit_hill_gw),
                    color = cbbPalette[i], 
                    size=1) 
    i <- i + 1L
  }
  
  curves_plot <- curves_plot +
    theme_bw() +
    geom_point(size = 4) + 
    scale_colour_manual(values=cbbPalette, guide = FALSE) +
    ylab("% Decreased Respiration") +
    xlab("log10 uM") +
    geom_text(aes(x = xval, y = yval, label = spid.ac50), size = 4, color = "#000000", show.legend = FALSE) +
    ggtitle(chemtitle) + theme (axis.title=element_text(size=16,face="bold"),
                                plot.title=element_text(size=18,face="bold", hjust = 0.5)) 
  
  ggsave(curves_plot, filename = "FENP.sc.png", path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2", 
         width = 7.5, height = 6.5, units = "in", device = "png")
}

# Hill fit by apid and calculate QC stats

apid.list <- NULL
spid.list <- NULL
AC50.um.list <- NULL
control.median.list <- NULL
control.mad.list <- NULL
vehicle.median.list <- NULL
vehicle.mad.list <- NULL
rZ.list <- NULL
rCV.list <- NULL

for (i in 1:length(unique(dt.plot$apid))){
  
  sub.dt <- NULL
  sub.dt <- dt.plot[apid == unique(dt.plot$apid)[i] & time == "T6",]
  
  dat_steve <- sub.dt
  setkey(dat_steve, spid)
  dat_steve <- dat_steve[,spid := as.factor(spid)]
  
  #Get hill parameters for each spid
  dat_hill <- data.table(spid = unique(dat_steve[,spid]))
  
  for(chem in dat_hill[,spid]){
    pipefit <- tcplFit_Lite(dat_steve[spid == chem,logc], dat_steve[spid == chem,resp], bmad)
    dat_hill[spid == chem, hill_tp := pipefit$hill_tp]
    dat_hill[spid == chem, hill_ga := pipefit$hill_ga]
    dat_hill[spid == chem, hill_gw := pipefit$hill_gw]
  }
  
  dat_hill$AC50.uM <- signif(10^dat_hill$hill_ga, digits = 3)
  
  apid.list <- c(apid.list, as.character(unique(dt.plot$apid)[i]))
  spid.list <- c(spid.list, as.character(unique(dt.plot$spid)))
  AC50.um.list <- c(AC50.um.list, dat_hill$AC50.uM)
  
  control.median <- sub.dt[logc == max(sub.dt$logc), median(rval)]
  control.mad <- sub.dt[logc == max(sub.dt$logc), mad(rval)]
  vehicle.median <- dt.control[apid == unique(sub.dt$apid) & time == "T6" & spid == "DMSO", median(rval)]
  vehicle.mad <- dt.control[apid == unique(sub.dt$apid) & time == "T6" & spid == "DMSO", mad(rval)]
  
  rZ <- 1 - 3 * (control.mad + vehicle.mad) / abs(control.median - vehicle.median)
  rZ.list <- c(rZ.list, rZ)
  
  rCV <- 100 * vehicle.mad / vehicle.median
  rCV.list <- c(rCV.list, rCV)
  
}

fenp.stat.table <- data.frame("apid" = apid.list, "spid" = spid.list, "ac50.um" = AC50.um.list, "rZ" = rZ.list, "rCV" = rCV.list)
write.csv(fenp.stat.table, file = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2/SC.FENP.table.csv", row.names = FALSE)

#Global Plot for DNP (all apid together)
dt.control <- dt[wllt %in% c("n", "p", "c"), ] #subset control data
dt.control[, bval := median(rval[time == "T6" & spid == "DMSO"]), by = apid]
dt.control$resp <- (dt.control$rval / dt.control$bval - 1) * 100

dt.plot <- dt.control[spid != "DMSO" & time == "T6", ]
dt.plot$logc <- as.numeric(dt.plot$logc)

for (i in 1:length(unique(dt.plot$spid))){
  
  chemtitle <- "2,4-Dinitrophenol"
  
  sub.dt <- dt.plot[spid == unique(dt.plot$spid)[i],]
  
  dat_steve <- sub.dt
  setkey(dat_steve, spid)
  dat_steve <- dat_steve[,spid := as.factor(spid)]
  
  #Get hill parameters for each spid
  dat_hill <- data.table(spid = unique(dat_steve[,spid]))
  
  for(chem in dat_hill[,spid]){
    pipefit <- tcplFit_Lite(dat_steve[spid == chem,logc], dat_steve[spid == chem,resp], bmad)
    dat_hill[spid == chem, hill_tp := pipefit$hill_tp]
    dat_hill[spid == chem, hill_ga := pipefit$hill_ga]
    dat_hill[spid == chem, hill_gw := pipefit$hill_gw]
  }
  
  spid.ac50 <- paste("AC50 = ", signif(10^dat_hill$hill_ga, digits = 3), "uM", sep = "")
  
  ymax <- max(0.5, max(sub.dt$resp)) + 0.1
  ymin <- min(-0.1, min(sub.dt$resp))
  xval <- min(dat_steve$logc) + 0.1*(max(dat_steve$logc) - min(dat_steve$logc))
  yval <- ymax - 0.1*(ymax-ymin)
  
  # The palette with black:
  cbbPalette <- c("#000000", "#ff0000")
  
  curves_plot <- ggplot(dat_steve, aes(x=logc, y=resp, color = spid))
  
  #plot the fitted curves, color by spid. This code may need to be modified if you have more colors than in cbbPalette
  i <- 1L
  for(chem in dat_hill[,spid]){
    chemical_fit_hill_tp <- dat_hill[spid == chem, hill_tp]
    chemical_fit_hill_ga <- dat_hill[spid == chem, hill_ga]
    chemical_fit_hill_gw <- dat_hill[spid == chem, hill_gw]
    curves_plot <- curves_plot +
      stat_function(fun = hill_curve, args=list(hill_tp = chemical_fit_hill_tp, 
                                                hill_ga = chemical_fit_hill_ga, 
                                                hill_gw = chemical_fit_hill_gw),
                    color = cbbPalette[i], 
                    size=1) 
    i <- i + 1L
  }
  
  curves_plot <- curves_plot +
    theme_bw() +
    geom_point(size = 4) + 
    #geom_hline(yintercept = threshold, , color = "#000000", linetype = 5) +
    #coord_cartesian(ylim=c(-10,175))+
    scale_colour_manual(values=cbbPalette, guide = FALSE) +
    ylab("% Increased Respiration") +
    xlab("log10 uM") +
    geom_text(aes(x = xval, y = yval, label = spid.ac50), size = 4, color = "#000000", show.legend = FALSE) +
    ggtitle(chemtitle) + theme (axis.title=element_text(size=16,face="bold"),
                                plot.title=element_text(size=18,face="bold", hjust = 0.5)) 
  
  ggsave(curves_plot, filename = "DNP.sc.png", path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2", 
         width = 7.5, height = 6.5, units = "in", device = "png")
}

# Hill fit by assay plate id and calculate QC stats

apid.list <- NULL
spid.list <- NULL
AC50.um.list <- NULL
control.median.list <- NULL
control.mad.list <- NULL
vehicle.median.list <- NULL
vehicle.mad.list <- NULL
rZ.list <- NULL
rCV.list <- NULL

for (i in 1:length(unique(dt.plot$apid))){
  
  sub.dt <- NULL
  sub.dt <- dt.plot[apid == unique(dt.plot$apid)[i] & time == "T6",]
  
  dat_steve <- sub.dt
  setkey(dat_steve, spid)
  dat_steve <- dat_steve[,spid := as.factor(spid)]
  
  #Get hill parameters for each spid
  dat_hill <- data.table(spid = unique(dat_steve[,spid]))
  
  for(chem in dat_hill[,spid]){
    pipefit <- tcplFit_Lite(dat_steve[spid == chem,logc], dat_steve[spid == chem,resp], bmad)
    dat_hill[spid == chem, hill_tp := pipefit$hill_tp]
    dat_hill[spid == chem, hill_ga := pipefit$hill_ga]
    dat_hill[spid == chem, hill_gw := pipefit$hill_gw]
  }
  
  dat_hill$AC50.uM <- signif(10^dat_hill$hill_ga, digits = 3)
  
  apid.list <- c(apid.list, as.character(unique(dt.plot$apid)[i]))
  spid.list <- c(spid.list, as.character(unique(dt.plot$spid)))
  AC50.um.list <- c(AC50.um.list, dat_hill$AC50.uM)
  
  control.median <- sub.dt[logc == max(sub.dt$logc), median(rval)]
  control.mad <- sub.dt[logc == max(sub.dt$logc), mad(rval)]
  vehicle.median <- dt.control[apid == unique(sub.dt$apid) & time == "T6" & spid == "DMSO", median(rval)]
  vehicle.mad <- dt.control[apid == unique(sub.dt$apid) & time == "T6" & spid == "DMSO", mad(rval)]
  
  rZ <- 1 - 3 * (control.mad + vehicle.mad) / abs(control.median - vehicle.median)
  rZ.list <- c(rZ.list, rZ)
  
  rCV <- 100 * vehicle.mad / vehicle.median
  rCV.list <- c(rCV.list, rCV)
  
}

dnp.stat.table <- data.frame("apid" = apid.list, "spid" = spid.list, "ac50.um" = AC50.um.list, "rZ" = rZ.list, "rCV" = rCV.list)
write.csv(dnp.stat.table, file = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2/SC.DNP.table.csv", row.names = FALSE)


####### Tier 2 Multi-Concentration Testing ########

### read-in data table
dt0 <- read.csv("L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/MC2/mc_seahorse.lvl0.merged.data.csv", stringsAsFactors = FALSE)
dt0 <- as.data.table(dt0)
dt0$time <- factor(dt0$time, levels = c("T1", "T2", "T3", "T4", "T5", "T6", "T7", "T8", "T9", "T10", "T11", "T12"))

### correct chem stock concnetrations and logc from ChemTracker info
samples <- tcplQuery("SELECT * FROM sample")
dt0$stck.correct <- ifelse(is.na(samples$stkc[match(dt0$spid[1:length(dt0$spid)], samples$spid)]), dt0$stock.conc.mM, samples$stkc[match(dt0$spid[1:length(dt0$spid)], samples$spid)])
dt0$conc.correction <- dt0$conc * dt0$stck.correct / dt0$stock.conc.mM
dt0$logc <- log10(dt0$conc.correction)
dt0$trim.conc <- signif(dt0$conc.correction, digits =3)

#remove wllq = 0, blank wells, convert rval to %norm, give DMSO a logc name, and define a dt w/o DMSO
dt <- dt0[wllq == 1 & spid != "blank", ]

### merge Rotenone (EPAPLT035K15) test with expanded range re-test
dt$spid[dt$spid == "EPAPLT0035K15.1"] <- "EPAPLT0035K15"

dt$norm <- dt$rval * 100

dt$trim.conc[dt$spid == "DMSO"] <- "DMSO"

dt2 <- dt[spid != "DMSO" & wllq == 1, ]
dt2$spid <- as.character(dt2$spid)

### match chmn to spid
dt.chem <- tcplLoadChem(field = "spid", val = unique(dt0$spid), include.spid = TRUE)
dt2$chmn <- ifelse(is.na(dt.chem$chnm[match(dt2$spid[1:length(dt2$spid)], dt.chem$spid)]), as.character(dt2$spid), dt.chem$chnm[match(dt2$spid[1:length(dt2$spid)], dt.chem$spid)])

### subset only spids active from SC2.1 analysis (filter out extra spids included in tier 2 from earlier analysis)
dt2 <- dt2[spid %in% c(SC.active.spids$x, "DNP", "Fenpyroximate"), ]

### Seahorse-style plots by test chemical/control plotted with DMSO

#set by chmn
dt2 <- setkey(dt2, chmn)

for (i in 1:length(unique(dt2$chmn))) {
  
  dt3 <- dt[spid == "DMSO"| spid %in% dt2$spid[dt2$chmn == unique(dt2$chmn)[i]], ]
  dt3 <- dt3[spid != "DNP", ]
  dt3 <- dt3[spid != "Fenpyroximate", ]
  concs <- sort(as.numeric(unique(dt3[spid != "DMSO", trim.conc])))
  concs.order <- c("DMSO", concs)
  dt3$trim.conc <- factor(dt3$trim.conc, levels = concs.order)
  dt4 <- dt3[ , .(median(norm), mad(norm)), by = .(trim.conc, time)]
  
  title <- unique(dt2$chmn)[i]
  title.font <- ifelse(nchar(title) < 50, 18, 14)
  
  palette <- if(length(concs) < 9){c("#000000", "#ff0000", "#ff8800", "#d6b811", "#00dd00", "#00aaff", "#1500ff", "#8000ff")}else{
    c("#000000", "#ff0000", "#ff1c00", "#ff3900", "#ff5500", "#d6b811", "#9cff00", "#55ff00", "#0eff00", "#00ff39",
      "#00ff7f", "#00ffc6", "#00f1ff", "#00aaff", "#0063ff", "#001cff", "#1500ff", "#3900ff", "#5c00ff", "#8000ff")}
  
  summplot <- ggplot(dt4, aes(x = time, y = V1, group = trim.conc, color = trim.conc)) + geom_line(aes(group = trim.conc)) +
    geom_errorbar(aes(ymin = dt4$V1 - dt4$V2, ymax = dt4$V1 + dt4$V2), width = 0.05) +
    geom_point(data = dt4, aes(x = time, y = V1, color = trim.conc), size = 2.5) +
    scale_color_manual(values = palette, labels = concs.order, name = "[Compound] uM") +
    scale_x_discrete(breaks = unique(dt4$time), labels= as.factor(seq(from = 6, to = 72, by = 6))) +
    coord_cartesian(ylim=c(0,275)) +
    geom_segment(aes(x = 3.2, y = 0, xend = 3.2, yend = 266), color = "#555555", linetype = 5) +
    geom_text(aes(x = 3.2, y = 275, label = "Test Compound"), size = 4, color = "#444444", show.legend = FALSE) +
    geom_segment(aes(x = 6.2, y = 0, xend = 6.2, yend = 266), color = "#555555", linetype = 5) +
    geom_text(aes(x = 6.2, y = 275, label = "FCCP"), size = 4, color = "#444444", show.legend = FALSE) +
    geom_segment(aes(x = 9.2, y = 0, xend = 9.2, yend = 266), color = "#555555", linetype = 5) +
    geom_text(aes(x = 9.2, y = 275, label = "ROT/AA"), size = 4, color = "#444444", show.legend = FALSE) +
    geom_text(aes(x = 1.7, y = 245, label = "Pre-Injection"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_text(aes(x = 4.7, y = 245, label = "Basal"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_text(aes(x = 7.7, y = 245, label = "Maximal"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_text(aes(x = 10.7, y = 245, label = "Inhibited"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_segment(aes(x = 0.8, y = 245, xend = 0.4, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 2.6, y = 245, xend = 3.1, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 4.2, y = 245, xend = 3.3, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 5.2, y = 245, xend = 6.1, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 7.2, y = 245, xend = 6.3, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 8.2, y = 245, xend = 9.1, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 10.2, y = 245, xend = 9.3, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    geom_segment(aes(x = 11.2, y = 245, xend = 12.2, yend = 245), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
    theme_bw() +
    ggtitle(title) +
    theme(panel.grid.major = element_line(colour = "grey"), panel.grid.major.x = element_blank()) +
    theme(panel.grid.minor = element_line(colour = "lightgrey")) +
    ylab("% Initial Respiration") +
    xlab("Time (mins)") +
    theme (axis.title=element_text(size=16,face="bold"),
           plot.title=element_text(size=title.font,face="bold", hjust = 0.5),
           legend.text=element_text(size=12), 
           legend.title=element_text(color="#000000"), 
           axis.text = element_text(size = 12)) 
  
  file.name <- paste("MC_chem", i, ".png", sep = "")
  
  ggsave(summplot, filename = file.name, path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/MC2/Seahorse_plots/hi-res", 
         width = 7.5, height = 6.5, units = "in", device = "png")
  
  ggsave(summplot, filename = file.name, path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/MC2/Seahorse_plots/lo-res", 
         width = 7.5, height = 6.5, units = "in", dpi = 150, device = "png")
}

### DMSO-only plot for mechanism figure
dt3 <- dt[spid == "DMSO", ]
concs <- sort(as.numeric(unique(dt3[spid != "DMSO", trim.conc])))
concs.order <- c("DMSO", concs)
dt3$trim.conc <- factor(dt3$trim.conc, levels = concs.order)
dt4 <- dt3[ , .(median(norm), mad(norm)), by = .(trim.conc, time)]

title <- "Mechanistic Determination"
title.font <- ifelse(nchar(title) < 50, 18, 14)

palette <- "#000000"

summplot <- ggplot(dt4, aes(x = time, y = V1, group = trim.conc, color = trim.conc)) + geom_line(aes(group = trim.conc)) +
  geom_errorbar(aes(ymin = dt4$V1 - dt4$V2, ymax = dt4$V1 + dt4$V2), width = 0.05) +
  geom_point(data = dt4, aes(x = time, y = V1, color = trim.conc), size = 2.5) +
  scale_color_manual(values = palette, labels = concs.order, name = "[Compound] uM") +
  scale_x_discrete(breaks = unique(dt4$time), labels= as.factor(seq(from = 6, to = 72, by = 6))) +
  coord_cartesian(ylim=c(0,275)) +
  geom_segment(aes(x = 3.2, y = 0, xend = 3.2, yend = 266), color = "#555555", linetype = 5) +
  geom_text(aes(x = 3.2, y = 275, label = "Test Compound"), size = 4, color = "#444444", show.legend = FALSE) +
  geom_segment(aes(x = 6.2, y = 0, xend = 6.2, yend = 266), color = "#555555", linetype = 5) +
  geom_text(aes(x = 6.2, y = 275, label = "FCCP"), size = 4, color = "#444444", show.legend = FALSE) +
  geom_segment(aes(x = 9.2, y = 0, xend = 9.2, yend = 266), color = "#555555", linetype = 5) +
  geom_text(aes(x = 9.2, y = 275, label = "ROT/AA"), size = 4, color = "#444444", show.legend = FALSE) +
  geom_text(aes(x = 1.7, y = 245, label = "Pre-Injection"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 4.7, y = 245, label = "Basal"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 7.7, y = 245, label = "Maximal"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 10.7, y = 245, label = "Inhibited"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_segment(aes(x = 0.8, y = 245, xend = 0.4, yend = 245), color = "#000000", size = 0.5, 
               arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
  geom_segment(aes(x = 2.6, y = 245, xend = 3.1, yend = 245), color = "#000000", size = 0.5, 
               arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
  geom_segment(aes(x = 4.2, y = 245, xend = 3.3, yend = 245), color = "#000000", size = 0.5, 
               arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
  geom_segment(aes(x = 5.2, y = 245, xend = 6.1, yend = 245), color = "#000000", size = 0.5, 
               arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
  geom_segment(aes(x = 7.2, y = 245, xend = 6.3, yend = 245), color = "#000000", size = 0.5, 
               arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
  geom_segment(aes(x = 8.2, y = 245, xend = 9.1, yend = 245), color = "#000000", size = 0.5, 
               arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
  geom_segment(aes(x = 10.2, y = 245, xend = 9.3, yend = 245), color = "#000000", size = 0.5, 
               arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
  geom_segment(aes(x = 11.2, y = 245, xend = 12.2, yend = 245), color = "#000000", size = 0.5, 
               arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=40)) +
  geom_rect(xmin = 3.4, xmax = 6, ymin = 0, ymax = dt4[time == "T6", V1] * 0.8,   fill = "#f3c982", show.legend = FALSE) +
  geom_rect(xmin = 3.4, xmax = 6, ymin = dt4[time == "T6", V1] * 1.2, ymax = 200,   fill = "#71daea", show.legend = FALSE) +
  geom_rect(xmin = 6.6, xmax = 9, ymin = 0, ymax = dt4[time == "T9", V1] * 0.8,   fill = "#89fbaa", show.legend = FALSE) +
  geom_rect(xmin = 6.6, xmax = 9, ymin = dt4[time == "T9", V1] * 1.2, ymax = 225,   fill = "#aaaaaa", show.legend = FALSE) +
  geom_rect(xmin = 10, xmax = 12.4, ymin = dt4[time == "T12", V1] * 1.3, ymax = 100,   fill = "#e49681", show.legend = FALSE) +
  geom_rect(xmin = 10, xmax = 12.4, ymin = 0, ymax = dt4[time == "T12", V1] * 0.7,   fill = "#aaaaaa", show.legend = FALSE) +
  geom_text(aes(x = 10.2, y = 90, label = "I"), size = 5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 3.6, y = 190, label = "II"), size = 5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 6.8, y = dt4[time == "T9", V1] * 0.8 - 10, label = "III"), size = 5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 3.6, y = dt4[time == "T6", V1] * 0.8 - 10, label = "IV"), size = 5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 4.7, y = 10, label = "aeid = 2442"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 4.7, y = dt4[time == "T6", V1] * 1.2 + 10, label = "aeid = 2443"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 7.8, y = 10, label = "aeid = 2444"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 7.8, y = dt4[time == "T9", V1] * 1.2 + 10, label = "aeid = 2445"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 11.2, y = 10, label = "aeid = 2446"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 11.2, y = dt4[time == "T12", V1] * 1.3 + 10, label = "aeid = 2447"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 4.7, y = (dt4[time == "T6", V1] * 1.2 + 200)/2, label = "Uncoupler"), size = 4.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 4.7, y = dt4[time == "T6", V1] * 0.8/2, label = "ATP synthase"), size = 4.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 7.8, y = dt4[time == "T9", V1] * 0.8/2, label = "ETC inhibitor"), size = 4.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 11.2, y = (dt4[time == "T12", V1] * 1.3 + 100)/2, label = "Redox cycler"), size = 4.5, color = "#000000", show.legend = FALSE) +
  theme_bw() +
  ggtitle(title) +
  theme(panel.grid.major = element_line(colour = "grey"), panel.grid.major.x = element_blank()) +
  theme(panel.grid.minor = element_line(colour = "lightgrey")) +
  ylab("% Initial Respiration") +
  xlab("Time (mins)") +
  theme (axis.title=element_text(size=16,face="bold"),
         plot.title=element_text(size=title.font,face="bold", hjust = 0.5),
         legend.text=element_text(size=12), 
         legend.title=element_text(color="#000000"), 
         axis.text = element_text(size = 12)) 

file.name <- "Mechanism.png"

ggsave(summplot, filename = file.name, path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/MC2/Seahorse_plots", 
       width = 7.5, height = 6.5, units = "in", device = "png")


### Tier 2 (multi-conc) QC ###

#remove wllq = 0, blank wells, convert rval to %norm, give DMSO a logc name, and define a dt w/o DMSO
dt <- dt0[wllq == 1 & spid != "blank", ]
dt$trim.conc[dt$spid == "DMSO"] <- "DMSO"


###plot control dose-reponse data

#Global Plot for Fenpyroximate (all apid together)
dt.control <- dt[wllt %in% c("n", "o", "m"), ] #subset control data
dt.control[, bval := median(rval[time == "T6" & spid == "DMSO"]), by = apid]
dt.control$resp <- (1 - dt.control$rval / dt.control$bval) * 100

dt.plot <- dt.control[spid != "DMSO" & time == "T6", ]
dt.plot$logc <- as.numeric(dt.plot$logc)

for (i in 1:length(unique(dt.plot$spid))){
  
  chemtitle <- unique(dt.plot$spid)[i]
  
  sub.dt <- dt.plot[spid == unique(dt.plot$spid)[i],]
  
  dat_steve <- sub.dt
  setkey(dat_steve, spid)
  dat_steve <- dat_steve[,spid := as.factor(spid)]
  
  #Get hill parameters for each spid
  dat_hill <- data.table(spid = unique(dat_steve[,spid]))
  
  for(chem in dat_hill[,spid]){
    pipefit <- tcplFit_Lite(dat_steve[spid == chem,logc], dat_steve[spid == chem,resp], bmad)
    dat_hill[spid == chem, hill_tp := pipefit$hill_tp]
    dat_hill[spid == chem, hill_ga := pipefit$hill_ga]
    dat_hill[spid == chem, hill_gw := pipefit$hill_gw]
  }
  
  spid.ac50 <- paste("AC50 = ", signif(10^dat_hill$hill_ga, digits = 3), "uM", sep = "")
  
  ymax <- max(0.5, max(sub.dt$resp)) + 0.1
  ymin <- min(-0.1, min(sub.dt$resp))
  xval <- min(dat_steve$logc) + 0.1*(max(dat_steve$logc) - min(dat_steve$logc))
  yval <- ymax - 0.1*(ymax-ymin)
  
  # The palette with black:
  cbbPalette <- c("#000000", "#ff0000")
  
  curves_plot <- ggplot(dat_steve, aes(x=logc, y=resp, color = spid))
  
  #plot the fitted curves, color by spid. This code may need to be modified if you have more colors than in cbbPalette
  i <- 1L
  for(chem in dat_hill[,spid]){
    chemical_fit_hill_tp <- dat_hill[spid == chem, hill_tp]
    chemical_fit_hill_ga <- dat_hill[spid == chem, hill_ga]
    chemical_fit_hill_gw <- dat_hill[spid == chem, hill_gw]
    curves_plot <- curves_plot +
      stat_function(fun = hill_curve, args=list(hill_tp = chemical_fit_hill_tp, 
                                                hill_ga = chemical_fit_hill_ga, 
                                                hill_gw = chemical_fit_hill_gw),
                    color = cbbPalette[i], 
                    size=1) 
    i <- i + 1L
  }
  
  curves_plot <- curves_plot +
    theme_bw() +
    geom_point(size = 4) + 
    #geom_hline(yintercept = threshold, , color = "#000000", linetype = 5) +
    #coord_cartesian(ylim=c(-10,150))+
    scale_colour_manual(values=cbbPalette, guide = FALSE) +
    ylab("% Decreased Respiration") +
    xlab("log10 uM") +
    geom_text(aes(x = xval, y = yval, label = spid.ac50), size = 4, color = "#000000", show.legend = FALSE) +
    ggtitle(chemtitle) + theme (axis.title=element_text(size=16,face="bold"),
                                plot.title=element_text(size=18,face="bold", hjust = 0.5)) 
  
  ggsave(curves_plot, filename = "FENP.mc.png", path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/MC2", 
         width = 7.5, height = 6.5, units = "in", device = "png")
}

# Hill fit by apid and calculate QC stats

apid.list <- NULL
spid.list <- NULL
AC50.um.list <- NULL
control.median.list <- NULL
control.mad.list <- NULL
vehicle.median.list <- NULL
vehicle.mad.list <- NULL
rZ.list <- NULL
rCV.list <- NULL

for (i in 1:length(unique(dt.plot$apid))){
  
  sub.dt <- NULL
  sub.dt <- dt.plot[apid == unique(dt.plot$apid)[i] & time == "T6",]
  
  dat_steve <- sub.dt
  setkey(dat_steve, spid)
  dat_steve <- dat_steve[,spid := as.factor(spid)]
  
  #Get hill parameters for each spid
  dat_hill <- data.table(spid = unique(dat_steve[,spid]))
  
  for(chem in dat_hill[,spid]){
    pipefit <- tcplFit_Lite(dat_steve[spid == chem,logc], dat_steve[spid == chem,resp], bmad)
    dat_hill[spid == chem, hill_tp := pipefit$hill_tp]
    dat_hill[spid == chem, hill_ga := pipefit$hill_ga]
    dat_hill[spid == chem, hill_gw := pipefit$hill_gw]
  }
  
  dat_hill$AC50.uM <- signif(10^dat_hill$hill_ga, digits = 3)
  
  apid.list <- c(apid.list, as.character(unique(dt.plot$apid)[i]))
  spid.list <- c(spid.list, as.character(unique(dt.plot$spid)))
  AC50.um.list <- c(AC50.um.list, dat_hill$AC50.uM)
  
  control.median <- sub.dt[logc == max(sub.dt$logc), median(rval)]
  control.mad <- sub.dt[logc == max(sub.dt$logc), mad(rval)]
  vehicle.median <- dt.control[apid == unique(sub.dt$apid) & time == "T6" & spid == "DMSO", median(rval)]
  vehicle.mad <- dt.control[apid == unique(sub.dt$apid) & time == "T6" & spid == "DMSO", mad(rval)]
  
  rZ <- 1 - 3 * (control.mad + vehicle.mad) / abs(control.median - vehicle.median)
  rZ.list <- c(rZ.list, rZ)
  
  rCV <- 100 * vehicle.mad / vehicle.median
  rCV.list <- c(rCV.list, rCV)
  
}

fenp.stat.table <- data.frame("apid" = apid.list, "spid" = spid.list, "ac50.um" = AC50.um.list, "rZ" = rZ.list, "rCV" = rCV.list)
write.csv(fenp.stat.table, file = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/MC2/MC.FENP.table.csv", row.names = FALSE)

#Global Plot for DNP (all apid together)
dt.control <- dt[wllt %in% c("n", "p", "c"), ] #subset control data
dt.control[, bval := median(rval[time == "T6" & spid == "DMSO"]), by = apid]
dt.control$resp <- (dt.control$rval / dt.control$bval - 1) * 100

dt.plot <- dt.control[spid != "DMSO" & time == "T6", ]
dt.plot$logc <- as.numeric(dt.plot$logc)

for (i in 1:length(unique(dt.plot$spid))){
  
  chemtitle <- "2,4-Dinitrophenol"
  
  sub.dt <- dt.plot[spid == unique(dt.plot$spid)[i],]
  
  dat_steve <- sub.dt
  setkey(dat_steve, spid)
  dat_steve <- dat_steve[,spid := as.factor(spid)]
  
  #Get hill parameters for each spid
  dat_hill <- data.table(spid = unique(dat_steve[,spid]))
  
  for(chem in dat_hill[,spid]){
    pipefit <- tcplFit_Lite(dat_steve[spid == chem,logc], dat_steve[spid == chem,resp], bmad)
    dat_hill[spid == chem, hill_tp := pipefit$hill_tp]
    dat_hill[spid == chem, hill_ga := pipefit$hill_ga]
    dat_hill[spid == chem, hill_gw := pipefit$hill_gw]
  }
  
  spid.ac50 <- paste("AC50 = ", signif(10^dat_hill$hill_ga, digits = 3), "uM", sep = "")
  
  ymax <- max(0.5, max(sub.dt$resp)) + 0.1
  ymin <- min(-0.1, min(sub.dt$resp))
  xval <- min(dat_steve$logc) + 0.1*(max(dat_steve$logc) - min(dat_steve$logc))
  yval <- ymax - 0.1*(ymax-ymin)
  
  # The palette with black:
  cbbPalette <- c("#000000", "#ff0000")
  
  curves_plot <- ggplot(dat_steve, aes(x=logc, y=resp, color = spid))
  
  #plot the fitted curves, color by spid. This code may need to be modified if you have more colors than in cbbPalette
  i <- 1L
  for(chem in dat_hill[,spid]){
    chemical_fit_hill_tp <- dat_hill[spid == chem, hill_tp]
    chemical_fit_hill_ga <- dat_hill[spid == chem, hill_ga]
    chemical_fit_hill_gw <- dat_hill[spid == chem, hill_gw]
    curves_plot <- curves_plot +
      stat_function(fun = hill_curve, args=list(hill_tp = chemical_fit_hill_tp, 
                                                hill_ga = chemical_fit_hill_ga, 
                                                hill_gw = chemical_fit_hill_gw),
                    color = cbbPalette[i], 
                    size=1) 
    i <- i + 1L
  }
  
  curves_plot <- curves_plot +
    theme_bw() +
    geom_point(size = 4) + 
    #geom_hline(yintercept = threshold, , color = "#000000", linetype = 5) +
    #coord_cartesian(ylim=c(-10,175))+
    scale_colour_manual(values=cbbPalette, guide = FALSE) +
    ylab("% Increased Respiration") +
    xlab("log10 uM") +
    geom_text(aes(x = xval, y = yval, label = spid.ac50), size = 4, color = "#000000", show.legend = FALSE) +
    ggtitle(chemtitle) + theme (axis.title=element_text(size=16,face="bold"),
                                plot.title=element_text(size=18,face="bold", hjust = 0.5)) 
  
  ggsave(curves_plot, filename = "DNP.mc.png", path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/MC2", 
         width = 7.5, height = 6.5, units = "in", device = "png")
}

# Hill fit by apid and calculate QC stats

apid.list <- NULL
spid.list <- NULL
AC50.um.list <- NULL
control.median.list <- NULL
control.mad.list <- NULL
vehicle.median.list <- NULL
vehicle.mad.list <- NULL
rZ.list <- NULL
rCV.list <- NULL

for (i in 1:length(unique(dt.plot$apid))){
  
  sub.dt <- NULL
  sub.dt <- dt.plot[apid == unique(dt.plot$apid)[i] & time == "T6",]
  
  dat_steve <- sub.dt
  setkey(dat_steve, spid)
  dat_steve <- dat_steve[,spid := as.factor(spid)]
  
  #Get hill parameters for each spid
  dat_hill <- data.table(spid = unique(dat_steve[,spid]))
  
  for(chem in dat_hill[,spid]){
    pipefit <- tcplFit_Lite(dat_steve[spid == chem,logc], dat_steve[spid == chem,resp], bmad)
    dat_hill[spid == chem, hill_tp := pipefit$hill_tp]
    dat_hill[spid == chem, hill_ga := pipefit$hill_ga]
    dat_hill[spid == chem, hill_gw := pipefit$hill_gw]
  }
  
  dat_hill$AC50.uM <- signif(10^dat_hill$hill_ga, digits = 3)
  
  apid.list <- c(apid.list, as.character(unique(dt.plot$apid)[i]))
  spid.list <- c(spid.list, as.character(unique(dt.plot$spid)))
  AC50.um.list <- c(AC50.um.list, dat_hill$AC50.uM)
  
  control.median <- sub.dt[logc == max(sub.dt$logc), median(rval)]
  control.mad <- sub.dt[logc == max(sub.dt$logc), mad(rval)]
  vehicle.median <- dt.control[apid == unique(sub.dt$apid) & time == "T6" & spid == "DMSO", median(rval)]
  vehicle.mad <- dt.control[apid == unique(sub.dt$apid) & time == "T6" & spid == "DMSO", mad(rval)]
  
  rZ <- 1 - 3 * (control.mad + vehicle.mad) / abs(control.median - vehicle.median)
  rZ.list <- c(rZ.list, rZ)
  
  rCV <- 100 * vehicle.mad / vehicle.median
  rCV.list <- c(rCV.list, rCV)
  
}

dnp.stat.table <- data.frame("apid" = apid.list, "spid" = spid.list, "ac50.um" = AC50.um.list, "rZ" = rZ.list, "rCV" = rCV.list)
write.csv(dnp.stat.table, file = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/MC2/MC.DNP.table.csv", row.names = FALSE)


### Mechanistic Binning and AUF calculations (post-tcpl analysis)###

#load lvl5/6 tcpl analysis
dt_l6 <- as.data.table(read.csv(file = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/MC2/tcpl_files/mc5_mc6_ncct_mito_nov2019.csv", stringsAsFactors = FALSE))

aeids <- 2442:2447

compiled.dt <- data.table(NULL)

for (i in 1:length(aeids)){
  
  temp.dt <- dt_l6[aeid == aeids[i], ]
  temp.dt <- temp.dt[!is.na(chnm), ]
  #subset to SC-actives only (several extra inactive spids were also run in Tier 2)
  temp.dt <- temp.dt[spid %in% SC.active.spids$x, ]
  setkey(temp.dt, spid)
  
  temp.dt[, modl_aic := min(cnst_aic, hill_aic, gnls_aic, na.rm = TRUE), by = spid]
  temp.dt$modl_tp <- ifelse(temp.dt$modl_aic == temp.dt$cnst_aic, 0, ifelse(temp.dt$modl_aic == temp.dt$hill_aic, temp.dt$hill_tp, temp.dt$gnls_tp))
  temp.dt$modl_ga <- ifelse(temp.dt$modl_aic == temp.dt$cnst_aic, 0, ifelse(temp.dt$modl_aic == temp.dt$hill_aic, temp.dt$hill_ga, temp.dt$gnls_ga))
  temp.dt$modl_gw <- ifelse(temp.dt$modl_aic == temp.dt$cnst_aic, 0, ifelse(temp.dt$modl_aic == temp.dt$hill_aic, temp.dt$hill_gw, temp.dt$gnls_gw))
  
  spid <- NULL
  chmn <- NULL
  AUFval <- NULL
  hitc <- NULL
  AC50val <- NULL
  
  min_logc <- min(temp.dt$logc_min, na.rm = TRUE)
  max_logc <- max(temp.dt$logc_max, na.rm = TRUE)
  
  for (j in 1:length(temp.dt$spid)){AUF <- integrate(hill_curve, lower = min_logc, upper = max_logc, 
                                                     hill_tp = temp.dt$modl_tp[temp.dt$spid == temp.dt$spid[j]], 
                                                     hill_ga = temp.dt$modl_ga[temp.dt$spid == temp.dt$spid[j]], 
                                                     hill_gw = temp.dt$modl_gw[temp.dt$spid == temp.dt$spid[j]])
  
  hitc.row <- ifelse(temp.dt$hitc == 1, "active", "inactive")
  
  spid <- c(spid, temp.dt$spid[j])
  chmn <- c(chmn, temp.dt$chnm[temp.dt$spid == temp.dt$spid[j]])
  hitc <- c(hitc, hitc.row)
  AUFval <- c(AUFval, AUF$value)
  AC50val <- c(AC50val, signif(10^temp.dt$modl_ga, digits = 3))
  
  }
  
  AUF.dt <- data.table(aeid = aeids[i], spid = spid, chmn = chmn, hitc = hitc, AUF = AUFval, AC50 = AC50val)
  AUF.summ.dt <- AUF.dt[, .(AUF = max(AUF, na.rm = TRUE), hitc = min(hitc, na.rm = TRUE), AC50 = min(AC50, na.rm = TRUE)), by = .(chmn, aeid)]
  AUF.summ.dt <- AUF.summ.dt[order(AUF.summ.dt$AUF, decreasing = TRUE), ]
  
  file.name <- paste("L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/MC2/", aeids[i], "_AUF", ".csv", sep = "")
  write.csv(AUF.summ.dt, file = file.name, row.names = FALSE)
  
  compiled.dt <- rbind(compiled.dt, AUF.summ.dt)
  
}

write.csv(compiled.dt, file = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/MC2/compiled_AUF_hitc_by_chmn.csv", row.names = FALSE)

# summarize mechanism for each test chemical using scheme higlighted in Figure 4H

summary.dt <- compiled.dt[, .(mechanism = ifelse(hitc[aeid == 2444] == "active" & chmn == "Gentian Violet", "Probe Interference",
                                                 ifelse(hitc[aeid == 2447] == "active", "Redox",
                                                        ifelse(hitc[aeid == 2443] == "active", "Uncoupler", 
                                                               ifelse(hitc[aeid == 2444] == "active", "ETC Inhibitor", 
                                                                      ifelse(hitc[aeid == 2442] == "active", "ATP Synthase Inhibitor", "Inactive"))))),
                              AUF = ifelse(hitc[aeid == 2447] == "active", AUF[aeid == 2447],
                                           ifelse(hitc[aeid == 2443] == "active", AUF[aeid == 2443], 
                                                  ifelse(hitc[aeid == 2444] == "active", AUF[aeid == 2444], 
                                                         ifelse(hitc[aeid == 2442] == "active", AUF[aeid == 2442], 0)))),
                              AC50 = ifelse(hitc[aeid == 2447] == "active", AC50[aeid == 2447],
                                            ifelse(hitc[aeid == 2443] == "active", AC50[aeid == 2443], 
                                                   ifelse(hitc[aeid == 2444] == "active", AC50[aeid == 2444], 
                                                          ifelse(hitc[aeid == 2442] == "active", AC50[aeid == 2442], 0))))), by = chmn]

write.csv(summary.dt, file = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/MC2/summary_mechanism_AUF_by_chmn.csv", row.names = FALSE)



####### Tier 3 Electron Flow Assay (EFA) Testing ########

### read-in data table
dt0 <- read.csv("L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/EFA/EFA.lvl0.merged.data.csv", stringsAsFactors = FALSE)
dt0 <- as.data.table(dt0)
dt0$time <- factor(dt0$time, levels = c("T1", "T2", "T3", "T4", "T5", "T6", "T7", "T8", "T9", "T10", "T11", "T12"))

### correct chem stock concnetrations and logc from ChemTracker info
samples <- tcplQuery("SELECT * FROM sample")
dt0$stck.correct <- ifelse(is.na(samples$stkc[match(dt0$spid[1:length(dt0$spid)], samples$spid)]), dt0$stock.mM, samples$stkc[match(dt0$spid[1:length(dt0$spid)], samples$spid)])
dt0$conc.correction <- dt0$final.uM * dt0$stck.correct / dt0$stock.mM
dt0$logc <- log10(dt0$conc.correction)
dt0$trim.conc <- signif(dt0$conc.correction, digits =3)

#remove wllq = 0, blank wells, give DMSO a logc name, and define a dt w/o DMSO

dt0[apid %in% c("EFA_2.2", "EFA_3.2"), "wllq"] <- 0  # failed plates due to software crash during assay
dt <- dt0[wllq == 1 & !spid %in% c("blank", "empty") ]

dt$trim.conc[dt$spid == "DMSO"] <- "DMSO"

#harmonize trim.conc for duplicate BPA spids for graphing (stck.concs differed slightly)
dt$trim.conc[dt$spid %in% c("EPAPLT0034D10", "EPAPLT0035K22") & dt$final.uM < 40] <- 32
dt$trim.conc[dt$spid %in% c("EPAPLT0034D10", "EPAPLT0035K22") & dt$final.uM > 90] <- 101
dt$trim.conc[dt$spid %in% c("EPAPLT0034D10", "EPAPLT0035K22") & dt$final.uM < 90 & dt$final.uM > 40] <- 65.5

dt2 <- dt[spid != "DMSO" & wllq == 1, ]

### match chmn to spid
dt.chem <- tcplLoadChem(field = "spid", val = unique(dt0$spid), include.spid = TRUE)
dt2$chmn <- ifelse(is.na(dt.chem$chnm[match(dt2$spid[1:length(dt2$spid)], dt.chem$spid)]), as.character(dt2$spid), dt.chem$chnm[match(dt2$spid[1:length(dt2$spid)], dt.chem$spid)])

### subset to only those ETCi spids from MC2 analysis and controls
dt2 <- dt2[spid %in% c("FENP", "MYXO") | chmn %in% as.character(unique(summary.dt$chmn[summary.dt$mechanism == "ETC Inhibitor"])), ]
setkey(dt2, chmn)

### Seahorse-style plots by test chemical/control plotted with DMSO

for (i in 1:length(unique(dt2$chmn))) {
  
  dt3 <- dt[spid == "DMSO" | spid == unique(dt2$spid[dt2$chmn == unique(dt2$chmn)[i]]), ]
  concs <- sort(as.numeric(unique(dt3[spid != "DMSO", trim.conc])))
  concs.order <- c("DMSO", concs)
  dt3$trim.conc <- factor(dt3$trim.conc, levels = concs.order)
  dt4 <- dt3[ , .(median(resp), mad(resp)), by = .(trim.conc, time)]
  
  title <- unique(dt2$chmn)[i]
  title.font <- ifelse(nchar(title) < 50, 18, 14)
  
  ymin <- 0
  ymax <- 150
  y.breaks <- seq(from = ymin, to = ymax, by = 25)
  
  c.shapes <- if(unique(dt3$spid[dt3$spid != "DMSO"])[1] %in% c("FENP", "MYXO")){c(23,0,1,2,19,17,15)}else{c(23,19,17,15)}
  c.palette <- if(unique(dt3$spid[dt3$spid != "DMSO"])[1] %in% c("FENP", "MYXO")){c("#000000", "#dcdcdc", "#c0c0c0", "#a9a9a9", "#959595", "#808080", "#696969")}else{c("#000000", "#ff0000", "#00dd00", "#1500ff")}
  
  summplot <- ggplot(dt4, aes(x = time, y = V1, group = trim.conc, color = trim.conc)) + geom_line(aes(group = trim.conc)) +
    geom_errorbar(aes(ymin = dt4$V1 - dt4$V2, ymax = dt4$V1 + dt4$V2), width = 0.05) +
    geom_point(data = dt4, aes(x = time, y = V1, color = trim.conc, shape = trim.conc), size = 2.5) +
    scale_color_manual(values = c.palette, labels = concs.order, name = "[Compound] uM") +
    scale_shape_manual(values = c.shapes, labels = concs.order, name = "[Compound] uM") +
    scale_y_continuous(limits = c(ymin, ymax), breaks = y.breaks) +
    scale_x_discrete(breaks = unique(dt4$time), labels= as.factor(seq(from = 3, to = 36, by = 3))) +
    geom_segment(aes(x = 2.2, y = ymin, xend = 2.2, yend = ymax - 5), color = "#555555", linetype = 5) +
    geom_text(aes(x = 2.2, y = ymax , label = "Test Compound"), size = 4, color = "#444444", show.legend = FALSE) +
    geom_segment(aes(x = 8.2, y = ymin, xend = 8.2, yend = ymax - 5), color = "#555555", linetype = 5) +
    geom_text(aes(x = 8.2, y = ymax, label = "ROT/SUCC"), size = 4, color = "#444444", show.legend = FALSE) +
    geom_segment(aes(x = 10.2, y = ymin, xend = 10.2, yend = ymax - 15), color = "#555555", linetype = 5) +
    geom_text(aes(x = 10.2, y = ymax - 10, label = "ASC/TMPD/AA"), size = 4, color = "#444444", show.legend = FALSE) +
    geom_text(aes(x = 1.3, y = ymax - 22, label = "Pre-Injection"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_text(aes(x = 5.2, y = ymax - 27, label = "Basal Respiration"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_text(aes(x = 9.2, y = ymax - 27, label = "Complex I"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_text(aes(x = 9.2, y = ymax - 31, label = "Bypass"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_text(aes(x = 11.2, y = ymax - 27, label = "Complex III"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_text(aes(x = 11.2, y = ymax - 31, label = "Bypass"), size = 3, color = "#000000", show.legend = FALSE) +
    geom_segment(aes(x = 3.9, y = ymax - 28, xend = 2.3, yend = ymax - 28), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=30)) +
    geom_segment(aes(x = 6.5, y = ymax - 28, xend = 8.1, yend = ymax - 28), color = "#000000", size = 0.5, 
                 arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=30)) +
    theme_bw() +
    ggtitle(title) +
    theme(panel.grid.major = element_line(colour = "grey"), panel.grid.major.x = element_blank()) +
    theme(panel.grid.minor = element_line(colour = "lightgrey")) +
    ylab("% Initial Respiration") +
    xlab("Time (mins)") +
    theme (axis.title=element_text(size=16,face="bold"),
           plot.title=element_text(size=title.font,face="bold", hjust = 0.5),
           legend.text=element_text(size=12), 
           legend.title=element_text(color="#000000"), 
           axis.text = element_text(size = 12)) 
  
  file.name <- paste("EFA_", i, ".png", sep = "")
  
  ggsave(summplot, filename = file.name, path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/EFA/Seahorse_plots/hi-res", 
         width = 7.5, height = 6.5, units = "in", dpi = 300 ,device = "png")
  ggsave(summplot, filename = file.name, path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/EFA/Seahorse_plots/lo-res", 
         width = 7.5, height = 6.5, units = "in", dpi = 150 ,device = "png")
}

### EFA Mechanism Plot (Figure 5A)

dt3 <- dt[spid == "DMSO", ]
concs <- sort(as.numeric(unique(dt3[spid != "DMSO", trim.conc])))
concs.order <- c("DMSO", concs)
dt3$trim.conc <- factor(dt3$trim.conc, levels = concs.order)
dt4 <- dt3[ , .(median(resp), mad(resp)), by = .(trim.conc, time)]

title <- "EFA Mechanisms"
title.font <- ifelse(nchar(title) < 50, 18, 14)

ymin <- 0
ymax <- 150
y.breaks <- seq(from = ymin, to = ymax, by = 25)

c.shapes <- if(unique(dt3$spid[dt3$spid != "DMSO"])[1] %in% c("FENP", "MYXO")){c(23,0,1,2,19,17,15)}else{c(23,19,17,15)}
c.palette <- if(unique(dt3$spid[dt3$spid != "DMSO"])[1] %in% c("FENP", "MYXO")){c("#000000", "#dcdcdc", "#c0c0c0", "#a9a9a9", "#959595", "#808080", "#696969")}else{c("#000000", "#ff0000", "#00dd00", "#1500ff")}

summplot <- ggplot(dt4, aes(x = time, y = V1, group = trim.conc, color = trim.conc)) + geom_line(aes(group = trim.conc)) +
  geom_errorbar(aes(ymin = dt4$V1 - dt4$V2, ymax = dt4$V1 + dt4$V2), width = 0.05) +
  geom_point(data = dt4, aes(x = time, y = V1, color = trim.conc, shape = trim.conc), size = 2.5) +
  scale_color_manual(values = c.palette, labels = concs.order, name = "Control") +
  scale_shape_manual(values = c.shapes, labels = concs.order, name = "Control") +
  scale_y_continuous(limits = c(ymin, ymax), breaks = y.breaks) +
  scale_x_discrete(breaks = unique(dt4$time), labels= as.factor(seq(from = 3, to = 36, by = 3))) +
  geom_segment(aes(x = 2.2, y = ymin, xend = 2.2, yend = ymax - 5), color = "#555555", linetype = 5) +
  geom_text(aes(x = 2.2, y = ymax , label = "Test Compound"), size = 4, color = "#444444", show.legend = FALSE) +
  geom_segment(aes(x = 8.2, y = ymin, xend = 8.2, yend = ymax - 5), color = "#555555", linetype = 5) +
  geom_text(aes(x = 8.2, y = ymax, label = "ROT/SUCC"), size = 4, color = "#444444", show.legend = FALSE) +
  geom_segment(aes(x = 10.2, y = ymin, xend = 10.2, yend = ymax - 15), color = "#555555", linetype = 5) +
  geom_text(aes(x = 10.2, y = ymax - 10, label = "ASC/TMPD/AA"), size = 4, color = "#444444", show.legend = FALSE) +
  geom_text(aes(x = 1.3, y = ymax - 22, label = "Pre-Injection"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 5.2, y = ymax - 27, label = "Basal Respiration"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 9.2, y = ymax - 27, label = "Complex I"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 9.2, y = ymax - 31, label = "Bypass"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 11.2, y = ymax - 27, label = "Complex III"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 11.2, y = ymax - 31, label = "Bypass"), size = 3, color = "#000000", show.legend = FALSE) +
  geom_segment(aes(x = 3.9, y = ymax - 28, xend = 2.3, yend = ymax - 28), color = "#000000", size = 0.5, 
               arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=30)) +
  geom_segment(aes(x = 6.5, y = ymax - 28, xend = 8.1, yend = ymax - 28), color = "#000000", size = 0.5, 
               arrow = arrow(length = unit(0.1, "inches"), type="closed", angle=30)) +
  geom_rect(mapping=aes(xmin=2.3, xmax=8.1, ymin=25.5, ymax=37), fill="#fae57c", color="black", alpha=0.5) +
  geom_rect(mapping=aes(xmin=8.3, xmax=10.1, ymin=25.5, ymax=49.5), fill="#90ce92", color="black", alpha=0.5) +
  geom_rect(mapping=aes(xmin=2.3, xmax=10.1, ymin=13, ymax=24.5), fill="#ffc0cb", color="black", alpha=0.5) +
  geom_rect(mapping=aes(xmin=2.3, xmax=12.1, ymin=0.5, ymax=12), fill="#b0e0e6", color="black", alpha=0.5) +
  geom_text(aes(x = 5.2, y = 31, label = "Complex I Inhibition"), size = 4, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 9.2, y = 40, label = "Complex II"), size = 4, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 9.2, y = 35, label = "Inhibition"), size = 4, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 6.2, y = 19, label = "Complex III Inhibition"), size = 4, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 7.2, y = 7, label = "Complex IV Inhibition"), size = 4, color = "#000000", show.legend = FALSE) +
  theme_bw() +
  ggtitle(title) +
  theme(panel.grid.major = element_line(colour = "grey"), panel.grid.major.x = element_blank()) +
  theme(panel.grid.minor = element_line(colour = "lightgrey")) +
  ylab("% Initial Respiration") +
  xlab("Time (mins)") +
  theme (axis.title=element_text(size=16,face="bold"),
         plot.title=element_text(size=title.font,face="bold", hjust = 0.5),
         legend.text=element_text(size=12), 
         legend.title=element_text(color="#000000"), 
         axis.text = element_text(size = 12)) 

file.name <- "Mechanism.png"

ggsave(summplot, filename = file.name, path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/EFA/Seahorse_plots", 
       width = 7.5, height = 6.5, units = "in", dpi = 300 ,device = "png")


######### Tox21/ToxCast Inter-assay analysis ########

### Load all Seahorse-tested spids (subset wllt == 't')
sc.lvl0 <- read.csv("L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/SC2/sc_seahorse.lvl0.merged.data.csv")
tested.spids <- unique(sc.lvl0[sc.lvl0$wllt == "t", "spid"])
tested.chems <- tcplLoadChem()[spid %in% tested.spids, ]
test.chem.list <- unique(tested.chems$chnm)

### Load MC2 tests chems by mechanism and determinative aeid
MC2.pos.chems <- summary.dt[summary.dt$mechanism %in% c("ATP Synthase Inhibitor", "ETC Inhibitor", "Uncoupler"), ]

### Pull aeids for mito and tox21_mmp assays (omit viability assays, Tox21_MMP channel-specific and APR mitotic and MEAN)
mito.assays <- tcplLoadAeid(fld = "aenm", val = tcplLoadAeid()$aenm[grep(tcplLoadAeid()$aenm, 
                                                                         pattern = c("MITO"), ignore.case = TRUE)])
mmp.assays <- tcplLoadAeid(fld = "aenm", val = tcplLoadAeid()$aenm[grep(tcplLoadAeid()$aenm, 
                                                                        pattern = "Tox21_MMP", ignore.case = TRUE)])
mito.aeids <- rbind(mito.assays, mmp.assays)
mito.aeids <- mito.aeids[!grep(mito.aeids$aenm, pattern = "mitotic", ignore.case = TRUE),]
mito.aeids <- mito.aeids[!grep(mito.aeids$aenm, pattern = "viability", ignore.case = TRUE),]
mito.aeids <- mito.aeids[!grep(mito.aeids$aenm, pattern = "fitc", ignore.case = TRUE),]
mito.aeids <- mito.aeids[!grep(mito.aeids$aenm, pattern = "rhodamine", ignore.case = TRUE),]
mito.aeids <- mito.aeids[!grep(mito.aeids$aenm, pattern = "mean", ignore.case = TRUE),]

### Non-seahorse aeids
toxcast.aeids <- mito.aeids[!grep(mito.aeids$aenm, pattern = "NCCT"),]

### Condense aeid to assay catagories (exclude Seahorse assays)
assay.catagories <- unique(ifelse(1:length(toxcast.aeids$aenm) %in% grep(x = toxcast.aeids$aenm, pattern = "down"), 
                                  substr(toxcast.aeids$aenm,1,nchar(toxcast.aeids$aenm)- 5) , 
                                  substr(toxcast.aeids$aenm,1,nchar(toxcast.aeids$aenm)- 3)))


### Condense assay catagories to assay technologies
assay.techs <- unique(ifelse(1:length(assay.catagories) %in% grep(x = assay.catagories, pattern = "1hr"), 
                             substr(assay.catagories,1,nchar(assay.catagories)- 4),
                             ifelse(1:length(assay.catagories) %in% grep(x = assay.catagories, pattern = "hr"),
                                    substr(assay.catagories,1,nchar(assay.catagories)- 5),
                                    ifelse(1:length(assay.catagories) %in% grep(x = assay.catagories, pattern = "1h"), 
                                           substr(assay.catagories,1,nchar(assay.catagories)- 3),
                                           ifelse(1:length(assay.catagories) %in% grep(x = assay.catagories, pattern = "h"), 
                                                  substr(assay.catagories,1,nchar(assay.catagories)- 4), assay.catagories)))))



### read in mc3 and mc5 ToxCast data
mc3.dat <- NULL
mc5.dat <- NULL

mc3.dat <-  as.data.table(read.csv("L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/RawMC3_ToxCast_by_aeid.csv",
                                   stringsAsFactors = FALSE))
mc5.dat <- as.data.table(read.csv("L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/RawMC5_ToxCast_by_aeid.csv",
                                  stringsAsFactors = FALSE))


### assay analysis

seahorse.actives <- MC2.pos.chems$chmn
seahorse.inactives <- test.chem.list[!test.chem.list %in% seahorse.actives]

### reference chemicals
ref.set <- read.csv("L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/ref.set.chems.csv", stringsAsFactors = FALSE)

#create seahorse dt
seahorse.dt <- data.table(chnm = ref.set$chnm)
seahorse.dt$assay <- ifelse(seahorse.dt$chnm %in% seahorse.actives, "active", "inactive")
seahorse.dt$assay <- as.factor(seahorse.dt$assay)
seahorse.dt$reference <- match.function(ref.set$activity, seahorse.dt$chnm, ref.set$chnm)
seahorse.dt$reference <- as.factor(seahorse.dt$reference)


### by ASSAY TECH

stat.dt <- NULL

for (i in 1:length(assay.techs)){
  
  temp.dat <- mc5.dat[grep(x = mc5.dat$aenm, pattern = unique(assay.techs)[i]),]
  temp.dat <- temp.dat[chnm %in% ref.set$chnm, ]
  
  aenm <- unique(assay.techs)[i]
  sample.size <- length(unique(temp.dat$chnm))
  
  hitc.chnm <- temp.dat[, .(hitc = ceiling(mean(hitc)), ga = if(sum(hitc) > 0){min(modl_ga[hitc ==1])}else{3}), by = chnm]
  hitc.chnm$ac50 <- ifelse(hitc.chnm$hitc == 1, 10^hitc.chnm$ga, NA)
  dt.hitc  <- hitc.chnm[, .(assay = ifelse(hitc > 0, "active", "inactive")), by = chnm]
  
  dt.hitc$reference <- match.function(ref.set$activity, dt.hitc$chnm, ref.set$chnm)
  dt.hitc$assay <- as.factor(dt.hitc$assay)
  dt.hitc$reference <- as.factor(dt.hitc$reference)
  
  
  assay.actives <- hitc.chnm[hitc == 1, chnm]
  assay.inactives <- hitc.chnm[hitc == 0, chnm]
  
  TP <- assay.actives[assay.actives %in% dt.hitc[reference == "active", chnm]]
  FP <- assay.actives[assay.actives %in% dt.hitc[reference == "inactive", chnm]]
  TN <- assay.inactives[assay.inactives %in% dt.hitc[reference == "inactive", chnm]]
  FN <- assay.inactives[assay.inactives %in% dt.hitc[reference == "active", chnm]]
  TPn <- as.numeric(length(TP))
  FPn <- as.numeric(length(FP))
  TNn <- as.numeric(length(TN))
  FNn <- as.numeric(length(FN))
  MCC <- (TPn * TNn - FPn * FNn) / sqrt((TPn + FNn) * (TNn + FPn) * (TPn + FPn) * (TNn + FNn))
  
  conf.matrix <- confusionMatrix(dt.hitc$assay, dt.hitc$reference)
  
  temp.row <- data.table(sample.size = sample.size, true.pos = TPn, true.neg = TNn, false.pos = FPn, false.neg = FNn, 
                         accuracy = conf.matrix$overall[1], 
                         kappa = conf.matrix$overall[2], 
                         no.information.rate = conf.matrix$overall[5], 
                         sensitivity = round(conf.matrix$byClass[1], digits = 3),
                         specificity = round(conf.matrix$byClass[2], digits = 3), 
                         pos.pred.value = conf.matrix$byClass[3],
                         neg.pred.value = conf.matrix$byClass[4],
                         precision = conf.matrix$byClass[5],
                         balanced.accuracy = round(conf.matrix$byClass[11], digits = 3),
                         MCC = round(MCC, digits = 3))
  
  temp.dt <- cbind(aenm, temp.row)
  stat.dt <- rbind(stat.dt, temp.dt)
  
}

tox21.ref.ac50 <- hitc.chnm

aenm <- "NCCT_Seahorse_RSA"

assay.actives <- seahorse.dt[assay == "active", chnm]
assay.inactives <- seahorse.dt[assay == "inactive", chnm]

TP <- assay.actives[assay.actives %in% dt.hitc[reference == "active", chnm]]
FP <- assay.actives[assay.actives %in% dt.hitc[reference == "inactive", chnm]]
TN <- assay.inactives[assay.inactives %in% dt.hitc[reference == "inactive", chnm]]
FN <- assay.inactives[assay.inactives %in% dt.hitc[reference == "active", chnm]]
TPn <- as.numeric(length(TP))
FPn <- as.numeric(length(FP))
TNn <- as.numeric(length(TN))
FNn <- as.numeric(length(FN))
MCC <- (TPn * TNn - FPn * FNn) / sqrt((TPn + FNn) * (TNn + FPn) * (TPn + FPn) * (TNn + FNn))

conf.matrix <- confusionMatrix(seahorse.dt$assay, seahorse.dt$reference)

temp.row <- data.table(sample.size = sample.size, true.pos = TPn, true.neg = TNn, false.pos = FPn, false.neg = FNn, 
                       accuracy = conf.matrix$overall[1], 
                       kappa = conf.matrix$overall[2], 
                       no.information.rate = conf.matrix$overall[5], 
                       sensitivity = round(conf.matrix$byClass[1], digits = 3),
                       specificity = round(conf.matrix$byClass[2], digits = 3), 
                       pos.pred.value = conf.matrix$byClass[3],
                       neg.pred.value = conf.matrix$byClass[4],
                       precision = conf.matrix$byClass[5],
                       balanced.accuracy = round(conf.matrix$byClass[11], digits = 3),
                       MCC = round(MCC, digits = 3))

temp.dt <- cbind(aenm, temp.row)
stat.dt <- rbind(stat.dt, temp.dt)

table5 <- stat.dt[, c("aenm", "sample.size", "true.pos", "true.neg", "false.pos", "false.neg", "sensitivity", "specificity", "balanced.accuracy", "MCC")]
table5 <- table5[order(balanced.accuracy, decreasing = TRUE), ]
setnames(table5, c("Assay Technology", "Reference Chemicals Tested", "True Pos.", "True Neg.", "False Pos.", "False Neg.", "Sensitivity", "Specificity", "Balanced Accuracy", "MCC"))
write.csv(table5, file = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/Table5.csv", row.names = FALSE)

### Tox21 v RFA AC50 table
AC50table <- tox21.ref.ac50[, c("chnm", "ac50")]
AC50table$seahorse <- match.function(MC2.pos.chems$AC50, AC50table$chnm, MC2.pos.chems$chmn)
table6 <- AC50table[chnm %in% ref.set[ref.set$activity == "active", ]$chnm, ]
table6 <- table6[!is.na(ac50) | !is.na(seahorse),]
table6$delta.log.ac50 <- ifelse(is.na(table6$ac50), 2.5 - log10(table6$seahorse), ifelse(is.na(table6$seahorse), log10(table6$ac50) - 2.5, log10(table6$ac50) - log10(table6$seahorse)))
table6$delta.log.abs <- abs(table6$delta.log.ac50)
table6 <- table6[order(delta.log.ac50, decreasing = TRUE), ]
names(table6) <- c("Refernce Chemical", "Tox21_MMP AC50 (uM)", "NCCT_RSA AC50 (uM)", "Delta Log AC50", "ABS Delta Log AC50")
write.csv(table6, file = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/Table6.csv", row.names = FALSE)

### Figure 6: Global (1,042 chems) Tox21 v RFA AC50
seahorse.pos <- MC2.pos.chems[, c("chmn" , "AC50")]
seahorse.neg <- data.frame(chmn = seahorse.inactives, AC50 = 316)
seahorse.all <- rbind(seahorse.pos, seahorse.neg)
seahorse.all$RSA.ga <- log10(seahorse.all$AC50)
seahorse.all.ga.only <- seahorse.all[, c("chmn", "RSA.ga")]

tox21.dat <- mc5.dat[grep(x = mc5.dat$aenm, pattern = "TOX21_MMP_ratio"),]
tox21.dat <- tox21.dat[chnm %in% seahorse.all$chmn, ]

tox21.ga <- tox21.dat[, .(hitc = ceiling(mean(hitc)), tox21.ga = if(sum(hitc) > 0){min(modl_ga[hitc ==1])}else{2.5}), by = chnm]
tox21.ga.only <- tox21.ga[, c("chnm", "tox21.ga")]

merge.ga <- merge.data.frame(seahorse.all.ga.only, tox21.ga.only, by.x = "chmn", by.y = "chnm")
merge.ga$ref <- ifelse(merge.ga$chmn %in% ref.set$chnm, "b", "a")
merge.ga$mech <- ifelse(merge.ga$chmn %in% MC2.pos.chems$chmn[MC2.pos.chems$mechanism == "ETC Inhibitor"], "ETC Inhibitor",
                        ifelse(merge.ga$chmn %in% MC2.pos.chems$chmn[MC2.pos.chems$mechanism == "Uncoupler"], "Uncoupler",
                               ifelse(merge.ga$chmn %in% MC2.pos.chems$chmn[MC2.pos.chems$mechanism == "ATP Synthase Inhibitor"], "ATP Synthase Inhibitor", 
                                      "RSA Inactive")))
merge.ga$mech <- factor(merge.ga$mech, levels = c("ETC Inhibitor", "Uncoupler", "ATP Synthase Inhibitor", "RSA Inactive"))

figure6a <- ggplot(merge.ga, aes(x = RSA.ga, y = tox21.ga, color = ref)) + geom_point(shape = 15) + 
  geom_abline(slope = 1, intercept = 0) + 
  geom_abline(slope = 1, intercept = 0.5, linetype = 5) + 
  geom_abline(slope = 1, intercept = -0.5, linetype = 5) + 
  scale_color_manual(values = c("#aaaaaa", "#000000")) +
  scale_x_continuous(breaks = c(-1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5)) +
  scale_y_continuous(breaks = c(-1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5)) +
  coord_cartesian(xlim=(c(-1.25,2.5)), ylim=(c(-1.25,2.5))) +
  theme_bw() +
  ggtitle("") +
  theme(panel.grid.major = element_line(colour = "grey"), panel.grid.major.x = element_blank(), panel.grid.major.y = element_blank()) +
  theme(panel.grid.minor = element_line(colour = "lightgrey")) +
  ylab("Tox21 MMP AC50 (log10 uM)") +
  xlab("RSA AC50 (log10 uM)") +
  theme (axis.title=element_text(size= 16,face="bold"),
         plot.title=element_text(size= title.font,face="bold", hjust = 0.5),
         legend.position = "none", 
         axis.text = element_text(size = 12)) +
  geom_text(aes(x = -1.24, y = 1.05, label = "1"), size = 3.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 2.45, y = 0.02, label = "2"), size = 3.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 0.62, y = 2.13, label = "3"), size = 3.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 1.52, y = 0.24, label = "4"), size = 3.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 0.8, y = 2.57, label = "5"), size = 3.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = -0.23, y = 1.04, label = "6"), size = 3.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 0.12, y = 1.35, label = "7"), size = 3.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = -0.32, y = 0.9, label = "8"), size = 3.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = -0.11, y = 1.06, label = "9"), size = 3.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 1.63, y = 0.49, label = "10"), size = 3.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 0.76, y = 0.05, label = "11"), size = 3.5, color = "#000000", show.legend = FALSE) +
  geom_text(aes(x = 0.21, y = 0.94, label = "12"), size = 3.5, color = "#000000", show.legend = FALSE)


ggsave(figure6a, filename = "figure6a.png", path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/ToxCast_plots", 
       width = 7.5, height = 6.5, units = "in", dpi = 300 ,device = "png")



figure6b <- ggplot(merge.ga, aes(x = RSA.ga, y = tox21.ga, color = mech, shape = mech)) + geom_point(size =2) + 
  geom_abline(slope = 1, intercept = 0) + 
  geom_abline(slope = 1, intercept = 0.5, linetype = 5) + 
  geom_abline(slope = 1, intercept = -0.5, linetype = 5) + 
  scale_color_manual(values = c("#ff0000", "#0000ff", "#00ff00","#000000")) +
  scale_shape_manual(values = c(18,17,16,15)) +
  scale_x_continuous(breaks = c(-1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5)) +
  scale_y_continuous(breaks = c(-1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5)) +
  coord_cartesian(xlim=(c(-1.25,2.5)), ylim=(c(-1.25,2.5))) +
  theme_bw() +
  ggtitle("") +
  theme(panel.grid.major = element_line(colour = "grey"), panel.grid.major.x = element_blank(), panel.grid.major.y = element_blank()) +
  theme(panel.grid.minor = element_line(colour = "lightgrey")) +
  ylab("Tox21 MMP AC50 (log10 uM)") +
  xlab("RSA AC50 (log10 uM)") +
  theme (axis.title=element_text(size= 16,face="bold"),
         plot.title=element_text(size= title.font,face="bold", hjust = 0.5),
         legend.text=element_text(size=10), 
         legend.title = element_blank(), 
         axis.text = element_text(size = 12))


ggsave(figure6b, filename = "figure6b.png", path = "L:/Lab/NCCT_Lab/Simmons_Lab/Projects/OXPHOS/ToxCast/ToxCast_plots", 
       width = 9.35, height = 6.5, units = "in", dpi = 300 ,device = "png")