From 329c63caefa70d3495ade7dcdcb6cc3c2cb8e710 Mon Sep 17 00:00:00 2001 From: Minh Date: Wed, 15 Jul 2015 02:49:57 -0500 Subject: [PATCH] Demo AUC --- R/GCAT/R/GCAT.main.R | 27 +- R/GCAT/R/class.well.R | 9 +- R/GCAT/R/fitted.calculations.R | 702 +++++++++--------- R/GCAT/R/plot.fit.R | 24 +- R/GCAT/man/GCAT.Rd | 20 + R/GCAT/man/gcat.analysis.main.Rd | 4 +- R/GCAT/man/gcat.fit.main.Rd | 4 +- R/GCAT/man/gcat.output.main.Rd | 8 +- R/GCAT/man/plot-well-missing-method.Rd | 4 +- R/GCAT/man/well-class.Rd | 2 + .../YPDAFEXglucoseTests_2-25-10__defaults.R | 2 +- Rails/app/assets/stylesheets/gcat.css.scss | 2 +- Rails/app/models/assay.rb | 23 +- Rails/app/views/assays/_form.html.erb | 9 + 14 files changed, 468 insertions(+), 372 deletions(-) create mode 100644 R/GCAT/man/GCAT.Rd diff --git a/R/GCAT/R/GCAT.main.R b/R/GCAT/R/GCAT.main.R index 8ba1cbe..3893a27 100644 --- a/R/GCAT/R/GCAT.main.R +++ b/R/GCAT/R/GCAT.main.R @@ -110,6 +110,7 @@ global.version.number = packageDescription(pkg="GCAT")$Version #' @param totalRange The heatmap specific range for the achieved growth on log scale. #' @param totalODRange The heatmap specific range for the achieved growth on linear (OD) scale. #' @param specRange The heatmap specific range for spec growth rate. +#' @param sranges The 2 boundaries for calculating the area under the curve. #' #' @return Depending on return.fit setting, an array of fitted well objects or a list of output files #' @@ -120,7 +121,7 @@ gcat.analysis.main = function(file.list, single.plate, layout.file = NULL, use.linear.param = F, use.loess = F, smooth.param=0.1, lagRange = NA, totalRange = NA, totalODRange = NA, specRange = NA, points.to.remove = 0, remove.jumps = F, time.input = NA, - plate.nrow = 8, plate.ncol = 12, input.skip.lines = 0, + plate.nrow = 8, plate.ncol = 12, sranges = NA, input.skip.lines = 0, multi.column.headers = c("Plate.ID", "Well", "OD", "Time"), single.column.headers = c("","A1"), layout.sheet.headers = c("Strain", "Media Definition"), silent = T, verbose = F, return.fit = F, overview.jpgs = T){ @@ -152,6 +153,16 @@ gcat.analysis.main = function(file.list, single.plate, layout.file = NULL, } } + if (!identical(sranges, NA)) { + sranges.string = paste(sranges, collapse = "-") + if (length(sranges) != 2) + exception("", paste("Should have 2 values for the bounds of integration. Bad range:", sranges.string)) + if (sranges[2] < sranges[1]) + exception("", paste("Bad range. The first range cannot be greater than the second:", sranges.string)) + if (sranges[1] < 0 | sranges[2] < 0) + exception("", paste("Bad range. Either one of the range is negative:", sranges.string)) + } + # MB: Now add.constant will always be 0. # No need to check. #if (add.constant < 0) @@ -179,7 +190,7 @@ gcat.analysis.main = function(file.list, single.plate, layout.file = NULL, use.linear.param=use.linear.param, use.loess=use.loess, smooth.param=smooth.param, plate.nrow = plate.nrow, plate.ncol = plate.ncol, multi.column.headers = multi.column.headers, single.column.headers = single.column.headers, layout.sheet.headers = layout.sheet.headers, - input.skip.lines = input.skip.lines, silent = silent, verbose = verbose), silent = T) + input.skip.lines = input.skip.lines, silent = silent, verbose = verbose, sranges = sranges), silent = T) # Return error message if the function fails. if(class(fitted.well.array) == "try-error") @@ -197,7 +208,7 @@ gcat.analysis.main = function(file.list, single.plate, layout.file = NULL, growth.cutoff = growth.cutoff, add.constant = add.constant, blank.value = blank.value, start.index = start.index, points.to.remove = points.to.remove, remove.jumps = remove.jumps, lagRange = lagRange, specRange = specRange, totalRange = totalRange, totalODRange = totalODRange, - out.dir = out.dir, graphic.dir = graphic.dir, overview.jpgs=overview.jpgs, + sranges = sranges, out.dir = out.dir, graphic.dir = graphic.dir, overview.jpgs=overview.jpgs, use.linear.param=use.linear.param, use.loess=use.loess, plate.ncol = plate.ncol, plate.nrow = plate.nrow, silent = silent, main.envir = main.envir), silent = T) @@ -261,6 +272,7 @@ gcat.analysis.main = function(file.list, single.plate, layout.file = NULL, #' @param layout.sheet.headers The headers of the layout file. #' @param growth.model What growth model should be used? #' @param backup.growth.model If the main growth model fails, the back up model will be used. +#' @param sranges The 2 boundaries for calculating the area under the curve. #' @param silent Surpress all messages. #' @param verbose Display all messages when analyzing each well. #' @@ -270,7 +282,7 @@ gcat.fit.main = function(file.name, input.data = NULL, load.type = "csv", layout normalize.method = "default", add.constant = 1, use.log = T, points.to.remove = 0, use.linear.param=F, use.loess=F, smooth.param=0.1, fall.cutoff = -0.0025, growth.cutoff = 0.05, remove.jumps = F, - plate.nrow = 8, plate.ncol = 12, input.skip.lines = 0, + plate.nrow = 8, plate.ncol = 12, sranges = NA, input.skip.lines = 0, multi.column.headers = c("Plate.ID", "Well", "OD", "Time"), single.column.headers = c("","A1"), layout.sheet.headers = c("Strain", "Media Definition"), growth.model = NA, backup.growth.model = NA, @@ -452,7 +464,7 @@ gcat.fit.main = function(file.name, input.data = NULL, load.type = "csv", layout # Return an error if there is a problem with model fitting if (class(well.array) == "try-error") stop("Error in : ", well.array) - + well.array = try(AUC_well(well.array, sranges[1], sranges[2], silent = silent)) if(!silent) cat("\ndone!\n") return(well.array) } @@ -495,11 +507,12 @@ gcat.fit.main = function(file.name, input.data = NULL, load.type = "csv", layout #' @param totalODRange The heatmap specific range for the achieved growth on linear (OD) scale. #' @param specRange The heatmap specific range for spec growth rate. #' @param main.envir starting environment of gcat.analysis.main(), captured as a list, printed out for debugging +#' @param sranges The 2 boundaries for calculating the area under the curve. #' #' @return A list of output files if success. gcat.output.main = function(fitted.well.array, out.prefix = "", source.file.list, upload.timestamp = NULL, add.constant, blank.value, start.index, growth.cutoff, points.to.remove, remove.jumps, - out.dir = getwd(), graphic.dir = paste(out.dir,"/pics",sep = ""), overview.jpgs = T, + out.dir = getwd(), graphic.dir = paste(out.dir,"/pics",sep = ""), overview.jpgs = T, sranges = NA, use.linear.param=F, use.loess=F, lagRange = NA, totalRange = NA, totalODRange = NA, specRange = NA, plate.nrow = 8, plate.ncol = 12, unlog = F, silent = T, main.envir){ @@ -571,7 +584,7 @@ gcat.output.main = function(fitted.well.array, out.prefix = "", source.file.list graphic.files = try(pdf.by.plate(fitted.well.array, out.prefix=out.prefix, upload.timestamp = upload.timestamp, unlog=unlog,constant.added=add.constant,overview.jpgs=overview.jpgs, lagRange = lagRange, specRange = specRange, totalRange = totalRange, - totalODRange = totalODRange, plate.ncol = plate.ncol, plate.nrow = plate.nrow),silent=silent) + totalODRange = totalODRange, sranges = sranges, plate.ncol = plate.ncol, plate.nrow = plate.nrow),silent=silent) if (class(graphic.files) == "try-error") stop("Error in : ", graphic.files) diff --git a/R/GCAT/R/class.well.R b/R/GCAT/R/class.well.R index 2819bdd..503b183 100644 --- a/R/GCAT/R/class.well.R +++ b/R/GCAT/R/class.well.R @@ -53,6 +53,7 @@ setOldClass("loess") #' @slot rss - residual sum of squares #' @slot loess - object returned by running loess on the normalized well data #' @slot nls - object returned by running nls on the normalized well data +#' @slot auc - The area under the curve. #' #' @export setClass("well", representation(position = "character", @@ -71,7 +72,8 @@ setClass("well", representation(position = "character", inflection.time = "numeric", rss = "numeric", loess = "loess", - nls = "nls")) + nls = "nls", + auc = "character")) #' Accessors for the well class #' @@ -304,13 +306,14 @@ setMethod("show", "well", #' @param show.calc draw lines that illustrate growth curve parameters #' @param draw.guess initial guess model. Drawn if specified #' @param well.number the number of the well in an array of wells +#' @param sranges The boundaries to calculate the AUC. #' @param ... additional arguments passed to the generic plot function #' #' @export setMethod("plot", signature(x = "well", y="missing"), function (x, y, constant.added = 1.0, xlim = NULL, ylim = NULL, - well.number = NULL, scale = 1, number.points = T, draw.symbols = F, show.text = T, show.calc = T, draw.guess = NULL, ...) + well.number = NULL, scale = 1, number.points = T, draw.symbols = F, show.text = T, show.calc = T, draw.guess = NULL, sranges = NA, ...) { # Determine the boundaries for the axes (if user did not specify them) if(is.null(ylim)){ @@ -347,7 +350,7 @@ setMethod("plot", # Show calculated parameters if specified. if (show.calc) - draw.calc.par(x, scale = scale * 0.5, constant.added = constant.added) + draw.calc.par(x, scale = scale * 0.5, constant.added = constant.added, sranges = sranges) # Draw initial guess if a model is specified. if (class(draw.guess) == "model"){ diff --git a/R/GCAT/R/fitted.calculations.R b/R/GCAT/R/fitted.calculations.R index 574d9be..366dd87 100644 --- a/R/GCAT/R/fitted.calculations.R +++ b/R/GCAT/R/fitted.calculations.R @@ -1,342 +1,360 @@ -#Copyright 2012 The Board of Regents of the University of Wisconsin System. -#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, -#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman -# -#This file is part of GCAT. -# -#GCAT is free software: you can redistribute it and/or modify -#it under the terms of the GNU Lesser General Public License as published by -#the Free Software Foundation, either version 3 of the License, or -#(at your option) any later version. -# -#GCAT is distributed in the hope that it will be useful, -#but WITHOUT ANY WARRANTY; without even the implied warranty of -#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -#GNU Lesser General Public License for more details. -# -#You should have received a copy of the GNU Lesser General Public License -#along with GCAT. If not, see . - -######################################################################## -# # -# Functions to calculate various things about wells based on fit model # -# # -######################################################################## - -# S3 generic for lag -lag <- function(fitted.well, ...) -{ - UseMethod("lag") -} - -# -# Common arguments: -# fitted.well - should be a well containing the results of , most functions will return NA if well has not been fit yet. -# unlog - should the value be returned on the linear scale as opposed to the log-transformed scale? -# constant.added - for returning values on the linear scale, what was the constant added before the log transform? -# digits - passed to the function, default is no rounding (infinity digits) - -# Transform values back to OD scale -unlog = function(x, constant.added) { - exp(x) - constant.added -} - -# Evaluate estimated OD at any timepoints using the fitted model -well.eval = function(fitted.well, Time = NULL){ - # If no timepoints are provided, use the ones collected in the experiment itself. - if(!is.numeric(Time)) - Time = data.from(fitted.well)$Time - - # Use of equation is deprecated. Use nls and loess models stored in the well object instead - # Attempt to use with the fitted equation and parameters to get estimates for OD at the given timepoints. - #output = try(eval(fitted.well@equation, fitted.well@fit.par), silent = T) - - # Predict log.OD value(s) using nls model if present. If no nls model, try using loess. - if (length(fitted.well@nls)>0) { - output = try(predict(fitted.well@nls,list(Time=Time)),silent=T) - } else if (length(fitted.well@loess)>0) { - output = try(predict(fitted.well@loess,Time),silent=T) - } else { - output = NA - } - - # Return values. If OD evaluation failed for any reason, return NULL. - if (is.numeric(output)){ - return(output) - } else { - return(NULL) - } -} - -# Evaluate model residuals using the measured vs. fitted log.OD values -model.residuals = function(fitted.well, unlog = F){ - measured.OD = data.from(fitted.well)[,2] - - # Use with no Time argument to get fitted OD values at measured timepoints. - predicted.OD = well.eval(fitted.well) - - # If all values are valid, return the differences - if (!is.numeric(predicted.OD)) - return(NA) - else - return(measured.OD - predicted.OD) - } - -# Evaluate deviations of log.OD values from the mean -dev.from.mean = function(fitted.well){ - measured.ODs = data.from(fitted.well,remove=T,na.rm=T)[,2] - - # Get the mean values of these measured ODs. - mean.ODs = mean(measured.ODs) - - if (!is.numeric(mean.ODs)) - return (NA) - else - return (measured.ODs - mean.ODs) -} - -# Get the residual sum of square. -rss = function(fitted.well){ - if (length(fitted.well@rss) == 0) - return (NA) - else - return (fitted.well@rss) -} - -# Calculate a metric for fit accuracy using squared residuals -model.good.fit = function(fitted.well, digits = Inf){ - # Sum of squared residuals - RSS = rss(fitted.well) - - # Total sum of squared - tot = sum(dev.from.mean(fitted.well)^2) - - # Coefficient of determination - return (1 - RSS/tot) - } - -# Output a string with values of fitted parameters -parameter.text = function(fitted.well){ - # Get a list of fitted parameters - fit.par = fitted.well@fit.par - - # Giving the parameter text descriptive names. - if (length(fitted.well@fit.par) != 0){ - names(fit.par)[1] = "A" - names(fit.par)[2] = "b" - names(fit.par)[3] = "lambda" - names(fit.par)[4] = "max.spec.growth.rate" - - if (fitted.well@model.name == "richards sigmoid"){ - names(fit.par)[5] = "shape.par" - } - - if (fitted.well@model.name == "richards sigmoid with linear par."){ - names(fit.par)[5] = "shape.param" - names(fit.par)[6] = "linear term" - } - - if (fitted.well@model.name == "logistic sigmoid with linear par.") - names(fit.par)[5] = "linear.term" - - # if loess, just show smoothing param - if(fitted.well@model.name == "local polynomial regression fit.") - fit.par = fitted.well@fit.par["smoothing parameter"] - } - - # Return nothing if the list is empty. Otherwise, concatenate the terms in the list with the parameter names. - if(!is.list(fit.par)) - return() - else{ - output = "" - i = 1 - while(i <= length(fit.par)){ - output = paste(output, names(fit.par)[i], "=", round(as.numeric(fit.par[i]),3), "; ", sep = "") - i = i + 1 - if (i %% 6 == 0) - output = paste(output, "\n") - } - output - } - } - -# Calculate maximum specific growth rate -max.spec.growth.rate = function(fitted.well, digits = Inf, ...){ - if(length(fitted.well@fit.par) == 0) - return(NA) - - round(fitted.well@fit.par$u,digits) -} - -# Calculate plateau log.OD from fitted parameters -plateau = function(fitted.well, digits = Inf){ - if(length(fitted.well@fit.par) == 0) - return(NA) - - plat = fitted.well@fit.par$A + fitted.well@fit.par$b - - if (!is.numeric(plat)) { - plat = NA - } else { - plat = round(plat, digits) - } - return(plat) -} - -# Calculate baseline log.OD from fitted parameters -baseline = function(fitted.well, digits = Inf){ - if(length(fitted.well@fit.par) == 0) - return(NA) - - base = fitted.well@fit.par$b - - # If A (plateau OD) is invalid, return NA. - if (!is.numeric(fitted.well@fit.par$A)) - base = NA - # If b (baseline OD) is invalid but plateau OD was valid, return zero. - else if (!is.numeric(base)) - base = 0 - else{ - base = round(base, digits) - } - return(base) - } - -# Calculate log.OD at inoculation from fitted parameters -inoc.log.OD = function(fitted.well, digits = Inf){ - # Evaluated the fitted model at the inoculation timepoint (should be zero from using from table2wells.R) - if (is.null(well.eval(fitted.well))) - return(NA) - else{ - inoc.time = fitted.well@screen.data$Time[fitted.well@start.index] - inoc.log.OD = well.eval(fitted.well, inoc.time) - if (is.na(inoc.log.OD)) inoc.log.OD = fitted.well@fit.par$b # need this in a special case: loess fits with start.index = 1 - return(round(inoc.log.OD, digits)) - } - } - -# Calculate max log.OD from model fit -max.log.OD = function(fitted.well, digits = Inf, ...){ - # Evaluated the fitted model at the final timepoint (just the last valid timepoint in the experiment) - if (is.null(well.eval(fitted.well))) - return(NA) - else{ - return(round(max(well.eval(fitted.well),na.rm=T), digits)) - } -} - -# Calculate projected growth: plateau minus the inoculated log.OD -projected.growth = function(fitted.well,digits=Inf) { - plateau(fitted.well,digits) - inoc.log.OD(fitted.well,digits) -} - - -# Calculate projected growth: plateau minus the inoculated log.OD -projected.growth.OD = function(fitted.well,constant.added,digits=Inf) { - value = unlog(plateau(fitted.well),constant.added) - unlog(inoc.log.OD(fitted.well),constant.added) - round(value,digits) -} - -# Calculate achieved growth: max.log.OD minus the inoculated log.OD -achieved.growth = function(fitted.well,digits=Inf) { - max.log.OD(fitted.well,digits) - inoc.log.OD(fitted.well,digits) -} - -# Calculate projected growth: plateau minus the inoculated log.OD -achieved.growth.OD = function(fitted.well,constant.added,digits=Inf) { - value = unlog(max.log.OD(fitted.well),constant.added) - unlog(inoc.log.OD(fitted.well),constant.added) - round(value,digits) -} - -# Did the curve come close to the plateau OD during the experiment? -reach.plateau = function(fitted.well, cutoff = 0.75){ - plat = plateau(fitted.well) - inoc = inoc.log.OD(fitted.well) - final = max.log.OD(fitted.well) - - if (!is.na(final)){ - # If the plateau is the same as the OD at inoculation, return TRUE - if ((plat - inoc) == 0) - return(T) - # If the difference between the final OD and inoculation OD is at least a certain proportion - # of the difference between the plateau and inoculated ODs, return TRUE. - else - return((final - inoc) / (plat - inoc) > cutoff) - } - else - return(T) - # If no final OD was calculated (if curve was not fit properly) just return T. - } - -# Calculate the lag time from the fitted OD -lag.time = function(fitted.well, digits = Inf, ...){ - if(length(fitted.well@fit.par) == 0) - return(NA) - - fitted.well@fit.par$lam -} - -# new params for GCAT 4.0 -# Get amplitude -amplitude = function(fitted.well){ - if(length(fitted.well@fit.par) == 0) - return(NA) - - return(fitted.well@fit.par$A) -} - -# Get shape parameter -shape.par = function(fitted.well){ - if(length(fitted.well@fit.par) == 0) - return(NA) - ifelse(is.null(fitted.well@fit.par$v), NA, fitted.well@fit.par$v) -} - -# Get standard error of the maximum specific growth rate value -max.spec.growth.rate.SE = function(fitted.well, ...){ - if(length(fitted.well@fit.std.err) == 0) - return(NA) - ifelse(is.null(fitted.well@fit.std.err$u), NA, fitted.well@fit.std.err$u) -} - -# Get standard error of the lag time value -lag.time.SE = function(fitted.well, ...){ - if(length(fitted.well@fit.std.err) == 0) - return(NA) - ifelse(is.null(fitted.well@fit.std.err$lam), NA, fitted.well@fit.std.err$lam) -} - -# Get standard error of the shape parameter -shape.par.SE = function(fitted.well){ - if(length(fitted.well@fit.std.err) == 0) - return(NA) - ifelse(is.null(fitted.well@fit.std.err$v), NA, fitted.well@fit.std.err$v) -} - -# Get standard error of the amplitude -amplitude.SE = function(fitted.well){ - if(length(fitted.well@fit.std.err) == 0) - return(NA) - ifelse(is.null(fitted.well@fit.std.err$A), NA, fitted.well@fit.std.err$A) -} - -# Get standard error of the baseline value -baseline.SE = function(fitted.well){ - if(length(fitted.well@fit.std.err) == 0) - return(NA) - ifelse(is.null(fitted.well@fit.std.err$b), NA, fitted.well@fit.std.err$b) -} - -# Calulate the inflection time value -inflection.time = function(well){ - if (length(well@loess) == 0 && length(well@nls) == 0) return(NA) # can' compute inflection time in the absence of a fit - data = data.from(well) - Time = data[,1] - t = seq(from = min(Time), to = max(Time), by = (max(Time)-min(Time))/1000) - y = well.eval(well,t) - if (is.null(y)) return(NA) - delta.t = diff(t) - dydt = diff(y)/delta.t - infl.index = which.max(dydt) - t[infl.index] -} +#Copyright 2012 The Board of Regents of the University of Wisconsin System. +#Contributors: Jason Shao, James McCurdy, Enhai Xie, Adam G.W. Halstead, +#Michael H. Whitney, Nathan DiPiazza, Trey K. Sato and Yury V. Bukhman +# +#This file is part of GCAT. +# +#GCAT is free software: you can redistribute it and/or modify +#it under the terms of the GNU Lesser General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#GCAT is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU Lesser General Public License for more details. +# +#You should have received a copy of the GNU Lesser General Public License +#along with GCAT. If not, see . + +######################################################################## +# # +# Functions to calculate various things about wells based on fit model # +# # +######################################################################## + +# S3 generic for lag +lag <- function(fitted.well, ...) +{ + UseMethod("lag") +} + +# +# Common arguments: +# fitted.well - should be a well containing the results of , most functions will return NA if well has not been fit yet. +# unlog - should the value be returned on the linear scale as opposed to the log-transformed scale? +# constant.added - for returning values on the linear scale, what was the constant added before the log transform? +# digits - passed to the function, default is no rounding (infinity digits) + +# Transform values back to OD scale +unlog = function(x, constant.added) { + exp(x) - constant.added +} + +# Evaluate estimated OD at any timepoints using the fitted model +well.eval = function(fitted.well, Time = NULL){ + # If no timepoints are provided, use the ones collected in the experiment itself. + if(!is.numeric(Time)) + Time = data.from(fitted.well)$Time + + # Use of equation is deprecated. Use nls and loess models stored in the well object instead + # Attempt to use with the fitted equation and parameters to get estimates for OD at the given timepoints. + #output = try(eval(fitted.well@equation, fitted.well@fit.par), silent = T) + + # Predict log.OD value(s) using nls model if present. If no nls model, try using loess. + if (length(fitted.well@nls)>0) { + output = try(predict(fitted.well@nls,list(Time=Time)),silent=T) + } else if (length(fitted.well@loess)>0) { + output = try(predict(fitted.well@loess,Time),silent=T) + } else { + output = NA + } + + # Return values. If OD evaluation failed for any reason, return NULL. + if (is.numeric(output)){ + return(output) + } else { + return(NULL) + } +} + +# Evaluate model residuals using the measured vs. fitted log.OD values +model.residuals = function(fitted.well, unlog = F){ + measured.OD = data.from(fitted.well)[,2] + + # Use with no Time argument to get fitted OD values at measured timepoints. + predicted.OD = well.eval(fitted.well) + + # If all values are valid, return the differences + if (!is.numeric(predicted.OD)) + return(NA) + else + return(measured.OD - predicted.OD) + } + +# Evaluate deviations of log.OD values from the mean +dev.from.mean = function(fitted.well){ + measured.ODs = data.from(fitted.well,remove=T,na.rm=T)[,2] + + # Get the mean values of these measured ODs. + mean.ODs = mean(measured.ODs) + + if (!is.numeric(mean.ODs)) + return (NA) + else + return (measured.ODs - mean.ODs) +} + +# Get the residual sum of square. +rss = function(fitted.well){ + if (length(fitted.well@rss) == 0) + return (NA) + else + return (fitted.well@rss) +} + +# Calculate a metric for fit accuracy using squared residuals +model.good.fit = function(fitted.well, digits = Inf){ + # Sum of squared residuals + RSS = rss(fitted.well) + + # Total sum of squared + tot = sum(dev.from.mean(fitted.well)^2) + + # Coefficient of determination + return (1 - RSS/tot) + } + +# Output a string with values of fitted parameters +parameter.text = function(fitted.well){ + # Get a list of fitted parameters + fit.par = fitted.well@fit.par + + # Giving the parameter text descriptive names. + if (length(fitted.well@fit.par) != 0){ + names(fit.par)[1] = "A" + names(fit.par)[2] = "b" + names(fit.par)[3] = "lambda" + names(fit.par)[4] = "max.spec.growth.rate" + + if (fitted.well@model.name == "richards sigmoid"){ + names(fit.par)[5] = "shape.par" + } + + if (fitted.well@model.name == "richards sigmoid with linear par."){ + names(fit.par)[5] = "shape.param" + names(fit.par)[6] = "linear term" + } + + if (fitted.well@model.name == "logistic sigmoid with linear par.") + names(fit.par)[5] = "linear.term" + + # if loess, just show smoothing param + if(fitted.well@model.name == "local polynomial regression fit.") + fit.par = fitted.well@fit.par["smoothing parameter"] + } + + # Return nothing if the list is empty. Otherwise, concatenate the terms in the list with the parameter names. + if(!is.list(fit.par)) + return() + else{ + output = "" + i = 1 + while(i <= length(fit.par)){ + output = paste(output, names(fit.par)[i], "=", round(as.numeric(fit.par[i]),3), "; ", sep = "") + i = i + 1 + if (i %% 6 == 0) + output = paste(output, "\n") + } + output + } + } + +# Calculate maximum specific growth rate +max.spec.growth.rate = function(fitted.well, digits = Inf, ...){ + if(length(fitted.well@fit.par) == 0) + return(NA) + + round(fitted.well@fit.par$u,digits) +} + +# Calculate plateau log.OD from fitted parameters +plateau = function(fitted.well, digits = Inf){ + if(length(fitted.well@fit.par) == 0) + return(NA) + + plat = fitted.well@fit.par$A + fitted.well@fit.par$b + + if (!is.numeric(plat)) { + plat = NA + } else { + plat = round(plat, digits) + } + return(plat) +} + +# Calculate baseline log.OD from fitted parameters +baseline = function(fitted.well, digits = Inf){ + if(length(fitted.well@fit.par) == 0) + return(NA) + + base = fitted.well@fit.par$b + + # If A (plateau OD) is invalid, return NA. + if (!is.numeric(fitted.well@fit.par$A)) + base = NA + # If b (baseline OD) is invalid but plateau OD was valid, return zero. + else if (!is.numeric(base)) + base = 0 + else{ + base = round(base, digits) + } + return(base) + } + +# Calculate log.OD at inoculation from fitted parameters +inoc.log.OD = function(fitted.well, digits = Inf){ + # Evaluated the fitted model at the inoculation timepoint (should be zero from using from table2wells.R) + if (is.null(well.eval(fitted.well))) + return(NA) + else{ + inoc.time = fitted.well@screen.data$Time[fitted.well@start.index] + inoc.log.OD = well.eval(fitted.well, inoc.time) + if (is.na(inoc.log.OD)) inoc.log.OD = fitted.well@fit.par$b # need this in a special case: loess fits with start.index = 1 + return(round(inoc.log.OD, digits)) + } + } + +# Calculate max log.OD from model fit +max.log.OD = function(fitted.well, digits = Inf, ...){ + # Evaluated the fitted model at the final timepoint (just the last valid timepoint in the experiment) + if (is.null(well.eval(fitted.well))) + return(NA) + else{ + return(round(max(well.eval(fitted.well),na.rm=T), digits)) + } +} + +# Calculate projected growth: plateau minus the inoculated log.OD +projected.growth = function(fitted.well,digits=Inf) { + plateau(fitted.well,digits) - inoc.log.OD(fitted.well,digits) +} + + +# Calculate projected growth: plateau minus the inoculated log.OD +projected.growth.OD = function(fitted.well,constant.added,digits=Inf) { + value = unlog(plateau(fitted.well),constant.added) - unlog(inoc.log.OD(fitted.well),constant.added) + round(value,digits) +} + +# Calculate achieved growth: max.log.OD minus the inoculated log.OD +achieved.growth = function(fitted.well,digits=Inf) { + max.log.OD(fitted.well,digits) - inoc.log.OD(fitted.well,digits) +} + +# Calculate projected growth: plateau minus the inoculated log.OD +achieved.growth.OD = function(fitted.well,constant.added,digits=Inf) { + value = unlog(max.log.OD(fitted.well),constant.added) - unlog(inoc.log.OD(fitted.well),constant.added) + round(value,digits) +} + +# Did the curve come close to the plateau OD during the experiment? +reach.plateau = function(fitted.well, cutoff = 0.75){ + plat = plateau(fitted.well) + inoc = inoc.log.OD(fitted.well) + final = max.log.OD(fitted.well) + + if (!is.na(final)){ + # If the plateau is the same as the OD at inoculation, return TRUE + if ((plat - inoc) == 0) + return(T) + # If the difference between the final OD and inoculation OD is at least a certain proportion + # of the difference between the plateau and inoculated ODs, return TRUE. + else + return((final - inoc) / (plat - inoc) > cutoff) + } + else + return(T) + # If no final OD was calculated (if curve was not fit properly) just return T. + } + +# Calculate the lag time from the fitted OD +lag.time = function(fitted.well, digits = Inf, ...){ + if(length(fitted.well@fit.par) == 0) + return(NA) + + fitted.well@fit.par$lam +} + +# new params for GCAT 4.0 +# Get amplitude +amplitude = function(fitted.well){ + if(length(fitted.well@fit.par) == 0) + return(NA) + + return(fitted.well@fit.par$A) +} + +# Get shape parameter +shape.par = function(fitted.well){ + if(length(fitted.well@fit.par) == 0) + return(NA) + ifelse(is.null(fitted.well@fit.par$v), NA, fitted.well@fit.par$v) +} + +# Get standard error of the maximum specific growth rate value +max.spec.growth.rate.SE = function(fitted.well, ...){ + if(length(fitted.well@fit.std.err) == 0) + return(NA) + ifelse(is.null(fitted.well@fit.std.err$u), NA, fitted.well@fit.std.err$u) +} + +# Get standard error of the lag time value +lag.time.SE = function(fitted.well, ...){ + if(length(fitted.well@fit.std.err) == 0) + return(NA) + ifelse(is.null(fitted.well@fit.std.err$lam), NA, fitted.well@fit.std.err$lam) +} + +# Get standard error of the shape parameter +shape.par.SE = function(fitted.well){ + if(length(fitted.well@fit.std.err) == 0) + return(NA) + ifelse(is.null(fitted.well@fit.std.err$v), NA, fitted.well@fit.std.err$v) +} + +# Get standard error of the amplitude +amplitude.SE = function(fitted.well){ + if(length(fitted.well@fit.std.err) == 0) + return(NA) + ifelse(is.null(fitted.well@fit.std.err$A), NA, fitted.well@fit.std.err$A) +} + +# Get standard error of the baseline value +baseline.SE = function(fitted.well){ + if(length(fitted.well@fit.std.err) == 0) + return(NA) + ifelse(is.null(fitted.well@fit.std.err$b), NA, fitted.well@fit.std.err$b) +} + +# Calulate the inflection time value +inflection.time = function(well){ + if (length(well@loess) == 0 && length(well@nls) == 0) return(NA) # can' compute inflection time in the absence of a fit + data = data.from(well) + Time = data[,1] + t = seq(from = min(Time), to = max(Time), by = (max(Time)-min(Time))/1000) + y = well.eval(well,t) + if (is.null(y)) return(NA) + delta.t = diff(t) + dydt = diff(y)/delta.t + infl.index = which.max(dydt) + t[infl.index] +} + +# Calculate the area under the curve with the specified time range. +AUC_well <- function(well.array, a = 0, b = 1, silent = TRUE) { + for (i in 1:length(well.array)) { + f <- function(x) {well.eval(well.array[[i]], x)} + auc.string = try(integrate(f, a, b), silent = silent) + if (class(auc.string) == "try-error") + { + well.array[[i]]@auc = "" + if (!silent) print(paste("Cannot calculate the area at well: ", i)) + } + else { + well.array[[i]]@auc = paste(abs(round(auc.string$value, digits = 3)), "with boundaries: (", a, ",", b, ")") + + } + } + return (well.array) +} \ No newline at end of file diff --git a/R/GCAT/R/plot.fit.R b/R/GCAT/R/plot.fit.R index e962f32..03c5868 100644 --- a/R/GCAT/R/plot.fit.R +++ b/R/GCAT/R/plot.fit.R @@ -140,7 +140,7 @@ draw.text = function(input.well, scale = 0.5, xlim = 0, ylim = 0,...){ # color = green if empty, blue if inoculated, red if inoculated but has no growth or empty but has growth. col2 = "blue" - text2 = paste(input.well@fit.info, input.well@model.name, input.well@add.info, "\n", parameter.text(input.well)) + text2 = paste(input.well@fit.info, input.well@model.name, input.well@add.info, "\n", parameter.text(input.well), "\n", "AUC: ", input.well@auc) if (length(input.well@fit.par) == 0) # no fit col2 = "red" @@ -196,7 +196,7 @@ draw.text = function(input.well, scale = 0.5, xlim = 0, ylim = 0,...){ # @details # \strong{show.num} - should curve parameters be labeled? # -draw.calc.par = function(input.well, scale = 0.5, unlog = F, constant.added, show.num = T){ +draw.calc.par = function(input.well, scale = 0.5, unlog = F, constant.added, show.num = T, sranges = NA){ # Don't do anything if well was not fit. if (is.null(well.eval(input.well))) @@ -221,8 +221,16 @@ draw.calc.par = function(input.well, scale = 0.5, unlog = F, constant.added, sho # ---- Specific growth rate ---- # lines(c(lag.x, inflection.time), c(lag.y, max.y), lty = 2, col = "red") - - + + # MB: Drawing boundaries for the AUC. + if (!is.na(sranges)) + { + #abline(v = sranges[1], lty = 1) + segments(sranges[1], baseline, sranges[1], well.eval(input.well, sranges[1]), lty = 1) + segments(sranges[2], baseline, sranges[2], well.eval(input.well, sranges[2]), lty = 1) + #abline(v = sranges[2], lty = 1) + } + # Blue dotted line at time of maximum growth, with text label for specific growth rate. abline(v = inflection.time, lty = 2, lw = (scale^2)*2, col = "blue") if(show.num) text(inflection.time, max.y, round(max.slope,3), col = "blue", cex = 1.5*scale, pos = 2) @@ -485,7 +493,7 @@ plate.overview = function(fitted.well.array, scale = 1, plate.ncol = 12, plate.n view.fit = function(fitted.data, indices = 1:length(fitted.data), unlog = F, constant.added, xlim = NULL, ylim = NULL, display.legend = T, show.text = T, show.calc = T, draw.guess = NULL, draw.symbols = F, number.points = T, - user.advance = T, show.residuals = F, scale = 1,...){ + user.advance = T, show.residuals = F, scale = 1, sranges = sranges,...){ if(!is.array(fitted.data)) fitted.data = list(fitted.data) @@ -530,7 +538,7 @@ view.fit = function(fitted.data, indices = 1:length(fitted.data), # plot the well fitted.well = fitted.data[[well.number]] plot(x=fitted.well, constant.added = constant.added, xlim = xlim, ylim = ylim, - unlog = unlog, well.number = well.number, scale = scale, number.points = T, draw.symbols = F, show.text = T, show.calc = T, draw.guess = NULL, ...) + unlog = unlog, well.number = well.number, scale = scale, number.points = T, draw.symbols = F, show.text = T, show.calc = T, draw.guess = NULL, sranges = sranges, ...) if(user.advance) cat("\n[", well.number, "] ", plate.name(fitted.well), " ", well.name(fitted.well), ".", sep = "") @@ -616,7 +624,7 @@ well.fit.legend = function(xlim, ylim, scale = 1, constant.added){ # Generate pdf files pdf.by.plate = function(fitted.data, out.prefix = "", upload.timestamp = NULL, out.dir = getwd(), unlog = F, constant.added, silent = T, overview.jpgs = T, plate.ncol = 12, plate.nrow = 8, - lagRange = NA, specRange = NA, totalRange = NA, totalODRange = NA, ...){ + lagRange = NA, specRange = NA, totalRange = NA, totalODRange = NA, sranges = sranges,...){ # Prepare timestamp for addition to output file names. filename.timestamp = strftime(upload.timestamp, format="_%Y-%m-%d_%H.%M.%S") @@ -678,7 +686,7 @@ pdf.by.plate = function(fitted.data, out.prefix = "", upload.timestamp = NULL, pdf(pdf.name, title = paste("R Graphics output for plate", plate.ID)) # Call to draw each well on the plate to the pdf. - view.fit.out = try(view.fit(fitted.data, indices = plate.indices, unlog=unlog, constant.added=constant.added, user.advance=F,...),silent=T) + view.fit.out = try(view.fit(fitted.data, indices = plate.indices, unlog=unlog, constant.added=constant.added, user.advance=F, sranges = sranges,...),silent=T) if(class(view.fit.out) == "try-error") stop("Error in : ", view.fit.out) diff --git a/R/GCAT/man/GCAT.Rd b/R/GCAT/man/GCAT.Rd new file mode 100644 index 0000000..1923532 --- /dev/null +++ b/R/GCAT/man/GCAT.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2 (4.1.1): do not edit by hand +% Please edit documentation in R/GCAT.main.R +\docType{package} +\name{GCAT} +\alias{GCAT} +\alias{GCAT-package} +\title{GCAT: Growth Curve Analysis Tool} +\description{ +Mathematical modeling and parameter estimation of high volume microbial growth data. +} +\details{ +GCAT input is in .csv format. GCAT analysis is accessed using \code{\link{gcat.analysis.main}} + + GCAT utilizes the \code{\link[stats]{nls}} function in the R stats package to fit logistic, Gompertz and Richards models to growth curve + data. Best model is selected automatically. Alternatively, the user may choose LOESS local regression fits, implemented using + \code{\link[stats]{loess}} function in the R stats package + + Internally, the data are stored in an array of \linkS4class{well} objects +} + diff --git a/R/GCAT/man/gcat.analysis.main.Rd b/R/GCAT/man/gcat.analysis.main.Rd index 80699c5..8da43bb 100644 --- a/R/GCAT/man/gcat.analysis.main.Rd +++ b/R/GCAT/man/gcat.analysis.main.Rd @@ -10,7 +10,7 @@ gcat.analysis.main(file.list, single.plate, layout.file = NULL, use.linear.param = F, use.loess = F, smooth.param = 0.1, lagRange = NA, totalRange = NA, totalODRange = NA, specRange = NA, points.to.remove = 0, remove.jumps = F, time.input = NA, - plate.nrow = 8, plate.ncol = 12, input.skip.lines = 0, + plate.nrow = 8, plate.ncol = 12, sranges = NA, input.skip.lines = 0, multi.column.headers = c("Plate.ID", "Well", "OD", "Time"), single.column.headers = c("", "A1"), layout.sheet.headers = c("Strain", "Media Definition"), silent = T, verbose = F, return.fit = F, @@ -59,6 +59,8 @@ gcat.analysis.main(file.list, single.plate, layout.file = NULL, \item{plate.ncol}{The number of columns in a plate.} +\item{sranges}{The 2 boundaries for calculating the area under the curve.} + \item{input.skip.lines}{If specified, this number of lines shall be skipped from the top when reading the input file with read.csv} \item{multi.column.headers}{The headers of the result tabular data when analyzing multiple plates at once.} diff --git a/R/GCAT/man/gcat.fit.main.Rd b/R/GCAT/man/gcat.fit.main.Rd index 2aec21f..e902c14 100644 --- a/R/GCAT/man/gcat.fit.main.Rd +++ b/R/GCAT/man/gcat.fit.main.Rd @@ -10,7 +10,7 @@ gcat.fit.main(file.name, input.data = NULL, load.type = "csv", add.constant = 1, use.log = T, points.to.remove = 0, use.linear.param = F, use.loess = F, smooth.param = 0.1, fall.cutoff = -0.0025, growth.cutoff = 0.05, remove.jumps = F, - plate.nrow = 8, plate.ncol = 12, input.skip.lines = 0, + plate.nrow = 8, plate.ncol = 12, sranges = NA, input.skip.lines = 0, multi.column.headers = c("Plate.ID", "Well", "OD", "Time"), single.column.headers = c("", "A1"), layout.sheet.headers = c("Strain", "Media Definition"), growth.model = NA, backup.growth.model = NA, @@ -63,6 +63,8 @@ By default(0) none are marked for removal.} \item{plate.ncol}{The number of columns in the input files.} +\item{sranges}{The 2 boundaries for calculating the area under the curve.} + \item{input.skip.lines}{If specified, this number of lines shall be skipped from the top when reading the input file with read.csv} \item{multi.column.headers}{The headers of the column when analyzing multiple plates.} diff --git a/R/GCAT/man/gcat.output.main.Rd b/R/GCAT/man/gcat.output.main.Rd index 96863d3..6774b00 100644 --- a/R/GCAT/man/gcat.output.main.Rd +++ b/R/GCAT/man/gcat.output.main.Rd @@ -8,9 +8,9 @@ gcat.output.main(fitted.well.array, out.prefix = "", source.file.list, upload.timestamp = NULL, add.constant, blank.value, start.index, growth.cutoff, points.to.remove, remove.jumps, out.dir = getwd(), graphic.dir = paste(out.dir, "/pics", sep = ""), overview.jpgs = T, - use.linear.param = F, use.loess = F, lagRange = NA, totalRange = NA, - totalODRange = NA, specRange = NA, plate.nrow = 8, plate.ncol = 12, - unlog = F, silent = T, main.envir) + sranges = NA, use.linear.param = F, use.loess = F, lagRange = NA, + totalRange = NA, totalODRange = NA, specRange = NA, plate.nrow = 8, + plate.ncol = 12, unlog = F, silent = T, main.envir) } \arguments{ \item{fitted.well.array}{A list of fitted well objects.} @@ -42,6 +42,8 @@ If NULL, defaults to the value of the first OD measurement of each well.} \item{overview.jpgs}{should jpgs be generated for each plate with the overview graphic? This is for backwards compatibility with the old web server.} +\item{sranges}{The 2 boundaries for calculating the area under the curve.} + \item{use.linear.param}{linear parameter is used or not?} \item{use.loess}{Is LOESS model going to be used?} diff --git a/R/GCAT/man/plot-well-missing-method.Rd b/R/GCAT/man/plot-well-missing-method.Rd index 6d7ec5a..533364a 100644 --- a/R/GCAT/man/plot-well-missing-method.Rd +++ b/R/GCAT/man/plot-well-missing-method.Rd @@ -8,7 +8,7 @@ \S4method{plot}{well,missing}(x, y, constant.added = 1, xlim = NULL, ylim = NULL, well.number = NULL, scale = 1, number.points = T, draw.symbols = F, show.text = T, show.calc = T, draw.guess = NULL, - ...) + sranges = NA, ...) } \arguments{ \item{x}{object of class well} @@ -35,6 +35,8 @@ \item{draw.guess}{initial guess model. Drawn if specified} +\item{sranges}{The boundaries to calculate the AUC.} + \item{...}{additional arguments passed to the generic plot function} } \description{ diff --git a/R/GCAT/man/well-class.Rd b/R/GCAT/man/well-class.Rd index 9572265..691a7cf 100644 --- a/R/GCAT/man/well-class.Rd +++ b/R/GCAT/man/well-class.Rd @@ -44,5 +44,7 @@ Class that contains well data \item{\code{loess}}{- object returned by running loess on the normalized well data} \item{\code{nls}}{- object returned by running nls on the normalized well data} + +\item{\code{auc}}{- The area under the curve.} }} diff --git a/R/GCAT/tests/YPDAFEXglucoseTests_2-25-10__defaults.R b/R/GCAT/tests/YPDAFEXglucoseTests_2-25-10__defaults.R index 21b2cfb..e146228 100755 --- a/R/GCAT/tests/YPDAFEXglucoseTests_2-25-10__defaults.R +++ b/R/GCAT/tests/YPDAFEXglucoseTests_2-25-10__defaults.R @@ -13,7 +13,7 @@ time.input=1/3600 out = gcat.analysis.main(file.list = INPUT.FILE, single.plate = T, layout.file = NULL, out.dir = OUTPUT.DIR, graphic.dir = OUTPUT.DIR, add.constant = 1, blank.value = NULL, start.index = 2, growth.cutoff = 0.05, - use.linear.param=F, use.loess=F, smooth.param=0.1, + use.linear.param=F, use.loess=F, smooth.param=0.1, sranges = c(5, 10), points.to.remove = integer(), remove.jumps = F, time.input=time.input, silent = F, verbose = T, return.fit = T, overview.jpgs = T) diff --git a/Rails/app/assets/stylesheets/gcat.css.scss b/Rails/app/assets/stylesheets/gcat.css.scss index 83a522d..4242abf 100644 --- a/Rails/app/assets/stylesheets/gcat.css.scss +++ b/Rails/app/assets/stylesheets/gcat.css.scss @@ -81,7 +81,7 @@ a:visited { color: #428bca !important;} /* MB: 1 option in GCAT. */ #options li { float: left; - width: 33%; + width: 45%; padding: 0 100px 0 0; } diff --git a/Rails/app/models/assay.rb b/Rails/app/models/assay.rb index 3802148..8c3c249 100644 --- a/Rails/app/models/assay.rb +++ b/Rails/app/models/assay.rb @@ -34,7 +34,7 @@ class Assay extend ActiveModel::Naming attr_accessor :input_file, :blank_value, :blank_value_input, :start_index, :remove_points, :remove_jumps, :plate_type, :plate_dimensions_row, :plate_dimensions_column, :timestamp_format, :growth_threshold, :layout_file,:filename,:content_type, :model, :loess_input, :console_out, :specg_min, - :specg_max, :totg_min, :totg_max, :totg_OD_min, :totg_OD_max, :lagT_min, :lagT_max,:transformation, :transformation_input + :specg_max, :totg_min, :totg_max, :totg_OD_min, :totg_OD_max, :lagT_min, :lagT_max,:transformation, :transformation_input, :range_a, :range_b # (1) Validation of input data file @@ -189,6 +189,12 @@ class Assay self.lagT_min = Float(self.lagT_min) end + # (8) Integration bound + if (self.range_a != '' && self.range_a != '') + self.range_a = Float(self.range_a) + self.range_b = Float(self.range_b) + end + ############################################################################################ @@ -380,12 +386,21 @@ class Assay else R.eval 'lagRange <- NA' end - + + # Integration bound. + if (self.range_a != '' && self.range_b != '') + R.assign 'range_a', self.range_a + R.assign 'range_b', self.range_b + R.eval "sranges <- c(range_a, range_b)" + else + R.eval "sranges <- NA" + end + # This block evaluates the files (csv or xlsx, single.plate or multiple.plate) R.eval ('R_file_return_value <- gcat.analysis.main(file, single.plate, layout.file, out.dir=out.dir, graphic.dir = out.dir, add.constant, blank.value, start.index, growth.cutoff, use.linear.param=use.linear.param, use.loess=use.loess, smooth.param=smooth.param, - lagRange = lagRange, totalRange = totalRange, totalODRange = totalODRange, specRange = specRange, - points.to.remove = points.to.remove, remove.jumps, time.input, plate.nrow = 8, + lagRange = lagRange, totalRange = totalRange, totalODRange = totalODRange, specRange = specRange, + points.to.remove = points.to.remove, remove.jumps, time.input, plate.nrow = 8, sranges = sranges, plate.ncol = 12, input.skip.lines = 0, multi.column.headers = c("Plate.ID", "Well", "OD", "Time"), single.column.headers = c("","A1"), layout.sheet.headers = c("Strain", "Media Definition"), silent = T, verbose = F, return.fit = F, overview.jpgs = T)') diff --git a/Rails/app/views/assays/_form.html.erb b/Rails/app/views/assays/_form.html.erb index 18d0271..ee94c0d 100644 --- a/Rails/app/views/assays/_form.html.erb +++ b/Rails/app/views/assays/_form.html.erb @@ -154,6 +154,15 @@ <%= f.text_field(:lagT_max, :size => '5') %> +
  • +

    Integration ranges

    +
    + <%= f.label "Enter the time range" %>
    + <%= f.text_field :range_a, :size => '5' %> + <%= f.label "-" %> + <%= f.text_field :range_b, :size => '5' %> +
    +
    • <%= f.submit %>