Yury V Bukhman
643 lines
38 KiB
R
643 lines
38 KiB
R
#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, Minh Bui, 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/>.
|
|
|
|
########################################################################
|
|
# #
|
|
# Package-level documentation and imports #
|
|
# #
|
|
########################################################################
|
|
#' GCAT: Growth Curve Analysis Tool
|
|
#'
|
|
#' 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
|
|
#'
|
|
#' @import pheatmap gplots methods
|
|
#' @docType package
|
|
#' @name GCAT
|
|
NULL
|
|
|
|
# Initialization
|
|
PLATE.LETTERS = paste(rep(c("", LETTERS), each = 26), rep(LETTERS, 26), sep="")
|
|
global.version.number = packageDescription(pkg="GCAT")$Version
|
|
|
|
########################################################################
|
|
# #
|
|
# Top-level functions for analysis of screening data from .csv files. #
|
|
# #
|
|
########################################################################
|
|
# This functions is called directly by the user interface.
|
|
# They in turn call the main function <gcat.fit.main> (below) multiple times for each data file provided in <file.list>.
|
|
|
|
# Arguments:
|
|
# file.list - a list of full paths to .csv files. all files must be in the same format (see <single.plate>)
|
|
# single.plate - are the file in the single plate (wide) format vs. the multi-plate (long) format?
|
|
# layout.file - (optional) provide full path to a layout file with strain and media definitions (applies to all files in list)
|
|
|
|
# out.dir - name a directory to output the table of curve parameters to (defaults to working directory)
|
|
# graphic.dir - name a directory to output the images of the fitted curves to (defaults to subdirectory "pics" of <out.dir> above)
|
|
|
|
|
|
# add.constant- should be a numeric constant that will be added to each curve before the log transform (defaults to 1)
|
|
# blank.value - user can enter a blank OD measurement for uninoculated wells. if NULL, defaults to the value of the first OD measurement of each well.
|
|
# start.index - which timepoint should be used as the first one after inoculation (defaults to the 2th one)
|
|
# growth.cutoff - minimum threshold for curve growth.
|
|
# points.to.remove - a list of numbers referring to troublesome points that should be removed across all wells.
|
|
# remove.jumps - should the slope checking function be on the lookout for large jumps in OD?
|
|
|
|
# silent - should messages be returned to the console?
|
|
# verbose - should sub-functions return messages to console? (when I say verbose, I mean it!)
|
|
|
|
# Returns:
|
|
# if <return.fit> = F (default), avector of full paths to all the files generated by the function.
|
|
# otherwise, the fitted array of well objects.
|
|
|
|
# Use this function to analyze any set of .csv files using the same plate layout info.
|
|
|
|
#' Analyze screening growth data from the given .csv files.
|
|
#'
|
|
#' Top-level GCAT function
|
|
#'
|
|
#' @param file.list A list of full paths to .csv files. all files must be in the same format (see <single.plate>)
|
|
#' @param single.plate The file in the single plate (wide) format vs. the multi-plate (long) format?
|
|
#' @param layout.file Full path to a layout file with strain and media definitions (applies to all files in list)
|
|
#' @param out.dir A directory to output the table of curve parameters to (defaults to working directory)
|
|
#' @param graphic.dir A directory to output the images of the fitted curves to (defaults to subdirectory "pics" of <out.dir> above)
|
|
#' @param use.linear.param Whether to use linear parameters or not?
|
|
#' @param use.loess Whether to use LOESS model or not?
|
|
#' @param smooth.param Smoothing parameter for LOESS model.
|
|
#' @param add.constant A numeric constant that will be added to each curve before the log transform (defaults to 1)
|
|
#' @param blank.value User can enter a blank OD measurement for uninoculated wells. if NULL, assumes the first OD measurement of each well to be blank.
|
|
#' @param start.index Which timepoint should be used as the first one after inoculation (defaults to the 2th one)
|
|
#' @param growth.cutoff Minimum threshold for curve growth.
|
|
#' @param points.to.remove A list of numbers referring to troublesome points that should be removed across all wells.
|
|
#' @param remove.jumps Should the slope checking function be on the lookout for large jumps in OD?
|
|
#' @param time.input The time setting in which the current system is running?
|
|
#' @param plate.nrow The number of rows in a plate.
|
|
#' @param plate.ncol The number of columns in a plate.
|
|
#' @param input.skip.lines If specified, this number of lines shall be skipped from the top when reading the input file with read.csv
|
|
#' @param multi.column.headers The headers of the result tabular data when analyzing multiple plates at once.
|
|
#' @param single.column.headers The headers of the result tebaular data when analyzaing a single plate.
|
|
#' @param layout.sheet.headers The headers of the layout file?
|
|
#' @param silent Shoulde messages be returned to the console?
|
|
#' @param verbose Should sub-functions return messages to console? (when I say verbose, I mean it!)
|
|
#' @param overview.jpgs Should GCAT enable an overview image?
|
|
#' @param return.fit Whether should a fit well object is returned or not.
|
|
#' @param lagRange The heatmap specific range for lag time.
|
|
#' @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.
|
|
#'
|
|
#' @return Depending on return.fit setting, an array of fitted well objects or a list of output files
|
|
#'
|
|
#' @export
|
|
gcat.analysis.main = function(file.list, single.plate, layout.file = NULL,
|
|
out.dir = getwd(), graphic.dir = paste(out.dir, "/pics", sep = ""),
|
|
add.constant = 0, blank.value, start.index, growth.cutoff = 0.05,
|
|
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,
|
|
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){
|
|
|
|
# Capture the starting environment for debugging
|
|
main.envir = c(as.list(environment()))
|
|
|
|
# Check blank value and start index
|
|
if (is.null(blank.value) && start.index==1) {
|
|
exception("", "If inoculation time point is 1, the user must specify a blank value")
|
|
}
|
|
|
|
# MB: Not the best solution.
|
|
if (is.na(time.input)) {
|
|
if (single.plate)
|
|
time.input = 1/3600
|
|
else
|
|
exception("", "time.input is NA.")
|
|
}
|
|
|
|
# MB: Now add.constant will always be 0.
|
|
# No need to check.
|
|
#if (add.constant < 0)
|
|
# exception("", "The constant r should not be negative.")
|
|
# End prototyping temporary solution.
|
|
|
|
# YB: seem to need this to avoid spurious discarding of some wells in example multiplate dataset: Trac ticket 1780.
|
|
# however, this causes an error in Rails
|
|
if(length(points.to.remove)==0) points.to.remove = 0
|
|
|
|
upload.timestamp = strftime(Sys.time(), format="%Y-%m-%d %H:%M:%S") # Get a timestamp for the time of upload.
|
|
fitted.well.array.master = list()
|
|
source.file.list = c()
|
|
|
|
dim(fitted.well.array.master) = c(plate.nrow,plate.ncol,0)
|
|
dimnames(fitted.well.array.master) = list(PLATE.LETTERS[1:plate.nrow], 1:plate.ncol, c())
|
|
|
|
for(file.name in file.list){
|
|
|
|
# Call <gcat.fit.main> on the file with single plate options
|
|
fitted.well.array = try(gcat.fit.main(file.name = file.name, load.type = "csv",
|
|
single.plate = single.plate, layout.file = layout.file, start.index = start.index,
|
|
time.input = time.input, add.constant = add.constant, blank.value = blank.value,
|
|
growth.cutoff = growth.cutoff, points.to.remove = points.to.remove, remove.jumps = remove.jumps,
|
|
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)
|
|
|
|
# Return error message if the function fails.
|
|
if(class(fitted.well.array) == "try-error")
|
|
return(as.character(fitted.well.array))
|
|
}
|
|
|
|
# Add fitted well array onto existing fitted wells
|
|
fitted.well.array.master = gcat.append.arrays(fitted.well.array.master, fitted.well.array, plate.ncol, plate.nrow)
|
|
|
|
# Remove the "processed_" tag from file names and add to the list of source files.
|
|
source.file.list = c(source.file.list, basename(paste(strsplit(file.name, "processed_")[[1]],collapse="/")))
|
|
|
|
out.files = try(gcat.output.main(fitted.well.array.master, out.prefix = "output",
|
|
source.file.list = source.file.list, upload.timestamp = upload.timestamp,
|
|
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,
|
|
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)
|
|
|
|
# Return file list or error message otherwise return "successful analysis" message (?)
|
|
|
|
# file.list = c("Data was successfully analyzed.", file.list) # <--- yet to be implemented. causes errors downstream right now
|
|
if(class(out.files) == "try-error") return(as.character(out.files))
|
|
|
|
if(return.fit) return(fitted.well.array.master)
|
|
else return(out.files)
|
|
}
|
|
|
|
########################################################################
|
|
# #
|
|
# Main function for analysis of screening data from input tables. #
|
|
# #
|
|
########################################################################
|
|
# This is the main function that handles all the analyses for files in both single and multiple plate formats.
|
|
# It is called by the top level function <analysis.main>
|
|
#
|
|
# It then calls the following functions on each member of the array:
|
|
# - curve normalization and standardization: <gcat.start.times>, <remove.points>, <normalize.ODs>, <transform.ODs>,
|
|
# - curve shape analysis before model fitting: <fill.slopes>, <check.curve>, <classify.curve>
|
|
# - to fit a nonlinear model to the growth data: <fit.model>
|
|
# Finally, it returns the fitted array of well objects.
|
|
|
|
#' Main analysis function for GCAT
|
|
#'
|
|
#' This is the main function that handles all the analyses for data files in both single and multiple plate formats.
|
|
#' It is called by the top level function \code{gcat.analysis.main} along with \code{gcat.output.main}.
|
|
#'
|
|
#' @param file.name Complete path and file name of a comma-separated values (.csv) file containing growth curve data
|
|
#' in the multiple-plate (long) format.
|
|
#' @param input.data A list of tables representing input files read with \code{read.table}. Used to save time in cases
|
|
#' of running multiple analyses on the same dataset. If used, the function will ignore \code{file.name} entirely.
|
|
#' @param load.type .csv by default.
|
|
#' @param layout.file Specifies the location of a layout file containing identifying information.
|
|
#' @param single.plate Whether the GCAT is analyzing a single plate or not.
|
|
#' @param blank.value Blank OD measurement for uninoculated wells. By default(NULL), the value of the first OD
|
|
#'measurement in each well is used.
|
|
#' @param start.index Which timepoint should be used as the first one after inoculation?
|
|
#' @param time.input Either a character describing the format used to convert timestamps in the input to numbers
|
|
#' representing number of seconds (see \code{strptime}), or a factor to divide entries in the Time column by to get the
|
|
#' numbers of hours.
|
|
#' @param normalize.method Describes the method used by \code{normalize.ODs} to normalize cell density values using blank reads.
|
|
#' @param add.constant A value for r in the log(OD + r) transformation.
|
|
#' @param use.log Should the analysis use log on all values.
|
|
#' @param points.to.remove A vector of integers specifying which timepoints should be removed across all wells.
|
|
#' By default(0) none are marked for removal.
|
|
#' @param use.linear.param Should the linear parameter be used or not.
|
|
#' @param use.loess Should the loess model be used or not.
|
|
#' @param smooth.param If loess model is used, this parameter define the smoothing parameter for the loess model.
|
|
#' @param fall.cutoff A cutoff used by \code{check.slopes} to decide on thresholds for jumps and tanking.
|
|
#' @param growth.cutoff A threshold used by check.growth to decide whether a well displays growth.
|
|
#' @param remove.jumps Should jumps in OD detected by the subfunction \code{check.slopes}?
|
|
#' @param plate.nrow The number of rows in the input files.
|
|
#' @param plate.ncol The number of columns in the input files.
|
|
#' @param input.skip.lines If specified, this number of lines shall be skipped from the top when reading the input file with read.csv
|
|
#' @param multi.column.headers The headers of the column when analyzing multiple plates.
|
|
#' @param single.column.headers The headers of the column when analyzing a single plate.
|
|
#' @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 silent Surpress all messages.
|
|
#' @param verbose Display all messages when analyzing each well.
|
|
#'
|
|
#' @return An array of well objects
|
|
gcat.fit.main = function(file.name, input.data = NULL, load.type = "csv", layout.file = NULL,
|
|
single.plate = F, blank.value = NULL, start.index = 2, time.input = NA,
|
|
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,
|
|
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,
|
|
silent = F, verbose = F){
|
|
|
|
# Explanation of arguments:
|
|
|
|
# ---File Handling---
|
|
# file.name - full path to an excel spreadsheet, .csv or tab-delimited text file, in either the single or multiple-plate format
|
|
# input.data - use pre-loaded data set (output from <gcat.load.data> function only). will override <file.name> if not NULL
|
|
# load.type - supports "csv."
|
|
# layout.file - full path to a file containing the plate layout in the same format as <file.name>. will not be used if <load.type> is "xlsx"
|
|
|
|
# ---Input file format---
|
|
# single.plate - true denotes data in single-plate format, i.e. simple OD output. false denotes multiple-plate robotic screening output.
|
|
# note: reading directly from excel to R results in timestamps being converted to days.
|
|
|
|
# ---Normalization and Transforms---
|
|
# blank.value - user can enter a blank OD measurement for uninoculated wells. if NULL, defaults to the value of the first OD measurement of each well.
|
|
# start.index - which timepoint should be used as the first one after inoculation (defaults to the 2th one)
|
|
# normalize.method - how should each growth curve be normalized? allowed values are:
|
|
# "first": subtracts the first OD, assumed to be the blank, from all ODs
|
|
# "none": does nothing, assumes no blank. highly recommend log(OD+1) transform in this case.
|
|
# "average.first": forces all filled wells on each plate to match the average value at <start.index> (after subtracting the first OD)
|
|
# add.constant - a numeric constant that will be added to each curve before the log transform (default is 0)
|
|
# use.log - should a log transform be applied to the data after normalization?
|
|
# points.to.remove - a list of numbers referring to troublesome points that should be removed across all wells.
|
|
|
|
# ---Pre-fitting processing---
|
|
# fall.cutoff - a cutoff value for determining whether OD falls significantly between two timepoints. see <check.slopes> in prefit.processing.R for details.
|
|
# growth.cutoff - a cutoff value for determining whether a well contains a successfully growing culture or not.
|
|
# remove.jumps - should the slope checking function be on the lookout for large jumps in OD?
|
|
|
|
# ---Model fitting---
|
|
# model - which parametrized growth model to use? can be richards, gompertz, or logistic. models are defined as objects of class model, see "model.class.R"
|
|
# backup.model - which model should be used if fitting using <model> fails? should ideally be simpler than the main model (less parameters)
|
|
|
|
# ---Miscellanous input/output preferences---
|
|
# silent - should messages be returned to the console?
|
|
# verbose - should sub-functions return messages to console? (when I say verbose, I mean it!)
|
|
# unlog - should exported graphics be transformed back to the OD scale?
|
|
# return.fit - should the function return an array of wells? if not, it will return a list of generated files.
|
|
|
|
|
|
########################################################################
|
|
# Read from .csv file #
|
|
########################################################################
|
|
#
|
|
# The functions used here are found in table2well.R
|
|
|
|
if(!silent) cat("\nReading input files...")
|
|
# Read from .csv or tab-delimited text file using <gcat.load.data> (in load.R)
|
|
# if <layout.file> is provided, it will be used here.
|
|
|
|
plate.layout = NULL
|
|
# Read layout file if it is specified.
|
|
if(!is.null(layout.file)){
|
|
if(load.type=="csv") plate.layout = read.csv(layout.file,header=T,stringsAsFactors=F)
|
|
else plate.layout = read.table(layout.file,header=T,sep="\t",stringsAsFactors=F)
|
|
if(!silent) cat("\n\tAdded plate layout information from", layout.file, "\n")
|
|
}
|
|
|
|
# Load the data
|
|
well.array = try(gcat.load.data(file.name = file.name, input.data = input.data,
|
|
plate.layout = plate.layout, plate.nrow = plate.nrow, plate.ncol = plate.ncol,
|
|
input.skip.lines = input.skip.lines, multi.column.headers = multi.column.headers,
|
|
single.column.headers = single.column.headers, layout.sheet.headers = layout.sheet.headers,
|
|
blank.value = blank.value, start.index = start.index, single.plate = single.plate,
|
|
load.type = load.type, silent=silent),silent=silent)
|
|
|
|
# Return an error if there is a problem with file loading.
|
|
if (class(well.array) == "try-error")
|
|
stop("Error in <gcat.load.data>: ", well.array)
|
|
|
|
|
|
# !---At this point, <well.array> is an array of well objects, each containing raw data and media/strain information if provided---
|
|
|
|
# Attempt to apply time formatting to all wells in array
|
|
well.array = try(aapply(well.array, gcat.start.times, start.index = start.index, time.input = time.input),silent=silent)
|
|
|
|
# Return an error if there is a problem with time formatting
|
|
if (class(well.array) == "try-error")
|
|
stop("Error in <gcat.start.times>: ", well.array)
|
|
|
|
########################################################################
|
|
# Perform normalization and transformation of raw data #
|
|
########################################################################
|
|
#
|
|
# The functions used here are found in normalize.and.transform.R
|
|
|
|
if(!silent) cat("\nProcessing raw data...")
|
|
|
|
# Set all timepoints to active for now using "points.to.remove=0" argument with <remove.points>
|
|
# adds an extra column to the "well.array" slot of each well specifying which points to remove when data is retrieved from the well
|
|
well.array = aapply(well.array, remove.points, points = 0)
|
|
|
|
# Normalize ODs using specified method and adding a constant if desired.
|
|
# sets the "norm" slot of each well to a value to be subtracted from OD values whenever data is retrieved from the well
|
|
well.array = try(normalize.ODs(well.array, normalize.method = normalize.method,
|
|
start.index = start.index, blank.value = blank.value, add.constant = add.constant),silent=silent)
|
|
|
|
# Return an error if there is a problem with normalization
|
|
if (class(well.array) == "try-error")
|
|
stop("Error in <normalize.ODs>: ", well.array)
|
|
#well.array = try(subtract.blank(well.array, blank.value), silent = silent)
|
|
|
|
# Transform ODs on the logarithmic scale, regardless of whether <use.log> is true
|
|
# an extra column of log-transformed values is added to the "well.array" slot of each well
|
|
# the "use.log" slot of each well is set instead to determine whether the transformed values will be returned when data is retrieved from the well.
|
|
well.array = try(aapply(well.array, transform.ODs, start.index = start.index, blank.value = blank.value, use.log = use.log, constant.added = add.constant),silent=silent)
|
|
|
|
# Return an error if there is a problem with transformation
|
|
if (class(well.array) == "try-error")
|
|
stop("Error in <transform.ODs>: ", well.array)
|
|
|
|
# Remove specified timepoints across wells (use "points.to.remove=NULL" if no points to remove)
|
|
well.array = try(aapply(well.array, remove.points, points = points.to.remove),silent=silent)
|
|
|
|
# Return an error if there is a problem with point removal
|
|
if (class(well.array) == "try-error")
|
|
stop("Error in <remove.points>: ", well.array)
|
|
|
|
|
|
|
|
########################################################################
|
|
# Pre-fitting data processing (analysis of curve shapes) #
|
|
########################################################################
|
|
#
|
|
# The functions used here are found in slope.analysis.R
|
|
|
|
# Estimate slope at each timepoint
|
|
# add a column to the "well.array" slot of each well with the local slope at each timepoint
|
|
|
|
well.array = try(aapply(well.array, calculate.slopes, silent=!verbose),silent=silent)
|
|
|
|
# Return an error if there is a problem with slope calculation
|
|
if (class(well.array) == "try-error")
|
|
stop("Error in <calculate.slopes>: ", well.array)
|
|
|
|
# Check slopes for tanking and/or jumping behavior
|
|
# fills the "curve.par" slot of each well with <tanking.start>, denoting the timepoint at which tanking starts (if none, value is NA)
|
|
# uses <remove.points> to remove all points after <tanking.start>
|
|
# It will also fill the "jump.error" slot with a status message, and try to use an automated process to remove the
|
|
# erroneous points if <remove.jumps> is true (default false).
|
|
|
|
well.array = try(aapply(well.array, check.slopes, fall.cutoff = fall.cutoff, remove.jumps = remove.jumps, silent=!verbose, draw = F),silent=silent)
|
|
|
|
# Return an error if there is a problem with slope analysis
|
|
if (class(well.array) == "try-error")
|
|
stop("Error in <check.slopes>: ", well.array)
|
|
|
|
|
|
# Check curves for growth above cutoff
|
|
# fills the "curve.par" slot of each well with <no.growth>, denoting whether the well has no detectable growth.
|
|
well.array = try(aapply(well.array, check.growth, growth.cutoff = growth.cutoff, start.index = start.index),silent=silent)
|
|
|
|
# Return an error if there is a problem with growth.check
|
|
if (class(well.array) == "try-error")
|
|
stop("Error in <check.growth>: ", well.array)
|
|
|
|
########################################################################
|
|
# Fit parameterized models to data #
|
|
########################################################################
|
|
#
|
|
# The functions used here are found in fit.model.R
|
|
|
|
# Fit each well with the selected model and attempt to catch failed fittings with the backup model
|
|
# skips wells designated as <no.growth> above
|
|
# fills the "fit.info" slot of each well with "success," "failed," or "skipped"
|
|
# if fit was successful:
|
|
# fills the "equation" and "model.name" slots with the relevant info for the successful model
|
|
# fills the "fit.par" slot with fitted parameters if fit is successful
|
|
|
|
if(!silent) cat("\nFitting models to data...")
|
|
well.array = aapply(well.array, fit.model, growth.model=growth.model,
|
|
backup.growth.model = backup.growth.model, use.linear.param=use.linear.param,
|
|
use.loess=use.loess, smooth.param=smooth.param, silent=!verbose)
|
|
|
|
# Return an error if there is a problem with model fitting
|
|
if (class(well.array) == "try-error")
|
|
stop("Error in <fit.model>: ", well.array)
|
|
|
|
if(!silent) cat("\ndone!\n")
|
|
return(well.array)
|
|
}
|
|
|
|
|
|
########################################################################
|
|
# #
|
|
# Output function for generating files from fitted data. #
|
|
# #
|
|
########################################################################
|
|
#' Output function for generating files from fitted data.
|
|
#'
|
|
#' Handles files and directories, calls \code{table.out}, \code{plate.overview} and \code{view.fit}
|
|
#' to generate output tables and graphics.
|
|
#'
|
|
#' @param fitted.well.array A list of fitted well objects.
|
|
#' @param out.prefix Prefix that is in the name of output files.
|
|
#' @param blank.value User can enter a blank OD measurement for uninoculated wells.
|
|
#' If NULL, defaults to the value of the first OD measurement of each well.
|
|
#' @param start.index Which timepoint should be used as the first one after inoculation (defaults to the 2th one)
|
|
#' @param growth.cutoff Minimum threshold for curve growth.
|
|
#' @param points.to.remove A list of numbers referring to troublesome points that should be removed across all wells.
|
|
#' @param remove.jumps Should the slope checking function be on the lookout for large jumps in OD?
|
|
#' @param out.dir name a directory to output the table of curve parameters to (defaults to working directory)
|
|
#' @param graphic.dir name a directory to output the images of the fitted curves to
|
|
#' (defaults to subdirectory "pics" of <out.dir> above)
|
|
#' @param overview.jpgs should jpgs be generated for each plate with the overview graphic?
|
|
#' This is for backwards compatibility with the old web server.
|
|
#' @param silent should messages be returned to the console?
|
|
#' @param unlog should exported graphics be transformed back to the OD scale?
|
|
#' @param source.file.list A list of the source files' names.
|
|
#' @param upload.timestamp The time format indicated by the user.
|
|
#' @param add.constant used to readjust for the constant added during the log transform when plotting ODs.
|
|
#' @param use.linear.param linear parameter is used or not?
|
|
#' @param use.loess Is LOESS model going to be used?
|
|
#' @param plate.nrow The number of rows for a plate.
|
|
#' @param plate.ncol The number of columns for a plate
|
|
#' @param lagRange The heatmap specific range for lag time.
|
|
#' @param totalRange The heatmap specific range for the achieved growth.
|
|
#' @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
|
|
#'
|
|
#' @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,
|
|
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){
|
|
|
|
# Prepare timestamp for addition to output file names.
|
|
filename.timestamp = strftime(upload.timestamp, format="_%Y-%m-%d_%H.%M.%S")
|
|
|
|
########################################################################
|
|
# Prepare to write to output files #
|
|
########################################################################
|
|
if(is.null(blank.value)) blank.value = "First timepoint in well"
|
|
|
|
if(!silent) cat("\nFinding/creating new output directories...")
|
|
|
|
old.wd = getwd()
|
|
# Create output directory if it doesn't exist
|
|
if(class(try(setwd(out.dir), silent = T)) == "try-error"){
|
|
if(!silent) cat("\ncreating new output directory")
|
|
if (class(try(dir.create(out.dir))) == "try-error")
|
|
stop("Error creating new output directory!")
|
|
}
|
|
|
|
# Create graphics directory if it doesn't exist
|
|
if(class(try(setwd(graphic.dir), silent = T)) == "try-error"){
|
|
if(!silent) cat("\ncreating new graphics directory")
|
|
if (class(try(dir.create(graphic.dir))) == "try-error")
|
|
stop("Error creating new graphics directory!")
|
|
}
|
|
|
|
########################################################################
|
|
# Populate a data table with fit results and write to file #
|
|
########################################################################
|
|
#
|
|
# The functions used here are found in table.output.R
|
|
|
|
# Creates a table with a row for each well and a column for each of various identifiers and fitted and calculated parameters.
|
|
|
|
if(!silent) cat("\nPopulating data table...")
|
|
table.fit = try(table.out(fitted.well.array, filename.timestamp=filename.timestamp,use.linear.param=use.linear.param, use.loess=use.loess, constant.added=add.constant))
|
|
|
|
# Return an error if there is a problem with returning the table
|
|
if (class(fitted.well.array) == "try-error")
|
|
stop("Error in <table.out>: ", fitted.well.array)
|
|
|
|
|
|
# Set working directory to <out.dir>
|
|
if (class(try(setwd(out.dir))) == "try-error")
|
|
stop("Error setting directory for table output")
|
|
|
|
# Write output table to file in <out.dir>
|
|
table.filename = paste(out.dir, "/", out.prefix, "_gcat.fit", filename.timestamp, ".txt", sep = "")
|
|
if (class(try(write.table(table.fit, table.filename, sep = "\t", row.names = F))) == "try-error")
|
|
stop("Error writing tabular output")
|
|
|
|
# ---If successfully written, add postscript and start a list of generated files.
|
|
generated.files = table.filename
|
|
|
|
########################################################################
|
|
# Write individual fit and overview graphics to file #
|
|
########################################################################
|
|
#
|
|
# The functions used here are found in graphic.output.R
|
|
|
|
if(!silent) cat("\nDrawing graphics...")
|
|
|
|
# Set working directory to <graphic.dir>
|
|
if (class(try(setwd(graphic.dir))) == "try-error")
|
|
stop("Error setting directory for graphic output")
|
|
# Use function <pdf.by.plate> to write fit graphics to file.
|
|
|
|
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)
|
|
|
|
if (class(graphic.files) == "try-error")
|
|
stop("Error in <pdf.by.plate>: ", graphic.files)
|
|
|
|
# If successfully written, add to the list of generated files.
|
|
generated.files = c(generated.files, graphic.files)
|
|
|
|
########################################################################
|
|
# Add a postscript to the output table with legend and file info. #
|
|
########################################################################
|
|
#
|
|
sink(table.filename, append = T)
|
|
analysis.timestamp = strftime(Sys.time(), format="%Y-%m-%d %H:%M:%S")
|
|
cat("\n# Raw OD values are adjusted and log-transformed before fitting a growth curve as follows: log.OD = log(OD - blank + const) where blank is OD of blank medium and const is specified by the user (1 by default)",
|
|
"\n# Values are reported on the above 'log.OD' scale unless otherwise specified.",
|
|
"\n# .SE columns report standard errors of those values that are estimated directly as parameters of global sigmoid models.",
|
|
"\n# .OD columns report values back-transformed to the linear 'OD - blank' scale.",
|
|
"\n")
|
|
|
|
cat("\n# -- Explanation of columns --",
|
|
"\n# - model: Name of the model the well was successfully fit with (if any)",
|
|
"\n# - lag.time: Lag time estimate inferred from the fitted model",
|
|
"\n# - inflection.time: inflection time point of the growth curve when drawn on the log scale",
|
|
"\n# - max.spec.growth.rate: maximum specific growth rate estimate inferred from the fitted model. Estimated as the first derivative of the growth curve at inflection time point",
|
|
"\n# - baseline: growth curve baseline. Global sigmoid model: baseline is parameter 'b' of the model. LOESS: baseline is the same as the lowest predicted log.OD value",
|
|
"\n# - amplitude: difference between upper plateau and baseline values. Global sigmoid model: amplitude is parameter 'A' of the model. LOESS: amplitude = max.log.OD - min.log.OD",
|
|
"\n# - plateau: upper asymptote value of the fitted model. Global sigmoid model: plateau = b + A. LOESS: plateau = max.log.OD",
|
|
"\n# - inoc.log.OD: log.OD value at inoculation. Estimated value from the fitted model is used, rather than the actual measurement",
|
|
"\n# - max.log.OD: maximal log.OD value reached during the experiment. Estimated value from the fitted model is used rather than the actual measurement",
|
|
"\n# - projected.growth: maximal projected growth over inoculation value. Global sigmoid model: projected.growth = plateau - inoc.log.OD. LOESS: not reported",
|
|
"\n# - achieved.growth: maximal growth over inoculation value actually achieved during the experiment. achieved.growth = max.log.OD - inoc.log.OD",
|
|
"\n# - shape.par: shape parameter of the Richard equation",
|
|
"\n# - R.squared: goodness of fit metric. Also known as coefficient of determination. R.squared is usually between 0 and 1. A value close to 1 indicates good fit.",
|
|
"\n# - RSS: residual sum of squares. Another goodness of fit metric. Smaller values indicate better fits.",
|
|
"\n# - empty: (Well indicator)",
|
|
"\n# - an 'E' indicates that the well was empty and no growth was detected. ",
|
|
"\n# - an 'I' indicates that the well was inoculated and growth was detected above the threshold. ",
|
|
"\n# - an 'E*' indicates that the well was empty and growth was detected (possible contamination). ",
|
|
"\n# - an '!' indicates that the well was inoculated and no growth was detected. ",
|
|
"\n# - asymp.not.reached: shows \"L\" if the bottom asymptote (baseline) was not reached and \"U\" if the upper asymptote (plateau) was not reached.",
|
|
"\n# - tank: (Tanking indicator) If a number is present then the growth trend was determined to tank at that timepoint index.",
|
|
"\n# - other: Additional flag column. Displays information about whether jumps in OD were detected and what was done about them.",
|
|
"\n# - pdf.file and page.no: location of the figure for this well in the output .pdf files."
|
|
)
|
|
|
|
# Analysis information
|
|
|
|
cat("\n#\n# -- Source file information--",
|
|
"\n# ", paste(source.file.list, collapse = "\n# "),
|
|
"\n# analyzed using GCAT v", global.version.number,
|
|
"\n# request sent: ", upload.timestamp,
|
|
"\n# completed: ", analysis.timestamp,
|
|
"\n#\n# -- Parameters used in current analysis --",
|
|
"\n# - Constant added to log(OD + n) transformation:", add.constant,
|
|
"\n# - Blank OD value: ", blank.value,
|
|
"\n# - Index of inoculation timepoint", start.index,
|
|
"\n# - Minimum growth threshold:", growth.cutoff,
|
|
"\n# - Removed points:", paste(points.to.remove, collapse = " "),
|
|
"\n# - Jump detection:", remove.jumps
|
|
)
|
|
|
|
# gcat.analysis.main() starting environment
|
|
cat("\n#\n# -- gcat.analysis.main() starting environment --\n#")
|
|
print(main.envir)
|
|
|
|
# Done with text file output
|
|
sink()
|
|
|
|
########################################################################
|
|
# Return values to R #
|
|
########################################################################
|
|
#
|
|
|
|
if(!silent) cat("\ndone!")
|
|
setwd(old.wd)
|
|
# Return list of generated files
|
|
return(generated.files)
|
|
}
|
|
|