minhdb/2015-GCAT
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Demo AUC

Browse Source
This commit is contained in:
Minh
2015-07-15 02:49:57 -05:00
parent 30e0658245
commit 329c63caef
14 changed files with 468 additions and 372 deletions
Download Patch File Download Diff File Expand all files Collapse all files
R/GCAT/R/GCAT.main.R
+20 -7
View File
@@ -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 <fit.model>: ", 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 <pdf.by.plate>: ", graphic.files)
R/GCAT/R/class.well.R
+6 -3
View File
@@ -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"){
R/GCAT/R/fitted.calculations.R
+360 -342
View File
@@ -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 <http://www.gnu.org/licenses/>.
########################################################################
# #
# 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 <fit.model>, 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 <round> 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 <eval> 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 <well.eval> 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 <start.times> 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
# <cutoff> 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 <http://www.gnu.org/licenses/>.
########################################################################
# #
# 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 <fit.model>, 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 <round> 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 <eval> 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 <well.eval> 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 <start.times> 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
# <cutoff> 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)
}
R/GCAT/R/plot.fit.R
+16 -8
View File
@@ -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 <view.fit> 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>: ", view.fit.out)
R/GCAT/man/GCAT.Rd
+20
View File
@@ -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
}
R/GCAT/man/gcat.analysis.main.Rd
+3 -1
View File
@@ -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.}
R/GCAT/man/gcat.fit.main.Rd
+3 -1
View File
@@ -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.}
R/GCAT/man/gcat.output.main.Rd
+5 -3
View File
@@ -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?}
R/GCAT/man/plot-well-missing-method.Rd
+3 -1
View File
@@ -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{
R/GCAT/man/well-class.Rd
+2
View File
@@ -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.}
}}
R/GCAT/tests/YPDAFEXglucoseTests_2-25-10__defaults.R
+1 -1
View File
@@ -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)
Rails/app/assets/stylesheets/gcat.css.scss
+1 -1
View File
@@ -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;
}
Rails/app/models/assay.rb
+19 -4
View File
@@ -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)')
Rails/app/views/assays/_form.html.erb
+9
View File
@@ -154,6 +154,15 @@
<%= f.text_field(:lagT_max, :size => '5') %>
</div>
</li>
<li>
<h3><span data-tooltip="Select value ranges to calculate the area under the curve">Integration ranges</span></h3>
<div class="field">
<%= f.label "Enter the time range" %><br />
<%= f.text_field :range_a, :size => '5' %>
<%= f.label "-" %>
<%= f.text_field :range_b, :size => '5' %>
</div>
</li>
<div class="actions" style="margin-top:30px; clear:both;">
<%= f.submit %>
</div>