## @importFrom grofit grofit.control gcFitModel gcFitSpline gcBootspline
## @importFrom methods is
### COMMON INTERFACES for growthrates/grofit packages
## TODO: more common functions for all packages
## use function with() to pass specific arguments
## * convertData with arguments grofit or growthrates
## * fitResults to get results
#' Add fits from various growth model interfaces to plate data object.
#'
#' @param fit result object from calls to batch model fitting functions
#' \code{\link{gcFit.2}}, any of \pkg{growthrates} batch functions (`all_'),
#' or \code{platexpress} interfaces \code{\link{segmented_plate}}, and
#' \code{\link{dpseg_plate}}.
#' @param data a platexpress data set, see \code{\link{readPlateData}}
#' @param ID ID for the new data set
#' @param ... arguments to \code{\link{addData}} (eg. \code{col} for
#' color selection)
#' @export
addModel <- function(fit, data, ID="model", ...) {
UseMethod("addModel", fit)
# if ( class(fit)=="gcFit" )
# return(addModel_gcFit(fit=fit, data=data, ID=ID, ...))
# else if ( class(fit)=="dpseg_list" )
# return(addModel_dpseg(fit=fit, data=data, ID=ID, ...))
# else if ( class(fit)=="segmented_list" )
# return(addModel_dpseg(fit=fit, data=data, ID=ID, ...))
# else if ( length(grep("^multiple_",class(fit)))==1 )
# return(addModel_growthrates(fit=fit, data=data, ID=ID, ...))
}
#' merge results to plate layout map
#'
#' a wrapper around R base function \code{\link[base:merge]{merge}}
#' for merging results per well into a \pkg{platexpress} layout map.
#' Usage is the same as for \code{\link[base:merge]{merge}} with
#' some defaults changed.
#' @param x a plate layout map or any data frame with a "well" column
#' specified by argument \code{by}
#' @param y results: any data frame with a "well" column
#' specified by argument \code{by}, see \code{\link[base:merge]{merge}}
#' @param ID prefix added to to the merged \code{y} data
#' @param by column name(s) by which to merge,
#' @param all.x keep all rows of \code{x}, see \code{\link[base:merge]{merge}}
#' @param sort sort the results on the \code{by} column? See
#' \code{\link[base:merge]{merge}}
#' @param ... further arguments to \code{\link[base:merge]{merge}}
#' @export
mergeResults <- function(x, y, ID, by = "well",
all.x=TRUE, sort=FALSE, ...) {
if ( !missing(ID) )
colnames(y)[2:ncol(y)] <- paste0(ID, "_",colnames(y)[2:ncol(y)])
merge(x=x, y=y, by="well", all.x=all.x, sort=sort, ...)
}
#### GROWTH PHASE SEGMENTATION
### SEGMENTED INTERFACE
#' Call \code{\link[segmented]{segmented}} for selected wells.
#'
#' The package \code{\link[segmented]{segmented}} splits curves
#' into linear segments, given a list of pre-defined breakpoints.
#' If the passed data is a biomass measure
#' (eg. OD), and option \code{log=TRUE} the slopes correspond to local
#' growth rate (per time unit). Missing (NA) or non-finite
#' y-values will be removed. \code{\link[segmented]{segmented}} has
#' a random initialization and sometimes fails. The function will
#' attempt \code{maxtry} of calls before giving up.
#' @param data \code{platexpress} data object
#' @param yid ID of the \code{platexpress} data to use
#' @param wells subset and plot order of wells
#' @param log use ln of the data
#' @param xid x-axis data ID in \code{platexpress} data
#' @param man apply a moving average \code{\link{ma}} with \code{n=man}
#' @param maxtry maximum number of attempts to run
#' \code{\link[segmented]{segmented}}
#' @param control a list of parameters for controlling the
#' \code{\link[segmented]{segmented}} fitting process.
#' See the documentation for \code{\link[segmented:seg.control]{seg.control}}
#' for details.
#' @param psis named list of breakpoints for wells, generated from
#' argument \code{psi} if missing
#' @param psi vector of breakpoints (x-values) to be passed to
#' \code{\link[segmented]{segmented}}, generated as equal spaced breakpoints
#' from argument \code{npsi} if missing
#' @param npsi number of equally spaced breakpoints, used if both
#' \code{psis} and \code{psi} are missing
#' @param plot plot each fit
#' @param verb progress messages
#' @param ... arguments passed to \code{\link[segmented]{segmented}}
#' @export
segmented_plate <- function(data, yid="OD", wells, log=TRUE, xid,
man=1, maxtry=5, control=segmented::seg.control(),
psis, psi, npsi=5, plot=FALSE, verb=0, ...) {
if ( missing(xid) ) xid <- data$xids[1]
X <- data[[xid]]
Y <- data[[yid]]$data
if ( man>1 )
Y <- apply(Y, 2, ma, man)
if ( log ) Y <- log(Y)
if ( missing(wells) )
wells <- colnames(Y)
if ( missing(psis) ) {
if ( missing(psi) )
psi <- seq(min(X),max(X),length.out=npsi+2)[1:npsi+1]
psis <- rep(list(psi), length(wells))
names(psis) <- wells
}
segments <- list()
for ( well in wells ) {
if ( verb>0 )
cat(paste("calculating segments for well", well, "\t"))
y <- Y[,well]
x <- X[!is.na(y)]
y <- y[!is.na(y)]
x <- x[is.finite(y)]
y <- y[is.finite(y)]
if ( plot ) {
plot(x,y, type="b",cex=.5,main=paste(well))
}
## linear model
out.lm <- lm(y~x)
## TODO: count and warn
test <- NA
class(test) <- "try-error"
mxtry <- maxtry
while( class(test)[1]=="try-error" & mxtry>0 ) {
test <- try(o <- segmented::segmented(out.lm,control=control,
psi=psis[[well]]))
if ( mxtry<5 )
warning(well, ": segmented failed in attempt: ",
maxtry-mxtry+1, " of ", maxtry)
mxtry <- mxtry -1
}
if ( verb>0 ) cat("\n")
if ( maxtry==0 ) {
warning(well, ": segmented failed")
next
}
segments[[well]] <- o
if ( plot ) {
segmented::plot.segmented(o,add=TRUE, col=2, lwd=3,
shade=TRUE, rug=FALSE)
segmented::lines.segmented(o,col=1, lwd=1)
abline(v=confint(o)$x[,1])
Sys.sleep(.4)
}
}
class(segments) <- "segmentedl"
invisible(segments)
}
#' Add results from \code{\link{segmented_plate}}
#' to \code{platexpress} object.
#'
#' Calls the \code{\link[segmented:predict.segmented]{predict.segmented}} method
#' to reconstruct the piecewise linear growth curve, and adds it to the
#' \code{platexpress} data object. Option \code{add.slopes}
#' allows to add a time-course of growth rates (slopes
#' of the linear segments) instead.
#' @param fit list of \code{\link[segmented:segmented]{segmented}} objects
#' @param data \code{platexpress} data object
#' @param ID data ID for the new object
#' @param add.slopes add slopes instead of reconstructed data
#' @param ... arguments passed to \code{\link{addModel}}
#'@export
addModel.segmentedl <- function(fit, data, ID="y", add.slopes=FALSE, ...) {
x <- data[[data$xids[1]]]
if ( add.slopes ) {
## get breakpoints and slopes
segs <- lapply(fit, function(o) c(min(x),
segmented::confint.segmented(o)$x[,1],
max(x)))
slps <- lapply(fit, function(o) segmented::slope(o)$x[,1])
## add mus as xy data
mus <- sapply(1:length(segs), function(i) {
sg <- segs[[i]]
sl <- slps[[i]]
yo <- rep(sl, each=10)
xo <- c(sapply(2:length(sg), function(j)
seq(sg[j-1]+1e-3, sg[j], length.out=10)))
approx(xo,yo, xout=x)$y
})
colnames(mus) <- names(fit)
## add to platexpress data to plot!
invisible(addData(data, ID=paste(ID,sep="_"),
dat=mus, ...))
} else {
## add modelled data
OD_segs <- lapply(fit, function(o) predict(o,newdata=data.frame(x=x)))
OD_segs <- do.call(cbind, OD_segs)
## add to platexpress data to plot!
invisible(addData(data, ID=paste(ID,sep="_"),dat=exp(OD_segs), ...))
}
}
### DPSEG INTERFACE
#' Call \code{\link[dpseg:dpseg]{dpseg}} for selected wells.
#'
#' The algorithm in \code{\link[dpseg:dpseg]{dpseg}} splits curves
#' into linear segments. If the passed data is a biomass measure
#' (eg. OD), and option \code{log=TRUE} the slopes correspond to local
#' growth rate (per time unit).
#' @param data \code{platexpress} data object
#' @param yid ID of the \code{platexpress} data to use
#' @param wells subset and plot order of wells
#' @param xid x-axis data ID in \code{platexpress} data
#' @param log use ln of the data
#' @param verb progress messages
#' @param ... arguments passed to \code{\link[dpseg:dpseg]{dpseg}}
#' @export
dpseg_plate <- function(data, yid="OD", wells, log=TRUE, xid, verb=0, ...) {
if ( missing(xid) ) xid <- data$xids[1]
x <- data[[xid]]
Y <- data[[yid]]$data
if ( log ) Y <- log(Y)
if ( missing(wells) )
wells <- colnames(Y)
## call dpseg - TODO: use parallel?
segments <- sapply(wells, function(well) {
if ( verb>0 )
cat(paste("calculating segmentation of", yid,
"from well", well, "\n"))
list(dpseg::dpseg(x=x, y=Y[,well], verb=verb, ...))
})
names(segments) <- wells
class(segments) <- "dpsegl"
invisible(segments)
}
#' Add results from \code{\link{dpseg_plate}}
#' to \code{platexpress} object.
#'
#' Calls the \code{\link[dpseg:predict.dpseg]{predict.dpseg}} method
#' to reconstruct the piecewise linear growth curve, and adds it to the
#' \code{platexpress} data object. Option \code{add.slopes}
#' allows to add a time-course of growth rates (slopes
#' of the linear segments) instead.
#' @param data \code{platexpress} data object
#' @param fit list of \code{\link[dpseg:dpseg]{dpseg}} objects
#' @param ID data ID for the new object
#' @param add.slopes add slopes instead of reconstructed data
#' @param ... arguments passed to \code{\link{addModel}}
#'@export
addModel.dpsegl <- function(fit, data, ID="y", add.slopes=FALSE, ...) {
if ( add.slopes ) {
## TODO: fix NA in last row!
## get breakpoints and slopes
segs <- lapply(fit, function(x) x$segments[,c("x1","x2")])
slps <- lapply(fit, function(x) x$segments[1:nrow(x$segments),"slope"])
x <- data[[data$xids[1]]]
## add mus as xy data
mus <- sapply(1:length(segs), function(i) {
sg <- segs[[i]]
sl <- slps[[i]]
yo <- rep(sl, each=2)
## TODO: *.999 looks very dirty
xo <- c(sapply(1:nrow(sg), function(j)
seq(sg[j,1], sg[j,2]*.999, length.out=2)))
approx(xo,yo, xout=x)$y
})
colnames(mus) <- names(fit)
## add to platexpress data to plot!
invisible(addData(data, ID=paste(ID,sep="_"), dat=mus, ...))
} else {
## add modelled data
OD_segs <- lapply(fit, function(o) predict(o)$y)
OD_segs <- do.call(cbind, OD_segs)
## add to platexpress data to plot!
invisible(addData(data, ID=paste(ID,sep="_"), dat=exp(OD_segs), ...))
}
}
#### GROWTH MODEL FITS
### GROFIT INTERFACE
### This would require to depend on grofit
## @examples
## data(ap12)
## grdat <- data2grofit(ap12data)
## fit <- gcFit.2(grdat$time, grdat$data)
### HIGH-LEVEL WRAPPER FOR GROFIT
#' high-level wrapper for \pkg{grofit}
#'
#' Calls \pkg{grofit} directly from \code{platexpress}
#' plate data and layout, using established settings. Please see
#' the documentation of package \pkg{grofit}
#' for details of the fitting procedure. This high-level wrapper uses
#' \code{platexpress} functions \code{\link{data2grofit}},
#' \code{\link{grofit.2.control}}, \code{\link{gcFit.2}} to convert
#' data, to set parameters and to call \pkg{grofit}.
#'
#' The function returns a list of fits, as returnted by
#' \code{platexpress}'s \code{\link{gcFit.2}}, a copy of
#' grofit's \code{\link[grofit:gcFit]{gcFit}} with modified
#' plotting and interactive behaviour but otherwise identical.
#' All run parameters of \code{\link[grofit:gcFit]{gcFit}}
#' can be set via option \code{control}, where \code{platexpress}'
#' \code{\link{grofit.2.control}} add the differing plot and interactive
#' settings.
#' @param data a platexpress data set, see \code{\link{readPlateData}}
#' @param plate plate layout map, see \code{\link{readPlateMap}}, columns
#' of this map can be converted to \pkg{grofit} data
#' annotation, and results will be returned in the same order of wells
#' @param yid data ID of the data to be converted; default
#' is to take the first of \code{data$dataIDs}
#' @param amount name of a column in \code{plate} that provides a numeric
#' 'amount', to be used as a 'dose' annotation in grofit
#' @param fields column IDs in the plate layout map to be used for
#' \pkg{grofit} data annotation; if \code{plate} is not
#' missing at least 1 column is required
#' @param control an object of class \code{grofit.control} with an additional
#' parameter \code{plot}, as provided by the \code{platexpress} function
#' \code{\link{grofit.2.control}}; if present arguments \code{model.type},
#' \code{nboot.gc} and \code{plot} will be ignored (unless the latter
#' is not present)
#' @param model.type models that \code{grofit} will attempt to fit,
#' see \code{\link[grofit:grofit.control]{grofit.control}}; will be igonred
#' if \code{control} is passed directly
#' @param nboot.gc numbers of bootstrap samples,
#' see \code{\link[grofit:grofit.control]{grofit.control}}; will be igonred
#' if \code{control} is passed directly
#' @param plot set to TRUE for plots even without \code{interactive} use
#' @param interactive set to TRUE for interactive growth curve fitting
#' @param col color of the model data in the plate data object
#' @param verb print messages
#' @param ... further arguments to \code{\link{grofit.2.control}}
#'@export
### do not add to model, just return list of results
grofit_plate <- function(data, plate, yid, amount,
fields=c("strain","medium","substance"),
control, model.type=c("richards","logistic",
"gompertz.exp","gompertz"),
nboot.gc=100, plot=TRUE, interactive=FALSE, verb=TRUE,
col="#0000FF", ...) {
## take the first data, if none was passed
if ( missing(yid) ) {
yid <- data$dataIDs[1]
if ( verb )
cat(paste("auto-selected data:", yid, "\n"))
}
## generate input data, incl. plate annotation
if ( missing(plate) )
gdat <- data2grofit(data, yid=yid, verb=verb)
else {
## if a plate map is available, use this for annotation!
## allow for full capabilities of grofit results: add dose
if ( missing(amount) )
gdat <- data2grofit(data, plate=plate, yid=yid,
wells=plate$well[!plate$blank],
eid=fields, verb=verb, ...)
else
gdat <- data2grofit(data, plate=plate, yid=yid,
wells=plate$well[!plate$blank],
dose=plate[[amount]][!plate$blank],
eid=fields, verb=verb, ...)
}
## set grofit paraameters
if ( missing(control) )
control <- grofit.2.control(plot=plot, interactive=interactive,
suppress.messages=TRUE,
nboot.gc=nboot.gc,
model.type=model.type)
if ( !"plot" %in% names(control) )
control$plot <- plot
## call grofit; redirect plots to detail pdf
odfit <- gcFit.2(gdat$time, gdat$data, control)
invisible(odfit)
}
### data2grofit: see AP12.R for example, TODO: fix example data and update file
#' interface to package \pkg{grofit}
#'
#' \code{\link{data2grofit}} : converts \code{\link{platexpress}} data to
#' \pkg{grofit} data format
#' @param data a platexpress data set, see \code{\link{readPlateData}}
#' @param yid data ID of the data to be converted for grofit, from \code{data$dataIDs}
#' @param min.time minimal time of the data to be used
#' @param max.time maximal time of the data to be used
#' @param wells column IDs of the data set to use, if missing all wells
#' are used
#' @param plate plate layout map, see \code{\link{readPlateMap}}, columns
#' of this map can be converted to \pkg{grofit} data
#' annotation
#' @param eid column IDs in the plate layout map to be used for
#' \pkg{grofit} data annotation; if missing but
#' \code{plate} is present, the columns 2 and 3 are used
#' @param dose vector of doses in each well, used as the third column of
#' \pkg{grofit} data annotation, where it can be used for
#' dose-response calculations
#' @param verb print messages
#' @details Returns a simple list with two entries \code{time} and \code{data},
#' as required for \pkg{grofit}.
#' @author Rainer Machne \email{raim@tbi.univie.ac.at}
#' @export
data2grofit <- function(data, yid, min.time, max.time, wells, plate, eid, dose, verb=TRUE) {
if ( missing(yid) ) {
yid <- data$dataIDs[1]
if ( verb )
cat(paste("auto-selected data:", yid, "\n"))
}
dat <- data[[yid]]$data
if ( missing(wells) )
wells <- colnames(dat)
else wells <- as.character(wells)
dat <- data.frame(dat[,wells])
## expand time to full matrix
## TODO: use internal time and non-interpolated data?
time <- data$Time
if ( !missing(max.time) ) {
dat <- dat[time<=max.time,]
time <- time[time<=max.time]
}
if ( !missing(min.time) ) {
dat <- dat[time>=min.time,]
time <- time[time>=min.time]
}
time <- t(matrix(rep(time, ncol(dat)), c(length(time), ncol(dat))))
## well annotation
## use well as TestId
if ( !missing(plate) ) {
if ( missing(eid) ) ## add second column as default
eid <- colnames(plate)[2]
## sub-selection
idx <- match(wells, as.character(plate[,"well"]))
annotation1 <- plate[idx,"well"]
## use both fields if two eids are given ..
annotation2 <- plate[idx,eid[1]]
if ( length(eid)>1 )
for ( k in 2:length(eid) )
annotation2 <- paste0(annotation2,"-",plate[idx,eid[k]])
annotation <- data.frame(annotation1, annotation2)
} else
annotation <- data.frame(cbind(colnames(dat),
rep("",ncol(dat))))
## dose information for grofit dose-response calculations
## TODO: this is ugly, do nicer!
found.dose <- !missing(dose)
if ( !found.dose )
if ( !missing(plate) ) ## get dose info from plate layout: TODO
if ( "amount" %in% colnames(plate) ) {
idx <- match(wells,as.character(plate[,"well"]))
dose <- as.numeric(plate[idx,"amount"])
found.dose <- TRUE
}
if ( !found.dose )
dose <- rep(0,ncol(dat))
## construct grofit data
grdat <- data.frame(annotation, dose, t(dat))
list(time=time, data=grdat)
}
#' parse \pkg{grofit} results
#'
#' parses the output of \code{\link{gcFit.2}} into a table
#' of the main model parameters, for each well, for both the
#' "best" model found by `grofit` as well the `spline` and `bootstrap`
#' fits. Parameter names are streamlined, ie. the ".model" suffix is removed.
#' @param fit grofit object, the result of a call to
#' \code{\link[grofit:gcFit]{gcFit}}
#' or the \code{platexpress} version \code{\link{gcFit.2}}
#' @param p the parameters to retrieve from \code{fit$gcTable},
#' the main results table of \code{\link[grofit:gcFit]{gcFit}}
#' @seealso \code{\link{data2grofit}}, \code{\link{growthratesResults}}
#' @export
### TODO: generic results.gcFit/.multiple_fits/.dpsegl/.segmented
### for the latter this may involve alignment of growth phases?
grofitResults <- function(fit, p=c("lambda.model","mu.model",
"A.model","used.model",
"lambda.spline","mu.spline","A.spline",
"lambda.bt","mu.bt","A.bt")) {
params <- fit$gcTable[,p]
colnames(params) <-sub("used","model",sub(".model","",colnames(params)))
rownames(params) <- as.character(fit$gcTable[,"TestId"])
data.frame(well=rownames(params), params)
}
#' Add results from \code{\link{grofit_plate}}
#' to \code{platexpress} object.
#'
#' Calls the \code{\link[stats:predict]{predict}} method
#' for the growth curves fits returned
#' by grofit's \code{\link[grofit:gcFit]{gcFit}} or via \code{platexpress}'s
#' \code{\link{gcFit.2}} or
#' \code{\link{grofit_plate}}, and adds it to the
#' \code{platexpress} data object.
#' @param fit a \code{\link[grofit:gcFit]{gcFit}} object
#' @param data \code{platexpress} data object
#' @param ID data ID for the new object
#' @param ... arguments passed to \code{\link{addData}}
#'@export
addModel.gcFit <- function(fit, data, ID="model", ... ) {
testid <- "TestId" # this requires wells being used as TestId in grofit
## copy first existing data set
newdat <- data[[data$dataIDs[[1]]]]$data
newdat[] <- NA
models <- rep(NA, ncol(newdat))
names(models) <- colnames(newdat)
## parse grofit results and use predict to add data
for ( i in 1:ncol(newdat) ) {
id <- colnames(newdat)[i]
idx <- which(as.character(fit$gcTable[,testid])==id)
if ( !length(idx) ) next
if ( !fit$gcFittedModels[[idx]]$fitFlag ) next
newdat[,i] <- stats::predict(fit$gcFittedModels[[idx]]$nls,
newdata=data.frame(time=data$Time))
models[i] <- fit$gcFittedModels[[idx]]$model
}
## TODO: use models info somehow?
## add and return
addData(data=data, ID=ID, dat= newdat,
processing="grofit prediction", ...)
}
#' wrapper of grofit function \code{\link[grofit:grofit.control]{grofit.control}}
#' which adds the new \code{plot} switch used in \code{\link{gcFit.2}}
#' @param ... parameters passed on to
#' \code{\link[grofit:grofit.control]{grofit.control}}
#' @param interactive set to TRUE for interactive growth curve fitting
#' @param plot set to TRUE for plots even without \code{interactive} use
#' @seealso \code{\link{gcFit.2}}
#' @export
grofit.2.control <- function(interactive=FALSE, plot=TRUE, ...) {
control <- grofit::grofit.control(...)
control$interactive <- interactive
control$plot <- plot
control
}
#' hack of grofit function \code{\link[grofit:gcFit]{gcFit}} to allow plotting
#' even when option \code{interative==FALSE}
#' @param time a matrix of measurment times for each well as required by
#' \pkg{grofit}, provided by \code{\link{data2grofit}}
#' @param data the data as produced by \code{\link{data2grofit}}
#' @param control the \pkg{grofit} control structure, see
#' \code{\link[grofit:grofit]{grofit.control}}, but with an additional
#' argument "plot" (TRUE or FALSE), allowing to plot individual fits
#' even if interactive is set to FALSE.
#' @details Adding the field "plot" allows to de-activate
#' \code{control$interactive}, so the whole fitting procedure runs
#' through all wells, yet results are plotted. \code{\link{platexpress}}
#' also offers a hack of \code{\link[grofit:grofit.control]{grofit.control}}
#' names \code{\link{grofit.2.control}} that adds these options.
#' @seealso \code{\link{gcFit.2}}
#' @export
gcFit.2 <- function (time, data, control = grofit.2.control()) {
requireNamespace("grofit")
if (methods::is(control) != "grofit.control")
stop("control must be of class grofit.control!")
if ((dim(time)[1]) != (dim(data)[1]))
stop("gcFit: Different number of datasets in data and time")
if (!is.element(control$fit.opt, c("s", "m", "b"))) {
warning("fit.opt must be set to 's', 'm' or 'b'. Changed to 'b'!")
control$fit.opt = "b"
}
## add gcFit.2 specific controls
if ( !"plot"%in%names(control) )
control$plot <- TRUE
out.table <- NULL
used.model <- NULL
fitpara.all <- list()
fitnonpara.all <- list()
boot.all <- list()
fitted.param <- NULL
fitted.nonparam <- NULL
bootstrap.param <- NULL
for (i in 1:dim(data)[1]) {
acttime <- as.numeric(as.matrix(time[i, ]))
actwell <- as.numeric(as.matrix((data[i, -1:-3])))
gcID <- as.matrix(data[i, 1:3])
if ((control$suppress.messages == FALSE)) {
cat("\n\n")
cat(paste("= ", as.character(i),
". growth curve =================================\n",
sep = ""))
cat("----------------------------------------------------\n")
}
if ((control$fit.opt == "m") || (control$fit.opt == "b")) {
fitpara <- grofit::gcFitModel(acttime, actwell, gcID, control)
fitpara.all[[i]] <- fitpara
}
else {
fitpara <- list(raw.x = acttime, raw.y = actwell,
gcID = gcID, fit.x = NA, fit.y = NA, parameters = list(A = NA,
mu = NA, lambda = NA, integral = NA), model = NA,
nls = NA, reliable = NULL, fitFlag = FALSE, control = control)
class(fitpara) <- "gcFitModel"
fitpara.all[[i]] <- fitpara
}
if ((control$fit.opt == "s") || (control$fit.opt == "b")) {
nonpara <- grofit::gcFitSpline(acttime, actwell, gcID, control)
fitnonpara.all[[i]] <- nonpara
}
else {
nonpara <- list(raw.x = acttime, raw.y = actwell,
gcID = gcID, fit.x = NA, fit.y = NA,
parameters = list(A=NA,mu=NA,lambda=NA,integral=NA),
parametersLowess = list(A=NA,mu=NA,lambda=NA),
spline = NA, reliable = NULL,
fitFlag = FALSE, control = control)
class(nonpara) <- "gcFitSpline"
fitnonpara.all[[i]] <- nonpara
}
wellname <- paste(as.character(data[i, 1]), as.character(data[i,
2]), as.character(data[i, 3]), sep = "-")
if ((control$plot == TRUE)) {
if (fitpara$fitFlag == TRUE) {
plot(fitpara, colData = "black", colModel = "blue", cex = 0.5)
plot(nonpara, add = TRUE, raw = FALSE, colData = 0,
colSpline = "red", cex = 1.5)
abline(h=summary(fitpara)$A,lwd=2,col="blue")
abline(v=summary(fitpara)$lambda,lwd=2,col="blue")
}
else {
plot(nonpara, add = FALSE, raw = TRUE, colData = "black",
colSpline = "red", cex = 0.5)
}
title(wellname)
if (control$fit.opt == "m")
legend(x = "bottomright", legend = fitpara$model,
col = "black", lty = 1)
if (control$fit.opt == "s")
legend(x = "bottomright", legend = "spline fit",
col = "red", lty = 1)
if (control$fit.opt == "b")
legend(x = "bottomright", legend = c(fitpara$model,
"spline fit"), col = c("blue", "red"), lty = c(1,
1))
}
reliability_tag <- NA
if (control$interactive == TRUE) {
answer <- readline("Are you satisfied (y/n)?")
if (substr(answer, 1, 1) == "n") {
cat("\n Tagged this well as unreliable !\n\n")
reliability_tag <- FALSE
fitpara.all[[i]]$reliable <- FALSE
fitnonpara.all[[i]]$reliable <- FALSE
}
else {
reliability_tag <- TRUE
fitpara.all[[i]]$reliable <- TRUE
fitnonpara.all[[i]]$reliable <- TRUE
cat("Well was (more ore less) o.k.\n")
}
}
else {
reliability_tag <- TRUE
}
if (control$interactive == TRUE)
graphics.off()
if ((control$nboot.gc > 0) && (reliability_tag == TRUE)) {
bt <- grofit::gcBootSpline(acttime, actwell, gcID, control)
boot.all[[i]] <- bt
}
else {
bt <- list(raw.x = acttime, raw.y = actwell, gcID = gcID,
boot.x = NA, boot.y = NA, boot.gcSpline = NA,
lambda = NA, mu = NA, A = NA, integral = NA,
bootFlag = FALSE, control = control)
class(bt) <- "gcBootSpline"
boot.all[[i]] <- bt
}
description <- data.frame(TestId = data[i, 1], AddId = data[i,
2], concentration = data[i, 3], reliability = reliability_tag,
used.model = fitpara$model, log.x = control$log.x.gc,
log.y = control$log.y.gc, nboot.gc = control$nboot.gc)
fitted <- cbind(description, summary(fitpara), summary(nonpara),
summary(bt))
out.table <- rbind(out.table, fitted)
}
gcFit <- list(raw.time = time, raw.data = data, gcTable = out.table,
gcFittedModels = fitpara.all, gcFittedSplines = fitnonpara.all,
gcBootSplines = boot.all, control = control)
class(gcFit) <- "gcFit"
## TODO: add fits to data
## wrapper: use platexpress data as input
## and add fits to data, aligned with cellGrowth::fitCellGrowth.2 wrapper
gcFit
}
### growthrates INTERFACE
## TODO: wrapper to call all_easylinearfits etc?
#' interface to package \pkg{growthrates}
#'
#' converts \code{platexpress} data for
#' use with \pkg{growthrates}
#' @param data platexpress data set, see \code{\link{readPlateData}}
#' @param yid data ID of the data to be converted for growthrates, from
#' \code{data$dataIDs}
#' @param wells column IDs of the data set to use, if missing all wells
#' are used
#' @param plate plate layout map, see \code{\link{readPlateMap}}, to skip
#' blanks wells
#' annotation
#' @seealso \code{\link{growthratesResults}}, \code{\link{data2grofit}}
#' @export
data2growthrates <- function(data, yid, wells, plate) {
if ( missing(wells) )
wells <- colnames(data[[yid]]$data)
## filter by plate & rm blanks
if ( !missing(plate) ) {
wells <- wells[wells%in%plate[,"well"]]
wells <- wells[!plate[match(wells, plate[,"well"]),"blank"]]
}
dat <- data[[yid]]$data[,wells,drop=FALSE]
value <- c(dat)
## NOTE: this also works for OD interpolated data, but name "time"
## is perhaps misleading
time <- rep(data[[data$xids[1]]], ncol(dat))
well <- rep(colnames(dat), each=nrow(dat))
df <- data.frame(time=time, value=value,
well=factor(well, levels=colnames(dat)))
df
}
## TODO: simple key word based switch
## for growthrates functions?
growthrates_plate <- function() {}
#' parse fitted parameters from package \pkg{growthrates}
#'
#' parses the output of \code{\link{gcFit.2}} into a table
#' of the main model parameters, for each well
#' @param fit growthrates object, the result of a call to fit functions
#' from package \pkg{growthrates}, such as
#' \code{\link[growthrates:all_easylinear]{all_easylinear}}
#' @param scale.richards if the \code{beta} parameter from Richard's model is present,
#' multiply growthrate \code{mu} with \code{beta} (original stored as \code{mumax})
#' @seealso \code{\link{data2growthrates}}, \code{\link{grofitResults}}
#' @export
growthratesResults <- function(fit, scale.richards=TRUE) {
res <- data.frame(growthrates::results(fit), stringsAsFactors = FALSE)
## replace names to match names in other packages
## lambda, mu; but keep K (A for some models, but K in Monod)
nms <- colnames(res)
nms <- sub("y0","X0",sub("lag","lambda",sub("mumax","mu",nms)))
colnames(res) <- nms
if ( scale.richards & "beta" %in% nms ) {
## NOTE/TODO: scaling mu from richards?
res[,"mumax"] <- res[,"mu"]
res[,"mu"] <- res[,"mu"] * res[,"beta"]
## TODO: also scale gompertz?
}
data.frame(res)
}
#' `predict' hack for \pkg{growthrates} fits
#'
#' \pkg{growthrates} provides several
#' functions to fit growth rates, but their result objects have currently
#' slighly different structures/interfaces. To obtain relevant data from
#' \code{\link[growthrates:fit_easylinear]{fit_easylinear}} and
#' \code{\link[growthrates:fit_spline]{fit_spline}} the fitted
#' functions and parameters have to be directly used instead of just
#' calling \code{\link[growthrates:predict]{predict}} methods.
#' @param fit fit object from \pkg{growthrates} fitting functions, either
#' a single fit (`fit_` functions) or a list of fits from batch functions
#' (`all_` functions)
#' @param time the time points at which growth data is to be predicted
#' @export
grpredict <- function(fit, time) {
## NOTE: predict not implemented for easylinear and returning
## the full spline fit for smooth.spline
## call recursively for lists of fits
if ( length(grep("^multiple_",class(fit)))==1 )
return(lapply(fit@fits, grpredict, time=time))
if ( class(fit)=="easylinear_fit" ) {
xy <- fit@FUN(time, fit@par)[,1:2]
## NOTE: easylinear requires to add the lag phase
xy[,1] <- xy[,1] + growthrates::coef(fit)["lag"]
## interpolate to requested time and convert to matrix
xy <- matrix(unlist(approx(x=xy[,1], y=xy[,2], xout=time)),
ncol=2, byrow=FALSE)
} else if ( class(fit)=="smooth.spline_fit" ) {
## NOTE: smooth.spline predict returns full spline fit as log(y)
xy <- fit@FUN(time, fit@par)[,1:2]
} else {
xy <- growthrates::predict(fit, newdata=list(time=time))[,c("time","y")]
}
xy
}
## adds predict data from growthrates to plate data; called from addModel
#' Add results from package \pkg{growthrates}
#' to \code{platexpress} object.
#'
#' Calls the appropriate \code{\link[stats:predict]{predict}} methods
#' for the growth curves fits returned by
#' by \pkg{growthrates}'s batch functions (`all_`),
#' and adds results to the \code{platexpress} data object.
#' @param fit a \code{\link[grofit:gcFit]{gcFit}} object
#' @param data \code{platexpress} data object
#' @param ID data ID for the new object
#' @param ... arguments passed to \code{\link{addData}}
#'@export
addModel.multiple_fits <- function(fit, data, ID="model", ... ) {
## global time and new data matrix
time <- data[[data$xids[1]]]
## newdat matrix, just copy first data set
newdat <- data[[data$dataIDs[1]]]$data
newdat[] <- NA
## get results
lst <- grpredict(fit, time=time)
## convert to matrix
lst <- lapply(lst, function(x) x[,2])
newdat[,names(lst)] <- matrix(unlist(lst),
ncol = length(lst), byrow = FALSE)
## add and return
addData(data=data, ID=ID, dat= newdat,
processing=paste("growthrates",class(fit),"prediction"), ...)
}
### cellGrowth INTERFACE
#' hack of the cellGrowth package function
#' \code{\link[cellGrowth:bandwidthCV]{bandwidthCV}}
#' to work with the \code{platexpress} data format, using adjusted default
#' values
#' TODO: select bandwidths from data$Time
#'
#' @param data platexpress data set, see \code{\link{readPlateData}}
#' @param yid data ID of the data to be converted for grofit, from
#' \code{data$dataIDs}
#' @param xid the x-axis to be used, defaults to the first available
#' (usually "Time")
#' @param wells column IDs of the data set to use, if missing all wells
#' are used
#' @param bandwidths see \code{\link[cellGrowth:bandwidthCV]{bandwidthCV}}
#' @param nFold see \code{\link[cellGrowth:bandwidthCV]{bandwidthCV}}
#' @param nWell see \code{\link[cellGrowth:bandwidthCV]{bandwidthCV}}
#' @param cutoff see \code{\link[cellGrowth:bandwidthCV]{bandwidthCV}}
#' @param scaleY see \code{\link[cellGrowth:bandwidthCV]{bandwidthCV}}
#' @seealso \code{\link[cellGrowth:bandwidthCV]{bandwidthCV}}
#'
#' @export
bandwidthCV.2 = function(data, yid, xid, wells,
bandwidths=seq(0.5,10, length.out=30), # hours - TODO: take 1/5th of data$Time
nFold=10,
nWell=50,
cutoff=0.95,
##calibration=identity,
scaleY=identity # log2
) {
if ( missing(xid) )
xid <- data$xids[1]
if ( missing(yid) ) # only use requested data
yid <- data$dataIDs[1]
## function to split sample in n folds
cv.folds = function (n ) {
split(sample(1:n), rep(1:nFold, length = n))
}
## function to predict max growth rate
cvpred_mu = function(cv, h, x, z) {
ls = lapply(
cv,
function(fo){
fit = tryCatch(
cellGrowth::fitCellGrowth(x = x[-fo], z = z[-fo],
model = "locfit", locfit.h = h),
warning=function(w){NA}, error=function(e){NA})
if( !is.null(attr(fit,"maxGrowthRate")) ) {
rv = list(pred=stats::predict( fit, x[fo]),
mu = attr(fit, "maxGrowthRate"))
}else{
rv = list(pred=rep(NA, length(fo)), mu = NA)
}
rv
})
list(pred = unlist(lapply(ls, function(x) x$pred)),
mu = unlist(lapply(ls, function(x) x$mu)))
}
## function to calculate error and std
err2_mustd_well = function(dat, w, hs) {
x = dat$Time
y = dat[[yid]]$data[ , match(w, wells)]
## negative ODs?
if ( !all(y>=0)) {
warning("Negative values ",yid," found. If you are using a calibration function, this might be the problem. Values are set to NAs")
y[y<0] = NA
}
z = scaleY( y )
cv = cv.folds(length(x))
pms = lapply(hs, function(h) cvpred_mu(cv,h,x,z))
pred = sapply(pms, function(x) x$pred)
mu = sapply(pms, function(x) x$mu)
err = pred - z[unlist(cv)]
list(err2 = colMeans(err^2),
err2std = apply(err^2, 2, sd),
mustd = apply(mu, 2, sd))
}
## do cv
ws = sample(length(wells),nWell) # wells to perform cv on
hs = bandwidths # possible bandwidths
## calculate mu error and std
err2_mustd = lapply(ws,
function(w) {
cat(paste("Treating well",w,colnames(data[[yid]]$data)[w],"\n"))
data <- err2_mustd_well(dat=data, w=wells[w], hs=hs)
})
err2 = sapply(err2_mustd, function(x) x$err2)
err2std = sapply(err2_mustd, function(x) x$err2std)
mustd = sapply(err2_mustd, function(x) x$mustd)
## which bandwidth are at one std of mini in average
onestd_of_mini = sapply(1:ncol(err2),
function(i) {
m = which.min(err2[,i])
err2[,i] <= err2[m,i] + err2std[m,i]
})
## onestd_of_mini returns a list
if( is.list(onestd_of_mini) )
onestd_of_mini = do.call(cbind,onestd_of_mini)
## calculate "optimal" bandwidth
bandwidth = hs[max( which(rowMeans(onestd_of_mini,na.rm=TRUE)>= cutoff))]
return(list(bandwidth=bandwidth,
wells=wells,
bandwidths=hs,
err2=err2,
err2std=err2std,
muStd=mustd,
oneStdOfMini=onestd_of_mini))
}
#' batch wrapper for fitCellGrowth
#' TODO: align in/out with gcFit.2 or optionally cellGrowth::fitCellGrowths
#'
#' @param data platexpress data set, see \code{\link{readPlateData}}
#' @param yid data ID of the data to be converted for grofit, from
#' \code{data$dataIDs}
#' @param xid the x-axis to be used, defaults to the first available
#' (usually "Time")
#' @param wells column IDs of the data set to use, if missing all wells
#' are used
#' @param ... further parameters to the \code{cellGrowth} package function
#' \code{\link[cellGrowth:fitCellGrowth]{fitCellGrowth}}
#' @seealso \code{\link[cellGrowth:fitCellGrowth]{fitCellGrowth}}
#'
#' @export
fitCellGrowths.2 = function(data, yid, xid, wells, ...) {
if ( missing(xid) )
xid <- data$xids[1]
if ( missing(yid) ) # only use requested data
yid <- data$dataIDs[1]
if ( missing(wells) )
wells <- colnames(data[[yid]]$data)
##fits <- matrix(NA,ncol=4,nrow=length(wells))
##rownames(fits) <- wells
fits <- NULL
x <- data[[xid]]
for ( well in wells ) {
cat(paste("well", well, "\n"))
y <- data[[yid]]$data[,well]
fit <- cellGrowth::fitCellGrowth(x, z=y, ...)
##fits[well,] <- unlist(attributes(fit)[c(3,4,5,6)])
fits <- append(fits, list(fit))
## TODO: align with grofit output and add to data
## calculate lag from pointOfMaxGrowthRate and MaxGrowthRate
## calculate X(0)
## TODO: add fits to data for plots
## data[[yid]]$fit <- fit...data
}
names(fits) <- wells
fits
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.