##----------------------------------------------------------------------------##
## R Script for Analyses comprising data for the manuscript
## "Mechanistic Toxicity Tests Based on an Adverse Outcome Pathway Network for Hepatic Steatosis"
##  
##
## Author -- Michelle Angrish
## Email  -- angrish.michelle@epa.gov
##
## January 19, 2017
##----------------------------------------------------------------------------##

library(data.table)
library(plyr)
library(dplyr)
library(stringr)

setwd("set-working-directory-to-folder-containig-supplementary-files")

lda <- read.csv(file.path("Supporting_file_1.csv"))
pa <- read.csv(file.path("Supporting_file_2.csv"))
fau <- read.csv(file.path("Supporting_file_3.csv"))

fae <- read.csv(file.path("Supporting_file_4.csv"))
fao <- read.csv(file.path("Supporting_file_5.csv"))
LDH <- read.csv(file.path("Supporting_file_6.csv"))

################################################################################
################################################################################
### sTART: CALCULATE MANUSCRIPT LDH DATA

###assign column names
colnames(LDH) <- c("chid", "stock.conc", "conc", "well", "LDH")

### Calculate DMSO med
DMSO <- LDH[LDH$chid == "DMSO",]
medDMSO <- median(DMSO$LDH)

###Calculate FC and log2 transform into positive scale
LDH$FCLDH <- abs((LDH$LDH)/(medDMSO))
LDH$logc <- log(LDH$conc)
LDH$log2FC <- log2(LDH$FCLDH)
LDH$abslog2FC <- abs(LDH$log2FC)
LDH$chid_conc <- str_c(LDH$chid, sep = "_", LDH$conc)

###calculate DMSO mad
dmso.mad <- LDH[LDH$chid == 'DMSO',]
madDMSO <- mad(dmso.mad$abslog2FC)

###convert to data.table
dat <- as.data.table(LDH)

###remove control rows
dat.samp <- dat[chid != 'DMSO',]

###calculate median response by sample
med.resp <- dat[,median(abslog2FC), by=chid_conc]
setnames(med.resp,c("chid_conc","resp"))
med.resp <- na.omit(med.resp)

### set BasalThreshold and calcualte hits
BMAD <-  4.025 * madDMSO # 3 * madDMSO
LDH.hits <- med.resp[resp > (BMAD),]

################################################################################
################################################################################
### sTART: CALCULATE MANUSCRIPT LIPID DROPLET ACCUMULATION DATA

###assign column names
colnames(lda) <- c("chid", "stock.conc", "conc", "well", "RFP")

### Calculate DMSO med
DMSO <- lda[lda$chid == "DMSO",]
medDMSO <- median(DMSO$RFP)

###Calculate FC and log2 transform into positive scale
lda$FCRFP <- abs((lda$RFP)/(medDMSO))
lda$logc <- log(lda$conc)
lda$log2FC <- log2(lda$FCRFP)
lda$abslog2FC <- abs(lda$log2FC)
lda$chid_conc <- str_c(lda$chid, sep = "_", lda$conc)

###calculate baseline DMSO mad
dmso.bmad <- lda[lda$chid == 'DMSO',]
bmadDMSO <- mad(dmso.bmad$abslog2FC)

###convert to data.table
dat <- as.data.table(lda)

###remove control rows
dat.samp <- dat[chid != 'DMSO',]

###calculate median response by sample
med.resp <- dat[,median(abslog2FC), by=chid_conc]
setnames(med.resp,c("chid_conc","resp"))
med.resp <- na.omit(med.resp)

### set BasalThreshold and calcualte hits
threshold <-  3 * bmadDMSO # 3 * madDMSO
lda.hits <- med.resp[resp > (threshold),]


################################################################################
################################################################################
### sTART: CALCULATE MANUSCRIPT PHOSPHOLIPID ACCUMULATION DATA

###assign column names
colnames(pa) <- c("chid", "stock.conc", "conc", "well", "RFP")

### Calculate DMSO med
DMSO <- pa[pa$chid == "DMSO",]
medDMSO <- median(DMSO$RFP)

###normalize, Calculate FC and log2 transform into positive scale
pa$FCRFP <- abs((pa$RFP)/(medDMSO))  
pa$logc <- log(pa$conc) 
pa$log2FC <- log2(pa$FCRFP) 
pa$abslog2FC <- abs(pa$log2FC) 
pa$chid_conc <- str_c(pa$chid, sep = "_", pa$conc) 

###calculate baseline DMSO mad
dmso.bmad <- pa[pa$chid == 'DMSO',] 
bmadDMSO <- mad(dmso.bmad$abslog2FC) 

###convert to data.table
dat <- as.data.table(pa)

###remove control rows
dat.samp <- dat[chid != 'DMSO',]

###calculate median response by sample
med.resp <- dat[,median(abslog2FC), by=chid_conc]
setnames(med.resp,c("chid_conc","resp"))
med.resp <- na.omit(med.resp)

### set BasalThreshold and calcualte hits
threshold <-  3 * bmadDMSO # 3 * madDMSO
pa.hits <- med.resp[resp > (threshold),]

################################################################################
################################################################################
### sTART: CALCULATE MANUSCRIPT FATTY ACID UPTAKE DATA

###assign column names
colnames(fau) <- c("chid", "stock.conc", "conc", "well", "FFA")

### Calculate DMSO med
DMSO <- fau[fau$chid == "DMSO",]
medDMSO <- median(DMSO$FFA)

###Calculate FC and log2 transform into positive scale
fau$FCFFA <- abs((fau$FFA)/(medDMSO))
fau$logc <- log(fau$conc)
fau$log2FC <- log2(fau$FCFFA)
fau$abslog2FC <- abs(fau$log2FC)
fau$chid_conc <- str_c(fau$chid, sep = "_", fau$conc)

