This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
# Housekeeping: removes all global environment data without removing functions
rm(list = setdiff(ls(), lsf.str()))
# ** Note: READ BEFORE RUNNING SCRIPT ** ####
# It is necessary to run the "Process Normalized Data" section below ONLY ONCE to generate up-to-date summary files:
# Once these final summary files have been generated within the working directory:
# All subsequent sections may be run WITHOUT RE-RUNNING the "process normalized data" section
# To collapse each section for ease in navigation and running separately:
# click the SINGLE arrow to the left of each section name such that it is pointing to the right
##############################################
# Set Working Directory/ Load All Packages/ Set Memory ####
# Defines the working directory as the folder in which all data files are located
# For this script to work, place the script file in the root folder and double click to open (DO NOT DRAG AND DROP INTO R STUDIO)
# For Recursive processing, the R.script must be placed in the root folder containing all other folders
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
library(xlsx)
library(xlsxjars)
library(rJava)
library(stringr)
library(gplots)
library(ggplot2)
library(plyr)
library(XLConnect)
library(XLConnectJars)
library(xtable)
library(tcltk)
library(tcltk2)
library(MESS)
library(gridExtra)
library(gridBase)
library(zoo)
library(grid)
# Change default Java memory to 4 Gigabytes
options(java.parameters = "-Xmx4g")
##########################################################
# Set all functions ####
argmax <- function(x, y, w=1, ...) {
require(zoo)
n <- length(y)
y.smooth <- loess(y ~ x, ...)$fitted
y.max <- rollapply(zoo(y.smooth), 2*w+1, max, align="center")
delta <- y.max - y.smooth[-c(1:w, n+1-1:w)]
i.max <- which(delta <= 0) + w
list(x=x[i.max], i=i.max, y.hat=y.smooth)
}
argmin <- function(x, y, w=1, ...) {
require(zoo)
n <- length(y)
y.smooth <- loess(y ~ x, ...)$fitted
y.min <- rollapply(zoo(y.smooth), 2*w+1, min, align="center")
delta <- y.min - y.smooth[-c(1:w, n+1-1:w)]
i.min <- which(delta >= 0) + w
list(x=x[i.min], i=i.min, y.hat=y.smooth)
}
# this is the core function responsible for smoothing normalized data and identifying peaks.
# The output is a table containing summary data pertaining to each dataset
test <- function(w, span, thres) {
peaks <- argmax(x, y, w=w, span=span )
nadirs <- argmin(x, y, w=w, span=span)
# # # # # # #
# define peak sizes as smoothed y maxima - smoothed y minima for all minima and maxima identified
sizes = sapply(peaks$i, function(peak.i) {
nadir.i = nadirs$i[max(which(nadirs$i < peak.i))]
size = peaks$y.hat[peak.i] - nadirs$y.hat[nadir.i]
# I Added this line as a "conditional argument for missing values": size is now always numeric to avoid problems when 'calling' "[sizes > thres]"
if(is.na(size)){size = 0}
return(size)
})
# # # # # # #
# Identify the first & last peak
# infinite (missing) values are converted to NA to assist in table cleanup and omission of non-responsive cells
first.peakx <- min(peaks$x[sizes > thres], na.rm=T)
if(!is.infinite(first.peakx)){first.peakx <- first.peakx}else{first.peakx = NA}
last.peakx <- max(peaks$x[sizes > thres], na.rm=T)
if(!is.infinite(last.peakx)){last.peakx <- last.peakx}else{last.peakx = NA}
# Identify the first local minima point immediately prior to the first identified peak
# if no peaks exist, this value will be designated NA
first.pre.troughX<-x[nadirs$i][max(which(x[nadirs$i]<first.peakx))]
if(!is.infinite(first.pre.troughX)){first.pre.troughX <- first.pre.troughX}else{first.pre.troughX = NA}
######
# Calculate and collect the Average Peak Durations & last.post.troughX:
# i.e.) the x value of local minima immediately following the final peak)
peakdur <- 0
preT.vector<-c()
if(length(sizes[sizes>thres])>0){
for(Q in 1:length(sizes[sizes>thres])){
if(Q==1){
preT <- x[nadirs$i][max(which(x[nadirs$i]<peaks$x[sizes > thres][Q]))]
posT <- x[nadirs$i][min(which(x[nadirs$i]>peaks$x[sizes > thres][Q]))]
preT.vector<-c(preT)
if(!is.na(preT) & !is.na(posT))
{
peakdur <- peakdur + (posT - preT)
}
}
if(Q>1){
preT<-x[nadirs$i][max(which(x[nadirs$i]<peaks$x[sizes > thres][Q]))]
posT<-x[nadirs$i][min(which(x[nadirs$i]>peaks$x[sizes > thres][Q]))]
if( preT < peaks$x[sizes > thres][Q-1])
{
peakdur<-peakdur
} else {
preT=preT
posT=posT
if(!is.na(preT) & !is.na(posT))
{
peakdur <- peakdur + (posT - preT)
}
preT.vector<-c(preT.vector, preT)
}
}
}
# the code above collects and adds each peaks' respsective duration based on the local minima points derived pre and post each peak
# derive the final identified local minima point immediately following the final peak
last.post.troughX<-x[nadirs$i][min(which(x[nadirs$i]>peaks$x[sizes > thres][length(sizes[sizes>thres])]))]
# if any peak is found near the end of a time lapse video such that a local minima was not defined after it (rare occasions):
# The last local minima point will be defined as the final x (time) value: ( i.e. last.post.troughX=x[length(x)] )
if(is.na(last.post.troughX)){
last.post.troughX=x[length(x)]
peakdur<-peakdur + (last.post.troughX - preT.vector[length(preT.vector)])
} else {
last.post.troughX<-last.post.troughX
peakdur<-peakdur
}
# calculate the average peak duration based on the number of peaks and the total sum peak duration for all peaks (derived above)
avgpeakdur <- peakdur/length(sizes[sizes>thres])
} else {
peakdur = NA
avgpeakdur = NA
last.post.troughX = NA
}
# Calculate First amplitude and Maximum Amplitudes respectively for future reference and analysis
# if the values returned are infinite (missing), these values are converted to "NA" for ease in omission and table organization down the road
first.amp <- sizes[sizes > thres][1]
max.amp <- max(sizes[sizes > thres])
if(!is.infinite(max.amp)){max.amp = max.amp}else{max.amp = NA}
# Ensure row names are numeric values
x = as.numeric(sub("^X","",rownames(data)))
# Set up an easily accessible table to use when calculating AUC (see xy.table below)
xy.table<-data.frame(x,y)
# Define the local minimum value x prior to first peak (per each cell; if non-existent = NA)
# this value stands for a specific time, and can also be used for raw data AUC calculation
x_initial<-as.numeric(first.pre.troughX)
# Define the y value on Normalized data which corresponds to this X(time)_initial
y_initial<-xy.table$y[xy.table$x==first.pre.troughX]
# Transform by subtraction: y = y-y_initial (i.e. initial y value becomes 0, subsequent values are y - y_initial)
y_sub<-y-y_initial
## If curve becomes negative, change value to 0
# (disclude from AUC integration' as alternate solution to adding negatives or adding absolute value of negatives)
y_sub_zero<-y_sub
y_sub_zero[y_sub_zero<0]<-0
# integrate newly formed y_sub vectors into xy.table
xy.table$y_sub<-y_sub
xy.table$y_sub_zero<-y_sub_zero
# make a base column of 0, helps with aesthetics when shading plots for AUC
xy.table$Base0<-0
# Incorporate column for Raw data y values
xy.table$Raw.y<-Raw.y
# Define the Y Value in Raw data which corresponds to this X_initial
Raw.y_initial<-xy.table$Raw.y[xy.table$x==x_initial]
# Transform Raw data (Raw.y - Raw.y_initial)
Raw.y_sub<-Raw.y - Raw.y_initial
# Change negative transformed Raw Data values to 0
Raw.y_sub_zero<-Raw.y_sub
Raw.y_sub_zero[Raw.y_sub_zero<0]<-0
# incorporate these columns in the xy.table
xy.table$Raw.y_sub<-Raw.y_sub
xy.table$Raw.y_sub_zero<-Raw.y_sub_zero
# subset responsive portion of the data for shading
shade.table<-xy.table[xy.table$x>=x_initial,]
shade.table<-shade.table[shade.table$x<=last.post.troughX,]
# Calculating number of peaks, frequency, AUC, average amplitude, duration, etc
if(any(sizes>thres, na.rm = TRUE)) {
numpeaks <- sum(sizes > thres, na.rm = TRUE)
average.amp<-(sum(sizes[sizes>thres]))/(numpeaks)
Duration<- last.post.troughX - first.pre.troughX
Total_Image_Time<-x[length(x)]
# AUC derivations
Norm.AUC<-auc(x,y_sub, from = x_initial, to = last.post.troughX)
Raw.AUC.with.negs<-auc(x,Raw.y_sub, from = x_initial, to = last.post.troughX)
Raw.AUC<-auc(x,Raw.y_sub_zero, from = x_initial, to = last.post.troughX)
# define frequency (mHz) and duty.cycle (deprecated?) for cells with greater than 1 peak
if(numpeaks>1) {
freq <- (numpeaks/(Duration))*1000
Duty.cycle <- avgpeakdur/Duration
} else {
freq = NA
Duty.cycle = NA
}
} else {
numpeaks = NA
average.amp = NA
Duration = NA
freq = NA
Duty.cycle = NA
Total_Image_Time = NA
Norm.AUC = NA
Raw.AUC.with.negs = NA
Raw.AUC = NA
}
#############################################
# Line Plot of smoothed y values over time(x)
# Resetting ymin and ymax to define plot window boundaries/scales for smoothed/normalized curve
ymin<-min(y)
ymax <-max(y)
# set up plot window bounds/scale
plot(x,y, cex=0.25, col="Gray", xlab = "Time (Seconds)", ylab = "Fold change from baseline", main=paste("Smoothed Curve of Normalized Data" ,"\n", "w = ", w, ", span = ", span, ", thres = ", thres, sep="")
,ylim = c(min(y)-1, max(y)+1))
# plot line for smoothed data
lines(x, peaks$y.hat, lwd=1)
# Show x,y values for local maxima, minima, and peaks (maxima with sizes>thres) on smoothed line plot
points(x[peaks$i], peaks$y.hat[peaks$i], col="Red", pch=19, cex=.5)
points(x[nadirs$i], nadirs$y.hat[nadirs$i], col="Green", pch=19, cex=.5)
points(x[peaks$i[sizes > thres]], nadirs$y.hat[peaks$i[sizes > thres]], col="Gold", pch=19, cex=.5)
# Show a red line from time(x) of first peak(y) to time(x) of final peak(y)
lines(c(first.peakx, last.peakx), c(max(peaks$y.hat)+.3,max(peaks$y.hat)+.3), col="red")
# Show Duration line from first.pre.troughX to last post.troughX:
# (refers to x(ie. time) value of minima before first (or only) peak, and x value of final post-peak minima
lines(c(first.pre.troughX, last.post.troughX ), c(min(peaks$y.hat)-.4,min(peaks$y.hat)-.4), col="Blue")
# Set legend
legend("topright", legend = paste("Total peaks = ", numpeaks, "\n freq = ", round(freq,digits=1), "mHz"), cex = 0.5)
########################
# AUC plot: for Raw Data (Inclusive of Negative values)
# setting ymin and ymax to define plot boundaries/ scale for the AUC plot window
ymin<-min(Raw.y_sub)
ymax <- max(Raw.y_sub)
# Create Area under the Curve plot (with soon to be shaded region)
# This will be the second plot: presented directly below the raw data plot called first in percellwrapper:
# I.e the first plot was generated by calling percellwrapper, which then calls this (test) function)
# making an initial plot window with all desired parameters and text
if(any(sizes>thres, na.rm = TRUE)){
plot(x, Raw.y_sub,
cex=0.25,
col="Gray",
xlab = "Time (Seconds)",
ylab = "Mean Grey Pixel Intensity (Arbitrary Units)",
ylim = c(ymin-1, ymax+1),
main=paste("Area Under the Curve: RAW DATA (With Negatives)",
"\n", "Blue = start ", " Purple = End")
)
# this line curve represents raw data with values for y - y_initial,
# This sets the initial baseline y=0 for t=0, and adjusts each subsequent y(t) accordingly
lines(x,Raw.y_sub, lwd =1)
points(x_initial,0, col="Cyan", pch=19, cex=1)
points(x[length(x)],Raw.y_sub[length(Raw.y_sub)], col = "Purple", pch=19, cex=1)
## Shading of the Area inclusive in AUC (per cell)
polygon(c(shade.table$x, rev(shade.table$x)), c(shade.table$Raw.y_sub, rev(shade.table$Base0)),
col = "gray70", border = NA)
}
# If the cell is unresponsive (no peaks){Skip AUC plot by instituting an empty placeholder slot}
else{plot.new()}
########################
# AUC plot: for RAW Data (Negatives = 0)
# setting ymin and ymax to define plot boundaries/ scale for the AUC plot window
ymin<-min(Raw.y_sub)
ymax<-max(Raw.y_sub)
# Create Area under the Curve plot (with soon to be shaded region)
# This will be the second plot: presented directly below the raw data plot called first in percellwrapper:
# I.e the first plot was generated by calling percellwrapper, which then calls this (test) function)
# making an initial plot window with all desired parameters and text
if(any(sizes>thres, na.rm = TRUE)){
plot(x, Raw.y_sub,
cex=0.25,
col="Gray",
xlab = "Time (Seconds)",
ylab = "Mean Grey Pixel Intensity (Arbitrary Units)",
ylim = c(ymin-1, ymax+1),
main=paste("Area Under the Curve: RAW DATA (Negatives = 0)",
"\n", "Blue = start ", " Purple = End")
)
# this line curve represents raw data with values for y - y_initial,
# This sets the initial baseline y=0 for t=0, and adjusts each subsequent y(t) accordingly
lines(x,Raw.y_sub, lwd =1)
points(x_initial,0, col="Cyan", pch=19, cex=1)
points(x[length(x)],Raw.y_sub[length(Raw.y_sub)], col = "Purple", pch=19, cex=1)
## Shading of the Area inclusive in AUC (per cell)
# change to last.post.troughX
polygon(c(shade.table$x, rev(shade.table$x)), c(shade.table$Raw.y_sub_zero, rev(shade.table$Base0)),
col = "gray70", border = NA)
}
# If the cell is unresponsive (no peaks){Skip AUC plot as stated above}
else{plot.new()}
# Return all relevant parameters to be compiled into a table (one row per each biological cell)
return(list(
first = c(first.peakx),
last = c(last.peakx),
first.pre.troughX = c(first.pre.troughX),
last.post.troughX = c(last.post.troughX),
Duration = c(Duration),
freq = c(freq),
Duty.cycle = c(Duty.cycle),
avg.peak.dur = c(avgpeakdur),
Total_Image_Time=c(Total_Image_Time),
fA = c(first.amp),
mA = c(max.amp),
numpeaks = c(numpeaks),
Average_Amplitude= c(average.amp),
Norm.AUC=c(Norm.AUC),
Raw.AUC.with.negs=c(Raw.AUC.with.negs),
Raw.AUC=c(Raw.AUC)
))
}
# this function allows for processing of all cells via the test function above, using parameters defined here as inputs for smoothing functions and threshold for peak identification
# this function also plots the normalized data prior to smoothing, before calling the test function (i.e. all analysis is iterated per cell per dataset, and 4 plots are generated with summary data for all cells in a dataset)
percellwrapper = function(cell, w = 3.5, span = 0.0145, thres = 0.35) {
Raw.y<<-as.numeric(Raw.data[,cell])
y <<- as.numeric(data[,cell])
x <<- as.numeric(sub("^X","",rownames(data)))
par(mfrow=c(4,1))
plot(x,y, type="l", ylim = c(min(y)-1, max(y)+1), xlab = "Time (Seconds)", ylab = "Fold change from baseline", main=paste("Cell",cell,"\n","\n", "Raw Data Normalized"))
results = test(w, span, thres)
return(results)
}
# Function for AOV standard error and statistical analysis/ comparison
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
## data: a data frame.
## measurevar: the name of a column that contains the variable to be summariezed
## groupvars: a vector containing names of columns that contain grouping variables
## na.rm: a boolean that indicates whether to ignore NA's
## conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
conf.interval=.95, .drop=TRUE) {
library(plyr)
# New version of length which can handle NA's: if na.rm==T, don't count them
length2 <- function (x, na.rm=FALSE) {
if (na.rm) sum(!is.na(x))
else length(x)
}
# This does the summary. For each group's data frame, return a vector with
# N, mean, and sd
datac <- ddply(data, groupvars, .drop=.drop,
.fun = function(xx, col) {
c(N = length2(xx[[col]], na.rm=na.rm),
mean = mean (xx[[col]], na.rm=na.rm),
sd = sd (xx[[col]], na.rm=na.rm)
)
},
measurevar
)
# Rename the "mean" column
datac <- rename(datac, c("mean" = measurevar))
datac$se <- datac$sd / sqrt(datac$N) # Calculate standard error of the mean
# Confidence interval multiplier for standard error
# Calculate t-statistic for confidence interval:
# e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
ciMult <- qt(conf.interval/2 + .5, datac$N-1)
datac$ci <- datac$se * ciMult
return(datac)
}
# Function: load summary data for each
# This function is called subsequent to the processing of all datasets and replicates via percellwrapper and test functions
# The output of the above functions are dataset and replicate specific summary tables, which are loaded and organized by replicate and dose by this function
load_combine <- function(dose_specific_files, experiment_drug_dose,
sheet_names = c("Summary_data","Response_Distributions","Monotonic_Amplitudes")) {
# One Final_summary excel WORKBOOK is generated for every position (field) for a given dose
# Each Final_summary WORKBOOK contains three sheets of summary data
# Load all worksheets for a given dose:
# If number of POSITION files for a given dose = 1:
if(length(dose_specific_files)==1)
{
# Iterates from the first sheet to the final sheet
# Loads the excel file and creates a data.frame object for the table within that sheet
# The data.frame objects are defined in the global environment: object names are derived from file names
for(i in 1:length(sheet_names))
{
# Placeholder which retrieves the current object's file name to be assigned during output
# Format: experiment_drug_dose_SheetName
variableN <- paste(experiment_drug_dose, sheet_names[i],sep="_")
# Begin at sheet 1, and continue to the final sheet (final sheet = length(sheet_names))
# Import the "i"th sheet:
temp_sheet <-read.xlsx(dose_specific_files, sheetIndex = i)
# Assigns output as a variable with relative identification
assign(eval(substitute(variableN)), temp_sheet, envir=.GlobalEnv)
}
rm(temp_sheet)
}
# Combination of summary data for a given dose, per position (imaging field). * Positions are designated first by base 0, and continued as integers (i.e p0,p1,p2,p3)
# if multiple positions (fields) were imaged in a well for a given dose, each will have its own excel summary workbook
# To Process cumulative data per each
if(length(dose_specific_files)>1){
# The code proceeds to combine the "i"th sheet from each position workbook corresponding to a given dose
# and combine them into one table per dose, per experiment - defined as a global environment object, per sheet
for(i in 1:length(sheet_names))
{
variableN <- paste(experiment_drug_dose, sheet_names[i],sep = "_")
# This retrieves sheet names/data for a given dose using a pattern which
data.list <- lapply(dose_specific_files, readWorksheetFromFile, sheet = i)
# Retrieves column names from data.frames
data.names <- lapply(data.list, names)
# Creates a "results data.frame" i.e.= begin an output dataframe for combined p0 & p1 data; starts by first inputing 'position 0 data'
result <- data.frame(Position = 0, data.list[[1]])
# Creates a working variable for new column name: position designation: formats according to data in position 0 dataframe
current.names <- c("Position", data.names[[1]])
# consistency in naming of rows, etc
for(i in seq_along(data.list)[-1]){
new.names <- setdiff(data.names[[i]], current.names)
current.names <- union(current.names, new.names)
for(nm in new.names)
result[[ nm ]] <- NA # Placeholder
temp <- as.data.frame(matrix(nrow = nrow(data.list[[i]]), ncol = ncol(result)))
names(temp) <- current.names
# I added -1 here to make 'sequenced data point#2 = '1';
# the line below ∆'s the position numbers to match base 0 and sequential integers, to dictate a sequence 0,1...n)
temp[[ "Position" ]] <- i-1
temp[ , data.names[[i]] ] <- data.list[[i]]
result <- rbind(result, temp) # results in an appended data.frame with position designations
rm(temp)
assign(eval(substitute(variableN)), result, envir=.GlobalEnv) # assigns iterative result as new variable in global.env
}
rm(result) #initial cleanup
}
}
}
########################
# Process Raw/ Normalized Data #####
# collect all file names for input of normalized data files
CSV_files<-list.files(pattern = "CSV")
# parse out dataset specific information for each of these data files
split_files<-strsplit(CSV_files, split = "_")
chopped_files<-strsplit(CSV_files, split = "_CSV_")
experimental.replicate<-sapply(chopped_files, (function(x) x[1]))
concentrations = sapply(split_files, (function(x) x[4]))
drugs = sapply(split_files, (function(x) x[3]))
experimental.replicate<-unique(experimental.replicate)
### Iterative processing of all normalized data files ###
# aesthetic progress bar for R script
CSV.count<-0
avg.time=0
progress.bar<-tkProgressBar("Smooth, Find Peaks, Summarize, and Plot","Running Normalized data \n Estimating Time Remaining",0,length(CSV_files),CSV.count, width = 420)
Sys.sleep(.5)
t00=as.numeric(proc.time()[3])
# first process by replicate
for(e in 1:length(experimental.replicate)){
cat("Processing Replicate:", experimental.replicate[e],"\n")
pattern.code <-paste(experimental.replicate[e],"CSV",sep="_")
experiment.specific.files <- as.character(list.files(pattern = pattern.code, recursive=T))
# parse out all files associated with each replicate (this accomodates the possibility for differing number of files per each replicate)
for(a in 1:length(experiment.specific.files))
{
t0=as.numeric(proc.time()[3])
cat("Processing Condition:", experiment.specific.files[a],"\n")
if(CSV.count<1){amt.done=paste0(0,"%", "\n", "Estimating Time Remaining")
} else if(CSV.count>0){
rem.time <- (round(avg.time*(length(CSV_files)-CSV.count)))
rem.time.str <- paste0(floor(rem.time/60), ":", rem.time-(floor(rem.time/60)*60))
amt.done<-paste0(round((as.numeric(CSV.count)/as.numeric(length(CSV_files))*100)),"% Time Remaining (MM:SS): ", rem.time.str)
}
setTkProgressBar(progress.bar,CSV.count,
label=paste0("Processing: ",
experiment.specific.files[a],
"\n", amt.done))
chopped_e_files<-strsplit(experiment.specific.files, split = "_CSV_")
file.prefix<-unlist(strsplit(chopped_e_files[[a]][1], split = ".csv"))
file.suffix<-unlist(strsplit(chopped_e_files[[a]][2], split = ".csv"))
file.name<-paste0(file.prefix, "_", file.suffix)
column.key<-unlist(strsplit(file.suffix, split = "_"))
column.key<-paste(column.key[2],column.key[3],sep = ".")
# conversion of characters to accomodate import of columns of an excel file which contains artifact frames to be removed prior to analysis
column.key<-gsub("[~]", ".", column.key)
column.key<-gsub("[-]", ".", column.key)
column.key<-gsub("[;]", ".", column.key)
data <- read.csv(experiment.specific.files[a],header=T,row.names=1)
data = t(data)
head(data)
x = as.numeric(sub("^X","",rownames(data)))
# Import raw (not normalized) data for AUC based calculations and reference further down the line
Raw.data.label<-gsub("CSV","RAW", experiment.specific.files[a])
Raw.data <-read.csv(Raw.data.label,header=T,row.names=1)
Raw.data = t(Raw.data)
head(Raw.data)
Raw.x = as.numeric(sub("^X","",rownames(Raw.data)))
#------------------------------------------------------------------------
# import artifact frame spreadsheet (void frames) to be excluded from analysis
# these artifacts come from improper opening of the shutter, resulting in a blank frame - careful inspection was supervised
void.file<-paste0(experimental.replicate[e],"_","voids",".xlsx")
void.sheet<-readWorksheetFromFile(void.file, sheet = "voids")
current.column.voids<-void.sheet[[column.key]]
current.column.voids<-current.column.voids[!is.na(current.column.voids)]
if(length(current.column.voids<1)){data=data}
# if the column contains anything apart from NA...
if(length(current.column.voids>0)){
# Process each column
for(Q in 1:length(current.column.voids)){
current.void<-current.column.voids[Q]
# if the void is a single frame number, remove it from the dataset
if(is.na(as.numeric(as.character(current.void)))==FALSE)
{
# Coerce data to character such that it can be coerced further (to numeric)
# (For some reason it is necessary to re-define as.character prior to as.numeric because it is actually a 'false character' upon import)
current.void<-as.numeric(as.character(current.void))
data = data[-(current.void), ]
}
# Otherwise (i.e. if the row contains a range of frames to be voided): remove these frames
else{
current.void.split<-unlist(strsplit(current.void, split = ","))
first<-as.numeric(as.character(current.void.split[1]))
last<-as.numeric(as.character(current.void.split[2]))
data = data[-(first:last), ]
}
} # ends the iteration for-loop for processing all rows within a column
} # ends the if statement for processing columns with voids
# Repeat the process of artifact frame removal for columns of raw mean intensity data (not normalized)
if(length(current.column.voids<1)){Raw.data=Raw.data}
# if the column contains anything apart from NA...
if(length(current.column.voids>0)){
# Process each column
for(Q in 1:length(current.column.voids)){
current.void<-current.column.voids[Q]
# if the void is a single frame number, remove it from the Raw.dataset
if(is.na(as.numeric(as.character(current.void)))==FALSE)
{
# Coerce Raw.data to character such that it can be coerced further (to numeric)
# (For some reason it is necessary to re-define as.character prior to as.numeric because it is actually a 'false character' upon import)
current.void<-as.numeric(as.character(current.void))
Raw.data = Raw.data[-(current.void), ]
}
# Otherwise (i.e. if the row contains a range of frames to be voided): remove these frames
else{
current.void.split<-unlist(strsplit(current.void, split = ","))
first<-as.numeric(as.character(current.void.split[1]))
last<-as.numeric(as.character(current.void.split[2]))
Raw.data = Raw.data[-(first:last), ]
}
} # ends the iteration for-loop for processing all rows within a column
} # ends the if statement for processing columns with voids
#-----------------------------------------------------------------------
### Begin processing and summarization of the current dataset iteration ###
# Begin encoding a pdf file which will contain the plots of each cell for the current dataset being iterated
pdf(paste0("plotted_results_of_",file.name,".pdf"))
# store the results of percellwrapper and test functions within a placeholder variable 'r'
r = sapply((1:dim(data)[2]), function(cell) {
percellwrapper(cell, w=3.5, span=.0145, thres = 0.35)
}, simplify="array")
# Completes PDF encoding of the file containing single cell plots for the current dataset (saves the file to disk)
dev.off()
# r is converted to a matrix with the summary results from each analysis
r = data.frame(t(as.matrix(r)))
head(r)
as.matrix(r)
# This generates an excel table for summary data, and eases the need to manually coerce data frame:
# i.e. This acts as a placeholder and assists in table organization by coercement upon immediate re-import
write.csv(as.matrix(r), "csv_raw_summary_of_.csv")
r = read.csv("csv_raw_summary_of_.csv")
# Calculate Response distributions etc
n = as.numeric(dim(data)[2])
n.mon <- as.numeric(sum(r$numpeaks<2,na.rm=T))
percent.Monotonic <- round((n.mon/ n *100), digits = 1)
n.osc <- as.numeric(sum(r$numpeaks>1,na.rm=T))
percent.Oscillatory <- round((n.osc/ n *100), digits = 1)
n.tot <- as.numeric(sum(r$numpeaks>0,na.rm=T))
percent.Total <- round((n.tot/ n *100), digits = 1)
n._no_resp <-as.numeric(n - n.tot)
percent._no_response <- round((n._no_resp/n *100), digits = 1)
# Retrieves maximum ampltidues from 'R matrix': Organizes max amplitude for all cells with >0 peaks
# Cell identification label is conserved here as well
cell.amps = r[!is.na(r$mA),]
cell.freqs = cell.amps[!is.na(cell.amps$freq),]
cell.freqs$Cell<-cell.freqs$X
cell.freqs$Frequency<-cell.freqs$freq
cell.freqs<-subset(cell.freqs, select = c(Cell,Frequency))
cell.amps$Cell<-cell.amps$X
cell.amps$Max_Amplitude<-cell.amps$mA
cell.amps<-subset(cell.amps, select = c(Cell, Max_Amplitude))
# Combined readout of amplitude & frequency data for all cells with peaks>0 (all responding cells)
# By providing all = true, NA's are transferred for subsequent derivation of monotonics
combo.test <- merge(cell.freqs,cell.amps, all= TRUE)
# Details for ONLY monotonics
monotonics <- combo.test[is.na(combo.test$Frequency),]
monotonics <- subset(monotonics, select = -c(Frequency))
# this is an add-on to ensure that the table is not empty if no monotonic cells are found within a dataset (lower doses mostly)
if(is.na(monotonics$Cell[1])){
monotonics<-monotonics[nrow(monotonics)+1, ]
monotonics$Cell<-"NA"
monotonics$Max_Amplitude<-"NA"
summary(monotonics$Max_Amplitude)
}
# Wrap up for quick reference while testing datasets
collective_summary <- as.data.frame(list(
c("Total n", "Correct Responses", "Monotonic", "Oscillatory","Unresponsive"),
c(n, n.tot, n.mon, n.osc, n._no_resp),
c(100, percent.Total, percent.Monotonic, percent.Oscillatory, percent._no_response)
))
# removes all non responsive cells; Additionally, rows with frequency values = to NA are preserved iff they have >0 peaks
r.table <- r[complete.cases(r[,(1:6)]),]
# Designate an output file name for the current datasets summary tables
previous.workbook<-paste0("Final_Summary_of_",file.name, ".xlsx")
# if a file exists in the current working directory of the same name, remove it first priror to writing it
# this speeds up processing when overwriting several summary sheets per each dataset
if(file.exists(previous.workbook))
{file.remove(previous.workbook)}
# Generates final cleaned up version RAW_summary_CSV: only includes responsive cells: monotonics have "freq = NA"
loadWorkbook(paste0("Final_Summary_of_",file.name, ".xlsx"), create = TRUE)
writeWorksheetToFile(paste0("Final_Summary_of_",file.name, ".xlsx"), data=r.table, sheet="Summary_data", startRow = 1, startCol = 1, header = TRUE)
# Generate a data.frame inclusive of final breakdown of response distributions, ready for export with correct titles
colnames(collective_summary) <- c(" Type"," # of cells"," % of population")
# Write additional sheets to the current final summary excel file for both response distributions and monotonic amplitudes for comparison
writeWorksheetToFile(paste0("Final_Summary_of_",file.name, ".xlsx"), data=collective_summary, sheet="Response_Distributions", startRow = 1, startCol = 1, header = TRUE)
writeWorksheetToFile(paste0("Final_Summary_of_",file.name, ".xlsx"), data=monotonics, sheet="Monotonic_Amplitudes", startRow = 1, startCol = 1, header = TRUE)
collective_summary
CSV.count<-CSV.count+1
avg.time <- (avg.time + (as.numeric(proc.time()[3])-t0))/CSV.count
}
}
# Add-on for aesthetic progress bar
total.runtime<-round(as.numeric(proc.time()[3])-t00, digits = 1)
total.runtime.str <- paste0(floor(total.runtime/60), ":", total.runtime-(floor(total.runtime/60)*60))
setTkProgressBar(progress.bar,CSV.count,label=paste0("Finished","\n","Run Time (MM:SS): ",total.runtime.str))
close(progress.bar)
# clean up WD folder by deleting the placeholder raw summary csv remaining after the final iteration
if(file.exists("csv_raw_summary_of_.csv")==TRUE){
file.remove("csv_raw_summary_of_.csv")}
##################################
# Combine positions and organize data for Statistical Analysis ####
# Once again, Remove all objects from global environment, but keep all functions
rm(list = setdiff(ls(), lsf.str()))
# download and install appropriate version of Java to match R
# i.e. 32-bit or 64-bit
# find the folder on your PC where Java JRE is located (deprecated?)
options(java.home="C:\\Program Files\\Java\\jre1.8.0_91\\")
options(java.parameters = "-Xmx4g")
# make a date variable for attachment to file names
full.date<-gsub(" ","-", as.character(date()))
full.date<-strsplit(full.date, split = "-")
current.date<-paste(full.date[[1]][2],full.date[[1]][3],full.date[[1]][5], sep = ".")
rm(full.date)
# define prefix for pattern based retrieval of all summary files to be imported
pre<-"Final_Summary_of"
# parse out parameters to act as inputs for all summary files identified above
all_files = list.files(pattern=pre, recursive=T)
exp = sapply(strsplit(all_files,"_"), function(x) x[4])
drug = sapply(strsplit(all_files,"_"), function(x) x[5])
conc = sapply(strsplit(all_files,"_"), function(x) x[6])
# start progress bar
Final.Summary.count<-0
Total_Output_Files<-as.numeric(length(unique(exp))*length(unique(drug))*length(unique(conc)))
t00=as.numeric(proc.time()[3])
progress.bar<-tkProgressBar("Combining Experimental Replicate Data by Dose","Preparing to Organize Replicate Data By Dose",0,Total_Output_Files,Final.Summary.count, width = 420)
Sys.sleep(.5)
# First loop through Experiments
for(e in unique(exp)){
cat("Processing experiment:",e,"\n")
# Next loop through Drug for that Experiment
for(d in unique(drug[exp == e])) {
cat("Processing drug:",d,"\n")
# Next loop through Concentrations for that Experiment and Drug
for(c in unique(conc[exp == e & drug == d])) {
cat("Processing drug dose:",c,"\n")
if(Final.Summary.count<1){amt.done=paste0(0,"%")
} else if(Final.Summary.count>0){amt.done<-paste0(round((as.numeric(Final.Summary.count)/as.numeric(Total_Output_Files)*100)),"%")
}
# update progress bar
setTkProgressBar(progress.bar,Final.Summary.count,
label=paste0("Combining Data Within ",
e," (",d,"): ",c, "\n", amt.done))
# dose_specific_files is used to call a list of positions for each dose called at the bottom level
pattern.code <-paste(pre,e,d,c,"p\\d*.xlsx",sep="_")
dose_specific_files <- as.character(list.files(pattern = pattern.code, recursive=T))
tryCatch(
# run function load combine
load_combine(dose_specific_files, paste(e,d,c,sep="_")), error=function(e) print(e))
Final.Summary.count<-Final.Summary.count+1
}
}
}
# end progress bar
total.runtime<-round(as.numeric(proc.time()[3])-t00, digits = 1)
total.runtime.str <- paste0(floor(total.runtime/60), ":", total.runtime-(floor(total.runtime/60)*60))
setTkProgressBar(progress.bar,Final.Summary.count,label=paste0("Finished","\n","Run Time (MM:SS): ",total.runtime.str))
Sys.sleep(.75)
close(progress.bar)
# organize relevant parameters to see inclusive categories per each in the environment (for reference and future processing)
Drugs<-unique(drug)
Experiments<-unique(exp)
Concentrations<-unique(conc)
# Cleanup of environment
rm(pre)
rm(pattern.code)
rm(drug)
rm(exp)
rm(conc)
rm(d)
rm(e)
rm(c)
# hard coded vector for iteration through sheets of the above combined and imported summary files
sheet.names = c("Summary_data","Response_Distributions")
# top of the hierarchy: experiment(i.e. replicate: run through all conditions/doses for replicate 1, then proceed to 2, 3, etc)
for(e in 1:length(Experiments)){
# defines all existing objects SPECIFIC TO ANY GIVEN EXPERIMENTAL REPLICATE[e]
exp.conc.sort <- objects(pattern = Experiments[e])
exp.conc.sort<-strsplit(exp.conc.sort,split="_")
# finds only the concentrations AVAILABLE FOR THAT EXPERIMENT
exp.conc = unique(sapply(exp.conc.sort, function(x) x[3]))
rm(exp.conc.sort)
# Below find the core loop which runs through 'Concentrations per experimental replicate:
# at this point, the 'Concentrations' vectore contains information for condition;dose 1,2,3, etc, and will be split later on.
for(c in 1:length(exp.conc))
{
for(s in 1:length(sheet.names))
{
if(sheet.names[s]=="Summary_data"){
if(e==1&c==1)
{
easy.name<-paste(Experiments[e],"oxo",exp.conc[c],"Summary_data", sep = "_")
new.name<-get(easy.name)
new.name$Group.Dose<-exp.conc[c]
new.name$Replicate<-unlist(strsplit(Experiments[e],split = ".Rep"))[2]
assign(eval(substitute(paste(easy.name))), new.name, envir=.GlobalEnv)
ANOVA.table<-new.name
}
else{
easy.name<-paste(Experiments[e],"oxo",exp.conc[c],"Summary_data", sep = "_")
new.name<-get(easy.name)
new.name$Group.Dose<-exp.conc[c]
new.name$Replicate<-unlist(strsplit(Experiments[e],split = ".Rep"))[2]
assign(eval(substitute(paste(easy.name))), new.name, envir=.GlobalEnv)
ANOVA.table<-rbind(ANOVA.table, new.name)
rm(easy.name)
rm(new.name)
} # Closes else condition for all concentrations other than e==1, c==1
} # Closes the "if sheet.names=Summary_Data" section
if(sheet.names[s]=="Response_Distributions"){
# if the current iteration = experiment 1, concentration 1:
# Creates a template table for responsive cells and total n:
if(e==1&c==1)
{
easy.name<-paste(Experiments[e],"oxo",exp.conc[c],"Response_Distributions", sep = "_")
new.name<-get(easy.name)
names(new.name)[2]<-"Type"
names(new.name)[3]<-"Number of Cells"
names(new.name)[4]<-"Percent of Population"
first<-new.name[new.name$Type=="Correct Responses",]
second<-new.name[new.name$Type=="Total n",]
new.name<-rbind(first,second)
new.name$Replicate<-unlist(strsplit(Experiments[e],split = ".Rep"))[2]
new.name$Group.Dose<-exp.conc[c]
new.name<-new.name[order(new.name$Position), ]
assign(eval(substitute(paste(easy.name))), new.name, envir=.GlobalEnv)
Response.table<-new.name
rm(first)
rm(second)
rm(easy.name)
rm(new.name)
}
# for each iterative concentration within experiment 1 and for all concentrations in subsequent experiments:
# create a new table of data for response distributions for the current iteration and concatenate to the first table
# each new iteration will concatenate the new table to the growing conglomerate
else{
easy.name<-paste(Experiments[e],"oxo",exp.conc[c],"Response_Distributions", sep = "_")
new.name<-get(easy.name)
names(new.name)[2]<-"Type"
names(new.name)[3]<-"Number of Cells"
names(new.name)[4]<-"Percent of Population"
first<-new.name[new.name$Type=="Correct Responses",]
second<-new.name[new.name$Type=="Total n",]
new.name<-rbind(first,second)
new.name$Replicate<-unlist(strsplit(Experiments[e],split = ".Rep"))[2]
new.name$Group.Dose<-exp.conc[c]
new.name<-new.name[order(new.name$Position), ]
assign(eval(substitute(paste(easy.name))), new.name, envir=.GlobalEnv)
Response.table<-rbind(Response.table, new.name)
rm(first)
rm(second)
rm(easy.name)
rm(new.name)
}
} # Closes "if sheet.names = "response distributions" section
} # closes the "sheet.names" loop
} # Closes the "Concentrations" loop
} # Closes the "Experiments" loop
# Cleanup
rm(e)
rm(c)
rm(s)
# Additional step: Separate the single "Group;Dose" column into separate Columns to be named "Condition" and "Dose"
# Placeholder variable J defines the "Group.Dose" Column in the conglomerate ANOVA.table as character vectors
J<-ANOVA.table$Group.Dose
# Split Group.Dose into a list containing 2 character vectors each;
# each [[list]] contains 2 character vectors: [1] = Group, [2] = Drug
Group.Dose.A<-strsplit(J, split=";")
rm(J)
# Defines the empty template table for Group and Dose outisde of the for loop
# This aids in adding data to these tables within the for loop
Group<-matrix(nrow=0, ncol=1)
Dose<-matrix(nrow=0, ncol=1)
# Generates iterative data and combines respsective Group and Dose data into separate single column tables
for(i in 1:length(Group.Dose.A)){
Group<-rbind(Group, c(Group.Dose.A[[i]][1]))
Dose<-rbind(Dose, c(Group.Dose.A[[i]][2]))
}
Group.A<-Group
Dose.A<-Dose
# cleanup
rm(i)
rm(Group)
rm(Dose)
# repeat this process for response table so that condition and dose vectors match the length of this table
# then incorporate thes as new columns as was completed above
K <-Response.table$Group.Dose
Group.Dose.B<-strsplit(K,split = ";")
rm(K)
Group<-matrix(nrow=0, ncol=1)
Dose<-matrix(nrow=0, ncol=1)
for(i in 1:length(Group.Dose.B)){
Group<-rbind(Group, c(Group.Dose.B[[i]][1]))
Dose<-rbind(Dose, c(Group.Dose.B[[i]][2]))
}
Group.B<-Group
Dose.B<-Dose
rm(Group)
rm(Dose)
rm(i)
# Column naming
colnames(Group.A)<-c("Condition")
colnames(Group.B)<-c("Condition")
colnames(Dose.A)<-c("Dose")
colnames(Dose.B)<-c("Dose")
# Remove the single column containing "Groups;Dose" in both ANOVA.table and Response.Table
ANOVA.table<-subset(ANOVA.table, select = -c(Group.Dose))
Response.table<-subset(Response.table, select = -c(Group.Dose))
rm(Group.Dose.A)
rm(Group.Dose.B)
# Insert new columns of "Condition" (i.e. previously was 'Group') and "Dose" into ANOVA.table and Response.table with appropriate vectors respectively
ANOVA.table$Condition<-Group.A
ANOVA.table$Dose<-Dose.A
Response.table$Condition<-Group.B
Response.table$Dose<-Dose.B
# Aesthetic Organization: Re-Order Both Tables by Condition
ANOVA.table<-ANOVA.table[order(ANOVA.table$Condition), ]
Response.table<-Response.table[order(Response.table$Condition),]
# Tidy up useful variables for later use
Dose<-unique(Dose.A)
Conditions <-unique(Group.A)
rm(Dose.A)
rm(Dose.B)
rm(Group.A)
rm(Group.B)
# Setup ANOVA.table with ordered factors for Dose (Molar) and numeric dose
# Define numeric parameters for converting text-based Dose to numeric Molar equivalent
fM <- .000000000000001
pM <- .000000000001
nM <- .000000001
µM <- .000001 ; uM = .000001 # Optional acceptable abbreviation for micromolar
mM <- .001
M <- 1
Split.Molar.Doses<-strsplit(Dose,split = "~")
for(x in 1:length(Split.Molar.Doses)){
if(x==1){
Numeric.Molar.Doses<-as.numeric(c((as.numeric(Split.Molar.Doses[[1]][1]))*(get(Split.Molar.Doses[[1]][2]))))
} else if(x>1){
next.Molar.Dose<-as.numeric((as.numeric(Split.Molar.Doses[[x]][1]))*(get(Split.Molar.Doses[[x]][2])))
Numeric.Molar.Doses<-as.numeric(c(Numeric.Molar.Doses,next.Molar.Dose))
rm(next.Molar.Dose)
}
}
rm(Split.Molar.Doses)
rm(x)
# Setup Tables with ordered factor for Dose (Molar) and numeric dose
ANOVA.table$named.dose<-ANOVA.table$Dose
ANOVA.table$Dose <- c(Numeric.Molar.Doses)[as.numeric(factor(ANOVA.table$Dose, levels = c(Dose)))]
Response.table$Molar.Dose<-c(Numeric.Molar.Doses)[as.numeric(factor(Response.table$Dose, levels = c(Dose)))]
# Set a variable for use in graphs later
Molar.Dose <- ANOVA.table$Dose
R.Molar.Dose<-Response.table$Molar.Dose
# Set a variable for graph representation on a log10[drug.dose] Molar scale
log10_Molar.Dose<- log10(Molar.Dose)
log10_R.Molar.Dose<-log10(R.Molar.Dose)
# Incorporate log converted Dose in matching columns of ANOVA.table
ANOVA.table$log10_Molar.Dose<-log10_Molar.Dose
Response.table$log10_Molar.Dose<-log10_R.Molar.Dose
# make a sorted vector of unique Molar Doses for the upcoming graph loops (see below)
Molar.Dose<-sort(unique(Molar.Dose))
R.Molar.Dose<-sort(unique(R.Molar.Dose))
# There now exists columns for numeric Dose in a Molar scale (ANOVA.table$Molar.Dose), as well as a column with matching log10(Dose) for graphs
# for easy observation and transfer to prism: create a table with only correct responses and their %
Response.table<-Response.table[order(Response.table$Molar.Dose),]
Correct.table<-Response.table[Response.table$Type=="Correct Responses",]
Totals.table<-Response.table[Response.table$Type=="Total n",]
Totals.table<-Totals.table[order(Totals.table$Replicate), ]
Experiments.Replicate<-Experiments
rm(Experiments)
# VISUALIZATION
# define experiment name for output
if(length(Experiments.Replicate)==1){
Experiment.name<-paste0(Experiments.Replicate[1])
} else if(length(Experiments.Replicate)>1){
Experiment.name<-unlist(strsplit(Experiments.Replicate[1], split = ".Rep"))[1]
} else {
Experiment.name=paste0("ERROR!!")
}
# prepare tables for statistics
ANOVA.table$Condition=factor(ANOVA.table$Condition)
ANOVA.table$Dose=factor(ANOVA.table$Dose)
# Subset ANOVA table to include only cells relevant for frequency (freq≠NA, i.e. >1 peak)
Freq.table <- ANOVA.table[ANOVA.table$numpeaks>1,]
# prepare tables for statistics
Response.table$Condition<-factor(Response.table$Condition)
Response.table$Dose<-factor(Response.table$Dose)
Response.table$Replicate<-factor(Response.table$Replicate)
# New parameters in ANOVA.table ( AUC/peak and AUC/sec (± negatives) ) <- NOW DEPRECATED, not used in subsequent analyses, was previously used to compare output
ANOVA.table$`Raw.AUC.per.sec.with.negs`<-ANOVA.table$Raw.AUC.with.negs/(ANOVA.table$Duration)
ANOVA.table$`Raw.AUC.per.sec`<-ANOVA.table$Raw.AUC/(ANOVA.table$Duration)
ANOVA.table$`Raw.AUC.per.peak.with.negs`<-ANOVA.table$Raw.AUC.with.negs/ANOVA.table$numpeaks
ANOVA.table$`Raw.AUC.per.peak`<-ANOVA.table$Raw.AUC/ANOVA.table$numpeaks
# Designate hard copy file name for Data Table output:
Anova.table.filename<-paste0("_",paste("Complete_Data_Tables_for",Experiment.name,"Experiment",current.date, sep = "_"), ".xlsx")
# If file exists in the working directory from previous run, delete file with same name before generating the new file
# this is protected by date, such that older versions are maintained for reference
if(file.exists(Anova.table.filename))
{file.remove(Anova.table.filename)}
# Create an excel file of relevant data tables summaries
wb<-loadWorkbook(Anova.table.filename, create = TRUE)
final.sheet.data<-c(paste("All Conditions","ANOVA.table",sep = "_"),
paste("Frequency table(>1 peaks)","Freq.table",sep = "_"),
paste("Total cells table","Totals.table", sep = "_"),
paste("Responsive cells table","Correct.table", sep = "_"))
# start progress bar
sheet.count<-0
t00=as.numeric(proc.time()[3])
progress.bar<-tkProgressBar("Processing Final Data into an Excel Workbook","Uploading Data Sheets",0,length(final.sheet.data),sheet.count, width = 400)
# For-loop: write all sheets to excel file as detailed above
for(s in 1:length(final.sheet.data)){
if(sheet.count<1){amt.done=paste0(0,"%")
} else if(sheet.count>0){
amt.done<-paste0(round((as.numeric(sheet.count)/as.numeric(length(final.sheet.data))*100)),"%")}
sheet.data.1<-strsplit(final.sheet.data,split = "_")[[s]][1]
sheet.data.2<-get(strsplit(final.sheet.data,split = "_")[[s]][2])
setTkProgressBar(progress.bar,sheet.count,
label=paste0("Processing Sheet: ",
sheet.data.1,
"\n", amt.done))
createSheet(wb,name = sheet.data.1)
writeWorksheet(wb,sheet.data.2,sheet.data.1)
sheet.count<-sheet.count+1
rm(sheet.data.1)
rm(sheet.data.2)
}
rm(s)
rm(amt.done)
# end progress bar
total.runtime<-round(as.numeric(proc.time()[3])-t00, digits = 1)
total.runtime.str <- paste0(floor(total.runtime/60), ":", total.runtime-(floor(total.runtime/60)*60))
setTkProgressBar(progress.bar,sheet.count,label=paste0("Finished","\n","Run Time (MM:SS): ",total.runtime.str))
Sys.sleep(.75)
close(progress.bar)
rm(progress.bar)
rm(final.sheet.data)
rm(sheet.count)
rm(sheet.names)
rm(total.runtime)
rm(total.runtime.str)
rm(t00)
##################################################################
# Statistics ####
# Collect Ordered Dose Names (mostly for aesthetic labeling of future tables)
Ordered.table<-ANOVA.table[order(ANOVA.table$Dose), ]
Ordered.table$named.dose<-gsub("~","",Ordered.table$named.dose)
Ordered.dose.names<-c(unique(Ordered.table$named.dose))
Ordered.log10_Molar.Dose<-c(unique(Ordered.table$log10_Molar.Dose))
# Factor by replicate prior to statistics
# NOTE: this has already been done for Response table above, but was held off for ANOVA or FREQ tables to allow calculations via column math, which is inhibited by factored columns
ANOVA.table$Replicate<-factor(ANOVA.table$Replicate)
Freq.table$Replicate<-factor(Freq.table$Replicate)
# Defininig placeholder variables for making statistics labels
A<-Conditions[Conditions=="CTRL"]
B<-Conditions[!Conditions=="CTRL"]
# Subsetting tables solely based on Control (for examination and reference, fewer lines of code)
Freq_CTRL.table<-as.data.frame(Freq.table[Freq.table$Condition=="CTRL",])
ANOVA_CTRL.table<-as.data.frame(ANOVA.table[ANOVA.table$Condition=="CTRL",])
# Create Frequency and ANOVA tables (as above), this time for the condition[!Control] - this is non-denominational to the experimental group!
Freq_Condition.table<-as.data.frame(Freq.table[Freq.table$Condition==B,])
condition.Freq.label<-paste0(paste("Freq",B, sep = "_"),".table")
assign(eval(substitute(paste(condition.Freq.label))), Freq_Condition.table, envir=.GlobalEnv)
rm(condition.Freq.label)
ANOVA_Condition.table<-as.data.frame(ANOVA.table[ANOVA.table$Condition==B,])
condition.Anova.label<-paste0(paste("ANOVA",B, sep = "_"),".table")
assign(eval(substitute(paste(condition.Anova.label))), ANOVA_Condition.table, envir=.GlobalEnv)
rm(condition.Anova.label)
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
### In each of the following sub-sections below:
## A two way ANOVA is run pertaining to the parameter compared to Conditions by dose (written to excel file)
## A One-Way anova is run per each control and experimental (condition) group per parameter (not written to excel file)
## Data is then re-organized into a table containing TUKEY HSD information pertaining to all conditions
## This data is then parsed (with labels) for each parameter based on significance (p<.05,.01,.001)
## Significance data is exported into the "Complete data tables" excel file which was generated at the end of the previous section
## Significance labels will also be incorporated into the correct locations on the plots below per condition, dose, and parameter!
# # # # # # # # # # # # # # #
# Frequency for cells with >2 peaks
freq.fit = aov(freq ~ Condition * Dose+Replicate, data=Freq.table)
summary(freq.fit)
stat_freq.AOV<-data.frame(summary(freq.fit)[[1]])
stat_freq.AOV$`Statistical Test`<-"Two-Way ANOVA"
sig.freq <- TukeyHSD(freq.fit)
sig.freq<-data.frame(sig.freq$`Condition:Dose`)
for(i in 1:length(Molar.Dose)){
if(i==1) {
dose.tag<-paste(paste0(B,":",Molar.Dose[1]),paste0(A,":",Molar.Dose[1]), sep = "-" )
sig.freq.final<-sig.freq[row.names(sig.freq)==dose.tag,]
sig.freq.final$Molar.Dose<-Molar.Dose[1]
sig.freq.final$named.dose<-Ordered.dose.names[1]
sig.freq.final$`Statistical Test`<-"TukeyHSD"
if((sig.freq.final$p.adj[1]<.001)==TRUE){
sig.freq.tag1<-"p < .001"
sig.freq.tag2<-"***"
sig.freq.final$significance<-sig.freq.tag1
sig.freq.final$`sig.label`<-sig.freq.tag2
} else if(((sig.freq.final$p.adj[1]>.001)==TRUE) & ((sig.freq.final$p.adj[1]<.01)==TRUE)){
sig.freq.tag1<-"p < .01 "
sig.freq.tag2<-"**"
sig.freq.final$significance<-sig.freq.tag1
sig.freq.final$`sig.label`<-sig.freq.tag2
} else if(((sig.freq.final$p.adj[1]>.01)==TRUE) & ((sig.freq.final$p.adj[1]<.05)==TRUE)){
sig.freq.tag1<-"p < .05"
sig.freq.tag2<-"*"
sig.freq.final$significance<-sig.freq.tag1
sig.freq.final$`sig.label`<-sig.freq.tag2
} else if(((sig.freq.final$p.adj[1]>=.05)==TRUE)){
sig.freq.tag1<-"Not Significant"
sig.freq.tag2<-" "
sig.freq.final$significance<-sig.freq.tag1
sig.freq.final$`sig.label`<-sig.freq.tag2
} else {
sig.freq.tag1<-paste("Error in calculating P value")
sig.freq.tag2<-"ERR0R!"
sig.freq.final$significance<-sig.freq.tag1
sig.freq.final$`sig.label`<-sig.freq.tag2
}
rm(sig.freq.tag1)
rm(sig.freq.tag2)
# for every iteration following the first:
} else if(i>1){
dose.tag<-paste(paste0(B,":",Molar.Dose[i]),paste0(A,":",Molar.Dose[i]), sep = "-" )
sig.freq.new<-sig.freq[row.names(sig.freq)==dose.tag,]
sig.freq.new$Molar.Dose<-Molar.Dose[i]
sig.freq.new$named.dose<-Ordered.dose.names[i]
sig.freq.new$`Statistical Test`<-"TukeyHSD"
if((sig.freq.new$p.adj[1]<.001)==TRUE){
sig.freq.tag1<-"p < .001"
sig.freq.tag2<-"***"
sig.freq.new$significance<-sig.freq.tag1
sig.freq.new$`sig.label`<-sig.freq.tag2
} else if(((sig.freq.new$p.adj[1]>.001)==TRUE)& ((sig.freq.new$p.adj[1]<.01)==TRUE)){
sig.freq.tag1<-"p < .01 "
sig.freq.tag2<-"**"
sig.freq.new$significance<-sig.freq.tag1
sig.freq.new$`sig.label`<-sig.freq.tag2
} else if(((sig.freq.new$p.adj[1]>.01)==TRUE)&((sig.freq.new$p.adj[1]<.05)==TRUE)){
sig.freq.tag1<-"p < .05"
sig.freq.tag2<-"*"
sig.freq.new$significance<-sig.freq.tag1
sig.freq.new$`sig.label`<-sig.freq.tag2
} else if((sig.freq.new$p.adj[1]>=.05)==TRUE){
sig.freq.tag1<-"Not Significant"
sig.freq.tag2<-" "
sig.freq.new$significance<-sig.freq.tag1
sig.freq.new$`sig.label`<-sig.freq.tag2
} else {
sig.freq.tag1<-paste("Error in calculating P value")
sig.freq.tag2<-"ERR0R!"
sig.freq.new$significance<-sig.freq.tag1
sig.freq.new$`sig.label`<-sig.freq.tag2
}
sig.freq.final<-rbind(sig.freq.new, sig.freq.final)
rm(sig.freq.new)
rm(sig.freq.tag1)
rm(sig.freq.tag2)
}
}
sig.freq.final<-sig.freq.final[order(sig.freq.final$Molar.Dose),]
stat_freq.TUK.sig_By_Dose<-sig.freq.final
freq.y.sig.labels<-c(stat_freq.TUK.sig_By_Dose$sig.label)
rm(sig.freq.final)
# One Way Anovas - Freq (one per condition)
Freq.CTRL.fit = aov(freq ~ Dose,data=Freq_CTRL.table)
summary(Freq.CTRL.fit)
TukeyHSD(Freq.CTRL.fit)
Freq.CTRL.TUK<-TukeyHSD(Freq.CTRL.fit)
Freq.condition.fit = aov(freq ~ Dose,data=Freq_Condition.table)
summary(Freq.condition.fit)
TukeyHSD(Freq.condition.fit)
Freq.condition.TUK<- TukeyHSD(Freq.condition.fit)
assign(eval(substitute(paste0("Freq.",B,".fit"))), Freq.condition.fit, envir=.GlobalEnv)
assign(eval(substitute(paste0("Freq.",B,".TUK"))), Freq.condition.TUK, envir=.GlobalEnv)
# Write Relevant statistics for the given parameter to the pre-existing 'completed analysis' excel file:
# Data includes 1) Two-way ANOVA [Parameter~Condition*Dose+Replicate]
# 2) TukeyHSD [Parameter vs Condition:Dose, for all matching pairs of doses between conditions]
createSheet(wb,name = "Frequency Stats")
writeWorksheet(wb,data=list((stat_freq.AOV),stat_freq.TUK.sig_By_Dose),"Frequency Stats",
startRow = c(2,13),
startCol = c(1,1),
header = c(TRUE,TRUE),
rownames = "Variable(s) Compared")
# # # # # # # # # # # # # # #
# Average Duration
Duration.fit = aov(Duration ~ Condition *Dose+Replicate, data=ANOVA.table)
summary(Duration.fit)
stat_Duration.AOV<-data.frame(summary(Duration.fit)[[1]])
stat_Duration.AOV$`Statistical Test`<-"Two-Way ANOVA"
sig.Duration<- TukeyHSD(Duration.fit)
sig.Duration<-data.frame(sig.Duration$`Condition:Dose`)
for(i in 1:length(Molar.Dose)){
if(i==1) {
dose.tag<-paste(paste0(B,":",Molar.Dose[1]),paste0(A,":",Molar.Dose[1]), sep = "-" )
sig.Duration.final<-sig.Duration[row.names(sig.Duration)==dose.tag,]
sig.Duration.final$Molar.Dose<-Molar.Dose[1]
sig.Duration.final$named.dose<-Ordered.dose.names[1]
sig.Duration.final$`Statistical Test`<-"TukeyHSD"
if((sig.Duration.final$p.adj[1]<.001)==TRUE){
sig.Duration.tag1<-"p < .001"
sig.Duration.tag2<-"***"
sig.Duration.final$significance<-sig.Duration.tag1
sig.Duration.final$`sig.label`<-sig.Duration.tag2
} else if(((sig.Duration.final$p.adj[1]>.001)==TRUE) & ((sig.Duration.final$p.adj[1]<.01)==TRUE)){
sig.Duration.tag1<-"p < .01 "
sig.Duration.tag2<-"**"
sig.Duration.final$significance<-sig.Duration.tag1
sig.Duration.final$`sig.label`<-sig.Duration.tag2
} else if(((sig.Duration.final$p.adj[1]>.01)==TRUE) & ((sig.Duration.final$p.adj[1]<.05)==TRUE)){
sig.Duration.tag1<-"p < .05"
sig.Duration.tag2<-"*"
sig.Duration.final$significance<-sig.Duration.tag1
sig.Duration.final$`sig.label`<-sig.Duration.tag2
} else if(((sig.Duration.final$p.adj[1]>=.05)==TRUE)){
sig.Duration.tag1<-"Not Significant"
sig.Duration.tag2<-" "
sig.Duration.final$significance<-sig.Duration.tag1
sig.Duration.final$`sig.label`<-sig.Duration.tag2
} else {
sig.Duration.tag1<-paste("Error in calculating P value")
sig.Duration.tag2<-"ERR0R!"
sig.Duration.final$significance<-sig.Duration.tag1
sig.Duration.final$`sig.label`<-sig.Duration.tag2
}
rm(sig.Duration.tag1)
rm(sig.Duration.tag2)
# for every iteration following the first:
} else if(i>1){
dose.tag<-paste(paste0(B,":",Molar.Dose[i]),paste0(A,":",Molar.Dose[i]), sep = "-" )
sig.Duration.new<-sig.Duration[row.names(sig.Duration)==dose.tag,]
sig.Duration.new$Molar.Dose<-Molar.Dose[i]
sig.Duration.new$named.dose<-Ordered.dose.names[i]
sig.Duration.new$`Statistical Test`<-"TukeyHSD"
if((sig.Duration.new$p.adj[1]<.001)==TRUE){
sig.Duration.tag1<-"p < .001"
sig.Duration.tag2<-"***"
sig.Duration.new$significance<-sig.Duration.tag1
sig.Duration.new$`sig.label`<-sig.Duration.tag2
} else if(((sig.Duration.new$p.adj[1]>.001)==TRUE)& ((sig.Duration.new$p.adj[1]<.01)==TRUE)){
sig.Duration.tag1<-"p < .01 "
sig.Duration.tag2<-"**"
sig.Duration.new$significance<-sig.Duration.tag1
sig.Duration.new$`sig.label`<-sig.Duration.tag2
} else if(((sig.Duration.new$p.adj[1]>.01)==TRUE)&((sig.Duration.new$p.adj[1]<.05)==TRUE)){
sig.Duration.tag1<-"p < .05"
sig.Duration.tag2<-"*"
sig.Duration.new$significance<-sig.Duration.tag1
sig.Duration.new$`sig.label`<-sig.Duration.tag2
} else if((sig.Duration.new$p.adj[1]>=.05)==TRUE){
sig.Duration.tag1<-"Not Significant"
sig.Duration.tag2<-" "
sig.Duration.new$significance<-sig.Duration.tag1
sig.Duration.new$`sig.label`<-sig.Duration.tag2
} else {
sig.Duration.tag1<-paste("Error in calculating P value")
sig.Duration.tag2<-"ERR0R!"
sig.Duration.new$significance<-sig.Duration.tag1
sig.Duration.new$`sig.label`<-sig.Duration.tag2
}
sig.Duration.final<-rbind(sig.Duration.new, sig.Duration.final)
rm(sig.Duration.new)
rm(sig.Duration.tag1)
rm(sig.Duration.tag2)
}
}
sig.Duration.final<-sig.Duration.final[order(sig.Duration.final$Molar.Dose),]
stat_Duration.TUK.sig_By_Dose<-sig.Duration.final
Duration.y.sig.labels<-c(stat_Duration.TUK.sig_By_Dose$sig.label)
rm(sig.Duration.final)
# One Way Anovas - Duration (one per condition)
Duration.CTRL.fit = aov(Duration ~ Dose,data=ANOVA_CTRL.table)
summary(Duration.CTRL.fit)
TukeyHSD(Duration.CTRL.fit)
Duration.CTRL.TUK<-TukeyHSD(Duration.CTRL.fit)
Duration.condition.fit = aov(Duration ~ Dose,data=ANOVA_Condition.table)
summary(Duration.condition.fit)
TukeyHSD(Duration.condition.fit)
Duration.condition.TUK<- TukeyHSD(Duration.condition.fit)
assign(eval(substitute(paste0("Duration.",B,".fit"))), Duration.condition.fit, envir=.GlobalEnv)
assign(eval(substitute(paste0("Duration.",B,".TUK"))), Duration.condition.TUK, envir=.GlobalEnv)
# Write Relevant statistics for the given parameter to the pre-existing 'completed analysis' excel file:
# Data includes 1) Two-way ANOVA [Parameter~Condition*Dose+Replicate]
# 2) TukeyHSD [Parameter vs Condition:Dose, for all matching pairs of doses between conditions]
createSheet(wb,name = "Avg Duration stats")
writeWorksheet(wb,data=list((stat_Duration.AOV),stat_Duration.TUK.sig_By_Dose),
sheet="Avg Duration stats",
startRow = c(2,13),
startCol = c(1,1), header = c(TRUE,TRUE), rownames = "Variable(s) Compared")
# # # # # # # # # # # # # # #
# Maximum Amplitude
mA.fit = aov(mA ~ Condition *Dose+Replicate, data=ANOVA.table)
summary(mA.fit)
stat_mA.AOV<-data.frame(summary(mA.fit)[[1]])
stat_mA.AOV$`Statistical Test`<-"Two-Way ANOVA"
sig.mA<- TukeyHSD(mA.fit)
sig.mA<-data.frame(sig.mA$`Condition:Dose`)
for(i in 1:length(Molar.Dose)){
if(i==1) {
dose.tag<-paste(paste0(B,":",Molar.Dose[1]),paste0(A,":",Molar.Dose[1]), sep = "-" )
sig.mA.final<-sig.mA[row.names(sig.mA)==dose.tag,]
sig.mA.final$Molar.Dose<-Molar.Dose[1]
sig.mA.final$named.dose<-Ordered.dose.names[1]
sig.mA.final$`Statistical Test`<-"TukeyHSD"
if((sig.mA.final$p.adj[1]<.001)==TRUE){
sig.mA.tag1<-"p < .001"
sig.mA.tag2<-"***"
sig.mA.final$significance<-sig.mA.tag1
sig.mA.final$`sig.label`<-sig.mA.tag2
} else if(((sig.mA.final$p.adj[1]>.001)==TRUE) & ((sig.mA.final$p.adj[1]<.01)==TRUE)){
sig.mA.tag1<-"p < .01 "
sig.mA.tag2<-"**"
sig.mA.final$significance<-sig.mA.tag1
sig.mA.final$`sig.label`<-sig.mA.tag2
} else if(((sig.mA.final$p.adj[1]>.01)==TRUE) & ((sig.mA.final$p.adj[1]<.05)==TRUE)){
sig.mA.tag1<-"p < .05"
sig.mA.tag2<-"*"
sig.mA.final$significance<-sig.mA.tag1
sig.mA.final$`sig.label`<-sig.mA.tag2
} else if(((sig.mA.final$p.adj[1]>=.05)==TRUE)){
sig.mA.tag1<-"Not Significant"
sig.mA.tag2<-" "
sig.mA.final$significance<-sig.mA.tag1
sig.mA.final$`sig.label`<-sig.mA.tag2
} else {
sig.mA.tag1<-paste("Error in calculating P value")
sig.mA.tag2<-"ERR0R!"
sig.mA.final$significance<-sig.mA.tag1
sig.mA.final$`sig.label`<-sig.mA.tag2
}
rm(sig.mA.tag1)
rm(sig.mA.tag2)
# for every iteration following the first:
} else if(i>1){
dose.tag<-paste(paste0(B,":",Molar.Dose[i]),paste0(A,":",Molar.Dose[i]), sep = "-" )
sig.mA.new<-sig.mA[row.names(sig.mA)==dose.tag,]
sig.mA.new$Molar.Dose<-Molar.Dose[i]
sig.mA.new$named.dose<-Ordered.dose.names[i]
sig.mA.new$`Statistical Test`<-"TukeyHSD"
if((sig.mA.new$p.adj[1]<.001)==TRUE){
sig.mA.tag1<-"p < .001"
sig.mA.tag2<-"***"
sig.mA.new$significance<-sig.mA.tag1
sig.mA.new$`sig.label`<-sig.mA.tag2
} else if(((sig.mA.new$p.adj[1]>.001)==TRUE)& ((sig.mA.new$p.adj[1]<.01)==TRUE)){
sig.mA.tag1<-"p < .01 "
sig.mA.tag2<-"**"
sig.mA.new$significance<-sig.mA.tag1
sig.mA.new$`sig.label`<-sig.mA.tag2
} else if(((sig.mA.new$p.adj[1]>.01)==TRUE)&((sig.mA.new$p.adj[1]<.05)==TRUE)){
sig.mA.tag1<-"p < .05"
sig.mA.tag2<-"*"
sig.mA.new$significance<-sig.mA.tag1
sig.mA.new$`sig.label`<-sig.mA.tag2
} else if((sig.mA.new$p.adj[1]>=.05)==TRUE){
sig.mA.tag1<-"Not Significant"
sig.mA.tag2<-" "
sig.mA.new$significance<-sig.mA.tag1
sig.mA.new$`sig.label`<-sig.mA.tag2
} else {
sig.mA.tag1<-paste("Error in calculating P value")
sig.mA.tag2<-"ERR0R!"
sig.mA.new$significance<-sig.mA.tag1
sig.mA.new$`sig.label`<-sig.mA.tag2
}
sig.mA.final<-rbind(sig.mA.new, sig.mA.final)
rm(sig.mA.new)
rm(sig.mA.tag1)
rm(sig.mA.tag2)
}
}
sig.mA.final<-sig.mA.final[order(sig.mA.final$Molar.Dose),]
stat_mA.TUK.sig_By_Dose<-sig.mA.final
mA.y.sig.labels<-c(stat_mA.TUK.sig_By_Dose$sig.label)
rm(sig.mA.final)
# One Way Anovas - mA (one per condition)
mA.CTRL.fit = aov(mA ~ Dose,data=ANOVA_CTRL.table)
summary(mA.CTRL.fit)
TukeyHSD(mA.CTRL.fit)
mA.CTRL.TUK<-TukeyHSD(mA.CTRL.fit)
mA.condition.fit = aov(mA ~ Dose,data=ANOVA_Condition.table)
summary(mA.condition.fit)
TukeyHSD(mA.condition.fit)
mA.condition.TUK<- TukeyHSD(mA.condition.fit)
assign(eval(substitute(paste0("mA.",B,".fit"))), mA.condition.fit, envir=.GlobalEnv)
assign(eval(substitute(paste0("mA.",B,".TUK"))), mA.condition.TUK, envir=.GlobalEnv)
# Write Relevant statistics for the given parameter to the pre-existing 'completed analysis' excel file:
# Data includes 1) Two-way ANOVA [Parameter~Condition*Dose+Replicate]
# 2) TukeyHSD [Parameter vs Condition:Dose, for all matching pairs of doses between conditions]
createSheet(wb,name = "Max Amplitude Stats")
writeWorksheet(wb, data=list((stat_mA.AOV),stat_mA.TUK.sig_By_Dose),
sheet="Max Amplitude Stats",
startRow = c(2,13),
startCol = c(1,1), header = c(TRUE,TRUE), rownames = "Variable(s) Compared")
# # # # # # # # # # # # # # #
# Average_Amplitude
Average_Amplitude.fit = aov(Average_Amplitude ~ Condition*Dose+Replicate,data=ANOVA.table)
summary(Average_Amplitude.fit)
stat_Average_Amplitude.AOV<-data.frame(summary(Average_Amplitude.fit)[[1]])
stat_Average_Amplitude.AOV$`Statistical Test`<-"Two-Way ANOVA"
sig.Average_Amplitude<- TukeyHSD(Average_Amplitude.fit)
sig.Average_Amplitude<-data.frame(sig.Average_Amplitude$`Condition:Dose`)
for(i in 1:length(Molar.Dose)){
if(i==1) {
dose.tag<-paste(paste0(B,":",Molar.Dose[1]),paste0(A,":",Molar.Dose[1]), sep = "-" )
sig.Average_Amplitude.final<-sig.Average_Amplitude[row.names(sig.Average_Amplitude)==dose.tag,]
sig.Average_Amplitude.final$Molar.Dose<-Molar.Dose[1]
sig.Average_Amplitude.final$named.dose<-Ordered.dose.names[1]
sig.Average_Amplitude.final$`Statistical Test`<-"TukeyHSD"
if((sig.Average_Amplitude.final$p.adj[1]<.001)==TRUE){
sig.Average_Amplitude.tag1<-"p < .001"
sig.Average_Amplitude.tag2<-"***"
sig.Average_Amplitude.final$significance<-sig.Average_Amplitude.tag1
sig.Average_Amplitude.final$`sig.label`<-sig.Average_Amplitude.tag2
} else if(((sig.Average_Amplitude.final$p.adj[1]>.001)==TRUE) & ((sig.Average_Amplitude.final$p.adj[1]<.01)==TRUE)){
sig.Average_Amplitude.tag1<-"p < .01 "
sig.Average_Amplitude.tag2<-"**"
sig.Average_Amplitude.final$significance<-sig.Average_Amplitude.tag1
sig.Average_Amplitude.final$`sig.label`<-sig.Average_Amplitude.tag2
} else if(((sig.Average_Amplitude.final$p.adj[1]>.01)==TRUE) & ((sig.Average_Amplitude.final$p.adj[1]<.05)==TRUE)){
sig.Average_Amplitude.tag1<-"p < .05"
sig.Average_Amplitude.tag2<-"*"
sig.Average_Amplitude.final$significance<-sig.Average_Amplitude.tag1
sig.Average_Amplitude.final$`sig.label`<-sig.Average_Amplitude.tag2
} else if(((sig.Average_Amplitude.final$p.adj[1]>=.05)==TRUE)){
sig.Average_Amplitude.tag1<-"Not Significant"
sig.Average_Amplitude.tag2<-" "
sig.Average_Amplitude.final$significance<-sig.Average_Amplitude.tag1
sig.Average_Amplitude.final$`sig.label`<-sig.Average_Amplitude.tag2
} else {
sig.Average_Amplitude.tag1<-paste("Error in calculating P value")
sig.Average_Amplitude.tag2<-"ERR0R!"
sig.Average_Amplitude.final$significance<-sig.Average_Amplitude.tag1
sig.Average_Amplitude.final$`sig.label`<-sig.Average_Amplitude.tag2
}
rm(sig.Average_Amplitude.tag1)
rm(sig.Average_Amplitude.tag2)
# for every iteration following the first:
} else if(i>1){
dose.tag<-paste(paste0(B,":",Molar.Dose[i]),paste0(A,":",Molar.Dose[i]), sep = "-" )
sig.Average_Amplitude.new<-sig.Average_Amplitude[row.names(sig.Average_Amplitude)==dose.tag,]
sig.Average_Amplitude.new$Molar.Dose<-Molar.Dose[i]
sig.Average_Amplitude.new$named.dose<-Ordered.dose.names[i]
sig.Average_Amplitude.new$`Statistical Test`<-"TukeyHSD"
if((sig.Average_Amplitude.new$p.adj[1]<.001)==TRUE){
sig.Average_Amplitude.tag1<-"p < .001"
sig.Average_Amplitude.tag2<-"***"
sig.Average_Amplitude.new$significance<-sig.Average_Amplitude.tag1
sig.Average_Amplitude.new$`sig.label`<-sig.Average_Amplitude.tag2
} else if(((sig.Average_Amplitude.new$p.adj[1]>.001)==TRUE)& ((sig.Average_Amplitude.new$p.adj[1]<.01)==TRUE)){
sig.Average_Amplitude.tag1<-"p < .01 "
sig.Average_Amplitude.tag2<-"**"
sig.Average_Amplitude.new$significance<-sig.Average_Amplitude.tag1
sig.Average_Amplitude.new$`sig.label`<-sig.Average_Amplitude.tag2
} else if(((sig.Average_Amplitude.new$p.adj[1]>.01)==TRUE)&((sig.Average_Amplitude.new$p.adj[1]<.05)==TRUE)){
sig.Average_Amplitude.tag1<-"p < .05"
sig.Average_Amplitude.tag2<-"*"
sig.Average_Amplitude.new$significance<-sig.Average_Amplitude.tag1
sig.Average_Amplitude.new$`sig.label`<-sig.Average_Amplitude.tag2
} else if((sig.Average_Amplitude.new$p.adj[1]>=.05)==TRUE){
sig.Average_Amplitude.tag1<-"Not Significant"
sig.Average_Amplitude.tag2<-" "
sig.Average_Amplitude.new$significance<-sig.Average_Amplitude.tag1
sig.Average_Amplitude.new$`sig.label`<-sig.Average_Amplitude.tag2
} else {
sig.Average_Amplitude.tag1<-paste("Error in calculating P value")
sig.Average_Amplitude.tag2<-"ERR0R!"
sig.Average_Amplitude.new$significance<-sig.Average_Amplitude.tag1
sig.Average_Amplitude.new$`sig.label`<-sig.Average_Amplitude.tag2
}
sig.Average_Amplitude.final<-rbind(sig.Average_Amplitude.new, sig.Average_Amplitude.final)
rm(sig.Average_Amplitude.new)
rm(sig.Average_Amplitude.tag1)
rm(sig.Average_Amplitude.tag2)
}
}
sig.Average_Amplitude.final<-sig.Average_Amplitude.final[order(sig.Average_Amplitude.final$Molar.Dose),]
stat_Average_Amplitude.TUK.sig_By_Dose<-sig.Average_Amplitude.final
Average_Amplitude.y.sig.labels<-c(stat_Average_Amplitude.TUK.sig_By_Dose$sig.label)
rm(sig.Average_Amplitude.final)
# One Way Anovas - Average Amp (one per condition)
Average_Amplitude.CTRL.fit = aov(Average_Amplitude ~ Dose,data=ANOVA_CTRL.table)
summary(Average_Amplitude.CTRL.fit)
TukeyHSD(Average_Amplitude.CTRL.fit)
Average_Amplitude.CTRL.TUK<-TukeyHSD(Average_Amplitude.CTRL.fit)
Average_Amplitude.condition.fit = aov(Average_Amplitude ~ Dose,data=ANOVA_Condition.table)
summary(Average_Amplitude.condition.fit)
TukeyHSD(Average_Amplitude.condition.fit)
Average_Amplitude.condition.TUK<-TukeyHSD(Average_Amplitude.condition.fit)
assign(eval(substitute(paste0("Average_Amplitude.",B,".fit"))), Average_Amplitude.condition.fit, envir=.GlobalEnv)
assign(eval(substitute(paste0("Average_Amplitude.",B,".TUK"))), Average_Amplitude.condition.TUK, envir=.GlobalEnv)
# Write Relevant statistics for the given parameter to the pre-existing 'completed analysis' excel file:
# Data includes 1) Two-way ANOVA [Parameter~Condition*Dose+Replicate]
# 2) TukeyHSD [Parameter vs Condition:Dose, for all matching pairs of doses between conditions]
createSheet(wb, name = "Average Amplitude Stats")
writeWorksheet(wb,
data=list((stat_Average_Amplitude.AOV),stat_Average_Amplitude.TUK.sig_By_Dose),
sheet="Average Amplitude Stats",
startRow = c(2,13),
startCol = c(1,1), header = c(TRUE,TRUE), rownames = "Variable(s) Compared")
# # # # # # # # # # # # # # #
# numpeaks
numpeaks.fit = aov(numpeaks ~ Condition*Dose+Replicate,data=ANOVA.table)
summary(numpeaks.fit)
stat_numpeaks.AOV<-data.frame(summary(numpeaks.fit)[[1]])
stat_numpeaks.AOV$`Statistical Test`<-"Two-Way ANOVA"
sig.numpeaks<- TukeyHSD(numpeaks.fit)
sig.numpeaks<-data.frame(sig.numpeaks$`Condition:Dose`)
for(i in 1:length(Molar.Dose)){
if(i==1) {
dose.tag<-paste(paste0(B,":",Molar.Dose[1]),paste0(A,":",Molar.Dose[1]), sep = "-" )
sig.numpeaks.final<-sig.numpeaks[row.names(sig.numpeaks)==dose.tag,]
sig.numpeaks.final$Molar.Dose<-Molar.Dose[1]
sig.numpeaks.final$named.dose<-Ordered.dose.names[1]
sig.numpeaks.final$`Statistical Test`<-"TukeyHSD"
if((sig.numpeaks.final$p.adj[1]<.001)==TRUE){
sig.numpeaks.tag1<-"p < .001"
sig.numpeaks.tag2<-"***"
sig.numpeaks.final$significance<-sig.numpeaks.tag1
sig.numpeaks.final$`sig.label`<-sig.numpeaks.tag2
} else if(((sig.numpeaks.final$p.adj[1]>.001)==TRUE) & ((sig.numpeaks.final$p.adj[1]<.01)==TRUE)){
sig.numpeaks.tag1<-"p < .01 "
sig.numpeaks.tag2<-"**"
sig.numpeaks.final$significance<-sig.numpeaks.tag1
sig.numpeaks.final$`sig.label`<-sig.numpeaks.tag2
} else if(((sig.numpeaks.final$p.adj[1]>.01)==TRUE) & ((sig.numpeaks.final$p.adj[1]<.05)==TRUE)){
sig.numpeaks.tag1<-"p < .05"
sig.numpeaks.tag2<-"*"
sig.numpeaks.final$significance<-sig.numpeaks.tag1
sig.numpeaks.final$`sig.label`<-sig.numpeaks.tag2
} else if(((sig.numpeaks.final$p.adj[1]>=.05)==TRUE)){
sig.numpeaks.tag1<-"Not Significant"
sig.numpeaks.tag2<-" "
sig.numpeaks.final$significance<-sig.numpeaks.tag1
sig.numpeaks.final$`sig.label`<-sig.numpeaks.tag2
} else {
sig.numpeaks.tag1<-paste("Error in calculating P value")
sig.numpeaks.tag2<-"ERR0R!"
sig.numpeaks.final$significance<-sig.numpeaks.tag1
sig.numpeaks.final$`sig.label`<-sig.numpeaks.tag2
}
rm(sig.numpeaks.tag1)
rm(sig.numpeaks.tag2)
# for every iteration following the first:
} else if(i>1){
dose.tag<-paste(paste0(B,":",Molar.Dose[i]),paste0(A,":",Molar.Dose[i]), sep = "-" )
sig.numpeaks.new<-sig.numpeaks[row.names(sig.numpeaks)==dose.tag,]
sig.numpeaks.new$Molar.Dose<-Molar.Dose[i]
sig.numpeaks.new$named.dose<-Ordered.dose.names[i]
sig.numpeaks.new$`Statistical Test`<-"TukeyHSD"
if((sig.numpeaks.new$p.adj[1]<.001)==TRUE){
sig.numpeaks.tag1<-"p < .001"
sig.numpeaks.tag2<-"***"
sig.numpeaks.new$significance<-sig.numpeaks.tag1
sig.numpeaks.new$`sig.label`<-sig.numpeaks.tag2
} else if(((sig.numpeaks.new$p.adj[1]>.001)==TRUE)& ((sig.numpeaks.new$p.adj[1]<.01)==TRUE)){
sig.numpeaks.tag1<-"p < .01 "
sig.numpeaks.tag2<-"**"
sig.numpeaks.new$significance<-sig.numpeaks.tag1
sig.numpeaks.new$`sig.label`<-sig.numpeaks.tag2
} else if(((sig.numpeaks.new$p.adj[1]>.01)==TRUE)&((sig.numpeaks.new$p.adj[1]<.05)==TRUE)){
sig.numpeaks.tag1<-"p < .05"
sig.numpeaks.tag2<-"*"
sig.numpeaks.new$significance<-sig.numpeaks.tag1
sig.numpeaks.new$`sig.label`<-sig.numpeaks.tag2
} else if((sig.numpeaks.new$p.adj[1]>=.05)==TRUE){
sig.numpeaks.tag1<-"Not Significant"
sig.numpeaks.tag2<-" "
sig.numpeaks.new$significance<-sig.numpeaks.tag1
sig.numpeaks.new$`sig.label`<-sig.numpeaks.tag2
} else {
sig.numpeaks.tag1<-paste("Error in calculating P value")
sig.numpeaks.tag2<-"ERR0R!"
sig.numpeaks.new$significance<-sig.numpeaks.tag1
sig.numpeaks.new$`sig.label`<-sig.numpeaks.tag2
}
sig.numpeaks.final<-rbind(sig.numpeaks.new, sig.numpeaks.final)
rm(sig.numpeaks.new)
rm(sig.numpeaks.tag1)
rm(sig.numpeaks.tag2)
}
}
sig.numpeaks.final<-sig.numpeaks.final[order(sig.numpeaks.final$Molar.Dose),]
stat_numpeaks.TUK.sig_By_Dose<-sig.numpeaks.final
numpeaks.y.sig.labels<-c(stat_numpeaks.TUK.sig_By_Dose$sig.label)
rm(sig.numpeaks.final)
# One Way Anovas - numpeaks (one per condition)
numpeaks.CTRL.fit = aov(numpeaks ~ Dose,data=ANOVA_CTRL.table)
summary(numpeaks.CTRL.fit)
TukeyHSD(numpeaks.CTRL.fit)
numpeaks.CTRL.TUK<-TukeyHSD(numpeaks.CTRL.fit)
numpeaks.condition.fit = aov(numpeaks ~ Dose,data=ANOVA_Condition.table)
summary(numpeaks.condition.fit)
TukeyHSD(numpeaks.condition.fit)
numpeaks.condition.TUK<-TukeyHSD(numpeaks.condition.fit)
assign(eval(substitute(paste0("numpeaks.",B,".fit"))), numpeaks.condition.fit, envir=.GlobalEnv)
assign(eval(substitute(paste0("numpeaks.",B,".TUK"))), numpeaks.condition.TUK, envir=.GlobalEnv)
# Write Relevant statistics for the given parameter to the pre-existing 'completed analysis' excel file:
# Data includes 1) Two-way ANOVA [Parameter~Condition*Dose+Replicate]
# 2) TukeyHSD [Parameter vs Condition:Dose, for all matching pairs of doses between conditions]
createSheet(wb, name = "Numpeaks Stats")
writeWorksheet(wb, data=list((stat_numpeaks.AOV),stat_numpeaks.TUK.sig_By_Dose),
sheet="Numpeaks Stats",
startRow = c(2,13),
startCol = c(1,1), header = c(TRUE,TRUE), rownames = "Variable(s) Compared")
# # # # # # # # # # # # # # #
# Normalized AUC
Norm.AUC.fit = aov(Norm.AUC ~ Condition*Dose+Replicate,data=ANOVA.table)
summary(Norm.AUC.fit)
stat_Norm.AUC.AOV<-data.frame(summary(Norm.AUC.fit)[[1]])
stat_Norm.AUC.AOV$`Statistical Test`<-"Two-Way ANOVA"
sig.Norm.AUC<- TukeyHSD(Norm.AUC.fit)
sig.Norm.AUC<-data.frame(sig.Norm.AUC$`Condition:Dose`)
for(i in 1:length(Molar.Dose)){
if(i==1) {
dose.tag<-paste(paste0(B,":",Molar.Dose[1]),paste0(A,":",Molar.Dose[1]), sep = "-" )
sig.Norm.AUC.final<-sig.Norm.AUC[row.names(sig.Norm.AUC)==dose.tag,]
sig.Norm.AUC.final$Molar.Dose<-Molar.Dose[1]
sig.Norm.AUC.final$named.dose<-Ordered.dose.names[1]
sig.Norm.AUC.final$`Statistical Test`<-"TukeyHSD"
if((sig.Norm.AUC.final$p.adj[1]<.001)==TRUE){
sig.Norm.AUC.tag1<-"p < .001"
sig.Norm.AUC.tag2<-"***"
sig.Norm.AUC.final$significance<-sig.Norm.AUC.tag1
sig.Norm.AUC.final$`sig.label`<-sig.Norm.AUC.tag2
} else if(((sig.Norm.AUC.final$p.adj[1]>.001)==TRUE) & ((sig.Norm.AUC.final$p.adj[1]<.01)==TRUE)){
sig.Norm.AUC.tag1<-"p < .01 "
sig.Norm.AUC.tag2<-"**"
sig.Norm.AUC.final$significance<-sig.Norm.AUC.tag1
sig.Norm.AUC.final$`sig.label`<-sig.Norm.AUC.tag2
} else if(((sig.Norm.AUC.final$p.adj[1]>.01)==TRUE) & ((sig.Norm.AUC.final$p.adj[1]<.05)==TRUE)){
sig.Norm.AUC.tag1<-"p < .05"
sig.Norm.AUC.tag2<-"*"
sig.Norm.AUC.final$significance<-sig.Norm.AUC.tag1
sig.Norm.AUC.final$`sig.label`<-sig.Norm.AUC.tag2
} else if(((sig.Norm.AUC.final$p.adj[1]>=.05)==TRUE)){
sig.Norm.AUC.tag1<-"Not Significant"
sig.Norm.AUC.tag2<-" "
sig.Norm.AUC.final$significance<-sig.Norm.AUC.tag1
sig.Norm.AUC.final$`sig.label`<-sig.Norm.AUC.tag2
} else {
sig.Norm.AUC.tag1<-paste("Error in calculating P value")
sig.Norm.AUC.tag2<-"ERR0R!"
sig.Norm.AUC.final$significance<-sig.Norm.AUC.tag1
sig.Norm.AUC.final$`sig.label`<-sig.Norm.AUC.tag2
}
rm(sig.Norm.AUC.tag1)
rm(sig.Norm.AUC.tag2)
# for every iteration following the first:
} else if(i>1){
dose.tag<-paste(paste0(B,":",Molar.Dose[i]),paste0(A,":",Molar.Dose[i]), sep = "-" )
sig.Norm.AUC.new<-sig.Norm.AUC[row.names(sig.Norm.AUC)==dose.tag,]
sig.Norm.AUC.new$Molar.Dose<-Molar.Dose[i]
sig.Norm.AUC.new$named.dose<-Ordered.dose.names[i]
sig.Norm.AUC.new$`Statistical Test`<-"TukeyHSD"
if((sig.Norm.AUC.new$p.adj[1]<.001)==TRUE){
sig.Norm.AUC.tag1<-"p < .001"
sig.Norm.AUC.tag2<-"***"
sig.Norm.AUC.new$significance<-sig.Norm.AUC.tag1
sig.Norm.AUC.new$`sig.label`<-sig.Norm.AUC.tag2
} else if(((sig.Norm.AUC.new$p.adj[1]>.001)==TRUE)& ((sig.Norm.AUC.new$p.adj[1]<.01)==TRUE)){
sig.Norm.AUC.tag1<-"p < .01 "
sig.Norm.AUC.tag2<-"**"
sig.Norm.AUC.new$significance<-sig.Norm.AUC.tag1
sig.Norm.AUC.new$`sig.label`<-sig.Norm.AUC.tag2
} else if(((sig.Norm.AUC.new$p.adj[1]>.01)==TRUE)&((sig.Norm.AUC.new$p.adj[1]<.05)==TRUE)){
sig.Norm.AUC.tag1<-"p < .05"
sig.Norm.AUC.tag2<-"*"
sig.Norm.AUC.new$significance<-sig.Norm.AUC.tag1
sig.Norm.AUC.new$`sig.label`<-sig.Norm.AUC.tag2
} else if((sig.Norm.AUC.new$p.adj[1]>=.05)==TRUE){
sig.Norm.AUC.tag1<-"Not Significant"
sig.Norm.AUC.tag2<-" "
sig.Norm.AUC.new$significance<-sig.Norm.AUC.tag1
sig.Norm.AUC.new$`sig.label`<-sig.Norm.AUC.tag2
} else {
sig.Norm.AUC.tag1<-paste("Error in calculating P value")
sig.Norm.AUC.tag2<-"ERR0R!"
sig.Norm.AUC.new$significance<-sig.Norm.AUC.tag1
sig.Norm.AUC.new$`sig.label`<-sig.Norm.AUC.tag2
}
sig.Norm.AUC.final<-rbind(sig.Norm.AUC.new, sig.Norm.AUC.final)
rm(sig.Norm.AUC.new)
rm(sig.Norm.AUC.tag1)
rm(sig.Norm.AUC.tag2)
}
}
sig.Norm.AUC.final<-sig.Norm.AUC.final[order(sig.Norm.AUC.final$Molar.Dose),]
stat_Norm.AUC.TUK.sig_By_Dose<-sig.Norm.AUC.final
Norm.AUC.y.sig.labels<-c(stat_Norm.AUC.TUK.sig_By_Dose$sig.label)
rm(sig.Norm.AUC.final)
# One Way Anovas - AUC (one per condition)
Norm.AUC.CTRL.fit = aov(Norm.AUC ~ Dose,data=ANOVA_CTRL.table)
summary(Norm.AUC.CTRL.fit)
TukeyHSD(Norm.AUC.CTRL.fit)
Norm.AUC.CTRL.TUK<-TukeyHSD(Norm.AUC.CTRL.fit)
Norm.AUC.condition.fit = aov(Norm.AUC ~ Dose,data=ANOVA_Condition.table)
summary(Norm.AUC.condition.fit)
TukeyHSD(Norm.AUC.condition.fit)
Norm.AUC.condition.TUK<-TukeyHSD(Norm.AUC.condition.fit)
assign(eval(substitute(paste0("Norm.AUC.",B,".fit"))), Norm.AUC.condition.fit, envir=.GlobalEnv)
assign(eval(substitute(paste0("Norm.AUC.",B,".TUK"))), Norm.AUC.condition.TUK, envir=.GlobalEnv)
# Write Relevant statistics for the given parameter to the pre-existing 'completed analysis' excel file:
# Data includes 1) Two-way ANOVA [Parameter~Condition*Dose+Replicate]
# 2) TukeyHSD [Parameter vs Condition:Dose, for all matching pairs of doses between conditions]
createSheet(wb, "Norm AUC Stats")
writeWorksheet(wb, data=list((stat_Norm.AUC.AOV),stat_Norm.AUC.TUK.sig_By_Dose),
sheet="Norm AUC Stats",
startRow = c(2,13),
startCol = c(1,1), header = c(TRUE,TRUE), rownames = "Variable(s) Compared")
# # # # # # # # # # # # # # #
# Raw.AUC (With Negs)
Raw.AUC.with.negs.fit = aov(Raw.AUC.with.negs ~ Condition *Dose+Replicate, data=ANOVA.table)
summary(Raw.AUC.with.negs.fit)
stat_Raw.AUC.with.negs.AOV<-data.frame(summary(Raw.AUC.with.negs.fit)[[1]])
stat_Raw.AUC.with.negs.AOV$`Statistical Test`<-"Two-Way ANOVA"
sig.Raw.AUC.with.negs<- TukeyHSD(Raw.AUC.with.negs.fit)
sig.Raw.AUC.with.negs<-data.frame(sig.Raw.AUC.with.negs$`Condition:Dose`)
for(i in 1:length(Molar.Dose)){
if(i==1) {
dose.tag<-paste(paste0(B,":",Molar.Dose[1]),paste0(A,":",Molar.Dose[1]), sep = "-" )
sig.Raw.AUC.with.negs.final<-sig.Raw.AUC.with.negs[row.names(sig.Raw.AUC.with.negs)==dose.tag,]
sig.Raw.AUC.with.negs.final$Molar.Dose<-Molar.Dose[1]
sig.Raw.AUC.with.negs.final$named.dose<-Ordered.dose.names[1]
sig.Raw.AUC.with.negs.final$`Statistical Test`<-"TukeyHSD"
if((sig.Raw.AUC.with.negs.final$p.adj[1]<.001)==TRUE){
sig.Raw.AUC.with.negs.tag1<-"p < .001"
sig.Raw.AUC.with.negs.tag2<-"***"
sig.Raw.AUC.with.negs.final$significance<-sig.Raw.AUC.with.negs.tag1
sig.Raw.AUC.with.negs.final$`sig.label`<-sig.Raw.AUC.with.negs.tag2
} else if(((sig.Raw.AUC.with.negs.final$p.adj[1]>.001)==TRUE) & ((sig.Raw.AUC.with.negs.final$p.adj[1]<.01)==TRUE)){
sig.Raw.AUC.with.negs.tag1<-"p < .01 "
sig.Raw.AUC.with.negs.tag2<-"**"
sig.Raw.AUC.with.negs.final$significance<-sig.Raw.AUC.with.negs.tag1
sig.Raw.AUC.with.negs.final$`sig.label`<-sig.Raw.AUC.with.negs.tag2
} else if(((sig.Raw.AUC.with.negs.final$p.adj[1]>.01)==TRUE) & ((sig.Raw.AUC.with.negs.final$p.adj[1]<.05)==TRUE)){
sig.Raw.AUC.with.negs.tag1<-"p < .05"
sig.Raw.AUC.with.negs.tag2<-"*"
sig.Raw.AUC.with.negs.final$significance<-sig.Raw.AUC.with.negs.tag1
sig.Raw.AUC.with.negs.final$`sig.label`<-sig.Raw.AUC.with.negs.tag2
} else if(((sig.Raw.AUC.with.negs.final$p.adj[1]>=.05)==TRUE)){
sig.Raw.AUC.with.negs.tag1<-"Not Significant"
sig.Raw.AUC.with.negs.tag2<-" "
sig.Raw.AUC.with.negs.final$significance<-sig.Raw.AUC.with.negs.tag1
sig.Raw.AUC.with.negs.final$`sig.label`<-sig.Raw.AUC.with.negs.tag2
} else {
sig.Raw.AUC.with.negs.tag1<-paste("Error in calculating P value")
sig.Raw.AUC.with.negs.tag2<-"ERR0R!"
sig.Raw.AUC.with.negs.final$significance<-sig.Raw.AUC.with.negs.tag1
sig.Raw.AUC.with.negs.final$`sig.label`<-sig.Raw.AUC.with.negs.tag2
}
rm(sig.Raw.AUC.with.negs.tag1)
rm(sig.Raw.AUC.with.negs.tag2)
# for every iteration following the first:
} else if(i>1){
dose.tag<-paste(paste0(B,":",Molar.Dose[i]),paste0(A,":",Molar.Dose[i]), sep = "-" )
sig.Raw.AUC.with.negs.new<-sig.Raw.AUC.with.negs[row.names(sig.Raw.AUC.with.negs)==dose.tag,]
sig.Raw.AUC.with.negs.new$Molar.Dose<-Molar.Dose[i]
sig.Raw.AUC.with.negs.new$named.dose<-Ordered.dose.names[i]
sig.Raw.AUC.with.negs.new$`Statistical Test`<-"TukeyHSD"
if((sig.Raw.AUC.with.negs.new$p.adj[1]<.001)==TRUE){
sig.Raw.AUC.with.negs.tag1<-"p < .001"
sig.Raw.AUC.with.negs.tag2<-"***"
sig.Raw.AUC.with.negs.new$significance<-sig.Raw.AUC.with.negs.tag1
sig.Raw.AUC.with.negs.new$`sig.label`<-sig.Raw.AUC.with.negs.tag2
} else if(((sig.Raw.AUC.with.negs.new$p.adj[1]>.001)==TRUE)& ((sig.Raw.AUC.with.negs.new$p.adj[1]<.01)==TRUE)){
sig.Raw.AUC.with.negs.tag1<-"p < .01 "
sig.Raw.AUC.with.negs.tag2<-"**"
sig.Raw.AUC.with.negs.new$significance<-sig.Raw.AUC.with.negs.tag1
sig.Raw.AUC.with.negs.new$`sig.label`<-sig.Raw.AUC.with.negs.tag2
} else if(((sig.Raw.AUC.with.negs.new$p.adj[1]>.01)==TRUE)&((sig.Raw.AUC.with.negs.new$p.adj[1]<.05)==TRUE)){
sig.Raw.AUC.with.negs.tag1<-"p < .05"
sig.Raw.AUC.with.negs.tag2<-"*"
sig.Raw.AUC.with.negs.new$significance<-sig.Raw.AUC.with.negs.tag1
sig.Raw.AUC.with.negs.new$`sig.label`<-sig.Raw.AUC.with.negs.tag2
} else if((sig.Raw.AUC.with.negs.new$p.adj[1]>=.05)==TRUE){
sig.Raw.AUC.with.negs.tag1<-"Not Significant"
sig.Raw.AUC.with.negs.tag2<-" "
sig.Raw.AUC.with.negs.new$significance<-sig.Raw.AUC.with.negs.tag1
sig.Raw.AUC.with.negs.new$`sig.label`<-sig.Raw.AUC.with.negs.tag2
} else {
sig.Raw.AUC.with.negs.tag1<-paste("Error in calculating P value")
sig.Raw.AUC.with.negs.tag2<-"ERR0R!"
sig.Raw.AUC.with.negs.new$significance<-sig.Raw.AUC.with.negs.tag1
sig.Raw.AUC.with.negs.new$`sig.label`<-sig.Raw.AUC.with.negs.tag2
}
sig.Raw.AUC.with.negs.final<-rbind(sig.Raw.AUC.with.negs.new, sig.Raw.AUC.with.negs.final)
rm(sig.Raw.AUC.with.negs.new)
rm(sig.Raw.AUC.with.negs.tag1)
rm(sig.Raw.AUC.with.negs.tag2)
}
}
sig.Raw.AUC.with.negs.final<-sig.Raw.AUC.with.negs.final[order(sig.Raw.AUC.with.negs.final$Molar.Dose),]
stat_Raw.AUC.with.negs.TUK.sig_By_Dose<-sig.Raw.AUC.with.negs.final
Raw.AUC.with.negs.y.sig.labels<-c(stat_Raw.AUC.with.negs.TUK.sig_By_Dose$sig.label)
rm(sig.Raw.AUC.with.negs.final)
# One Way Anovas - Raw.AUC.with.negs (one per condition)
Raw.AUC.with.negs.CTRL.fit = aov(Raw.AUC.with.negs ~ Dose,data=ANOVA_CTRL.table)
summary(Raw.AUC.with.negs.CTRL.fit)
TukeyHSD(Raw.AUC.with.negs.CTRL.fit)
Raw.AUC.with.negs.CTRL.TUK<-TukeyHSD(Raw.AUC.with.negs.CTRL.fit)
Raw.AUC.with.negs.condition.fit = aov(Raw.AUC.with.negs ~ Dose,data=ANOVA_Condition.table)
summary(Raw.AUC.with.negs.condition.fit)
TukeyHSD(Raw.AUC.with.negs.condition.fit)
Raw.AUC.with.negs.condition.TUK<- TukeyHSD(Raw.AUC.with.negs.condition.fit)
assign(eval(substitute(paste0("Raw.AUC.with.negs.",B,".fit"))), Raw.AUC.with.negs.condition.fit, envir=.GlobalEnv)
assign(eval(substitute(paste0("Raw.AUC.with.negs.",B,".TUK"))), Raw.AUC.with.negs.condition.TUK, envir=.GlobalEnv)
# Write Relevant statistics for the given parameter to the pre-existing 'completed analysis' excel file:
# Data includes 1) Two-way ANOVA [Parameter~Condition*Dose+Replicate]
# 2) TukeyHSD [Parameter vs Condition:Dose, for all matching pairs of doses between conditions]
createSheet(wb, "AUC (Raw, negs) Stats" )
writeWorksheet(wb, data=list((stat_Raw.AUC.with.negs.AOV),stat_Raw.AUC.with.negs.TUK.sig_By_Dose),
sheet="AUC (Raw, negs) Stats",
startRow = c(2,13),
startCol = c(1,1), header = c(TRUE,TRUE), rownames = "Variable(s) Compared")
# # # # # # # # # # # # # # #
# Raw Data AUC (negatives = 0)
Raw.AUC.fit = aov(Raw.AUC ~ Condition*Dose+Replicate,data=ANOVA.table)
summary(Raw.AUC.fit)
stat_Raw.AUC.AOV<-data.frame(summary(Raw.AUC.fit)[[1]])
stat_Raw.AUC.AOV$`Statistical Test`<-"Two-Way ANOVA"
sig.Raw.AUC<- TukeyHSD(Raw.AUC.fit)
sig.Raw.AUC<-data.frame(sig.Raw.AUC$`Condition:Dose`)
for(i in 1:length(Molar.Dose)){
if(i==1) {
dose.tag<-paste(paste0(B,":",Molar.Dose[1]),paste0(A,":",Molar.Dose[1]), sep = "-" )
sig.Raw.AUC.final<-sig.Raw.AUC[row.names(sig.Raw.AUC)==dose.tag,]
sig.Raw.AUC.final$Molar.Dose<-Molar.Dose[1]
sig.Raw.AUC.final$named.dose<-Ordered.dose.names[1]
sig.Raw.AUC.final$`Statistical Test`<-"TukeyHSD"
if((sig.Raw.AUC.final$p.adj[1]<.001)==TRUE){
sig.Raw.AUC.tag1<-"p < .001"
sig.Raw.AUC.tag2<-"***"
sig.Raw.AUC.final$significance<-sig.Raw.AUC.tag1
sig.Raw.AUC.final$`sig.label`<-sig.Raw.AUC.tag2
} else if(((sig.Raw.AUC.final$p.adj[1]>.001)==TRUE) & ((sig.Raw.AUC.final$p.adj[1]<.01)==TRUE)){
sig.Raw.AUC.tag1<-"p < .01 "
sig.Raw.AUC.tag2<-"**"
sig.Raw.AUC.final$significance<-sig.Raw.AUC.tag1
sig.Raw.AUC.final$`sig.label`<-sig.Raw.AUC.tag2
} else if(((sig.Raw.AUC.final$p.adj[1]>.01)==TRUE) & ((sig.Raw.AUC.final$p.adj[1]<.05)==TRUE)){
sig.Raw.AUC.tag1<-"p < .05"
sig.Raw.AUC.tag2<-"*"
sig.Raw.AUC.final$significance<-sig.Raw.AUC.tag1
sig.Raw.AUC.final$`sig.label`<-sig.Raw.AUC.tag2
} else if(((sig.Raw.AUC.final$p.adj[1]>=.05)==TRUE)){
sig.Raw.AUC.tag1<-"Not Significant"
sig.Raw.AUC.tag2<-" "
sig.Raw.AUC.final$significance<-sig.Raw.AUC.tag1
sig.Raw.AUC.final$`sig.label`<-sig.Raw.AUC.tag2
} else {
sig.Raw.AUC.tag1<-paste("Error in calculating P value")
sig.Raw.AUC.tag2<-"ERR0R!"
sig.Raw.AUC.final$significance<-sig.Raw.AUC.tag1
sig.Raw.AUC.final$`sig.label`<-sig.Raw.AUC.tag2
}
rm(sig.Raw.AUC.tag1)
rm(sig.Raw.AUC.tag2)
# for every iteration following the first:
} else if(i>1){
dose.tag<-paste(paste0(B,":",Molar.Dose[i]),paste0(A,":",Molar.Dose[i]), sep = "-" )
sig.Raw.AUC.new<-sig.Raw.AUC[row.names(sig.Raw.AUC)==dose.tag,]
sig.Raw.AUC.new$Molar.Dose<-Molar.Dose[i]
sig.Raw.AUC.new$named.dose<-Ordered.dose.names[i]
sig.Raw.AUC.new$`Statistical Test`<-"TukeyHSD"
if((sig.Raw.AUC.new$p.adj[1]<.001)==TRUE){
sig.Raw.AUC.tag1<-"p < .001"
sig.Raw.AUC.tag2<-"***"
sig.Raw.AUC.new$significance<-sig.Raw.AUC.tag1
sig.Raw.AUC.new$`sig.label`<-sig.Raw.AUC.tag2
} else if(((sig.Raw.AUC.new$p.adj[1]>.001)==TRUE)& ((sig.Raw.AUC.new$p.adj[1]<.01)==TRUE)){
sig.Raw.AUC.tag1<-"p < .01 "
sig.Raw.AUC.tag2<-"**"
sig.Raw.AUC.new$significance<-sig.Raw.AUC.tag1
sig.Raw.AUC.new$`sig.label`<-sig.Raw.AUC.tag2
} else if(((sig.Raw.AUC.new$p.adj[1]>.01)==TRUE)&((sig.Raw.AUC.new$p.adj[1]<.05)==TRUE)){
sig.Raw.AUC.tag1<-"p < .05"
sig.Raw.AUC.tag2<-"*"
sig.Raw.AUC.new$significance<-sig.Raw.AUC.tag1
sig.Raw.AUC.new$`sig.label`<-sig.Raw.AUC.tag2
} else if((sig.Raw.AUC.new$p.adj[1]>=.05)==TRUE){
sig.Raw.AUC.tag1<-"Not Significant"
sig.Raw.AUC.tag2<-" "
sig.Raw.AUC.new$significance<-sig.Raw.AUC.tag1
sig.Raw.AUC.new$`sig.label`<-sig.Raw.AUC.tag2
} else {
sig.Raw.AUC.tag1<-paste("Error in calculating P value")
sig.Raw.AUC.tag2<-"ERR0R!"
sig.Raw.AUC.new$significance<-sig.Raw.AUC.tag1
sig.Raw.AUC.new$`sig.label`<-sig.Raw.AUC.tag2
}
sig.Raw.AUC.final<-rbind(sig.Raw.AUC.new, sig.Raw.AUC.final)
rm(sig.Raw.AUC.new)
rm(sig.Raw.AUC.tag1)
rm(sig.Raw.AUC.tag2)
}
}
sig.Raw.AUC.final<-sig.Raw.AUC.final[order(sig.Raw.AUC.final$Molar.Dose),]
stat_Raw.AUC.TUK.sig_By_Dose<-sig.Raw.AUC.final
Raw.AUC.y.sig.labels<-c(stat_Raw.AUC.TUK.sig_By_Dose$sig.label)
rm(sig.Raw.AUC.final)
# One Way Anovas - Raw.AUC (one per condition)
Raw.AUC.CTRL.fit = aov(Raw.AUC ~ Dose,data=ANOVA_CTRL.table)
summary(Raw.AUC.CTRL.fit)
TukeyHSD(Raw.AUC.CTRL.fit)
Raw.AUC.CTRL.TUK<-TukeyHSD(Raw.AUC.CTRL.fit)
Raw.AUC.condition.fit = aov(Raw.AUC ~ Dose,data=ANOVA_Condition.table)
summary(Raw.AUC.condition.fit)
TukeyHSD(Raw.AUC.condition.fit)
Raw.AUC.condition.TUK<-TukeyHSD(Raw.AUC.condition.fit)
assign(eval(substitute(paste0("Raw.AUC.",B,".fit"))), Raw.AUC.condition.fit, envir=.GlobalEnv)
assign(eval(substitute(paste0("Raw.AUC.",B,".TUK"))), Raw.AUC.condition.TUK, envir=.GlobalEnv)
# Write Relevant statistics for the given parameter to the pre-existing 'completed analysis' excel file:
# Data includes 1) Two-way ANOVA [Parameter~Condition*Dose+Replicate]
# 2) TukeyHSD [Parameter vs Condition:Dose, for all matching pairs of doses between conditions]
createSheet(wb,"AUC (Raw,0=negs)" )
writeWorksheet(wb, data=list((stat_Raw.AUC.AOV),stat_Raw.AUC.TUK.sig_By_Dose),
sheet="AUC (Raw,0=negs)",
startRow = c(2,13),
startCol = c(1,1), header = c(TRUE,TRUE), rownames = "Variable(s) Compared")
# # # # # # # # # # # # # # #
# Average Peak Duration
avg.peak.dur.fit = aov(avg.peak.dur ~ Condition *Dose+Replicate, data=ANOVA.table)
summary(avg.peak.dur.fit)
stat_avg.peak.dur.AOV<-data.frame(summary(avg.peak.dur.fit)[[1]])
stat_avg.peak.dur.AOV$`Statistical Test`<-"Two-Way ANOVA"
sig.avg.peak.dur<- TukeyHSD(avg.peak.dur.fit)
sig.avg.peak.dur<-data.frame(sig.avg.peak.dur$`Condition:Dose`)
for(i in 1:length(Molar.Dose)){
if(i==1) {
dose.tag<-paste(paste0(B,":",Molar.Dose[1]),paste0(A,":",Molar.Dose[1]), sep = "-" )
sig.avg.peak.dur.final<-sig.avg.peak.dur[row.names(sig.avg.peak.dur)==dose.tag,]
sig.avg.peak.dur.final$Molar.Dose<-Molar.Dose[1]
sig.avg.peak.dur.final$named.dose<-Ordered.dose.names[1]
sig.avg.peak.dur.final$`Statistical Test`<-"TukeyHSD"
if((sig.avg.peak.dur.final$p.adj[1]<.001)==TRUE){
sig.avg.peak.dur.tag1<-"p < .001"
sig.avg.peak.dur.tag2<-"***"
sig.avg.peak.dur.final$significance<-sig.avg.peak.dur.tag1
sig.avg.peak.dur.final$`sig.label`<-sig.avg.peak.dur.tag2
} else if(((sig.avg.peak.dur.final$p.adj[1]>.001)==TRUE) & ((sig.avg.peak.dur.final$p.adj[1]<.01)==TRUE)){
sig.avg.peak.dur.tag1<-"p < .01 "
sig.avg.peak.dur.tag2<-"**"
sig.avg.peak.dur.final$significance<-sig.avg.peak.dur.tag1
sig.avg.peak.dur.final$`sig.label`<-sig.avg.peak.dur.tag2
} else if(((sig.avg.peak.dur.final$p.adj[1]>.01)==TRUE) & ((sig.avg.peak.dur.final$p.adj[1]<.05)==TRUE)){
sig.avg.peak.dur.tag1<-"p < .05"
sig.avg.peak.dur.tag2<-"*"
sig.avg.peak.dur.final$significance<-sig.avg.peak.dur.tag1
sig.avg.peak.dur.final$`sig.label`<-sig.avg.peak.dur.tag2
} else if(((sig.avg.peak.dur.final$p.adj[1]>=.05)==TRUE)){
sig.avg.peak.dur.tag1<-"Not Significant"
sig.avg.peak.dur.tag2<-" "
sig.avg.peak.dur.final$significance<-sig.avg.peak.dur.tag1
sig.avg.peak.dur.final$`sig.label`<-sig.avg.peak.dur.tag2
} else {
sig.avg.peak.dur.tag1<-paste("Error in calculating P value")
sig.avg.peak.dur.tag2<-"ERR0R!"
sig.avg.peak.dur.final$significance<-sig.avg.peak.dur.tag1
sig.avg.peak.dur.final$`sig.label`<-sig.avg.peak.dur.tag2
}
rm(sig.avg.peak.dur.tag1)
rm(sig.avg.peak.dur.tag2)
# for every iteration following the first:
} else if(i>1){
dose.tag<-paste(paste0(B,":",Molar.Dose[i]),paste0(A,":",Molar.Dose[i]), sep = "-" )
sig.avg.peak.dur.new<-sig.avg.peak.dur[row.names(sig.avg.peak.dur)==dose.tag,]
sig.avg.peak.dur.new$Molar.Dose<-Molar.Dose[i]
sig.avg.peak.dur.new$named.dose<-Ordered.dose.names[i]
sig.avg.peak.dur.new$`Statistical Test`<-"TukeyHSD"
if((sig.avg.peak.dur.new$p.adj[1]<.001)==TRUE){
sig.avg.peak.dur.tag1<-"p < .001"
sig.avg.peak.dur.tag2<-"***"
sig.avg.peak.dur.new$significance<-sig.avg.peak.dur.tag1
sig.avg.peak.dur.new$`sig.label`<-sig.avg.peak.dur.tag2
} else if(((sig.avg.peak.dur.new$p.adj[1]>.001)==TRUE)& ((sig.avg.peak.dur.new$p.adj[1]<.01)==TRUE)){
sig.avg.peak.dur.tag1<-"p < .01 "
sig.avg.peak.dur.tag2<-"**"
sig.avg.peak.dur.new$significance<-sig.avg.peak.dur.tag1
sig.avg.peak.dur.new$`sig.label`<-sig.avg.peak.dur.tag2
} else if(((sig.avg.peak.dur.new$p.adj[1]>.01)==TRUE)&((sig.avg.peak.dur.new$p.adj[1]<.05)==TRUE)){
sig.avg.peak.dur.tag1<-"p < .05"
sig.avg.peak.dur.tag2<-"*"
sig.avg.peak.dur.new$significance<-sig.avg.peak.dur.tag1
sig.avg.peak.dur.new$`sig.label`<-sig.avg.peak.dur.tag2
} else if((sig.avg.peak.dur.new$p.adj[1]>=.05)==TRUE){
sig.avg.peak.dur.tag1<-"Not Significant"
sig.avg.peak.dur.tag2<-" "
sig.avg.peak.dur.new$significance<-sig.avg.peak.dur.tag1
sig.avg.peak.dur.new$`sig.label`<-sig.avg.peak.dur.tag2
} else {
sig.avg.peak.dur.tag1<-paste("Error in calculating P value")
sig.avg.peak.dur.tag2<-"ERR0R!"
sig.avg.peak.dur.new$significance<-sig.avg.peak.dur.tag1
sig.avg.peak.dur.new$`sig.label`<-sig.avg.peak.dur.tag2
}
sig.avg.peak.dur.final<-rbind(sig.avg.peak.dur.new, sig.avg.peak.dur.final)
rm(sig.avg.peak.dur.new)
rm(sig.avg.peak.dur.tag1)
rm(sig.avg.peak.dur.tag2)
}
}
sig.avg.peak.dur.final<-sig.avg.peak.dur.final[order(sig.avg.peak.dur.final$Molar.Dose),]
stat_avg.peak.dur.TUK.sig_By_Dose<-sig.avg.peak.dur.final
avg.peak.dur.y.sig.labels<-c(stat_avg.peak.dur.TUK.sig_By_Dose$sig.label)
rm(sig.avg.peak.dur.final)
# One Way Anovas - avg.peak.dur (one per condition)
avg.peak.dur.CTRL.fit = aov(avg.peak.dur ~ Dose,data=ANOVA_CTRL.table)
summary(avg.peak.dur.CTRL.fit)
TukeyHSD(avg.peak.dur.CTRL.fit)
avg.peak.dur.CTRL.TUK<-TukeyHSD(avg.peak.dur.CTRL.fit)
avg.peak.dur.condition.fit = aov(avg.peak.dur ~ Dose,data=ANOVA_Condition.table)
summary(avg.peak.dur.condition.fit)
TukeyHSD(avg.peak.dur.condition.fit)
avg.peak.dur.condition.TUK<- TukeyHSD(avg.peak.dur.condition.fit)
assign(eval(substitute(paste0("avg.peak.dur.",B,".fit"))), avg.peak.dur.condition.fit, envir=.GlobalEnv)
assign(eval(substitute(paste0("avg.peak.dur.",B,".TUK"))), avg.peak.dur.condition.TUK, envir=.GlobalEnv)
# Write Relevant statistics for the given parameter to the pre-existing 'completed analysis' excel file:
# Data includes 1) Two-way ANOVA [Parameter~Condition*Dose+Replicate]
# 2) TukeyHSD [Parameter vs Condition:Dose, for all matching pairs of doses between conditions]
createSheet(wb, "avg peak duration stats")
writeWorksheet(wb, data=list((stat_avg.peak.dur.AOV),stat_avg.peak.dur.TUK.sig_By_Dose),
sheet="avg peak duration stats",
startRow = c(2,13),
startCol = c(1,1), header = c(TRUE,TRUE), rownames = "Variable(s) Compared")
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Lastly, SAVE WORKBOOK to disk: without this, no file will be created within the working directory folder
# ALSO NOTE: The workbook has not yet been saved since its creation at the end of the previous 'Combine data section'
saveWorkbook(wb)
rm(wb)
################
# Plots (Please Note: "Statistics" section MUST be run prior to running this section) ####
Dose.name<- c(Ordered.dose.names)
Parameter<-c("Frequencies","Maximum Amplitudes", "Ca2+ Oscillations")
# start progress bar
graph.A.collection = 55
graph.A.count<-0
t00=as.numeric(proc.time()[3])
progress.bar<-tkProgressBar("Encoding Plots into a PDF file","Preparing to Encode Plots into a PDF file",0,graph.A.collection,graph.A.count, width = 400)
Sys.sleep(.5)
# designate an output name for the PDF file which will contain all plots generated below
pdf.name<-paste0("_",Experiment.name, current.date,".pdf")
# Begin encoding the data below into a PDF file: Encoding terminates at "dev.off()"
pdf(pdf.name)
# Loop for generation of Distribution & Density plots for each dose and experimental condition
for(p in 1:length(Parameter)){
for(d in 1: length(Dose.name))
{
if(graph.A.count<1){amt.done=paste0(0,"%")
} else if(graph.A.count>0){
amt.done<-paste0(round((as.numeric(graph.A.count)/as.numeric(graph.A.collection)*100)),"%")}
dense.title<-paste("Geometric Density:", Parameter[p], paste0("(", Dose.name[d], ")"))
distr.title<-paste("Distribution of", Parameter[p],"at", Dose.name[d], "Oxo-M", sep = " ")
setTkProgressBar(progress.bar,graph.A.count,
label=paste0("Generating Frequency Plots for: ",
Dose.name[d],
"\n", amt.done))
if(Parameter[p]=="Frequencies"){
setTkProgressBar(progress.bar,
graph.A.count,
label=paste0("Generating Frequency plots for: ",
Dose.name[d],
"\n", amt.done))
freq.dense <- ggplot(Freq.table[Freq.table$Dose==Molar.Dose[d],], aes(x=freq,fill=Condition)) +
geom_density(alpha=0.3) +
ggtitle(dense.title) +
xlab("Frequency (mHz)") +
scale_color_manual(values=c('green','red')) +
scale_fill_manual(values=c('green','red'))
graph.A.count<-graph.A.count+1
cdat <- ddply(Freq.table[Freq.table$Dose==Molar.Dose[d],], "Condition", summarize, freq.mean=mean(freq))
freq.distr<-ggplot(Freq.table[Freq.table$Dose==Molar.Dose[d],], aes(x=freq,colour=Dose)) +
geom_histogram(binwidth=4,colour='black',fill='white') +
facet_grid(Condition ~ . ) +
ggtitle(distr.title) +
xlab("Frequency (mHz)") +
geom_vline(data=cdat, aes(xintercept=freq.mean),linetype="dashed",size=1,colour='red') +
theme_bw()
print(freq.dense)
print(freq.distr)
graph.A.count<-graph.A.count+1
}
if(Parameter[p]=="Maximum Amplitudes"){
setTkProgressBar(progress.bar,graph.A.count,
label=paste0("Generating Amplitude plots for: ",
Dose.name[d],
"\n", amt.done))
mA.dense <- ggplot(ANOVA.table[ANOVA.table$Dose==Molar.Dose[d],], aes(x=mA,fill=Condition)) +
geom_density(alpha=0.3) +
ggtitle(dense.title) +
xlab("Fold change from baseline (pixel intensity)") +
scale_color_manual(values=c('green','red')) +
scale_fill_manual(values=c('green','red'))
graph.A.count<-graph.A.count+1
cdat <- ddply(ANOVA.table[ANOVA.table$Dose==Molar.Dose[d],], "Condition", summarize, mA.mean=mean(mA))
mA.distr<-ggplot(ANOVA.table[ANOVA.table$Dose==Molar.Dose[d],], aes(x=mA,colour=Dose)) +
geom_histogram(binwidth=.25,colour='black',fill='white') +
facet_grid(Condition ~ . ) +
ggtitle(distr.title) +
xlab("Fold change from baseline (pixel intensity)") +
geom_vline(data=cdat, aes(xintercept=mA.mean),linetype="dashed",size=1,colour='red') +
theme_bw()
print(mA.dense)
print(mA.distr)
graph.A.count<-graph.A.count+1
}
if(p==3){
setTkProgressBar(progress.bar,graph.A.count,
label=paste0("Generating Oscillation plots for: ",
Dose.name[d],
"\n", amt.done))
numpeaks.dense <- ggplot(ANOVA.table[ANOVA.table$Dose==Molar.Dose[d],], aes(x=numpeaks,fill=Condition)) +
geom_density(alpha=0.3) +
ggtitle(dense.title) +
xlab("Peaks") +
scale_color_manual(values=c('green','red')) +
scale_fill_manual(values=c('green','red'))
graph.A.count<-graph.A.count+1
cdat <- ddply(ANOVA.table[ANOVA.table$Dose==Molar.Dose[d],], "Condition", summarize, numpeaks.mean=mean(numpeaks))
numpeaks.distr<-ggplot(ANOVA.table[ANOVA.table$Dose==Molar.Dose[d],], aes(x=numpeaks,colour=Dose)) +
geom_histogram(binwidth=1,colour='black',fill='white') +
facet_grid(Condition ~ .) +
ggtitle(distr.title) +
xlab("Peaks") +
geom_vline(data=cdat, aes(xintercept=numpeaks.mean),linetype="dashed",size=1,colour='red') +
theme_bw()
print(numpeaks.dense)
print(numpeaks.distr)
graph.A.count<-graph.A.count+1
}
}
}
# Cleanup
rm(p)
rm(d)
rm(dense.title)
rm(distr.title)
rm(freq.dense)
rm(mA.dense)
rm(numpeaks.dense)
rm(freq.distr)
rm(mA.distr)
rm(numpeaks.distr)
rm(cdat)
# update for progress bar
amt.done<-paste0(round((as.numeric(graph.A.count)/as.numeric(graph.A.collection)*100)),"%")
# # # # # # # # # # # # # # # #
# Duration
# Cumulative Plot line
# # # # # # # # # # # # # # # #
tgc <- summarySE(ANOVA.table, measurevar="Duration",groupvars=c("Condition","log10_Molar.Dose"))
# Obtain a vector of y values for placement of significance labels, wherein:
# y coordinate = y(data.point.mean) + se(data.point.mean) + ( se(data.point.mean)/10 )
# obtain y coordinate for first data.point attach subsequent y coordinates to vector with each iteration (for all log10_Molar dose values)
for(p in 1:length(Ordered.log10_Molar.Dose)){
if(p==1){
Duration.y.sig.positions = tgc$Duration[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)])/10)}
else if(p>1){
new.Duration.y.sig.positions<-tgc$Duration[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)])/10)
Duration.y.sig.positions<-c(Duration.y.sig.positions,new.Duration.y.sig.positions)
rm(new.Duration.y.sig.positions)
}
}
rm(p)
# significance (star) labels (defined in statistics section, are added here for reference):
# Corresponds to TukeyHSD$Condition:Dose output for experimental-condition-matched-doses
Duration.y.sig.labels<-c(stat_Duration.TUK.sig_By_Dose$sig.label)
# # # # # # # # #
ggplot(tgc, aes(x=log10_Molar.Dose, y=Duration, colour=Condition)) +
annotate("text", x = c(Ordered.log10_Molar.Dose), y = Duration.y.sig.positions, label = Duration.y.sig.labels, size = 7.5) +
geom_errorbar(aes(ymin=Duration-se,ymax=Duration+se), width=.035, colour="black") +
geom_line() +
geom_point(size=3, shape=21,aes(fill=Condition)) +
xlab("log[Oxo] M") + ylab("Duration (sec)") +
ggtitle(paste0("Effect of ", Experiment.name, " on Response Duration")) +
theme_bw() +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# Duration
# Cumulative Scatter/Jitter plot
# # # # # # # # # # # # # # # #
ggplot(ANOVA.table, aes(x=log10_Molar.Dose, y=Duration, colour=Condition)) +
geom_jitter(size=.75, shape=1, width = .15, height = .8) +
xlab("log[Oxo] M") +
ylab("Duration (sec)") +
facet_grid(. ~ Condition) +
ggtitle(paste0("Effect of ", Experiment.name, " on Response Duration")) +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
theme_gray() + theme(panel.margin = unit(2, "lines")) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreeen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# Frequency
# Cumulative Bar Chart
# # # # # # # # # # # # # # # #
setTkProgressBar(progress.bar,graph.A.count,
label=paste0("Generating Cumulative Frequency plots ",
"\n", amt.done))
fdat <- ddply(Freq.table[Freq.table$Dose==Freq.table$Dose,], "Condition", summarize, freq.mean=mean(freq))
ggplot(Freq.table[Freq.table$Dose==Freq.table$Dose,], aes(x=freq,colour=Dose)) +
geom_histogram(binwidth=4,colour='black',fill='white') +
facet_grid(Condition ~ . ) +
ggtitle("Distribution of Ca2+ Response Frequency") +
xlab("Frequency (mHz)") +
geom_vline(data=fdat, aes(xintercept=freq.mean),linetype="dashed",size=1,colour='red') +
theme_bw()
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# Frequency
# Cumulative Density plots
# # # # # # # # # # # # # # # #
ggplot(Freq.table[Freq.table$Dose==Freq.table$Dose,], aes(x=freq,fill=Condition)) +
geom_density(alpha=.6) +
ggtitle("Geometric Density: Frequency") +
xlab("Frequency (mHz)") +
xlim(NA,150) +
scale_color_manual(values=c('green','firebrick2')) +
scale_fill_manual(values=c('green','firebrick2'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# Frequency
# Cumulative Plot line
# # # # # # # # # # # # # # # #
tgc <- summarySE(Freq.table, measurevar="freq",groupvars=c("Condition","log10_Molar.Dose"))
# Obtain a vector of y values for placement of significance labels
# y coordinate = y(data.point.mean) + se(data.point.mean) + ( se(data.point.mean)/10 )
# obtain y coordinate for first data.point attach subsequent y coordinates to vector with each iteration (for all log10_Molar dose values)
for(p in 1:length(Ordered.log10_Molar.Dose)){
if(p==1){
freq.y.sig.positions = tgc$freq[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)])/10)}
else if(p>1){
new.freq.y.sig.positions<-tgc$freq[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)])/10)
freq.y.sig.positions<-c(freq.y.sig.positions,new.freq.y.sig.positions)
rm(new.freq.y.sig.positions)
}
}
rm(p)
# significance (star) labels (defined in statistics section, added here for reference):
# Corresponds to TukeyHSD$Condition:Dose for condition-matched-doses
freq.y.sig.labels<-c(stat_freq.TUK.sig_By_Dose$sig.label)
# # # # # # # # #
ggplot(tgc, aes(x=log10_Molar.Dose, y=freq, colour=Condition)) +
annotate("text", x = c(Ordered.log10_Molar.Dose), y = freq.y.sig.positions, label = freq.y.sig.labels, size = 7.5) +
geom_errorbar(aes(ymin=freq-se,ymax=freq+se), width=.035, colour="black") +
geom_line() +
geom_point(size=3, shape=21,aes(fill=Condition)) +
xlab("log[Oxo] M") +
ylab("Frequency (mHz)") +
ggtitle(paste0("Effect of ", Experiment.name," on Ca2+ frequency (>1 peaks)")) +
theme_bw() +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreen','tomato'))
graph.A.count<-graph.A.count+1
amt.done<-paste0(round((as.numeric(graph.A.count)/as.numeric(graph.A.collection)*100)),"%")
setTkProgressBar(progress.bar,graph.A.count,
label=paste0("Generating Cumulative Amplitude plots ",
"\n", amt.done))
# # # # # # # # # # # # # # # #
# Duty cycle
# Cumulative Plot line
# # # # # # # # # # # # # # # #
tgc <- summarySE(Freq.table, measurevar="Duty.cycle",groupvars=c("Condition","log10_Molar.Dose"))
# Obtain a vector of y values for placement of significance labels
# y coordinate = y(data.point.mean) + se(data.point.mean) + ( se(data.point.mean)/10 )
# obtain y coordinate for first data.point attach subsequent y coordinates to vector with each iteration (for all log10_Molar dose values)
for(p in 1:length(Ordered.log10_Molar.Dose)){
if(p==1){
Duty.cycle.y.sig.positions = tgc$Duty.cycle[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)])/10)}
else if(p>1){
new.Duty.cycle.y.sig.positions<-tgc$Duty.cycle[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)])/10)
Duty.cycle.y.sig.positions<-c(Duty.cycle.y.sig.positions,new.Duty.cycle.y.sig.positions)
rm(new.Duty.cycle.y.sig.positions)
}
}
rm(p)
# significance (star) labels (defined in statistics section, added here for reference):
# Corresponds to TukeyHSD$Condition:Dose output for condition-matched-doses
Duty.cycle.y.sig.labels<-c(stat_Duty.cycle.TUK.sig_By_Dose$sig.label)
# # # # # # # # #
ggplot(tgc, aes(x=log10_Molar.Dose, y=Duty.cycle, colour=Condition)) +
annotate("text", x = c(Ordered.log10_Molar.Dose), y = Duty.cycle.y.sig.positions, label = Duty.cycle.y.sig.labels, size = 7.5) +
geom_errorbar(aes(ymin=Duty.cycle-se,ymax=Duty.cycle+se), width=.035, colour="black") +
geom_line() +
geom_point(size=3, shape=21,aes(fill=Condition)) +
xlab("log[Oxo] M") + ylab("Duty Ratio") +
ggtitle(paste0("Effect of ", Experiment.name, " on Duty Cycle")) +
theme_bw() +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# Duty Cycle
# Cumulative Scatter/Jitter plot
# # # # # # # # # # # # # # # #
ggplot(Freq.table, aes(x=log10_Molar.Dose, y=Duty.cycle, colour=Condition)) +
geom_jitter(size=.75, shape=1, width = .15, height = .8) +
xlab("log[Oxo] M") +
ylab("Duty Ratio") +
facet_grid(. ~ Condition) +
ggtitle(paste0("Effect of ", Experiment.name, " on Duty Cycle")) +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
theme_gray() + theme(panel.margin = unit(2, "lines")) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreeen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# Maxiumum Amplitude
# Cumulative Bar Chart
# # # # # # # # # # # # # # # #
mdat <- ddply(ANOVA.table[ANOVA.table$Dose==ANOVA.table$Dose,], "Condition", summarize, mA.mean=mean(mA))
ggplot(ANOVA.table[ANOVA.table$Dose==ANOVA.table$Dose,], aes(x=mA,colour=Dose)) +
geom_histogram(binwidth=.25,colour='black',fill='white') +
facet_grid(Condition ~ . ) +
xlab("Fold change from baseline (pixel intensity)") +
ggtitle("Distribution of Maximum Response Amplitudes") +
geom_vline(data=mdat, aes(xintercept=mA.mean),linetype="dashed",size=1,colour='red') +
theme_bw()
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# Max Amplitude
# Cumulative Density Distribution
# # # # # # # # # # # # # # # #
ggplot(ANOVA.table[ANOVA.table$Dose==ANOVA.table$Dose,], aes(x=mA,fill=Condition)) +
geom_density(alpha=0.6) +
ggtitle("Geometric Density: Maximum Amplitude") +
xlab("Fold change from baseline (pixel intensity)") +
scale_color_manual(values=c('green','firebrick2')) +
scale_fill_manual(values=c('green','firebrick2'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# Max Amplitude
# Cumulative Plot line
# # # # # # # # # # # # # # # #
tgc <- summarySE(ANOVA.table, measurevar="mA",groupvars=c("Condition","log10_Molar.Dose"))
# Obtain a vector of y values for placement of significance labels
# y coordinate = y(data.point.mean) + se(data.point.mean) + ( se(data.point.mean)/10 )
# obtain y coordinate for first data.point attach subsequent y coordinates to vector with each iteration (for all log10_Molar dose values)
for(p in 1:length(Ordered.log10_Molar.Dose)){
if(p==1){
mA.y.sig.positions = tgc$mA[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)])/10)}
else if(p>1){
new.mA.y.sig.positions<-tgc$mA[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)])/10)
mA.y.sig.positions<-c(mA.y.sig.positions,new.mA.y.sig.positions)
rm(new.mA.y.sig.positions)
}
}
rm(p)
# significance (star) labels (defined in statistics section, added here for reference):
# Corresponds to TukeyHSD$Condition:Dose output for condition-matched-doses
mA.y.sig.labels<-c(stat_mA.TUK.sig_By_Dose$sig.label)
# # # # # # # # #
ggplot(tgc, aes(x=log10_Molar.Dose, y=mA, colour=Condition)) +
annotate("text", x = c(Ordered.log10_Molar.Dose), y = mA.y.sig.positions, label = mA.y.sig.labels, size = 7.5) +
geom_errorbar(aes(ymin=mA-se,ymax=mA+se), width=.035, colour="black") +
geom_line() +
geom_point(size=3, shape=21,aes(fill=Condition)) +
xlab("log[Oxo] M") + ylab("Maximum Amplitude") +
ggtitle(paste0("Effect of ", Experiment.name, " on Ca2+ max amplitude")) +
theme_bw() +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# Max Amplitude
# Cumulative Scatter/Jitter plot
# # # # # # # # # # # # # # # #
ggplot(ANOVA.table, aes(x=log10_Molar.Dose, y=mA, colour=Condition)) +
geom_jitter(size=.75, shape=1, width = .15, height = .8) +
xlab("log[Oxo] M") +
ylab("Maximum Amplitude") +
facet_grid(. ~ Condition) +
ggtitle(paste0("Effect of ", Experiment.name, " on Ca2+ Maximum Amplitude")) +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
theme_gray() + theme(panel.margin = unit(2, "lines")) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreeen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# Average Amplitude
# Cumulative Plot line
# # # # # # # # # # # # # # # #
tgc <- summarySE(ANOVA.table, measurevar="Average_Amplitude",groupvars=c("Condition","log10_Molar.Dose"))
# Obtain a vector of y values for placement of significance labels
# y coordinate = y(data.point.mean) + se(data.point.mean) + ( se(data.point.mean)/10 )
# obtain y coordinate for first data.point attach subsequent y coordinates to vector with each iteration (for all log10_Molar dose values)
for(p in 1:length(Ordered.log10_Molar.Dose)){
if(p==1){
Average_Amplitude.y.sig.positions = tgc$Average_Amplitude[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)])/10)}
else if(p>1){
new.Average_Amplitude.y.sig.positions<-tgc$Average_Amplitude[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)])/10)
Average_Amplitude.y.sig.positions<-c(Average_Amplitude.y.sig.positions,new.Average_Amplitude.y.sig.positions)
rm(new.Average_Amplitude.y.sig.positions)
}
}
rm(p)
# significance (star) labels (defined in statistics section, added here for reference):
# Corresponds to TukeyHSD$Condition:Dose output for condition-matched-doses
Average_Amplitude.y.sig.labels<-c(stat_Average_Amplitude.TUK.sig_By_Dose$sig.label)
# # # # # # # # #
ggplot(tgc, aes(x=log10_Molar.Dose, y=Average_Amplitude, colour=Condition)) +
annotate("text", x = c(Ordered.log10_Molar.Dose), y = Average_Amplitude.y.sig.positions, label = Average_Amplitude.y.sig.labels, size = 7.5) +
geom_errorbar(aes(ymin=Average_Amplitude-se,ymax=Average_Amplitude+se), width=.035, colour="black") +
geom_line() +
geom_point(size=3, shape=21,aes(fill=Condition)) +
xlab("log[Oxo] M") + ylab("Average Amplitude (Fold Change from Baseline") +
ggtitle(paste0("Effect of ", Experiment.name, " on Average Amplitude")) +
theme_bw() +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreen','tomato'))
graph.A.count<-graph.A.count+1
amt.done<-paste0(round((as.numeric(graph.A.count)/as.numeric(graph.A.collection)*100)),"%")
setTkProgressBar(progress.bar,graph.A.count,
label=paste0("Generating Cumulative Oscillation plots ",
"\n", amt.done))
# # # # # # # # # # # # # # # #
# Average Amplitude
# Cumulative Scatter/Jitter plot
# # # # # # # # # # # # # # # #
ggplot(ANOVA.table, aes(x=log10_Molar.Dose, y=Average_Amplitude, colour=Condition)) +
geom_jitter(size=.75, shape=1, width = .15, height = .8) +
xlab("log[Oxo] M") +
ylab("Average Peak Amplitude (AU)") +
facet_grid(. ~ Condition) +
ggtitle(paste0("Effect of ", Experiment.name, " on Ca2+ Average Peak Amplitude")) +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
theme_gray() + theme(panel.margin = unit(2, "lines")) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreeen','tomato'))
# # # # # # # # # # # # # # # #
# numpeaks
# Cumulative Bar Chart
# # # # # # # # # # # # # # # #
pdat <- ddply(ANOVA.table[ANOVA.table$Dose==ANOVA.table$Dose,], "Condition", summarize, numpeaks.mean=mean(numpeaks))
ggplot(ANOVA.table[ANOVA.table$Dose==ANOVA.table$Dose,], aes(x=numpeaks,colour=Dose)) +
geom_histogram(binwidth=1,colour='black',fill='white') +
facet_grid(Condition ~ . ) +
xlab("Peaks") +
ggtitle("Distribution of Total Ca2+ Oscillations") +
geom_vline(data=pdat, aes(xintercept=numpeaks.mean),linetype="dashed",size=1,colour='red') +
theme_bw()
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# numpeaks
# Cumulative Density Distribution
# # # # # # # # # # # # # # # #
ggplot(ANOVA.table[ANOVA.table$Dose==ANOVA.table$Dose,], aes(x=numpeaks,fill=Condition)) +
geom_density(alpha=0.6) +
ggtitle("Geometric Density: Total Ca2+ Oscillations") +
xlab("Peaks") +
scale_color_manual(values=c('green','firebrick2')) +
scale_fill_manual(values=c('green','firebrick2'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# numpeaks
# Cumulative Plot line
# # # # # # # # # # # # # # # #
tgc <- summarySE(ANOVA.table, measurevar="numpeaks",groupvars=c("Condition","log10_Molar.Dose"))
# Obtaining a vector of y values for placement of significance labels
# y coordinate = y(data.point.mean) + se(data.point.mean) + ( se(data.point.mean)/10 )
# obtain y coordinate for first data.point attach subsequent y coordinates to vector with each iteration (for all log10_Molar dose values)
for(p in 1:length(Ordered.log10_Molar.Dose)){
if(p==1){
numpeaks.y.sig.positions = tgc$numpeaks[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)])/10)}
else if(p>1){
new.numpeaks.y.sig.positions<-tgc$numpeaks[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)])/10)
numpeaks.y.sig.positions<-c(numpeaks.y.sig.positions,new.numpeaks.y.sig.positions)
rm(new.numpeaks.y.sig.positions)
}
}
rm(p)
# significance (star) labels (defined in statistics section, added here for reference):
# Correspond to TukeyHSD$Condition:Dose for condition-matched-doses
numpeaks.y.sig.labels<-c(stat_numpeaks.TUK.sig_By_Dose$sig.label)
# # # # # # # # #
ggplot(tgc, aes(x=log10_Molar.Dose, y=numpeaks, colour=Condition)) +
annotate("text", x = c(Ordered.log10_Molar.Dose), y = numpeaks.y.sig.positions, label = numpeaks.y.sig.labels, size = 7.5) +
geom_errorbar(aes(ymin=numpeaks-se,ymax=numpeaks+se), width=.035, colour="black") +
geom_line() +
geom_point(size=3, shape=21,aes(fill=Condition)) +
xlab("log[Oxo] M") + ylab("Peaks") +
ggtitle(paste0("Effect of ", Experiment.name, " on Number of Ca2+ Oscillations")) +
theme_bw() +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # #
# numpeaks
# Cumulative Scatter/Jitter plot
# # # # # # # # # # # # # # # #
ggplot(ANOVA.table, aes(x=log10_Molar.Dose, y=numpeaks, colour=Condition)) +
geom_jitter(size=.75, shape=1, width = .15, height = 1.1) +
xlab("log[Oxo] M") +
ylab("Peaks") +
facet_grid(. ~ Condition) +
ggtitle(paste0("Effect of ", Experiment.name, " on Number of Ca2+ Oscillations")) +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
theme_gray() + theme(panel.margin = unit(2, "lines")) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreeen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # # # # # #
# Normalized Area under the curve (AUC)
# Cumulative Plot line
# # # # # # # # # # # # # # # # # # # #
tgc <- summarySE(ANOVA.table, measurevar="Norm.AUC",groupvars=c("Condition","log10_Molar.Dose"))
# Obtain a vector of y values for placement of significance labels
# y coordinate = y(data.point.mean) + se(data.point.mean) + ( se(data.point.mean)/10 )
# obtain y coordinate for first data.point attach subsequent y coordinates to vector with each iteration (for all log10_Molar dose values)
for(p in 1:length(Ordered.log10_Molar.Dose)){
if(p==1){
Norm.AUC.y.sig.positions = tgc$Norm.AUC[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)])/10)}
else if(p>1){
new.Norm.AUC.y.sig.positions<-tgc$Norm.AUC[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)])/10)
Norm.AUC.y.sig.positions<-c(Norm.AUC.y.sig.positions,new.Norm.AUC.y.sig.positions)
rm(new.Norm.AUC.y.sig.positions)
}
}
rm(p)
# significance (star) labels (defined in statistics section, added here for reference):
# Corresponds to TukeyHSD$Condition:Dose output for condition-matched-doses
Norm.AUC.y.sig.labels<-c(stat_Norm.AUC.TUK.sig_By_Dose$sig.label)
# # # # # # # # #
ggplot(tgc, aes(x=log10_Molar.Dose, y=Norm.AUC, colour=Condition)) +
annotate("text", x = c(Ordered.log10_Molar.Dose), y = Norm.AUC.y.sig.positions, label = Norm.AUC.y.sig.labels, size = 7.5) +
geom_errorbar(aes(ymin=Norm.AUC-se,ymax=Norm.AUC+se), width=.035, colour="black") +
geom_line() +
geom_point(size=3, shape=21,aes(fill=Condition)) +
xlab("log[Oxo] M") + ylab("Normalized Area Under the Curve") +
ggtitle(paste0("Effect of ", Experiment.name, " on Area Under the Curve: Normalized Data")) +
theme_bw() +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # #
# AUC (Raw, WITH NEGATIVES)
# Cumulative Plot line
# # # # # # # # # # # # # #
tgc <- summarySE(ANOVA.table, measurevar="Raw.AUC.with.negs",groupvars=c("Condition","log10_Molar.Dose"))
# Obtain a vector of y values for placement of significance labels
# y coordinate = y(data.point.mean) + se(data.point.mean) + ( se(data.point.mean)/10 )
# obtain y coordinate for first data.point attach subsequent y coordinates to vector with each iteration (for all log10_Molar dose values)
for(p in 1:length(Ordered.log10_Molar.Dose)){
if(p==1){
Raw.AUC.with.negs.y.sig.positions = tgc$Raw.AUC.with.negs[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)])/10)}
else if(p>1){
new.Raw.AUC.with.negs.y.sig.positions<-tgc$Raw.AUC.with.negs[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)])/10)
Raw.AUC.with.negs.y.sig.positions<-c(Raw.AUC.with.negs.y.sig.positions,new.Raw.AUC.with.negs.y.sig.positions)
rm(new.Raw.AUC.with.negs.y.sig.positions)
}
}
rm(p)
# significance (star) labels (defined in statistics section, added here for reference):
# Corresponds to TukeyHSD$Condition:Dose output for condition-matched-doses
Raw.AUC.with.negs.y.sig.labels<-c(stat_Raw.AUC.with.negs.TUK.sig_By_Dose$sig.label)
# # # # # # # # #
ggplot(tgc, aes(x=log10_Molar.Dose, y=Raw.AUC.with.negs, colour=Condition)) +
annotate("text", x = c(Ordered.log10_Molar.Dose), y = Raw.AUC.with.negs.y.sig.positions, label = Raw.AUC.with.negs.y.sig.labels, size = 7.5) +
geom_errorbar(aes(ymin=Raw.AUC.with.negs-se,ymax=Raw.AUC.with.negs+se), width=.035, colour="black") +
geom_line() +
geom_point(size=3, shape=21,aes(fill=Condition)) +
xlab("log[Oxo] M") + ylab("Area Under the Curve (Arbitrary Units)") +
ggtitle(paste0("Effect of ", Experiment.name, " on Area Under the Curve: Raw Data (With Negatives)")) +
theme_bw() +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # #
# AUC (Raw, NEGATIVES = 0)
# Cumulative Plot line
# # # # # # # # # # # # # #
tgc <- summarySE(ANOVA.table, measurevar="Raw.AUC",groupvars=c("Condition","log10_Molar.Dose"))
# Obtain a vector of y values for placement of significance labels
# y coordinate = y(data.point.mean) + se(data.point.mean) + ( se(data.point.mean)/10 )
# obtain y coordinate for first data.point attach subsequent y coordinates to vector with each iteration (for all log10_Molar dose values)
for(p in 1:length(Ordered.log10_Molar.Dose)){
if(p==1){
Raw.AUC.y.sig.positions = tgc$Raw.AUC[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)])/10)}
else if(p>1){
new.Raw.AUC.y.sig.positions<-tgc$Raw.AUC[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)])/10)
Raw.AUC.y.sig.positions<-c(Raw.AUC.y.sig.positions,new.Raw.AUC.y.sig.positions)
rm(new.Raw.AUC.y.sig.positions)
}
}
rm(p)
# significance (star) labels (defined in statistics section, added here for reference):
# Corresponds to TukeyHSD$Condition:Dose output for condition-matched-doses
Raw.AUC.y.sig.labels<-c(stat_Raw.AUC.TUK.sig_By_Dose$sig.label)
# # # # # # # # #
ggplot(tgc, aes(x=log10_Molar.Dose, y=Raw.AUC, colour=Condition)) +
annotate("text", x = c(Ordered.log10_Molar.Dose), y = Raw.AUC.y.sig.positions, label = Raw.AUC.y.sig.labels, size = 7.5) +
geom_errorbar(aes(ymin=Raw.AUC-se,ymax=Raw.AUC+se), width=.035, colour="black") +
geom_line() +
geom_point(size=3, shape=21,aes(fill=Condition)) +
xlab("log[Oxo] M") + ylab("Area Under the Curve (Arbitrary Units)") +
ggtitle(paste0("Effect of ", Experiment.name, " on Area Under the Curve: Raw Data (Negatives = 0)")) +
theme_bw() +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # # #
# AUC, (Raw, WITH NEGATIVES)
# Cumulative Scatter/jitter plot
# # # # # # # # # # # # # # # # #
ggplot(ANOVA.table, aes(x=log10_Molar.Dose, y=Raw.AUC.with.negs, colour=Condition)) +
geom_jitter(size=.75, shape=1, width = .15, height = .8) +
xlab("log[Oxo] M") +
ylab("Area Under the Curve (Arbitrary Units)") +
facet_grid(. ~ Condition) +
ggtitle(paste0("Effect of ", Experiment.name, " on Area Under the Curve: Raw Data (With Negatives)")) +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
theme_gray() + theme(panel.margin = unit(2, "lines")) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreeen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # # #
# AUC (Raw, NEGATIVES = 0)
# Cumulative Scatter/jitter plot
# # # # # # # # # # # # # # # # #
ggplot(ANOVA.table, aes(x=log10_Molar.Dose, y=Raw.AUC, colour=Condition)) +
geom_jitter(size=.75, shape=1, width = .15, height = .8) +
xlab("log[Oxo] M") +
ylab("Area Under the Curve (Arbitrary Units)") +
facet_grid(. ~ Condition) +
ggtitle(paste0("Effect of ", Experiment.name, " on Area Under the Curve: Raw Data (Negatives = 0)")) +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
theme_gray() + theme(panel.margin = unit(2, "lines")) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreeen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # #
# Average Peak Duration
# Cumulative Plot line
# # # # # # # # # # # # #
tgc <- summarySE(ANOVA.table, measurevar="avg.peak.dur",groupvars=c("Condition","log10_Molar.Dose"))
# Obtain a vector of y values for placement of significance labels
# y coordinate = y(data.point.mean) + se(data.point.mean) + ( se(data.point.mean)/10 )
# obtain y coordinate for first data.point attach subsequent y coordinates to vector with each iteration (for all log10_Molar dose values)
for(p in 1:length(Ordered.log10_Molar.Dose)){
if(p==1){
avg.peak.dur.y.sig.positions = tgc$avg.peak.dur[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[1] & (tgc$Condition== B)])/10)}
else if(p>1){
new.avg.peak.dur.y.sig.positions<-tgc$avg.peak.dur[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)] +
(tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)]) +
((tgc$se[tgc$log10_Molar.Dose==Ordered.log10_Molar.Dose[p] & (tgc$Condition== B)])/10)
avg.peak.dur.y.sig.positions<-c(avg.peak.dur.y.sig.positions,new.avg.peak.dur.y.sig.positions)
rm(new.avg.peak.dur.y.sig.positions)
}
}
rm(p)
# significance (star) labels (defined in statistics section, added here for reference):
# Corresponds to TukeyHSD$Condition:Dose output for condition-matched-doses
avg.peak.dur.y.sig.labels<-c(stat_avg.peak.dur.TUK.sig_By_Dose$sig.label)
# # # # # # # # #
ggplot(tgc, aes(x=log10_Molar.Dose, y=avg.peak.dur, colour=Condition)) +
annotate("text", x = c(Ordered.log10_Molar.Dose), y = avg.peak.dur.y.sig.positions, label = avg.peak.dur.y.sig.labels, size = 7.5) +
geom_errorbar(aes(ymin=avg.peak.dur-se,ymax=avg.peak.dur+se), width=.035, colour="black") +
geom_line() +
geom_point(size=3, shape=21,aes(fill=Condition)) +
xlab("log[Oxo] M") +
ylab("Average Peak Duration (sec)") +
ggtitle(paste0("Effect of ", Experiment.name, " on Average Peak Duration")) +
theme_bw() +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreen','tomato'))
graph.A.count<-graph.A.count+1
# # # # # # # # # # # # # # # # #
# Average Peak Duration
# Cumulative Scatter/Jitter plot
# # # # # # # # # # # # # # # # #
ggplot(ANOVA.table, aes(x=log10_Molar.Dose, y=avg.peak.dur, colour=Condition)) +
geom_jitter(size=.75, shape=1, width = .15, height = .8) +
xlab("log[Oxo] M") +
ylab("Average Peak Duration (sec)") +
facet_grid(. ~ Condition) +
ggtitle(paste0("Effect of ", Experiment.name, " on Average Peak Duration")) +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
theme_gray() +
theme(panel.margin = unit(2, "lines")) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreeen','tomato'))
graph.A.count<-graph.A.count+1
amt.done<-paste0(round((as.numeric(graph.A.count)/as.numeric(graph.A.collection)*100)),"%")
setTkProgressBar(progress.bar,graph.A.count,
label=paste0("Generating Cumulative Oscillation plots ",
"\n", amt.done))
Sys.sleep(.5)
# Finish encoding and saving PDF file containing all plots generated above
dev.off()
total.runtime<-round(as.numeric(proc.time()[3])-t00, digits = 1)
total.runtime.str <- paste0(floor(total.runtime/60), ":", total.runtime-(floor(total.runtime/60)*60))
setTkProgressBar(progress.bar,graph.A.count,label=paste0("Graph Upload Complete: Opening PDF","\n","Run Time (MM:SS): ",total.runtime.str))
Sys.sleep(.75)
close(progress.bar)
# Locates and opens the pdf containing all final summary plots generated above
system2('open', args = pdf.name, wait = FALSE)
print("Fin")
#########################################################################################
# ICC Information ####
# download and install appropriate version of Java to match R
# Collect all ICC files based on simple pattern
all.ICC.files<-list.files(pattern = "ICC", recursive = TRUE)
# quick look to ensure everything is in order
all.ICC.files
# 1st Run-through = CTRL -> replicate 1 -> collect dose -> combine dose -> outputs REP1_CTRL table (***all results [1-4] are 'technical replicates' and should be combined as such)
#
# REPEAT FOR:
#
# replicate 2 -> collect dose -> combine dose -> outputs REP2_CTRL.table (***all results [1-4] are 'technical replicates' and should be combined as such)
# 2nd Run-through = "Condition tested"-> replicate 1 -> collect dose -> combine dose -> outputs REP1_Condition.table (***all results [1-4] are 'technical replicates' and should be combined as such)
#
# REPEAT FOR:
#
# replicate 2 -> collect dose -> combine dose -> outputs REP2_Condition.table (***all results [1-4] are 'technical replicates' and should be combined as such)
for(c in 1:length(Conditions)){
Condition.files<-list.files(pattern = "ICC")
Condition.files<-Condition.files[grep(Conditions[c], Condition.files)]
# Simply Condition.files, but split
designation_index<-strsplit(Condition.files, split ="_")
# Due to condition files drawing only relevant files, the 2nd file of any list should be the same i.e. "Null"
Condition.name<-unlist(strsplit(designation_index[[1]][4],split = ".xlsx"))
## designate correct file name per iteration
file.name<-paste(Condition.name, "Count","table", sep="_")
for(i in 1:length(Condition.files)){
# make the table first
if(i==1){
replicate.split<-strsplit(designation_index[[i]][3], split = ".xlsx")
replicate.name<-replicate.split[[1]][1]
tech.rep.split<-strsplit(designation_index[[i]][1], split = "ICC")
tech.rep.name<-tech.rep.split[[1]][2]
Dose.name<-designation_index[[1]][2]
big.table <- readWorksheetFromFile(Condition.files[1],sheet=1)
big.table$Replicate<-replicate.name
big.table$Position<-tech.rep.name
big.table$Condition<-Condition.name
big.table$Dose<-Dose.name
}
# then generate each new table, and bind to the first by row: results will ultimately be included in one conjoined table
# column data is preserved per cell label
if(i>1){
replicate.split<-strsplit(designation_index[[i]][3], split = ".xlsx")
replicate.name<-replicate.split[[1]][1]
tech.rep.split<-strsplit(designation_index[[i]][1], split = "ICC")
tech.rep.name<-tech.rep.split[[1]][2]
Dose.name<-designation_index[[i]][2]
next.block<-readWorksheetFromFile(Condition.files[i],sheet=1)
next.block$Replicate<-replicate.name
next.block$Position<-tech.rep.name
next.block$Condition<-Condition.name
next.block$Dose<-Dose.name
big.table<-rbind(big.table,next.block)
}}
big.table<-subset(big.table, select = -c(Slice,Type.1,Type.2,Type.3))
big.table<-rename(big.table, c("Type.4"="DAPI_Count", "Type.5"="Olig2_Count"))
big.table<-big.table[big.table$DAPI_Count>0,] # no more 0
# some data seems to have single slice information containing "total counts"
# FIJI results which have many slices (which all = 0) contain a blank row followed by TOTAL SUM
# this blank row gets translated as a row of NA values in R, leaving each table with a counted slice total and a duplicate row of the same values
# the code below ensures that when this happens, only even rows are reported so as to eliminate these duplicates
if(any(is.na(big.table))==TRUE)
{
big.table<-big.table[complete.cases(big.table),] # no more NA
}
# aesthetic reordering of big.table (and subsequent output)
big.table<-big.table[c("Condition","Replicate","Olig2_Count","DAPI_Count","Position","Dose")]
# # # # # # # # #
assign(eval(substitute(file.name)), big.table, envir=.GlobalEnv)
}
tables<-objects(pattern = "_table" )
# concatenate all ICC counts tables from the global environment into one consolidated table
for(t in 1:length(tables)){
if(t==1){
first<-get(tables[1])
}
if(t>1){
subsequent<-get(tables[t])
first<-rbind(first,subsequent)
}
assign(eval(substitute("Counts.table")), first, envir=.GlobalEnv)
}
# Aesthetic re-ordering of data in ICC counts table
Counts.table<-Counts.table[order(Counts.table$Replicate), ]
# Cleanup
rm(first)
rm(subsequent)
rm(next.block)
rm(big.table)
rm(t)
rm(i)
rm(c)
rm(Condition.name)
rm(Dose.name)
rm(file.name)
rm(replicate.name)
rm(tech.rep.name)
rm(replicate.split)
rm(tech.rep.split)
rm(Condition.files)
rm(designation_index)
# Calculation of ICC counts as relative percentages
Total_DAPI_Count<-as.numeric(sum(Counts.table$DAPI_Count,na.rm = TRUE))
Total_Olig2_Count<-as.numeric(sum(Counts.table$Olig2_Count,na.rm = TRUE))
Total.Olig2.Percent<-paste0(round((Total_Olig2_Count/Total_DAPI_Count)*100, digits = 0)," %")
ICC.Replicates<-unique(Counts.table$Replicate)
REP1<-subset(Counts.table, Replicate==ICC.Replicates[1])
rep1.Total_Olig2_Count<-sum(REP1$Olig2_Count)
rep1.Total_DAPI_Count<-sum(REP1$DAPI_Count)
rep1.Olig2.Percent<-paste0("Rep1 Olig2/DAPI: ", round((rep1.Total_Olig2_Count/rep1.Total_DAPI_Count)*100, digits = 0), " %")
REP2<-subset(Counts.table, Replicate==ICC.Replicates[2])
rep2.Total_Olig2_Count<-sum(REP2$Olig2_Count)
rep2.Total_DAPI_Count<-sum(REP2$DAPI_Count)
rep2.Olig2.Percent<-paste0("Rep2 Olig2/DAPI: ",round((rep2.Total_Olig2_Count/rep2.Total_DAPI_Count)*100, digits = 0), " %")
ICC_All.table<-rbind(REP1,REP2)
ICC_All.table<-ICC_All.table[order(ICC_All.table$Condition), ]
# # # # # # # # # # # #
## Total Counts
print.noquote(c("Total_Olig2_Count:", Total_Olig2_Count))
print.noquote(c("Total_DAPI_Count:", Total_DAPI_Count))
print.noquote(c("Total.Olig2.Percent:", Total.Olig2.Percent))
# # # # # # # # # # # #
## Total Counts By Replicate
# # # # # # # # #
# Replicate 1 #
# # # # # # # # #
print.noquote(c("rep1.Total_Olig2_Count:", rep1.Total_Olig2_Count))
print.noquote(c("rep1.Total_DAPI_Count:", rep1.Total_DAPI_Count))
print.noquote(c("rep1.Olig2.Percent:", rep1.Olig2.Percent))
# # # # # # # # #
# Replicate 2 #
# # # # # # # # #
print.noquote(c("rep2.Total_Olig2_Count:", rep2.Total_Olig2_Count))
print.noquote(c("rep2.Total_DAPI_Count:", rep2.Total_DAPI_Count))
print.noquote(c("rep2.Olig2.Percent:", rep2.Total_DAPI_Count))
Olig2.avg.rep.percent<-paste0(round(((
((rep1.Total_Olig2_Count/rep1.Total_DAPI_Count)*100)
+((rep2.Total_Olig2_Count/rep2.Total_DAPI_Count)*100))/2),
digits = 1), " %")
print.noquote(c("Olig2 Replicate %'s averaged:", Olig2.avg.rep.percent))
ICC.file.label<-paste0("_ICC_Dtaining_Data_For_", Experiment.name, "_Experiment_", current.date,".xlsx")
if(file.exists(ICC.file.label))
{file.remove(ICC.file.label)}
wb<-loadWorkbook(ICC.file.label, create=TRUE)
createSheet(wb,name = "ICC_counts")
writeWorksheet(wb, data=Counts.table,sheet = "ICC_counts")
saveWorkbook(wb)
rm(wb)
# # # # # # # # # # # #
######################
# Organization of Additional Results ####
### FREQ Control #
# Total cells with >1 peak in control
nrow(Freq_CTRL.table)
# Quick look: Average across Samples
((nrow(Freq_CTRL.table[Freq_CTRL.table$Replicate==1,]))+(nrow(Freq_CTRL.table[Freq_CTRL.table$Replicate==2,])))/2
### FREQ M3R #
# Hard Coded for quickly checking M3R results
# Total M3R cells with >1 peak
nrow(Freq_M3R.table)
# Average across Samples
((nrow(Freq_M3R.table[Freq_M3R.table$Replicate==1,]))+(nrow(Freq_M3R.table[Freq_M3R.table$Replicate==2,])))/2
# ANOVA (numpeaks,amps,etc) cell numbers as well as average across both samples
### Control #
# Total
nrow(ANOVA_CTRL.table)
# Average across Samples
((nrow(ANOVA_CTRL.table[ANOVA_CTRL.table$Replicate==1,]))+(nrow(ANOVA_CTRL.table[ANOVA_CTRL.table$Replicate==2,])))/2
### M3R #
# Total
nrow(ANOVA_M3R.table)
# Average across Samples
((nrow(ANOVA_M3R.table[ANOVA_M3R.table$Replicate==1,]))+(nrow(ANOVA_M3R.table[ANOVA_M3R.table$Replicate==2,])))/2
## FREQ AS A % OF TOTAL RESPONSIVE (I.E. [SUM OF ALL IN FREQ.TABLE/ SUM OF ALL IN ANOVA.TABLE] * 100)
# ALL
(nrow(Freq.table)/nrow(ANOVA.table))*100
# CTRL
(nrow(Freq_CTRL.table)/nrow(ANOVA_CTRL.table))*100
# M3R
(nrow(Freq_M3R.table)/nrow(ANOVA_M3R.table))*100
# # # # # #
Totals.table<-Response.table[Response.table$Type=="Total n",]
Correct.table<-Response.table[Response.table$Type=="Correct Responses",]
Control.table<-Totals.table[Totals.table$Condition=="CTRL",]
M3R.table<-Totals.table[Totals.table$Condition=="M3R",]
# Total Cells Tested
Total.cells.tested<-sum(Totals.table$`Number of Cells`)
Total.cells.tested
# Total cells with Correct Responses
Total.cells.correct<-sum(Correct.table$`Number of Cells`)
Total.cells.correct
## Per Condition:
# % Control responsive
(sum(Correct.table$`Number of Cells`[Correct.table$Condition=="CTRL"])/Total.cells.correct)*100
# % M3R responsive
(sum(Correct.table$`Number of Cells`[Correct.table$Condition=="M3R"])/Total.cells.correct)*100
## By Dose
# ( Conducted in graphpad )
### adding these numbers to table (for organization)
Control.table$Sum.all.conditions.Tested<-Total.cells.tested
Control.table$Total.per.condition<-sum(Totals.table$`Number of Cells`[Totals.table$Condition=="CTRL"])
Total.CTRL.tested<-sum(Totals.table$`Number of Cells`[Totals.table$Condition=="CTRL"])
Control.table$Sum.all.correct<-Total.cells.correct
Control.table$Correct.per.condition<-sum(Correct.table$`Number of Cells`[Correct.table$Condition=="CTRL"])
Total.CTRL.correct<-sum(Correct.table$`Number of Cells`[Correct.table$Condition=="CTRL"])
M3R.table$Sum.all.conditions.Tested<-Total.cells.tested
M3R.table$Total.per.condition<-sum(Totals.table$`Number of Cells`[Totals.table$Condition=="M3R"])
Total.M3R.tested<-sum(Totals.table$`Number of Cells`[Totals.table$Condition=="M3R"])
M3R.table$Sum.all.correct<-Total.cells.correct
M3R.table$Correct.per.condition<-sum(Correct.table$`Number of Cells`[Correct.table$Condition=="M3R"])
Total.M3R.Correct<-sum(Correct.table$`Number of Cells`[Correct.table$Condition=="M3R"])
A_Final.Breakdown<-rbind(Control.table,M3R.table)
A_Final.Breakdown<-subset(A_Final.Breakdown, select = -c(Type,`Number of Cells`, `Percent of Population`))
A_Final.Breakdown$Avg.per.image.field.whole<-Avg.cells.per.image.field<-sum(Totals.table$`Number of Cells`/76)
# Additional Breakdown Info
# Total cells per technicical replicate (technical replicate denotes all cells imaged in a single well)
# i.e.: For each condition:: (2 positions imaged per well; 1 well per plate (tech.rep); 2 plates per sample; 2 samples tested)
# and sample replicate (n = 2)
#
totals.rep1<-Totals.table[Totals.table$Replicate==1,]
totals.rep2<-Totals.table[Totals.table$Replicate==2,]
## averages per sample
Total.cells.rep1<-sum(totals.rep1$`Number of Cells`)
Total.cells.rep2<-sum(totals.rep2$`Number of Cells`)
# CTRL
CTRL.rep1<-sum(totals.rep1$`Number of Cells`[totals.rep1$Condition=="CTRL"])
CTRL.rep2<-sum(totals.rep2$`Number of Cells`[totals.rep1$Condition=="CTRL"])
CTRL.AVG.samples<-(CTRL.rep1 + CTRL.rep2)/2
# M3R
M3R.rep1<-sum(totals.rep1$`Number of Cells`[totals.rep1$Condition=="M3R"])
M3R.rep2<-sum(totals.rep2$`Number of Cells`[totals.rep2$Condition=="M3R"])
M3R.AVG.samples<-(M3R.rep1+M3R.rep2)/2
# # # # # # #
# TECH REPS
totals.rep1.1<-totals.rep1[totals.rep1$Position<2,]
totals.rep1.2<-totals.rep1[totals.rep1$Position>1,]
sum(totals.rep1.1$`Number of Cells`)
sum(totals.rep1.2$`Number of Cells`)
totals.rep2.1<-totals.rep2[totals.rep2$Position<2,]
totals.rep2.2<-totals.rep2[totals.rep2$Position>1,]
#Total Number of Cells replicate 2
sum(totals.rep2.1$`Number of Cells`)
sum(totals.rep2.2$`Number of Cells`)
correct.rep1<-Correct.table[Correct.table$Replicate==1,]
correct.rep2<-Correct.table[Correct.table$Replicate==2,]
correct.rep1.1<-correct.rep1[correct.rep1$Position<2,]
correct.rep1.2<-correct.rep1[correct.rep1$Position>1,]
sum(correct.rep1.1$`Number of Cells`)
sum(correct.rep1.2$`Number of Cells`)
correct.rep2.1<-correct.rep2[correct.rep2$Position<2,]
correct.rep2.2<-correct.rep2[correct.rep2$Position>1,]
sum(correct.rep2.1$`Number of Cells`)
sum(correct.rep2.2$`Number of Cells`)
sum(correct.rep1$`Number of Cells`)
sum(correct.rep2$`Number of Cells`)
###########################################
# Quick export of additional statistical information for analysis in prism AND/OR Graphical elements for figure add-ons ########
# # # # # # # # # # # # # # #
# Significance by parameter:#
# # # # # # # # # # # # # # #
parameters<-c()
for(P in 7:ncol(ANOVA.table)){
if(P==7){
parameters<-c(colnames(ANOVA.table)[7])
} else if(P>7){
parameters<-c(parameters,colnames(ANOVA.table)[P])
}
}
ANOVA.ctrl.table<-ANOVA.table[ANOVA.table$Condition=="CTRL",]
Average.amp<-summarySE(ANOVA.table, measurevar="Average_Amplitude",groupvars=c("Condition","log10_Molar.Dose"))
Average.amp$`_____`<-" "
Average.amp$`Average Amplitude`<-" "
Average.amp<-Average.amp[,c(ncol(Average.amp),2:ncol(Average.amp)-1)]
Maximum.amp<-summarySE(ANOVA.table, measurevar="mA",groupvars=c("Condition","log10_Molar.Dose"))
Maximum.amp$`_____`<-" "
Maximum.amp$`Maximum Amplitude`<-" "
Maximum.amp<-Maximum.amp[,c(ncol(Maximum.amp),2:ncol(Maximum.amp)-1)]
Number.of.peaks<-summarySE(ANOVA.table, measurevar="numpeaks",groupvars=c("Condition","log10_Molar.Dose"))
Number.of.peaks$`_____`<-" "
Number.of.peaks$`Number of peaks`<-" "
Number.of.peaks<-Number.of.peaks[,c(ncol(Number.of.peaks),2:ncol(Number.of.peaks)-1)]
Normalized.AUC.with.negs<-summarySE(ANOVA.table, measurevar="Norm.AUC",groupvars=c("Condition","log10_Molar.Dose"))
Normalized.AUC.with.negs$`_____`<-" "
Normalized.AUC.with.negs$`Normalized AUC with negs`<-" "
Normalized.AUC.with.negs<-Normalized.AUC.with.negs[,c(ncol(Normalized.AUC.with.negs),2:ncol(Normalized.AUC.with.negs)-1)]
Stats.table<-cbind(Average.amp,Maximum.amp,Number.of.peaks,Normalized.AUC.with.negs)
stat.file.label<-paste("_",Experiment.name,"Average_stats_for_significant_parameters",current.date,".xlsx")
if(file.exists(stat.file.label))
{file.remove(stat.file.label)}
wb<-loadWorkbook(stat.file.label,create = TRUE)
createSheet(wb, name="significance ± SEM")
writeWorksheet(wb,Stats.table,sheet = "significance ± SEM")
saveWorkbook(wb)
rm(wb)
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # #
# Total N & cumulative number of peaks by dose: #
# # # # # # # # # # # # # # # # # # # # # # # # #
Control_stat_table<-data.frame(Molar.Dose)
GNB4_stat_table<-data.frame(Molar.Dose)
Control_nCell.dose<-c()
GNB4_nCell.dose<-c()
for(W in 1:length(Molar.Dose)){
if(W==1){
ANOVA.table.search<-ANOVA.table[ANOVA.table$Dose==Molar.Dose[1]&ANOVA.table$Condition==A,]
Control_nCell_by.dose<-nrow(ANOVA.table.search)
Control_cum.numpeaks_by.dose<-sum(ANOVA.table.search$numpeaks)
}
if(W>1)
{
ANOVA.table.search<-ANOVA.table[ANOVA.table$Dose==Molar.Dose[W]&ANOVA.table$Condition==A,]
new.N<-nrow(ANOVA.table.search)
Control_nCell_by.dose<-c(Control_nCell_by.dose,new.N)
new.numpeaks<-sum(ANOVA.table.search$numpeaks)
Control_cum.numpeaks_by.dose<-c(Control_cum.numpeaks_by.dose,new.numpeaks)
}
}
for(W in 1:length(Molar.Dose)){
if(W==1){
ANOVA.table.search<-ANOVA.table[ANOVA.table$Dose==Molar.Dose[1]&ANOVA.table$Condition==B,]
GNB4_nCell_by.dose<-nrow(ANOVA.table.search)
GNB4_cum.numpeaks_by.dose<-sum(ANOVA.table.search$numpeaks)
}
if(W>1)
{
ANOVA.table.search<-ANOVA.table[ANOVA.table$Dose==Molar.Dose[W]&ANOVA.table$Condition==B,]
new.N<-nrow(ANOVA.table.search)
GNB4_nCell_by.dose<-c(GNB4_nCell_by.dose,new.N)
new.numpeaks<-sum(ANOVA.table.search$numpeaks)
GNB4_cum.numpeaks_by.dose<-c(GNB4_cum.numpeaks_by.dose,new.numpeaks)
}
}
rm(W)
Control_stat_table$N<-Control_nCell_by.dose
Control_stat_table$`cumulative peaks`<-Control_cum.numpeaks_by.dose
GNB4_stat_table$N<-GNB4_nCell.by.dose
GNB4_stat_table$`cumulative peaks`<-GNB4_cum.numpeaks_by.dose
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# setting up individual tables for by dose as variables of name: dose_X.table, where X = dose (1,2, etc)
for(I in 1:length(Molar.Dose)){
dose.table<-ANOVA.table[ANOVA.table$Dose==Molar.Dose[I],]
dose.table.label<-paste0("Dose_",I,".table___",Ordered.dose.names[I])
assign(eval(substitute(paste(dose.table.label))), dose.table, envir=.GlobalEnv)
rm(dose.table)
rm(dose.table.label)
}
rm(I)
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Multi parameter dotplots for GNB4 experiments #
# # # # # # # # # # # # # # # # # ## # # # # # # # # # # # # #
relevant.doses<-c(Molar.Dose[2],Molar.Dose[4],Molar.Dose[6])
Relevant.table<-rbind(Dose_2.table___250nM,Dose_4.table___2.5µM,Dose_6.table___25µM)
Relevant.table$Dose=factor(Relevant.table$Dose)
Relevant.table<-Relevant.table[order(Relevant.table$Dose), ]
relevant.doses.big<-Relevant.table$Dose
relevant.log10_Molar.Dose<-Relevant.table$log10_Molar.Dose
# pdf("_Three parameter graphs.pdf")
# # # # # # # # # #
### Average Amp ###
# # # # # # # # # #
# figure graph (different styles are commented out, uncomment as per user preference)
A<-ggplot(Relevant.table, aes(x=Relevant.table$Dose, y=Average_Amplitude, colour=Condition, shape=Condition)) +
# geom_jitter(size=.5, position = position_jitterdodge(jitter.width =0, dodge.width = .1)) +
geom_dotplot(binaxis='y', stackdir='center', dotsize=1,binwidth = .02,position = position_dodge(width = .8)) +
# geom_boxplot( outlier.shape = NA) +
xlab("log[Oxo] M") +
ylab("Average Fold Change from Baseline") +
ggtitle(paste0("Average Amplitude")) +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
theme_bw() +
scale_y_continuous(limits = c(0,max(Relevant.table$Average_Amplitude)+.05)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) +
theme(axis.text.x=element_text(angle = -90, hjust = 0)) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreeen','tomato'))
## Average Amp ~ Figure add-on: error bars ± sem; same scale ##
avg.amp.summary <- summarySE(Relevant.table, measurevar="Average_Amplitude",groupvars=c("Condition","Dose"))
A_mean_sem <- ggplot(avg.amp.summary, aes(x=avg.amp.summary$Dose, y=Average_Amplitude, colour=Condition)) +
geom_boxplot(lwd = .15) +
geom_errorbar(aes(ymin=Average_Amplitude-se,ymax=Average_Amplitude+se), width=.05) +
xlab("log[Oxo] M") +
ylab("Average Fold Change from Baseline") +
ggtitle(paste0("Average Amplitude")) +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
theme_bw() +
scale_y_continuous(limits = c(0,max(Relevant.table$Average_Amplitude)+.05)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) +
theme(axis.text.x=element_text(angle = -90, hjust = 0)) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreeen','tomato'))
# # # # # # # # # # #
# number of peaks #
# # # # # # # # # # #
# figure graph (different styles are commented out, uncomment as per user preference)
B<-ggplot(Relevant.table, aes(x=Relevant.table$Dose, y=numpeaks, colour=Condition, shape = Condition)) +
# geom_jitter(size=, position = position_jitterdodge(jitter.width =0, dodge.width = .1)) +
geom_dotplot(binaxis='y', stackdir='center', dotsize=1,binwidth = .01,position = position_dodge(width = 1.1)) +
# geom_boxplot(outlier.shape = NA) +
xlab("log[Oxo] M") +
ylab("Peaks") +
ggtitle(paste0("Number of Peaks")) +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
theme_bw() +
theme(panel.margin = unit(2, "lines")) +
theme(axis.text.x=element_text(angle = -90, hjust = 0)) +
scale_y_continuous(limits = c(0,max(Relevant.table$numpeaks)+.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreen','tomato'))
## Peaks ~ Figure add-on: error bars ± sem; same scale ##
peaks.summary <- summarySE(Relevant.table, measurevar="numpeaks",groupvars=c("Condition","Dose"))
B_mean_sem<-ggplot(peaks.summary, aes(x=peaks.summary$Dose, y=numpeaks, colour=Condition)) +
geom_boxplot(lwd = .15) +
geom_errorbar(aes(ymin=numpeaks-se,ymax=numpeaks+se), width=.05) +
xlab("log[Oxo] M") +
ylab("Peaks") +
ggtitle(paste0("Number of Peaks")) +
theme(legend.justification=c(0,1),legend.position=c(0,1)) +
theme_bw() +
theme(panel.margin = unit(2, "lines")) +
theme(axis.text.x=element_text(angle = -90, hjust = 0)) +
scale_y_continuous(limits = c(0,max(Relevant.table$numpeaks)+.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreen','tomato'))
# # # # # # #
# FREQ #
# # # # # # #
Relevant.freq.table<-Relevant.table[Relevant.table$numpeaks>1,]
relevant.freq.doses.big<-Relevant.freq.table$Dose
# figure graph (different styles are commented out, uncomment as per user preference)
C<-ggplot(Relevant.freq.table, aes(x=Relevant.freq.table$Dose, y=freq, colour=Condition,shape=Condition)) +
# geom_jitter(size=.5, position = position_jitterdodge(jitter.width =0, dodge.width = .1)) +
geom_dotplot(binaxis='y', stackdir='center', dotsize=1,binwidth = .3, position = position_dodge(width = 1.2)) +
# geom_boxplot(outlier.shape = NA) +
xlab("log[Oxo] M") +
ylab("Frequency (mHz)") +
ggtitle(paste0("Frequency")) +
theme(legend.position ="none") +
theme_bw() +
theme(axis.text.x=element_text(angle = -90, hjust = 0)) +
scale_y_continuous(limits = c(0,max(Relevant.freq.table$freq)+5)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreeen','tomato'))
## Frequency ~ Figure add-on: error bars ± sem; same scale ##
freq.summary <- summarySE(Relevant.freq.table, measurevar="freq",groupvars=c("Condition","Dose"))
C_mean_sem<-ggplot(freq.summary, aes(x=freq.summary$Dose, y=freq, colour=Condition)) +
geom_boxplot(lwd = .15) +
geom_errorbar(aes(ymin=freq-se,ymax=freq+se), width=.05) +
xlab("log[Oxo] M") +
ylab("Frequency (mHz)") +
ggtitle(paste0("Frequency")) +
theme(legend.position ="none") +
theme_bw() +
theme(axis.text.x=element_text(angle = -90, hjust = 0)) +
scale_y_continuous(limits = c(0,max(Relevant.freq.table$freq)+5)) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) +
scale_color_manual(values=c('limegreen','tomato')) +
scale_fill_manual(values=c('limegreeen','tomato'))
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
pdf("_Mean ± SEM (lines)_.pdf")
grid.arrange(A_mean_sem + theme(legend.position="none") + theme(legend.position="none"),
B_mean_sem + theme(legend.position="none") + theme(legend.position="none"),
C_mean_sem + theme(legend.position="none") + theme(legend.position="none"),
nrow=1)
A_mean_sem + theme(legend.position="none")
B_mean_sem + theme(legend.position="none")
C_mean_sem + theme(legend.position="none")
dev.off()
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
monotonic.table<-ANOVA.table[ANOVA.table$numpeaks==1,]
monotonic.CTRL.table<-ANOVA.table[ANOVA.table$numpeaks==1&ANOVA.table$Condition=="CTRL",]
monotonic.Condition.table<-ANOVA.table[ANOVA.table$numpeaks==1&ANOVA.table$Condition!="CTRL",]
PEAK <- summarySE(Relevant.table, measurevar="numpeaks" ,groupvars=c("Condition","log10_Molar.Dose"))
AMP <- summarySE(Relevant.table, measurevar="numpeaks" ,groupvars=c("Condition","log10_Molar.Dose"))
FREQ<- summarySE(Relevant.freq.table, measurevar="numpeaks" ,groupvars=c("Condition","log10_Molar.Dose"))
###########################################################################################################################