###calculate DMSO mad
dmso.bmad <- fau[fau$chid == 'DMSO',]
bmadDMSO <- mad(dmso.bmad$abslog2FC)

###convert to data.table
dat <- as.data.table(fau)

###remove control rows
dat.samp <- dat[chid != 'DMSO',]

###calculate median response by sample
med.resp <- dat[,median(abslog2FC), by=chid_conc]
setnames(med.resp,c("chid_conc","resp"))
med.resp <- na.omit(med.resp)


### set BasalThreshold and calcualte hits
threshold <-  3 * bmadDMSO # 3 * madDMSO
fau.hits <- med.resp[resp > (threshold),]


################################################################################
################################################################################
### sTART: CALCULATE MANUSCRIPT FATTY ACID EFFLUX DATA


###assign column names
colnames(fae) <- c("chid", "stock.conc", "conc", "well", "ApoB")

### Calculate DMSO med
DMSO <- fae[fae$chid == "DMSO",]
medDMSO <- median(DMSO$ApoB)

###Calculate FC and log2 transform into positive scale
fae$FCApoB <- abs((fae$ApoB)/(medDMSO))
fae$logc <- log(fae$conc)
fae$log2FC <- log2(fae$FCApoB)
fae$abslog2FC <- abs(fae$log2FC)
fae$chid_conc <- str_c(fae$chid, sep = "_", fae$conc)


###calculate DMSO mad
dmso.bmad <- fae[fae$chid == 'DMSO',]
bmadDMSO <- mad(dmso.bmad$abslog2FC)

###convert to data.table
dat <- as.data.table(fae)

###remove control rows
dat.samp <- dat[chid != 'DMSO',]

###calculate median response by sample
med.resp <- dat[,median(abslog2FC), by=chid_conc]
setnames(med.resp,c("chid_conc","resp"))
med.resp <- na.omit(med.resp)

### set BasalThreshold and calcualte hits
threshold <-  3 * bmadDMSO # 3 * madDMSO
fae.hits <- med.resp[resp > (threshold),]


################################################################################
################################################################################
### sTART: CALCULATE MANUSCRIPT FATTY ACID OXIDATION DATA


FAOXdat <- fao

###assign column names
colnames(FAOXdat) <- c("chid","pretreatment", "stock.conc", "conc", "well", "T1", "T2", "T3", "T4", "T5", "T6", "T7", "T8", "T9", "T10", "T11", "T12", "T13", "T14", "T15")

### Extract BSA and PALM data
BSA <- FAOXdat[FAOXdat$pretreatment == "BSA",]
PALM <- FAOXdat[FAOXdat$pretreatment == "PALM",]

####subset data by BSA and PALM pretreatment
BSAdat <- BSA[, c("chid", "pretreatment", "T1","T2", "T3", "T4", "T5", "T6", "T7", "T8", "T9", "T10", "T11", "T12", 
                  "T13", "T14", "T15")]
PALMdat <- PALM[, c("chid", "pretreatment", "T1","T2", "T3", "T4", "T5", "T6", "T7", "T8", "T9", "T10", "T11", "T12", 
                    "T13", "T14", "T15")]

### calcualte the median DMSO of ROT/AA (T15) for BSA and PALM separately.

BSAdat.table <- as.data.table(BSAdat)
PALMdat.table <- as.data.table(PALMdat)

###Calculate median for all bsa samples
med.BSA <- BSAdat.table[,median(T15), by=chid]
setnames(med.BSA,c("chid","med.BSAT15"))

###Calculate median for all PALM samples
med.PALM <- PALMdat.table[,median(T15), by=chid]
setnames(med.PALM,c("chid","med.PALMT15"))

####subset value for bsa and palm med DMSO at T15
medDMSO.BSA <- med.BSA[which(med.BSA$chid == "DMSO")]
medDMSO.palm <- med.PALM[which(med.PALM$chid == "DMSO")]

### baseline to median ROT/AA (T15) of DMSO control wells and scale 0 to 100
### Do this for the BSA and PALM separately
BSA.DMSOT15med <- medDMSO.BSA$med.BSAT15
PALM.DMSOT15med <- medDMSO.palm$med.PALMT15

BSAdat$T1scale <- (BSAdat$T1-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T2scale <- (BSAdat$T2-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T3scale <- (BSAdat$T3-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T4scale <- (BSAdat$T4-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T5scale <- (BSAdat$T5-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T6scale <- (BSAdat$T6-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T7scale <- (BSAdat$T7-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T8scale <- (BSAdat$T8-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T9scale <- (BSAdat$T9-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T10scale <- (BSAdat$T10-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T11scale <- (BSAdat$T11-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T12scale <- (BSAdat$T12-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T13scale <- (BSAdat$T13-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T14scale <- (BSAdat$T14-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100
BSAdat$T15scale <- (BSAdat$T15-BSA.DMSOT15med)/(100-BSA.DMSOT15med) *100

PALMdat$T1scale <- (PALMdat$T1-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T2scale <- (PALMdat$T2-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T3scale <- (PALMdat$T3-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T4scale <- (PALMdat$T4-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T5scale <- (PALMdat$T5-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T6scale <- (PALMdat$T6-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T7scale <- (PALMdat$T7-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T8scale <- (PALMdat$T8-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T9scale <- (PALMdat$T9-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T10scale <- (PALMdat$T10-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T11scale <- (PALMdat$T11-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T12scale <- (PALMdat$T12-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T13scale <- (PALMdat$T13-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T14scale <- (PALMdat$T14-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100
PALMdat$T15scale <- (PALMdat$T15-PALM.DMSOT15med)/(100-PALM.DMSOT15med) *100

### calcualte the median DMSO of ROT/AA (T15) for BSA and PALM separately.

BSAdat.table <- as.data.table(BSAdat)
PALMdat.table <- as.data.table(PALMdat)

###calculate DMSO mad
dmso.madPALM <- PALMdat[PALMdat$chid == 'DMSO',]
madDMSO.PALM <- mad(dmso.madPALM$T6scale)

dmso.madBSA <- BSAdat[BSAdat$chid == 'DMSO',]
madDMSO.BSA <- mad(dmso.madBSA$T6scale)

###remove control rows
#dat.samp <- dat[chid != 'blank',]
#dat.samp <- dat.samp[chid != 'DMSO',]
#dat.samp <- dat.samp[chid != 'DNP',]
#dat.samp <- dat.samp[chid != 'Fenpyroximate',]

###calculate median response by sample
medocr.palm <- PALMdat.table[,median(T6scale), by=chid]
setnames(medocr.palm,c("chid","med.OCR"))

medocr.BSA <- BSAdat.table[,median(T6scale), by=chid]
setnames(medocr.BSA,c("chid","med.OCR"))

### calculate BasalThreshold
Basal.madPALM <- PALMdat.table[chid == 'DMSO',mad(T6scale)]
BasalThreshold.palm <- 3 * madDMSO.PALM #10 * Basal.mad
ocr.hits.uppalm <- medocr.palm[med.OCR > (100 + BasalThreshold.palm),]
ocr.hits.downpalm <- medocr.palm[med.OCR < (100 - BasalThreshold.palm),]

Basal.madBSA <- BSAdat.table[chid == 'DMSO',mad(T6scale)]
BasalThreshold.BSA <- 3 * madDMSO.BSA #10 * Basal.mad
ocr.hits.upBSA <- medocr.BSA[med.OCR > (100 + BasalThreshold.BSA),]
ocr.hits.downBSA <- medocr.BSA[med.OCR < (100 - BasalThreshold.BSA),]

###calculate MaxResp.median and MaxRespThreshold
MaxResp.median.palm <- PALMdat.table[chid == 'DMSO',median(T12scale)]
MaxResp.mad.palm <- PALMdat.table[chid == 'DMSO',mad(T12scale)]
MaxRespThreshold.palm <- 1.5 * MaxResp.mad.palm

MaxResp.median.BSA <- BSAdat.table[chid == 'DMSO',median(T12scale)]
MaxResp.mad.BSA <- BSAdat.table[chid == 'DMSO',mad(T12scale)]
MaxRespThreshold.BSA <- 1.5 * MaxResp.mad.BSA

### calculate hits using MaxRespThreshold
med.MR.palm <- PALMdat.table[,median(T12scale), by=chid]
setnames(med.MR.palm,c("chid","med.MR"))
MR.hits.up.palm <- med.MR.palm[med.MR > (MaxResp.median.palm + MaxRespThreshold.palm),]
MR.hits.down.palm <- med.MR.palm[med.MR < (MaxResp.median.palm - MaxRespThreshold.palm),]

med.MR.BSA <- BSAdat.table[,median(T12scale), by=chid]
setnames(med.MR.BSA,c("chid","med.MR"))
MR.hits.up.BSA <- med.MR.BSA[med.MR > (MaxResp.median.BSA + MaxRespThreshold.BSA),]
MR.hits.down.BSA <- med.MR.BSA[med.MR < (MaxResp.median.BSA - MaxRespThreshold.BSA),]


################################################################################
################################################################################
### sTART: CALCULATE MANUSCRIPT GENE EXPRESSOIN DATA

#####delta delta Ct method#####
## delta Ct1 = Ct(TargetA_treated) - Ct(RefB-treated)
## delta Ct2 = Ct(TargetA_control) - Ct(RefB_control)
## delta delta Ct = delta Ct1(treated) - delta Ct2(control)
## Normalized target gene expression level = 2^delta delta Ct

## Create a sample annotation file and read it in
### rbind Ct data

ActbFasn <- fread(file.path("Supporting_file_7.csv"))
Cd36Gapdh <- fread(file.path("Supporting_file_8.csv"))
Dgat1Dgat2 <- fread(file.path("Supporting_file_9.csv"))
HprtSrebf1 <- fread(file.path("Supporting_file_10.csv"))
MttpPp1a <- fread(file.path("Supporting_file_11.csv"))
TbpNrl13 <- fread(file.path("Supporting_file_12.csv"))
ApoBCpt1a <- fread(file.path("Supporting_file_13.csv"))
AcacaPparg <- fread(file.path("Supporting_file_14.csv"))

PCRdata <- rbind(ActbFasn, Cd36Gapdh, Dgat1Dgat2, HprtSrebf1,  MttpPp1a, TbpNrl13, ApoBCpt1a, AcacaPparg)


###assign column names
colnames(PCRdata) <- c("well", "Fluor", "Target", "Sample", "Ct", "Ctmean", "Conc", "Chid")

###Use Gapdh as house keeping gene
###Calculate the delta Ct1 for all target genes
Fasn <- PCRdata[PCRdata$Target == "Fasn",]##extract FASN data
Gapdh <- PCRdata[PCRdata$Target == "Gapdh",]
### calculate delta Ct1
Fasn$FasndeltaCt1 <- (Fasn$Ct - Gapdh$Ct)

###Calculate delta Ct2
Gapdhcontrol <- Gapdh[Gapdh$Sample == "0.1 DMSO",]  # extract Gapdh (ref) control (DMSO) data
Fasncontrol <- Fasn[Fasn$Sample == "0.1 DMSO",] # extract Fasn (target) control data
Fasn$FasndeltaCt2 <- (Fasncontrol$Ct - Gapdhcontrol$Ct)  ##calcualte delta Ct2
### deltaCT2 <- Gapdhcontrol  ##create a data.frame to store delta Ct2 data
FasndeltaCt2.med <- median(Fasn$FasndeltaCt2)

Fasn$deltadeltaCt <- (Fasn$FasndeltaCt1) - (FasndeltaCt2.med)  #calculate the delta delta Ct
Fasn$FC <- 2^abs(Fasn$deltadeltaCt) # Calculate the fold change
Fasn$log2FC <- abs(log2(Fasn$FC))

#extract DMSO data and calculate the median
DMSOFasn <- Fasn[Fasn$Sample == "0.1 DMSO",]
medDMSOFASN <- median(DMSOFasn$log2FC)

################################################################################

####Calcualte hits
###calculate DMSO mad
bmadDMSO.fasn <- mad(DMSOFasn$log2FC)

###convert to data.table
dat.fasn <- as.data.table(Fasn)

###assign column names
colnames(dat.fasn) <- c("well", "Fluor", "Target", "Sample", "Ct", "Ctmean", "Conc", "Chid", "FasndeltaCt1", "deltadeltaCt", "FC", "abslog2FC", "log2FC")

###remove control rows
dat.fasn.1 <- dat.fasn[dat.fasn$Chid != 'DMSO',]
dat.fasn.2 <- dat.fasn.1[dat.fasn.1$Chid != 'untreated',]

###calculate median response by sample
med.resp.fasn <- dat.fasn.2[,median(log2FC), by=Sample]
setnames(med.resp.fasn,c("Sample","resp"))
med.resp.fasn <- na.omit(med.resp.fasn)

### set Basal median absolute deviation and calcualte hits
threshold.fasn <-  3 * bmadDMSO.fasn # 3 * madDMSO
Fasn.hits <- med.resp.fasn[resp > (threshold.fasn),]

#################################################################################
################################################################################
Cd36 <- PCRdata[PCRdata$Target == "Cd36",]##extract Cd36 data

### calculate delta Ct1
Cd36$Cd36deltaCt1 <- (Cd36$Ct - Gapdh$Ct)

###Calculate delta Ct2
Gapdhcontrol <- Gapdh[Gapdh$Sample == "0.1 DMSO",]  # extract Gapdh (ref) control (DMSO) data
Cd36control <- Cd36[Cd36$Sample == "0.1 DMSO",] # extract Cd36 (target) control data
Cd36$Cd36deltaCt2 <- (Cd36control$Ct - Gapdhcontrol$Ct)  ##calcualte delta Ct2
### deltaCT2 <- Gapdhcontrol  ##create a data.frame to store delta Ct2 data
Cd36deltaCt2.med <- median(Cd36$Cd36deltaCt2)

Cd36$deltadeltaCt <- (Cd36$Cd36deltaCt1) - (Cd36deltaCt2.med)  #calculate the delta delta Ct
Cd36$FC <- 2^abs(Cd36$deltadeltaCt) # Calculate the fold change
Cd36$log2FC <- abs(log2(Cd36$FC))

#extract DMSO data and calculate the median
DMSOCd36 <- Cd36[Cd36$Sample == "0.1 DMSO",]
medDMSOCd36 <- median(DMSOCd36$log2FC)

################################################################################

####Calcualte hits
###calculate DMSO mad
bmadDMSO.Cd36 <- mad(DMSOCd36$log2FC)

###convert to data.table
dat.Cd36 <- as.data.table(Cd36)

###assign column names
colnames(dat.Cd36) <- c("well", "Fluor", "Target", "Sample", "Ct", "Ctmean", "Conc", "Chid", "Cd36deltaCt1", "deltadeltaCt", "FC", "abslog2FC", "log2FC")

###remove control rows
dat.Cd36.1 <- dat.Cd36[dat.Cd36$Chid != 'DMSO',]
dat.Cd36.2 <- dat.Cd36.1[dat.Cd36.1$Chid != 'untreated',]

###calculate median response by sample
med.resp.Cd36 <- dat.Cd36.2[,median(log2FC), by=Sample]
setnames(med.resp.Cd36,c("Sample","resp"))
med.resp.Cd36 <- na.omit(med.resp.Cd36)

### set Basal median absolute deviation and calcualte hits
threshold.Cd36 <-  3 * bmadDMSO.Cd36 # 3 * madDMSO
Cd36.hits <- med.resp.Cd36[resp > (threshold.Cd36),]

###################################################################################
################################################################################
####Dgat1

Dgat1 <- PCRdata[PCRdata$Target == "Dgat1",]##extract Dgat1 data

### calculate delta Ct1
Dgat1$Dgat1deltaCt1 <- (Dgat1$Ct - Gapdh$Ct)

###Calculate delta Ct2
Gapdhcontrol <- Gapdh[Gapdh$Sample == "0.1 DMSO",]  # extract Gapdh (ref) control (DMSO) data
Dgat1control <- Dgat1[Dgat1$Sample == "0.1 DMSO",] # extract Dgat1 (target) control data
Dgat1$Dgat1deltaCt2 <- (Dgat1control$Ct - Gapdhcontrol$Ct)  ##calcualte delta Ct2
### deltaCT2 <- Gapdhcontrol  ##create a data.frame to store delta Ct2 data
Dgat1deltaCt2.med <- median(Dgat1$Dgat1deltaCt2)

Dgat1$deltadeltaCt <- (Dgat1$Dgat1deltaCt1) - (Dgat1deltaCt2.med)  #calculate the delta delta Ct
Dgat1$FC <- 2^abs(Dgat1$deltadeltaCt) # Calculate the fold change
Dgat1$log2FC <- abs(log2(Dgat1$FC))

#extract DMSO data and calculate the median
DMSODgat1 <- Dgat1[Dgat1$Sample == "0.1 DMSO",]
medDMSODgat1 <- median(DMSODgat1$log2FC)

################################################################################

####Calcualte hits
###calculate DMSO mad
bmadDMSO.Dgat1 <- mad(DMSODgat1$log2FC)

###convert to data.table
dat.Dgat1 <- as.data.table(Dgat1)

###assign column names
colnames(dat.Dgat1) <- c("well", "Fluor", "Target", "Sample", "Ct", "Ctmean", "Conc", "Chid", "Dgat1deltaCt1", "deltadeltaCt", "FC", "abslog2FC", "log2FC")

###remove control rows
dat.Dgat1.1 <- dat.Dgat1[dat.Dgat1$Chid != 'DMSO',]
dat.Dgat1.2 <- dat.Dgat1.1[dat.Dgat1.1$Chid != 'untreated',]

###calculate median response by sample
med.resp.Dgat1 <- dat.Dgat1.2[,median(log2FC), by=Sample]
setnames(med.resp.Dgat1,c("Sample","resp"))
med.resp.Dgat1 <- na.omit(med.resp.Dgat1)

### set BasalThreshold and calcualte hits
threshold.Dgat1 <-  3 * bmadDMSO.Dgat1 # 3 * madDMSO
Dgat1.hits <- med.resp.Dgat1[resp > (threshold.Dgat1),]

##################################################################################
################################################################################

####Dgat2

Dgat2 <- PCRdata[PCRdata$Target == "Dgat2",]##extract Dgat2 data

### calculate delta Ct1
Dgat2$Dgat2deltaCt1 <- (Dgat2$Ct - Gapdh$Ct)

###Calculate delta Ct2
Gapdhcontrol <- Gapdh[Gapdh$Sample == "0.1 DMSO",]  # extract Gapdh (ref) control (DMSO) data
Dgat2control <- Dgat2[Dgat2$Sample == "0.1 DMSO",] # extract Dgat2 (target) control data
Dgat2$Dgat2deltaCt2 <- (Dgat2control$Ct - Gapdhcontrol$Ct)  ##calcualte delta Ct2
### deltaCT2 <- Gapdhcontrol  ##create a data.frame to store delta Ct2 data
Dgat2deltaCt2.med <- median(Dgat2$Dgat2deltaCt2)

Dgat2$deltadeltaCt <- (Dgat2$Dgat2deltaCt1) - (Dgat2deltaCt2.med)  #calculate the delta delta Ct
Dgat2$FC <- 2^abs(Dgat2$deltadeltaCt) # Calculate the fold change
Dgat2$log2FC <- abs(log2(Dgat2$FC))

#extract DMSO data and calculate the median
DMSODgat2 <- Dgat2[Dgat2$Sample == "0.1 DMSO",]
medDMSODgat2 <- median(DMSODgat2$log2FC)

################################################################################

####Calcualte hits
###calculate DMSO mad
bmadDMSO.Dgat2 <- mad(DMSODgat2$log2FC)

###convert to data.table
dat.Dgat2 <- as.data.table(Dgat2)

###assign column names
colnames(dat.Dgat2) <- c("well", "Fluor", "Target", "Sample", "Ct", "Ctmean", "Conc", "Chid", "Dgat2deltaCt1", "deltadeltaCt", "FC", "abslog2FC", "log2FC")

###remove control rows
dat.Dgat2.1 <- dat.Dgat2[dat.Dgat2$Chid != 'DMSO',]
dat.Dgat2.2 <- dat.Dgat2.1[dat.Dgat2.1$Chid != 'untreated',]

###calculate median response by sample
med.resp.Dgat2 <- dat.Dgat2.2[,median(log2FC), by=Sample]
setnames(med.resp.Dgat2,c("Sample","resp"))
med.resp.Dgat2 <- na.omit(med.resp.Dgat2)

### set BasalThreshold and calcualte hits
threshold.Dgat2 <-  3 * bmadDMSO.Dgat2  # 3 * madDMSO
Dgat2.hits <- med.resp.Dgat2[resp > (threshold.Dgat2),]

################################################################################

####Srebf1

Srebf1 <- PCRdata[PCRdata$Target == "Srebf1",]##extract Srebf1 data

### calculate delta Ct1
Srebf1$Srebf1deltaCt1 <- (Srebf1$Ct - Gapdh$Ct)

###Calculate delta Ct2
Gapdhcontrol <- Gapdh[Gapdh$Sample == "0.1 DMSO",]  # extract Gapdh (ref) control (DMSO) data
Srebf1control <- Srebf1[Srebf1$Sample == "0.1 DMSO",] # extract Srebf1 (target) control data
Srebf1$Srebf1deltaCt2 <- (Srebf1control$Ct - Gapdhcontrol$Ct)  ##calcualte delta Ct2
### deltaCT2 <- Gapdhcontrol  ##create a data.frame to store delta Ct2 data
Srebf1deltaCt2.med <- median(Srebf1$Srebf1deltaCt2)

Srebf1$deltadeltaCt <- (Srebf1$Srebf1deltaCt1) - (Srebf1deltaCt2.med)  #calculate the delta delta Ct
Srebf1$FC <- 2^abs(Srebf1$deltadeltaCt) # Calculate the fold change
Srebf1$log2FC <- abs(log2(Srebf1$FC))

#extract DMSO data and calculate the median
DMSOSrebf1 <- Srebf1[Srebf1$Sample == "0.1 DMSO",]
medDMSOSrebf1 <- median(DMSOSrebf1$log2FC)

################################################################################

####Calcualte hits
###calculate DMSO mad
bmadDMSO.Srebf1 <- mad(DMSOSrebf1$log2FC)

###convert to data.table
dat.Srebf1 <- as.data.table(Srebf1)

###assign column names
colnames(dat.Srebf1) <- c("well", "Fluor", "Target", "Sample", "Ct", "Ctmean", "Conc", "Chid", "Srebf1deltaCt1", "deltadeltaCt", "FC", "abslog2FC", "log2FC")

###remove control rows
dat.Srebf1.1 <- dat.Srebf1[dat.Srebf1$Chid != 'DMSO',]
dat.Srebf1.2 <- dat.Srebf1.1[dat.Srebf1.1$Chid != 'untreated',]

###calculate median response by sample
med.resp.Srebf1 <- dat.Srebf1.2[,median(log2FC), by=Sample]
setnames(med.resp.Srebf1,c("Sample","resp"))
med.resp.Srebf1 <- na.omit(med.resp.Srebf1)

### set BasalThreshold and calcualte hits
threshold.Srebf1 <-  3 * bmadDMSO.Srebf1# 3 * madDMSO

################################################################################
####Tbp

Tbp <- PCRdata[PCRdata$Target == "TBP",]##extract Tbp data

### calculate delta Ct1
Tbp$TbpdeltaCt1 <- (Tbp$Ct - Gapdh$Ct)

###Calculate delta Ct2
Gapdhcontrol <- Gapdh[Gapdh$Sample == "0.1 DMSO",]  # extract Gapdh (ref) control (DMSO) data
Tbpcontrol <- Tbp[Tbp$Sample == "0.1 DMSO",] # extract Tbp (target) control data
Tbp$TbpdeltaCt2 <- (Tbpcontrol$Ct - Gapdhcontrol$Ct)  ##calcualte delta Ct2
### deltaCT2 <- Gapdhcontrol  ##create a data.frame to store delta Ct2 data
TbpdeltaCt2.med <- median(Tbp$TbpdeltaCt2)

Tbp$deltadeltaCt <- (Tbp$TbpdeltaCt1) - (TbpdeltaCt2.med)  #calculate the delta delta Ct
Tbp$FC <- 2^abs(Tbp$deltadeltaCt) # Calculate the fold change
Tbp$log2FC <- abs(log2(Tbp$FC))

#extract DMSO data and calculate the median
DMSOTbp <- Tbp[Tbp$Sample == "0.1 DMSO",]
medDMSOTbp <- median(DMSOTbp$log2FC)

################################################################################

####Calcualte hits
###calculate DMSO mad
bmadDMSO.Tbp <- mad(DMSOTbp$log2FC)

###convert to data.table
dat.Tbp <- as.data.table(Tbp)

###assign column names
colnames(dat.Tbp) <- c("well", "Fluor", "Target", "Sample", "Ct", "Ctmean", "Conc", "Chid", "TbpdeltaCt1", "deltadeltaCt", "FC", "abslog2FC", "log2FC")

###remove control rows
dat.Tbp.1 <- dat.Tbp[dat.Tbp$Chid != 'DMSO',]
dat.Tbp.2 <- dat.Tbp.1[dat.Tbp.1$Chid != 'untreated',]

###calculate median response by sample
med.resp.Tbp <- dat.Tbp.2[,median(log2FC), by=Sample]
setnames(med.resp.Tbp,c("Sample","resp"))
med.resp.Tbp <- na.omit(med.resp.Tbp)

### set BasalThreshold and calcualte hits
threshold.Tbp <-  3 * bmadDMSO.Tbp # 3 * madDMSO
Tbp.hits <- med.resp.Tbp[resp > (threshold.Tbp),]

################################################################################

####Nrl13

Nrl13 <- PCRdata[PCRdata$Target == "NR1H3",]##extract Nrl13 data

### calculate delta Ct1
Nrl13$Nrl13deltaCt1 <- (Nrl13$Ct - Gapdh$Ct)

###Calculate delta Ct2
Gapdhcontrol <- Gapdh[Gapdh$Sample == "0.1 DMSO",]  # extract Gapdh (ref) control (DMSO) data
Nrl13control <- Nrl13[Nrl13$Sample == "0.1 DMSO",] # extract Nrl13 (target) control data
Nrl13$Nrl13deltaCt2 <- (Nrl13control$Ct - Gapdhcontrol$Ct)  ##calcualte delta Ct2
### deltaCT2 <- Gapdhcontrol  ##create a data.frame to store delta Ct2 data
Nrl13deltaCt2.med <- median(Nrl13$Nrl13deltaCt2)

Nrl13$deltadeltaCt <- (Nrl13$Nrl13deltaCt1) - (Nrl13deltaCt2.med)  #calculate the delta delta Ct
Nrl13$FC <- 2^abs(Nrl13$deltadeltaCt) # Calculate the fold change
Nrl13$log2FC <- abs(log2(Nrl13$FC))

#extract DMSO data and calculate the median
DMSONrl13 <- Nrl13[Nrl13$Sample == "0.1 DMSO",]
medDMSONrl13 <- median(DMSONrl13$log2FC)

#######################################################################################
####Calcualte hits
###calculate DMSO mad
bmadDMSO.Nrl13 <- mad(DMSONrl13$log2FC)

###convert to data.table
dat.Nrl13 <- as.data.table(Nrl13)

###assign column names
colnames(dat.Nrl13) <- c("well", "Fluor", "Target", "Sample", "Ct", "Ctmean", "Conc", "Chid", "Nrl13deltaCt1", "deltadeltaCt", "FC", "abslog2FC", "log2FC")

###remove control rows
dat.Nrl13.1 <- dat.Nrl13[dat.Nrl13$Chid != 'DMSO',]
dat.Nrl13.2 <- dat.Nrl13.1[dat.Nrl13.1$Chid != 'untreated',]

###calculate median response by sample
med.resp.Nrl13 <- dat.Nrl13.2[,median(log2FC), by=Sample]
setnames(med.resp.Nrl13,c("Sample","resp"))
med.resp.Nrl13 <- na.omit(med.resp.Nrl13)

### set BasalThreshold and calcualte hits
threshold.Nrl13 <-  3 * bmadDMSO.Nrl13 # 3 * madDMSO
Nrl13.hits <- med.resp.Nrl13[resp> (threshold.Nrl13),]

################################################################################

###Acaca
Acaca <- PCRdata[PCRdata$Target == "ACACA",]##extract Acaca data
### calculate delta Ct1
Acaca$AcacadeltaCt1 <- (Acaca$Ct - Gapdh$Ct)

###Calculate delta Ct2
Gapdhcontrol <- Gapdh[Gapdh$Sample == "0.1 DMSO",]  # extract Gapdh (ref) control (DMSO) data
Acacacontrol <- Acaca[Acaca$Sample == "0.1 DMSO",] # extract Acaca (target) control data
Acaca$AcacadeltaCt2 <- (Acacacontrol$Ct - Gapdhcontrol$Ct)  ##calcualte delta Ct2
### deltaCT2 <- Gapdhcontrol  ##create a data.frame to store delta Ct2 data
AcacadeltaCt2.med <- median(Acaca$AcacadeltaCt2)

Acaca$deltadeltaCt <- (Acaca$AcacadeltaCt1) - (AcacadeltaCt2.med)  #calculate the delta delta Ct
Acaca$FC <- 2^abs(Acaca$deltadeltaCt) # Calculate the fold change
Acaca$log2FC <- abs(log2(Acaca$FC))

#extract DMSO data and calculate the median
DMSOAcaca <- Acaca[Acaca$Sample == "0.1 DMSO",]
medDMSOAcaca <- median(DMSOAcaca$log2FC)

################################################################################

####Calcualte hits
###calculate DMSO mad
bmadDMSO.Acaca <- mad(DMSOAcaca$log2FC)

###convert to data.table
dat.Acaca <- as.data.table(Acaca)

###assign column names
colnames(dat.Acaca) <- c("well", "Fluor", "Target", "Sample", "Ct", "Ctmean", "Conc", "Chid", "AcacadeltaCt1", "deltadeltaCt", "FC", "abslog2FC", "log2FC")

###remove control rows
dat.Acaca.1 <- dat.Acaca[dat.Acaca$Chid != 'DMSO',]
dat.Acaca.2 <- dat.Acaca.1[dat.Acaca.1$Chid != 'untreated',]

###calculate median response by sample
med.resp.Acaca <- dat.Acaca.2[,median(log2FC), by=Sample]
setnames(med.resp.Acaca,c("Sample","resp"))
med.resp.Acaca <- na.omit(med.resp.Acaca)

### set BasalThreshold and calcualte hits
threshold.Acaca <-  3 * bmadDMSO.Acaca# 3 * madDMSO
Acaca.hits <- med.resp.Acaca[resp > (threshold.Acaca),]

#################################################################################
################################################################################

###Pparg
Pparg <- PCRdata[PCRdata$Target == "PPARG",]##extract Pparg data
### calculate delta Ct1
Pparg$PpargdeltaCt1 <- (Pparg$Ct - Gapdh$Ct)

###Calculate delta Ct2
Gapdhcontrol <- Gapdh[Gapdh$Sample == "0.1 DMSO",]  # extract Gapdh (ref) control (DMSO) data
Ppargcontrol <- Pparg[Pparg$Sample == "0.1 DMSO",] # extract Pparg (target) control data
Pparg$PpargdeltaCt2 <- (Ppargcontrol$Ct - Gapdhcontrol$Ct)  ##calcualte delta Ct2
### deltaCT2 <- Gapdhcontrol  ##create a data.frame to store delta Ct2 data
PpargdeltaCt2.med <- median(Pparg$PpargdeltaCt2)

Pparg$deltadeltaCt <- (Pparg$PpargdeltaCt1) - (PpargdeltaCt2.med)  #calculate the delta delta Ct
Pparg$FC <- 2^abs(Pparg$deltadeltaCt) # Calculate the fold change
Pparg$log2FC <- abs(log2(Pparg$FC))

#extract DMSO data and calculate the median
DMSOPparg <- Pparg[Pparg$Sample == "0.1 DMSO",]
medDMSOPparg <- median(DMSOPparg$log2FC)

################################################################################

####Calcualte hits
###calculate DMSO mad
bmadDMSO.Pparg <- mad(DMSOPparg$log2FC)

###convert to data.table
dat.Pparg <- as.data.table(Pparg)

###assign column names
colnames(dat.Pparg) <- c("well", "Fluor", "Target", "Sample", "Ct", "Ctmean", "Conc", "Chid", "PpargdeltaCt1", "deltadeltaCt", "FC", "abslog2FC", "log2FC")

###remove control rows
dat.Pparg.1 <- dat.Pparg[dat.Pparg$Chid != 'DMSO',]
dat.Pparg.2 <- dat.Pparg.1[dat.Pparg.1$Chid != 'untreated',]

###calculate median response by sample
med.resp.Pparg <- dat.Pparg.2[,median(log2FC), by=Sample]
setnames(med.resp.Pparg,c("Sample","resp"))
med.resp.Pparg <- na.omit(med.resp.Pparg)

### set BasalThreshold and calcualte hits
threshold.Pparg <-  3 * bmadDMSO.Pparg # 3 * madDMSO
Pparg.hits <- med.resp.Pparg[resp > (threshold.Pparg),]

#################################################################################
################################################################################


###Apob
Apob <- PCRdata[PCRdata$Target == "ApoB100",]##extract Apob data
### calculate delta Ct1
Apob$ApobdeltaCt1 <- (Apob$Ct - Gapdh$Ct)

###Calculate delta Ct2
Gapdhcontrol <- Gapdh[Gapdh$Sample == "0.1 DMSO",]  # extract Gapdh (ref) control (DMSO) data
Apobcontrol <- Apob[Apob$Sample == "0.1 DMSO",] # extract Apob (target) control data
Apob$ApobdeltaCt2 <- (Apobcontrol$Ct - Gapdhcontrol$Ct)  ##calcualte delta Ct2
### deltaCT2 <- Gapdhcontrol  ##create a data.frame to store delta Ct2 data
ApobdeltaCt2.med <- median(Apob$ApobdeltaCt2)

Apob$deltadeltaCt <- (Apob$ApobdeltaCt1) - (ApobdeltaCt2.med)  #calculate the delta delta Ct
Apob$FC <- 2^abs(Apob$deltadeltaCt) # Calculate the fold change
Apob$log2FC <- abs(log2(Apob$FC))

#extract DMSO data and calculate the median
DMSOApob <- Apob[Apob$Sample == "0.1 DMSO",]
medDMSOApob <- median(DMSOApob$log2FC)

################################################################################

####Calcualte hits
###calculate DMSO mad
bmadDMSO.Apob <- mad(DMSOApob$log2FC)

###convert to data.table
dat.Apob <- as.data.table(Apob)

###assign column names
colnames(dat.Apob) <- c("well", "Fluor", "Target", "Sample", "Ct", "Ctmean", "Conc", "Chid", "ApobdeltaCt1", "deltadeltaCt", "FC", "abslog2FC", "log2FC")

###remove control rows
dat.Apob.1 <- dat.Apob[dat.Apob$Chid != 'DMSO',]
dat.Apob.2 <- dat.Apob.1[dat.Apob.1$Chid != 'untreated',]

###calculate median response by sample
med.resp.Apob <- dat.Apob.2[,median(log2FC), by=Sample]
setnames(med.resp.Apob,c("Sample","resp"))
med.resp.Apob <- na.omit(med.resp.Apob)

### set BasalThreshold and calcualte hits
threshold.Apob <-  3 * bmadDMSO.Apob# 3 * madDMSO
Apob.hits <- med.resp.Apob[resp > (threshold.Apob),]

##############################################################################################################
##############################################################################################################

###Cpt1a
Cpt1a <- PCRdata[PCRdata$Target == "Cpt1a",]##extract Cpt1a data
### calculate delta Ct1
Cpt1a$Cpt1adeltaCt1 <- (Cpt1a$Ct - Gapdh$Ct)

###Calculate delta Ct2
Gapdhcontrol <- Gapdh[Gapdh$Sample == "0.1 DMSO",]  # extract Gapdh (ref) control (DMSO) data
Cpt1acontrol <- Cpt1a[Cpt1a$Sample == "0.1 DMSO",] # extract Cpt1a (target) control data
Cpt1a$Cpt1adeltaCt2 <- (Cpt1acontrol$Ct - Gapdhcontrol$Ct)  ##calcualte delta Ct2
### deltaCT2 <- Gapdhcontrol  ##create a data.frame to store delta Ct2 data
Cpt1adeltaCt2.med <- median(Cpt1a$Cpt1adeltaCt2)

Cpt1a$deltadeltaCt <- (Cpt1a$Cpt1adeltaCt1) - (Cpt1adeltaCt2.med)  #calculate the delta delta Ct
Cpt1a$FC <- 2^abs(Cpt1a$deltadeltaCt) # Calculate the fold change
Cpt1a$log2FC <- abs(log2(Cpt1a$FC))

#extract DMSO data and calculate the median
DMSOCpt1a <- Cpt1a[Cpt1a$Sample == "0.1 DMSO",]
medDMSOCpt1a <- median(DMSOCpt1a$log2FC)

################################################################################

####Calcualte hits
###calculate DMSO mad
bmadDMSO.Cpt1a <- mad(DMSOCpt1a$log2FC)

###convert to data.table
dat.Cpt1a <- as.data.table(Cpt1a)

###assign column names
colnames(dat.Cpt1a) <- c("well", "Fluor", "Target", "Sample", "Ct", "Ctmean", "Conc", "Chid", "Cpt1adeltaCt1", "deltadeltaCt", "FC", "abslog2FC", "log2FC")

###remove control rows
dat.Cpt1a.1 <- dat.Cpt1a[dat.Cpt1a$Chid != 'DMSO',]
dat.Cpt1a.2 <- dat.Cpt1a.1[dat.Cpt1a.1$Chid != 'untreated',]

###calculate median response by sample
med.resp.Cpt1a <- dat.Cpt1a.2[,median(log2FC), by=Sample]
setnames(med.resp.Cpt1a,c("Sample","resp"))
med.resp.Cpt1a <- na.omit(med.resp.Cpt1a)

### set BasalThreshold and calcualte hits
threshold.Cpt1a <-  3 * bmadDMSO.Cpt1a  # 3 * madDMSO
Cpt1a.hits <- med.resp.Cpt1a[resp > (threshold.Cpt1a),]

############################################################################################################
###########################################################################################################