Nothing
R/fluorescence_workflows.R
In QurvE: Robust and User-Friendly Analysis of Growth and Fluorescence Curves
Defines functions
fl.workflow
flFit
Documented in
flFit
fl.workflow
#' Perform a fluorescence curve analysis on all samples in the provided dataset.
#'
#' \code{flFit} performs all computational fluorescence fitting operations based on the user input.
#'
#' @param fl_data Either... \enumerate{ \item a \code{grodata} object created with \code{\link{read_data}} or \code{\link{parse_data}},
#' \item a list containing a \code{'time'} matrix (for x_type == "time") or \code{'growth'} dataframe (for x_type == "growth") and a \code{'fluorescence'} dataframes,
#' or \item a dataframe containing (normalized) fluorescence values (if a \code{time} matrix or \code{growth} dataframe is provided as separate argument).}
#' @param time (optional) A matrix containing time values for each sample.
#' @param growth (optional) A dataframe containing growth values for each sample and sample identifiers in the first three columns.
#' @param control A \code{fl.control} object created with \code{\link{fl.control}}, defining relevant fitting options.
#' @param parallelize Run linear fits and bootstrapping operations in parallel using all but one available processor cores
#' @param ... Further arguments passed to the shiny app.
#'
#' @return An \code{flFit} object that contains all fluorescence fitting results, compatible with
#' various plotting functions of the QurvE package.
#' \item{raw.x}{Raw x matrix passed to the function as \code{time} (for x_type = 'time') or \code{growth} (for x_type = 'growth').}
#' \item{raw.fl}{Raw growth dataframe passed to the function as \code{data}.}
#' \item{flTable}{Table with fluorescence parameters and related statistics for each fluorescence curve evaluation performed by the function. This table, which is also returned by the generic \code{summary.flFit} method applied to a \code{flFit} object, is used as an input for \code{\link{fl.drFit}}.}
#' \item{flFittedLinear}{List of all \code{flFitLinear} objects, generated by the call of \code{\link{flFitLinear}}. Note: access to each object in the list via double brace: flFittedLinear\[\[#n\]\].}
#' \item{flFittedSplines}{List of all \code{flFitSpline} objects, generated by the call of \code{\link{flFitSpline}}. Note: access to each object via double brace: flFittedSplines\[\[#n\]\].}
#' \item{flBootSplines}{List of all \code{flBootSpline} objects, generated by the call of \code{\link{flBootSpline}}. Note: access to each object via double brace: flFittedSplines\[\[#n\]\].}
#' \item{control}{Object of class \code{fl.control} containing list of options passed to the function as \code{control}.}
#'
#' @details
#' Common response parameters used in dose-response analysis:<br><br><b>Linear fit:</b><br>- max_slope.linfit: Fluorescence increase rate<br>- lambda.linfit: Lag time<br>- dY.linfit: Maximum Fluorescence - Minimum Fluorescence<br>- A.linfit: Maximum fluorescence<br><br><b>Spline fit:</b><br>- max_slope.spline: Fluorescence increase rate<br>- lambda.spline: Lag time<br>- dY.spline: Maximum Fluorescence - Minimum Fluorescence<br>- A.spline: Maximum fluorescence<br>- integral.spline: Integral<br><br><b>Parametric fit:</b><br>- max_slope.model: Fluorescence increase rate<br>- lambda.model: Lag time<br>- dY.model: Maximum Fluorescence - Minimum Fluorescence<br>- A.model: Maximum fluorescence<br>- integral.model: Integral'
#'
#' @family workflows
#' @family fluorescence fitting functions
#' @family dose-response analysis functions
#'
#' @export
#'
#' @importFrom foreach %dopar%
#'
#' @examples
#' # load example dataset
#' input <- read_data(data.growth = system.file("lac_promoters_growth.txt", package = "QurvE"),
#' data.fl = system.file("lac_promoters_fluorescence.txt", package = "QurvE"),
#' csvsep = "\t",
#' csvsep.fl = "\t" )
#'
#' # Define fit controls
#' control <- fl.control(fit.opt = "s",
#' x_type = "time", norm_fl = TRUE,
#' dr.parameter = "max_slope.spline",
#' dr.method = "model",
#' suppress.messages = TRUE)
#'
#' # Run curve fitting workflow
#' res <- flFit(fl_data = input$norm.fluorescence,
#' time = input$time,
#' control = control,
#' parallelize = FALSE)
#'
#' summary(res)
#'
flFit <- function(fl_data, time = NULL, growth = NULL, control= fl.control(), parallelize = TRUE, ...)
{
old.options <- options()
on.exit(options(old.options))
# Define objects based on additional function calls
call <- match.call()
## remove strictly defined arguments
call$time <- call$growth <- call$fl_data <- call$control<- call$parallelize <- NULL
arglist <- sapply(call, function(x) x)
arglist <- unlist(arglist)[-1]
## Assign additional arguments (...) as R objects
if(length(arglist) > 0){
for(i in 1:length(arglist)){
assign(names(arglist)[i], arglist[[i]])
}
}
x_type <- control$x_type
if(!(class(fl_data) %in% c("list", "grodata"))){
if (x_type == "time" && is.numeric(as.matrix(time)) == FALSE)
stop("Need a numeric matrix for 'time' (for x_type = 'time') or a grodata object created with read_data() or parse_data() in the 'fl_data' argument.")
if (x_type == "growth" && is.numeric(as.matrix(growth[-1:-3])) == FALSE)
stop("Need a dataframe for 'growth' (for x_type = 'growth') or a grodata object created with read_data() or parse_data() in the 'fl_data' argument.")
if (is.numeric(as.matrix(fl_data[-1:-3])) == FALSE)
stop("Need a dataframe for 'fl_data' or a grodata object created with read_data() or parse_data() in the 'fl_data' argument.")
} else if(!is.null(fl_data)){
time <- fl_data$time
growth <- fl_data$growth
fl_data <- fl_data$fluorescence
}
# /// check if start growth values are above min.growth in all samples
if(!is.null(growth) && length(growth) > 1){
max.growth <- unlist(lapply(1:nrow(growth), function (x) max(as.numeric(as.matrix(growth[x,-1:-3]))[!is.na(as.numeric(as.matrix(growth[x,-1:-3])))])))
if(is.numeric(control$min.growth) && control$min.growth != 0){
if(!is.na(control$min.growth) && all(as.numeric(max.growth) < control$min.growth)){
stop(paste0("The chosen global start growth value (min.growth) is larger than every value in your dataset.\nThe maximum value in your dataset is: ",
max(as.numeric(max.growth))))
}
}
}
# /// check input parameters
if (is(control)!="fl.control") stop("control must be of class fl.control!")
# Check presence of data for chosen fits
if(x_type == "growth" && is.null(growth))
stop("To perform a fits on fluorescence vs. growth data, please provide a 'growth' data matrix of the same dimensions as 'fl_data'.")
if(x_type == "time" && is.null(time))
stop("To perform a fits on fluorescence vs. time data, please provide a 'time' data matrix of the same dimensions as 'fl_data'.")
# /// check number of datasets
if(control$x_type == "growth"){
if ( (dim(growth)[1])!=(dim(fl_data)[1]) ) stop("flFit: Different number of datasets in fl_data and growth")
x <- growth[,-(1:3)]
}
if(control$x_type == "time"){
if ( (dim(time)[1])!=(dim(fl_data)[1]) ) stop("flFit: Different number of datasets in fl_data and time")
x <- time
}
# /// check fitting options
if (!all(control$fit.opt %in% c("s", "l"))){
options(warn=1)
if(control$suppress.messages==F) message("fit.opt must contain 's', and/or 'l'. Changed to c('s', 'l') (both fit methods)!")
fit.opt=c('s', 'l')
options(warn=0)
}
# /// Initialize some parameters
out.table <- NULL
fitnonpara.all <- list()
fitlinear.all <- list()
boot.all <- list()
fitted.param <- NULL
fitted.nonparam <- NULL
bootstrap.param <- NULL
reliability_tag_linear <- NA
reliability_tag_nonpara <- NA
if(control$interactive == FALSE && parallelize == TRUE &&
dim(fl_data)[1] > 30 &&
(
("l" %in% control$fit.opt) || ("a" %in% control$fit.opt) ||
("s" %in% control$fit.opt && control$nboot.fl > 0)
)
){
x.ls <- lapply(1:nrow(x), function(j) x[j, ][!is.na(x[j, ])][!is.na(fl_data[j, -1:-3])])
wells.ls <- lapply(1:nrow(fl_data), function(j) as.numeric(fl_data[j, -1:-3][!is.na(x[j, ])][!is.na(fl_data[j, -1:-3])]))
IDs.ls <- lapply(1:nrow(fl_data), function(j) as.matrix(fl_data[j, 1:3]))
wellnames.ls <- lapply(1:nrow(fl_data), function(j) paste(as.character(fl_data[j,1]), as.character(fl_data[j,2]),as.character(fl_data[j,3]), sep=" | "))
# Set up computing clusters (all available processor cores - 1)
cl <- parallel::makeCluster(parallel::detectCores(all.tests = FALSE, logical = TRUE)-1)
doParallel::registerDoParallel(cl)
# Perform linear fits in parallel
if (("l" %in% control$fit.opt) || ("a" %in% control$fit.opt)){
fitlinear.all <- foreach::foreach(i = 1:dim(fl_data)[1]
) %dopar% {
if(control$x_type == "growth"){
QurvE::flFitLinear(growth = x.ls[[i]], fl_data = wells.ls[[i]], ID = IDs.ls[[i]], control = control)
} else {
QurvE::flFitLinear(time = x.ls[[i]], fl_data = wells.ls[[i]], ID = IDs.ls[[i]], control = control)
}
}
} else {
# /// generate list with empty objects
fitlinear.all <- lapply(1:nrow(fl_data), function(j) list(raw.x = x.ls[[j]],
raw.fl_data = wells.ls[[j]],
filt.x = NA,
filt.fl_data = NA,
log.fl_data = NA,
ID = IDs.ls[[j]],
FUN = NA,
fit = NA,
par = c(y0 = NA, y0_lm = NA, mumax = 0, mu.se = NA, lag = NA, tmax_start = NA, tmax_end = NA,
t_turn = NA, mumax2 = NA, y0_lm2 = NA, lag2 = NA, tmax2_start = NA,
tmax2_end = NA),
ndx = NA, ndx2 = NA,
quota = NA,
rsquared = NA, rsquared2 = NA,
control = control,
fitFlag = FALSE, fitFlag2 = FALSE)
)
}
# Perform spline bootstrappings in parallel
if ((("s" %in% control$fit.opt) || ("a" %in% control$fit.opt) ) &&
(control$nboot.fl > 10) ){
boot.all <- foreach::foreach(i = 1:dim(fl_data)[1]
) %dopar% {
if(control$x_type == "growth"){
QurvE::flBootSpline(growth = x.ls[[i]], fl_data = wells.ls[[i]], ID = IDs.ls[[i]], control = control)
} else {
QurvE::flBootSpline(time = x.ls[[i]], fl_data = wells.ls[[i]], ID = IDs.ls[[i]], control = control)
}
}
}
else{
# /// create empty gcBootSpline object
boot.all <- lapply(1:nrow(fl_data), function(j) list(raw.x=x.ls[[j]],
raw.fl_data=wells.ls[[j]],
ID =IDs.ls[[j]],
boot.x=NA,
boot.y=NA,
boot.gcSpline=NA,
lambda=NA, mu=NA, A=NA, integral=NA,
bootFlag=FALSE, control=control
)
)
}
parallel::stopCluster(cl = cl)
# Assign classes to list elements
for(i in 1:length(fitlinear.all)){
class(fitlinear.all[[i]]) <- "flFitLinear"
}
# for(i in 1:length(fitpara.all)){
# class(fitpara.all[[i]]) <- "gcFitModel"
# }
# for(i in 1:length(fitnonpara.all)){
# class(fitnonpara.all[[i]]) <- "gcFitSpline"
# }
for(i in 1:length(boot.all)){
class(boot.all[[i]]) <- "flBootSpline"
}
}
reliability_tag <- c()
# /// loop over all wells
for (i in 1:dim(fl_data)[1]){
# Progress indicator for shiny app
if(exists("shiny") && shiny == TRUE){
shiny::incProgress(
amount = 1/(dim(fl_data)[1]),
message = "Computations completed")
}
# /// conversion, to handle even data.frame inputs
actx <-
as.numeric(as.matrix(x[i, ]))[!is.na(as.numeric(as.matrix(x[i, ])))][!is.na(as.numeric(as.matrix((fl_data[i, -1:-3]))))]
actwell <-
as.numeric(as.matrix((fl_data[i, -1:-3])))[!is.na(as.numeric(as.matrix(x[i, ])))][!is.na(as.numeric(as.matrix((fl_data[i, -1:-3]))))]
ID <- as.matrix(fl_data[i,1:3])
wellname <- paste(as.character(fl_data[i,1]), as.character(fl_data[i,2]),as.character(fl_data[i,3]), sep=" | ")
if(control$suppress.messages==FALSE){
cat("\n\n")
cat(paste("=== ", as.character(i), ". [", wellname, "] fluorescence curve =================================\n", sep=""))
cat("----------------------------------------------------\n")
}
if(parallelize == FALSE || control$interactive == TRUE ||
dim(fl_data)[1] <= 30 ||
!("l" %in% control$fit.opt || "a" %in% control$fit.opt || ("s" %in% control$fit.opt && control$nboot.fl > 10))
){
# /// Linear regression fl_data
if ("l" %in% control$fit.opt){
if(control$x_type == "growth"){
fitlinear <- flFitLinear(growth = actx, fl_data = actwell, ID = ID, control = control)
} else {
fitlinear <- flFitLinear(time = actx, fl_data = actwell, ID = ID, control = control)
}
fitlinear.all[[i]] <- fitlinear
}
else{
# /// generate empty object
fitlinear <- list(x.in = actx, fl.in = actwell,
raw.x = actx, raw.fl = actwell,
filt.x = actx, filt.fl = actwell,
ID = ID, FUN = grow_exponential, fit = NA, par = c(
y0 = NA, dY= NA, A = NA, y0_lm = NA, max_slope = 0, tD = NA, slope.se = NA, lag = NA, x.max_start = NA, x.max_end = NA,
x.turn = NA, max_slope2 = NA, tD2 = NA, y0_lm2 = NA, lag2 = NA, x.max2_start = NA,
x.max2_end = NA), ndx = NA, ndx.in = NA, ndx2 = NA, ndx2.in = NA, quota = 0.95, rsquared = NA, rsquared2 = NA, control = control, fitFlag = FALSE, fitFlag2 = FALSE)
class(fitlinear) <- "flFitLinear"
fitlinear.all[[i]] <- fitlinear
}
# /// plot linear fit
if ((control$interactive == TRUE)) {
if (("l" %in% control$fit.opt) || ("a" %in% control$fit.opt)) {
answer_satisfied <- "n"
reliability_tag_linear <- NA
while ("n" %in% answer_satisfied) {
try(plot(fitlinear, log = ""))
mtext(side = 3, line = 0, adj = 0,
outer = FALSE,
cex = 1,
wellname)
answer_satisfied <- readline("Are you satisfied with the linear fit (y/n)?\n\n")
if ("n" %in% answer_satisfied) {
test_answer <- readline("Enter: t0, h, quota, min.growth, R2, RSD, tmax, growth.max >>>>\n\n [Skip (enter 'n'), or adjust fit parameters (see ?flFitLinear).\n Leave {blank} at a given position if standard parameters are desired.]\n\n")
if ("n" %in% test_answer) {
if(control$suppress.messages==FALSE){
cat("\n Tagged the linear fit of this sample as unreliable !\n\n")
}
reliability_tag_linear <- FALSE
fitlinear$reliable <- FALSE
fitlinear.all[[i]]$reliable <- FALSE
answer_satisfied <- "y"
} # end if ("n" %in% test_answer)
else {
new_params <- unlist(strsplit(test_answer, split = ","))
t0_new <- ifelse(!is.na(as.numeric(new_params[1])), as.numeric(new_params[1]), control$t0)
h_new <- if(!is.na(as.numeric(new_params[2]))){
as.numeric(new_params[2])
} else {
control$lin.h
}
quota_new <- ifelse(!is.na(as.numeric(new_params[3])), as.numeric(new_params[3]), 0.95)
min.growth_new <- ifelse(!is.na(as.numeric(new_params[4])), as.numeric(new_params[4]), control$min.growth)
R2_new <- ifelse(!is.na(as.numeric(new_params[5])), as.numeric(new_params[5]), control$lin.R2)
RSD_new <- ifelse(!is.na(as.numeric(new_params[6])), as.numeric(new_params[6]), control$lin.RSD)
tmax_new <- ifelse(!is.na(as.numeric(new_params[7])), as.numeric(new_params[7]), control$tmax)
max.growth_new <- ifelse(!is.na(as.numeric(new_params[8])), as.numeric(new_params[8]), control$max.growth)
control_new <- control
control_new$t0 <- t0_new
control_new$lin.h <- h_new
control_new$lin.R2 <- R2_new
control_new$lin.RSD <- RSD_new
control_new$tmax <- tmax_new
control_new$max.growth <- max.growth_new
if(is.numeric(min.growth_new)){
if(!is.na(min.growth_new) && all(as.vector(actwell) < min.growth_new)){
message(paste0("Start growth values need to be greater than 'min.growth'.\nThe minimum start value in your dataset is: ",
min(as.vector(actwell)),". 'min.growth' was not adjusted."), call. = FALSE)
} else if(!is.na(min.growth_new)){
control_new$min.growth <- min.growth_new
}
}
if ("l" %in% control$fit.opt){
if(control$x_type == "growth"){
fitlinear <- flFitLinear(growth = actx, fl_data = actwell, ID = ID, control = control_new, quota = quota_new)
} else {
fitlinear <- flFitLinear(time = actx, fl_data = actwell, ID = ID, control = control_new, quota = quota_new)
}
fitlinear.all[[i]] <- fitlinear
}
fitlinear.all[[i]] <- fitlinear
} #end else
} # end if ("n" %in% test_answer)
else{
reliability_tag_linear <- TRUE
fitlinear$reliable <- TRUE
fitlinear.all[[i]]$reliable <- TRUE
if(control$suppress.messages==FALSE){
cat("Sample was (more or less) o.k.\n")
}
} # end else
} # end while ("n" %in% answer_satisfied)
} # end if (("l" %in% control$fit.opt) || ("a" %in% control$fit.opt))
} # end if ((control$interactive == TRUE))
else {
reliability_tag_linear <- TRUE
fitlinear$reliable <- TRUE
fitlinear.all[[i]]$reliable <- TRUE
}
} # # control$interactive == TRUE || dim(fl_data)[1] <= 30
# /// Non parametric fit
if ("s" %in% control$fit.opt){
if(control$x_type == "growth"){
nonpara <- flFitSpline(growth = actx, fl_data = actwell, ID = ID, control = control)
} else {
nonpara <- flFitSpline(time = actx, fl_data = actwell, ID = ID, control = control)
}
fitnonpara.all[[i]] <- nonpara
}
else{
# /// generate empty object
nonpara <- list(raw.x = actx, raw.fl = actwell,
fit.x = rep(NA, length(actx)), fit.fl = rep(NA, length(actwell)),
parameters = list(A = NA, dY = NA, max_slope = NA, x.max = NA, lambda = NA, b.tangent = NA, max_slope2 = NA, x.max2 = NA,
lambda2 = NA, b.tangent2 = NA, integral = NA),
spline = NA, reliable = NULL, fitFlag = FALSE, fitFlag2 = FALSE,
control = control)
class(nonpara) <- "flFitSpline"
fitnonpara.all[[i]] <- nonpara
}
# /// plotting parametric fit
if ((control$interactive == TRUE)) {
# /// plotting nonparametric fit
if (("s" %in% control$fit.opt) || ("a" %in% control$fit.opt)) {
if (nonpara$fitFlag == TRUE) {
answer_satisfied <- "n"
reliability_tag_nonpara <- NA
while ("n" %in% answer_satisfied) {
plot(nonpara, add=FALSE, raw=TRUE,slope = TRUE, colData=1, cex=1, plot=T, export=F)
answer_satisfied <- readline("Are you satisfied with the spline fit (y/n)?\n\n")
if ("n" %in% answer_satisfied) {
test_answer <- readline("Enter: smooth.fl, t0, min.growth, tmax, max.growth >>>> \n\n [Skip (enter 'n'), or smooth.fl, t0, and min.growth (see ?fl.control).\n Leave {blank} at a given position if standard parameters are desired.]\n\n ")
if ("n" %in% test_answer) {
if(control$suppress.messages==FALSE){
cat("\n Tagged the linear fit of this sample as unreliable !\n\n")
}
reliability_tag_nonpara <- FALSE
nonpara$reliable <- FALSE
fitnonpara.all[[i]]$reliable <- FALSE
fitnonpara.all[[i]]$FitFlag <- FALSE
answer_satisfied <- "y"
} # end if ("n" %in% test_answer)
else{
new_params <- unlist(strsplit(test_answer, split = ","))
if(!is.na(as.numeric(new_params[2])) && as.numeric(new_params[2]) != ""){
t0_new <- as.numeric(new_params[2])
} else {
t0_new <- control$t0
}
smooth.fl_new <- as.numeric(new_params[1])
control_new <- control
if(!is.na(smooth.fl_new) && smooth.fl_new != ""){
control_new$smooth.fl <- smooth.fl_new
}
control_new$t0 <- t0_new
min.growth_new <- as.numeric(new_params[3])
if(!is.na(min.growth_new)){
if(is.numeric(min.growth_new) && min.growth_new != 0 && all(as.vector(actx) < min.growth_new)){
message(paste0("Start growth values need to be below 'min.growth'.\nThe minimum start value in your dataset is: ",
min(as.vector(growth[,4])),". 'min.growth' was not adjusted."), call. = FALSE)
} else if(!is.na(min.growth_new)){
control_new$min.growth <- min.growth_new
}
}
tmax_new <- as.numeric(new_params[4])
if(!is.na(tmax_new) && tmax_new != ""){
control_new$tmax <- tmax_new
}
max.growth_new <- as.numeric(new_params[5])
if(!is.na(max.growth_new) && max.growth_new != ""){
control_new$max.growth <- max.growth_new
}
if(control$x_type == "growth"){
nonpara <- flFitSpline(growth = actx, fl_data = actwell, ID = ID, control = control_new)
} else {
nonpara <- flFitSpline(time = actx, fl_data = actwell, ID = ID, control = control_new)
}
fitnonpara.all[[i]] <- nonpara
} #end else
} # end if ("n" %in% answer_satisfied)
else{
reliability_tag_nonpara <- TRUE
nonpara$reliable <- TRUE
fitnonpara.all[[i]]$reliable <- TRUE
fitnonpara.all[[i]]$FitFlag <- TRUE
if(control$suppress.messages==FALSE){
cat("Sample was (more or less) o.k.\n")
}
} # end else
} # end while ("n" %in% answer_satisfied)
} # end if (nonpara$fitFlag == TRUE)
} # end if (("s" %in% control$fit.opt) || ("a" %in% control$fit.opt) )
} # end of if((control$interactive == TRUE))
else{
reliability_tag_nonpara <- TRUE
}
if(parallelize == FALSE || control$interactive == TRUE ||
dim(fl_data)[1] <= 30 ||
!("l" %in% control$fit.opt || "a" %in% control$fit.opt || ("s" %in% control$fit.opt && control$nboot.fl > 10))
){
# /// Beginn Bootstrap
if ((("s" %in% control$fit.opt) ) &&
(control$nboot.fl > 0) && (reliability_tag_nonpara ==TRUE) && nonpara$fitFlag == TRUE){
if(control$x_type == "growth") bt <- flBootSpline(growth = actx, fl_data = actwell, ID = ID, control = control)
if(control$x_type == "time") bt <- flBootSpline(time = actx, fl_data = actwell, ID = ID, control = control)
boot.all[[i]] <- bt
} # /// end of if (control$nboot.fl ...)
else{
# /// create empty flBootSpline object
bt <- list(raw.x=actx, raw.fl=actwell, ID =ID, boot.x=NA, boot.y=NA, boot.flSpline=NA,
lambda=NA, mu=NA, A=NA, integral=NA, bootFlag=FALSE, control=control)
class(bt) <- "flBootSpline"
boot.all[[i]] <- bt
}
} # if(interactive == TRUE || dim(fl_data)[1] <= 30 ||
reliability_tag <- c(reliability_tag, any(reliability_tag_linear, reliability_tag_nonpara))
# create output table
# description <- data.frame(TestId=fl_data[i,1], AddId=fl_data[i,2],concentration=fl_data[i,3],
# reliability_tag=reliability_tag,
# log.x.spline=control$log.x.spline, log.y.spline=control$log.y.spline,
# log.x.lin=control$log.x.lin, log.y.spline=control$log.y.lin, nboot.fl=control$nboot.fl)
#
# fitted <- cbind(description, summary.flFitLinear(fitlinear), summary.flFitSpline(nonpara), summary.flBootSpline(bt))
#
# out.table <- rbind(out.table, fitted)
# class(out.table) <- c("data.frame", "flTable")
} # /// end of for (i in 1:dim(fl_data)[1])
# Assign names to list elements
names(fitlinear.all) <- names(fitnonpara.all) <- names(boot.all) <- paste0(as.character(fl_data[,1]), " | ", as.character(fl_data[,2]), " | ", as.character(fl_data[,3]))
# create output table
description <- lapply(1:nrow(fl_data), function(x) data.frame(TestId = fl_data[x,1], AddId = fl_data[x,2],concentration = fl_data[x,3],
reliability_tag = reliability_tag[x],
log.x.spline = control$log.x.spline,
log.y.spline = control$log.y.spline,
log.x.lin = control$log.x.lin,
log.y.lin =control$log.y.lin,
nboot.fl = control$nboot.fl
)
)
fitted <- lapply(1:length(fitlinear.all), function(x) cbind(description[[x]],
summary.flFitLinear(fitlinear.all[[x]]),
summary.flFitSpline(fitnonpara.all[[x]]),
summary.flBootSpline(boot.all[[x]])
)
)
df <- data.frame()
out.table <- do.call(rbind, fitted)
class(out.table) <- c("data.frame", "flTable")
flFit <- list(raw.x = x, raw.fl = fl_data, flTable = out.table, flFittedLinear = fitlinear.all, flFittedSplines = fitnonpara.all, flBootSplines = boot.all, control=control)
class(flFit) <- "flFit"
invisible(flFit)
}
#' Run a complete fluorescence curve analysis and dose-reponse analysis workflow.
#'
#' \code{fl.workflow} runs \code{\link{fl.control}} to create a \code{fl.control} object and then performs all computational fitting operations based on the user input. Finally, if desired, a final report is created in PDF or HTML format that summarizes all results obtained.
#'
#' @param grodata A \code{grodata} object created with \code{\link{read_data}} or \code{\link{parse_data}}, containing fluorescence data and data for the independent variable (i.e., time or growth).
#' @param time (optional) A matrix containing time values for each sample (if a \code{fl_data} dataframe is provided as separate argument).
#' @param growth (optional) A dataframe containing growth data (if a \code{fl_data} matrix is provided as separate argument).
#' @param fl_data (optional) A dataframe containing fluorescence data (if a \code{time} matrix or \code{growth} dataframe is provided as separate argument).
#' @param ec50 (Logical) Perform dose-response analysis (\code{TRUE}) or not (\code{FALSE}).
#' @param mean.grp (\code{"all"}, a string vector, or a list of string vectors) Define groups to combine into common plots in the final report based on sample identifiers (if \code{report == TRUE}). Partial matches with sample/group names are accepted. Note: The maximum number of sample groups (with unique condition/concentration indicators) is 50. If you have more than 50 groups, option \code{"all"} will produce the error \code{! Insufficient values in manual scale. [Number] needed but only 50 provided}.
#' @param mean.conc (A numeric vector, or a list of numeric vectors) Define concentrations to combine into common plots in the final report (if \code{report == TRUE}).
#' @param fit.opt (Character or character vector) Indicates whether the program should perform a linear regression (\code{"l"}), model fit (\code{"m"}), spline fit (\code{"s"}), or all (\code{"a"}). Combinations can be freely chosen by providing a character vector, e.g. \code{fit.opt = c("l", "s")} Default: \code{fit.opt = c("l", "s")}.
#' @param x_type (Character) Which data type shall be used as independent variable? Options are \code{'growth'} and \code{'time'}.
#' @param norm_fl (Logical) use normalized (to growth) fluorescence data in fits. Has an effect only when \code{x_type = 'time'}
#' @param t0 (Numeric) Minimum time value considered for linear and spline fits (if \code{x_type = 'time'}).
#' @param tmax (Numeric) Maximum time value considered for linear and spline fits (if \code{x_type = 'time'})..
#' @param min.growth (Numeric) Indicate whether only values above a certain threshold should be considered for linear regressions or spline fits (if \code{x_type = 'growth'}).
#' @param max.growth (Numeric) Indicate whether only growth values below a certain threshold should be considered for linear regressions or spline fits (if \code{x_type = 'growth'}).
#' @param log.x.lin (Logical) Indicates whether _ln(x+1)_ should be applied to the independent variable for _linear_ fits. Default: \code{FALSE}.
#' @param log.x.spline (Logical) Indicates whether _ln(x+1)_ should be applied to the independent variable for _spline_ fits. Default: \code{FALSE}.
#' @param log.y.lin (Logical) Indicates whether _ln(y/y0)_ should be applied to the fluorescence data for _linear_ fits. Default: \code{FALSE}
#' @param log.y.spline (Logical) Indicates whether _ln(y/y0)_ should be applied to the fluorescence data for _spline_ fits. Default: \code{FALSE}
#' @param lin.h (Numeric) Manually define the size of the sliding window used in \code{\link{flFitLinear}}. If \code{NULL}, h is calculated for each samples based on the number of measurements in the fluorescence increase phase of the plot.
#' @param lin.R2 (Numeric) \ifelse{html}{\out{R<sup>2</sup>}}{\eqn{R^2}} threshold for \code{\link{flFitLinear}}.
#' @param lin.RSD (Numeric) Relative standard deviation (RSD) threshold for the calculated slope in \code{\link{flFitLinear}}.
#' @param lin.dY (Numeric) Threshold for the minimum fraction of growth increase a linear regression window should cover. Default: 0.05 (5%).
#' @param biphasic (Logical) Shall \code{\link{flFitLinear}} and \code{\link{flFitSpline}} try to extract fluorescence parameters for two different phases (as observed with, e.g., regulator-promoter systems with varying response in different growth stages) (\code{TRUE}) or not (\code{FALSE})?
#' @param interactive (Logical) Controls whether the fit for each sample and method is controlled manually by the user. If \code{TRUE}, each fit is visualized in the _Plots_ pane and the user can adjust fitting parameters and confirm the reliability of each fit per sample. Default: \code{TRUE}.
#' @param dr.parameter (Character or numeric) The response parameter in the output table to be used for creating a dose response curve. See \code{\link{fl.drFit}} for further details. Default: \code{"max_slope.spline"}, which represents the maximum slope of the spline fit Typical options include: \code{"max_slope.linfit"}, \code{"dY.linfit"}, \code{"max_slope.spline"}, and \code{"dY.spline"}.
#' @param dr.method (Character) Perform either a smooth spline fit on response parameter vs. concentration data (\code{"spline"}) or fit a biosensor response model (proposed by Meyer et al., 2019).
#' @param dr.have.atleast (Numeric) Minimum number of different values for the response parameter one should have for estimating a dose response curve. Note: All fit procedures require at least six unique values. Default: \code{6}.
#' @param smooth.dr (Numeric) Smoothing parameter used in the spline fit by smooth.spline during dose response curve estimation. Usually (not necessesary) in (0; 1]. See \code{\link{smooth.spline}} for further details. Default: \code{NULL}.
#' @param log.x.dr (Logical) Indicates whether \code{ln(x+1)} should be applied to the concentration data of the dose response curves. Default: \code{FALSE}.
#' @param log.y.dr (Logical) Indicates whether \code{ln(y+1)} should be applied to the response data of the dose response curves. Default: \code{FALSE}.
#' @param nboot.dr (Numeric) Defines the number of bootstrap samples for EC50 estimation. Use \code{nboot.dr = 0} to disable bootstrapping. Default: \code{0}.
#' @param nboot.fl (Numeric) Number of bootstrap samples used for nonparametric curve fitting with \code{\link{flBootSpline}}. Use \code{nboot.fl = 0} to disable the bootstrap. Default: \code{0}
#' @param smooth.fl (Numeric) Parameter describing the smoothness of the spline fit; usually (not necessary) within (0;1]. \code{smooth.gc=NULL} causes the program to query an optimal value via cross validation techniques. Especially for datasets with few data points the option \code{NULL} might cause a too small smoothing parameter. This can result a too tight fit that is susceptible to measurement errors (thus overestimating slopes) or produce an error in \code{\link{smooth.spline}} or lead to overfitting. The usage of a fixed value is recommended for reproducible results across samples. See \code{\link{smooth.spline}} for further details. Default: \code{0.55}
#' @param growth.thresh (Numeric) Define a threshold for growth. Only if any growth value in a sample is greater than \code{growth.thresh} (default: 1.5) times the start growth, further computations are performed. Else, a message is returned.
#' @param suppress.messages (Logical) Indicates whether messages (information about current fluorescence curve, EC50 values etc.) should be displayed (\code{FALSE}) or not (\code{TRUE}). This option is meant to speed up the high-throughput processing data. Note: warnings are still displayed. Default: \code{FALSE}.
#' @param neg.nan.act (Logical) Indicates whether the program should stop when negative fluorescence values or NA values appear (\code{TRUE}). Otherwise, the program removes these values silently (\code{FALSE}). Improper values may be caused by incorrect data or input errors. Default: \code{FALSE}.
#' @param clean.bootstrap (Logical) Determines if negative values which occur during bootstrap should be removed (\code{TRUE}) or kept (\code{FALSE}). Note: Infinite values are always removed. Default: \code{TRUE}.
#' @param report (Character or NULL) Create a PDF (\code{'pdf'}) and/or HTML (\code{'html'}) report after running all computations. Define \code{NULL} if no report should be created. Default: (\code{c('pdf', 'html')})
#' @param out.dir {Character or \code{NULL}} Define the name of a folder in which all result files (tables and reports) are stored. If \code{NULL}, the folder will be named with a combination of "FluorescenceResults_" and the current date and time.
#' @param out.nm {Character or \code{NULL}} Define the name of the report files. If \code{NULL}, the files will be named with a combination of "FluorescenceReport_" and the current date and time.
#' @param export.fig (Logical) Export all figures created in the report as separate PNG and PDF files (\code{TRUE}) or not (\code{FALSE}). Only effective if \code{report = TRUE}.
#' @param export.res (Logical) Create tab-separated TXT files containing calculated parameters and dose-response analysis results as well as an .RData file for the resulting `flFitRes` object.
#' @param parallelize Run linear fits and bootstrapping operations in parallel using all but one available processor cores
#' @param ... Further arguments passed to the shiny app.
#'
#' @return A \code{flFitRes} object that contains all computation results, compatible with various plotting functions of the QurvE package and with \code{\link{fl.report}}.
#' \item{time}{Raw time matrix passed to the function as \code{time} (if no \code{grofit} object is provided. Else, extracted from \code{grofit}).}
#' \item{data}{Raw data dataframe passed to the function as \code{grodata}.}
#' \item{flFit}{\code{flFit} object created with the call of \code{\link{flFit}} on fluorescence data.}
#' \item{drFit}{\code{drFit} or \code{drFitfl} object created with the call of \code{\link{growth.drFit}} or \code{\link{fl.drFit}} for fluorescence data (based on the \code{dr.method} argument in \code{control}; see \code{\link{fl.control}}).}
#' \item{expdesign}{Experimental design table inherited from \code{grodata} or created from the identifier columns (columns 1-3) in \code{data}.}
#' \item{control}{Object of class \code{fl.control} created with the call of \code{\link{fl.control}}.}
#'
#' @export
#'
#' @examples
#' # load example dataset
#' input <- read_data(data.growth = system.file("lac_promoters_growth.txt", package = "QurvE"),
#' data.fl = system.file("lac_promoters_fluorescence.txt", package = "QurvE"),
#' csvsep = "\t",
#' csvsep.fl = "\t")
#'
#' # Run workflow
#' res <- fl.workflow(grodata = input, ec50 = FALSE, fit.opt = "s",
#' x_type = "time", norm_fl = TRUE,
#' dr.parameter = "max_slope.spline",
#' suppress.messages = TRUE,
#' parallelize = FALSE)
#'
#' plot(res, data.type = "raw", legend.ncol = 3, basesize = 15)
#'
fl.workflow <- function(grodata = NULL,
time = NULL,
growth = NULL,
fl_data = NULL,
ec50 = TRUE,
mean.grp = NA,
mean.conc = NA,
fit.opt = c("l", "s"),
x_type = c("growth", "time"),
norm_fl = TRUE,
t0 = 0,
tmax = NA,
min.growth = 0,
max.growth = NA,
log.x.lin = FALSE,
log.x.spline = FALSE,
log.y.lin = FALSE,
log.y.spline = FALSE,
lin.h = NULL,
lin.R2 = 0.97,
lin.RSD = 0.07,
lin.dY = 0.05,
biphasic = FALSE,
interactive = FALSE,
dr.parameter = "max_slope.spline",
dr.method = c("model", "spline"),
dr.have.atleast = 5,
smooth.dr = NULL,
log.x.dr = FALSE,
log.y.dr = FALSE,
nboot.dr = 0,
nboot.fl = 0,
smooth.fl = 0.75,
growth.thresh = 1.5,
suppress.messages = FALSE,
neg.nan.act = FALSE,
clean.bootstrap = TRUE,
report = NULL,
out.dir = NULL,
out.nm = NULL,
export.fig = FALSE,
export.res = FALSE,
parallelize = TRUE,
...)
{
if(ec50 == TRUE){
dr.parameter.fit.method <- gsub(".+\\.", "", dr.parameter)
if((dr.parameter.fit.method == "spline" && !any(fit.opt %in% c("s"))) ||
(dr.parameter.fit.method == "linfit" && !any(fit.opt %in% c("l")))
)
message("The chosen 'dr.parameter' is not compatible with the selected fitting options ('fit.opt'). Dose-response analysis will not be performed.")
}
if(exists("lin.h") && !is.null(lin.h) && (is.na(lin.h) || lin.h == "")) lin.h <- NULL
if ( isTRUE(export.fig) && is.null(report) ){
message(
"The export of plots as separate files (`export.fig = TRUE`) is only valid if `report != NULL`."
)
}
# Define objects based on additional function calls
call <- match.call()
## remove strictly defined arguments
call$grodata <- call$time <- call$growth <- call$fl_data <- call$ec50 <- call$mean.grp <- call$mean.conc <- call$neg.nan.act <- call$clean.bootstrap <- call$suppress.messages <- call$export.res <-
call$fit.opt <- call$t0 <- call$min.growth <- call$log.x.lin <- call$log.x.spline <- call$log.y.spline <- call$log.y.lin <- call$biphasic <- call$norm_fl <- call$x_type <- call$tmax <- call$max.growth <-
call$lin.h <- call$lin.R2 <- call$lin.RSD <- call$lin.dY <- call$interactive <- call$nboot.fl <- call$smooth.fl <- call$dr.method <- call$growth.thresh <- call$parallelize <-
call$dr.have.atleast <- call$dr.parameter <- call$smooth.dr <- call$log.x.dr <- call$log.y.dr <- call$nboot.dr <- call$report <- call$out.dir <- call$out.nm <- call$export.fig <- NULL
arglist <- sapply(call, function(x) x)
arglist <- unlist(arglist)[-1]
## Assign additional arguments (...) as R objects
if(length(arglist) > 0){
for(i in 1:length(arglist)){
assign(names(arglist)[i], arglist[[i]])
}
}
x_type <- match.arg(x_type)
dr.method <- match.arg(dr.method)
if(!is.null(grodata) && !(is(grodata)=="list") && !(is(grodata)=="grodata")){
if (is.numeric(as.matrix(time)) == FALSE)
stop("Need a numeric matrix for 'time' or a grodata object created with read_data() or parse_data().")
if (is.numeric(as.matrix(fl_data[-1:-3])) == FALSE)
stop("Need a numeric matrix for 'fl_data' or a grodata object created with read_data() or parse_data().")
if (is.logical(ec50) == FALSE)
stop("Need a logical value for 'ec50'")
} else {
if(!is.null(grodata$time)) time <- grodata$time
if(!is.null(grodata$growth)) growth <- grodata$growth
if(!is.null(grodata$expdesign)) expdesign <- grodata$expdesign
if(!is.null(grodata$fluorescence)) fluorescence <- grodata$fluorescence
# if(!is.null(grodata$fluorescence2)) fluorescence2 <- grodata$fluorescence2
if(!is.null(grodata$norm.fluorescence)) norm.fluorescence <- grodata$norm.fluorescence
# if(!is.null(grodata$norm.fluorescence2)) norm.fluorescence2 <- grodata$norm.fluorescence2
if(!is.null(time)) time <- time
if(!is.null(growth)) growth <- growth
if(!is.null(fl_data)) fluorescence <- fl_data
}
control <- fl.control(fit.opt = fit.opt, norm_fl = norm_fl, x_type = x_type, t0 = t0, min.growth = min.growth, log.x.lin = log.x.lin,
log.x.spline = log.x.spline, log.y.lin = log.y.lin, log.y.spline = log.y.spline, tmax = tmax, max.growth = max.growth,
lin.h = lin.h, lin.R2 = lin.R2, lin.RSD = lin.RSD, lin.dY = lin.dY, dr.have.atleast = dr.have.atleast,
smooth.dr = smooth.dr, log.x.dr = log.x.dr, log.y.dr = log.y.dr, nboot.dr = nboot.dr,
biphasic = biphasic, interactive = interactive, nboot.fl = nboot.fl, dr.parameter = dr.parameter, dr.method = dr.method, clean.bootstrap = clean.bootstrap,
smooth.fl = smooth.fl, growth.thresh = growth.thresh, suppress.messages = suppress.messages, neg.nan.act = neg.nan.act)
nboot.fl <- control$nboot.fl
nboot.dr <- control$nboot.dr
out.flFit <- NA
out.drFit <- NA
# /// fit of fluorescence curves -----------------------------------
if(norm_fl == TRUE && x_type == "time" && (!is.null(norm.fluorescence) && length(norm.fluorescence) > 1 && !all(is.na(norm.fluorescence)))){
if ((control$suppress.messages==FALSE)){
cat("\n\n")
cat(paste("=== Performing Fits for fluorescence =================================\n"))
cat("----------------------------------------------------\n")
}
if(exists("shiny") && shiny == TRUE){
out.flFit <- flFit(time = time, growth = growth, fl_data = norm.fluorescence, control = control, shiny = TRUE, parallelize = parallelize)
} else {
out.flFit <- flFit(time = time, growth = growth, fl_data = norm.fluorescence, control = control, shiny = FALSE, parallelize = parallelize)
}
} else if (!is.null(fluorescence) && length(fluorescence) > 1 && !all(is.na(fluorescence))){
if(exists("shiny") && shiny == TRUE){
out.flFit <- flFit(time = time, growth = growth, fl_data = fluorescence, control = control, shiny = TRUE, parallelize = parallelize)
} else {
out.flFit <- flFit(time = time, growth = growth, fl_data = fluorescence, control = control, shiny = FALSE, parallelize = parallelize)
}
}
# if(norm_fl == TRUE && x_type == "time" && (!is.null(norm.fluorescence2) && length(norm.fluorescence2) > 1 && !all(is.na(norm.fluorescence2)))){
# if ((control$suppress.messages==FALSE)){
# cat("\n\n")
# cat(paste("=== Performing Fits for Fluorescence 2 =================================\n"))
# cat("----------------------------------------------------\n")
# }
# if(exists("shiny") && shiny == TRUE){
# out.flFit2 <- flFit(time = time, growth = growth, fl_data = norm.fluorescence2, control = control, shiny = TRUE)
# } else {
# out.flFit2 <- flFit(time = time, growth = growth, fl_data = norm.fluorescence2, control = control, shiny = FALSE)
# }
# } else if (!is.null(fluorescence2) && length(fluorescence2) > 1 && !all(is.na(fluorescence2))){
# if ((control$suppress.messages==FALSE)){
# cat("\n\n")
# cat(paste("=== Performing Fits for Fluorescence 2 =================================\n"))
# cat("----------------------------------------------------\n")
# }
# if(exists("shiny") && shiny == TRUE){
# out.flFit2 <- flFit(time = time, growth = growth, fl_data = fluorescence2, control = control, shiny = TRUE)
# } else {
# out.flFit2 <- flFit(time = time, growth = growth, fl_data = fluorescence2, control = control, shiny = FALSE)
# }
# }
# /// Estimate EC50 values
if (ec50 == TRUE &&
!((dr.parameter.fit.method == "spline" && !any(fit.opt %in% c("s"))) ||
(dr.parameter.fit.method == "linfit" && !any(fit.opt %in% c("l")))
)
) {
if (!is.null(fluorescence) && length(fluorescence) > 1 && !all(is.na(fluorescence))){
if ((control$suppress.messages==FALSE)){
cat("\n\n")
cat(paste("=== Performing Dose-Response Analysis for fluorescence =================================\n"))
cat("----------------------------------------------------\n")
}
if(control$dr.method == "spline"){
out.drFit <- growth.drFit(summary.flFit(out.flFit), control)
boot.ec1 <- out.drFit$boot.ec
} else {
out.drFit <- fl.drFit(summary.flFit(out.flFit), control)
boot.ec1 <- NA
}
EC50.table1 <- out.drFit$drTable
}
# if (!is.null(fluorescence2) && length(fluorescence2) > 1 && !all(is.na(fluorescence2))){
# if ((control$suppress.messages==FALSE)){
# cat("\n\n")
# cat(paste("=== Performing Dose-Response Analysis for Fluorescence 2 =================================\n"))
# cat("----------------------------------------------------\n")
# }
# if(control$dr.method == "spline"){
# out.drFit2 <- growth.drFit(summary.flFit(out.flFit2), control)
# boot.ec2 <- out.drFit2$boot.ec
# } else {
# out.drFit2 <- fl.drFit(summary.flFit(out.flFit2), control)
# boot.ec2 <- NA
# }
# EC50.table2 <- out.drFit2$drTable
# }
}
# ///
na.obj <- NA
flFitRes <- list(time = time, data = grodata, flFit = get(ifelse(exists("out.flFit"), "out.flFit", "na.obj")),
# flFit2 = get(ifelse(exists("out.flFit2"), "out.flFit2", "na.obj")),
drFit = get(ifelse(exists("out.drFit"), "out.drFit", "na.obj")),
# drFit2 = get(ifelse(exists("out.drFit2"), "out.drFit2", "na.obj")),
expdesign = expdesign, control = control)
class(flFitRes) <- "flFitRes"
if(!exists("shiny") || shiny != TRUE){
if(!is.null(out.dir)){
wd <- paste0(out.dir)
} else {
wd <- paste(getwd(), "/FluorescenceResults_", format(Sys.time(),
"%Y%m%d_%H%M%S"), sep = "")
}
if(export.res)
dir.create(wd, showWarnings = FALSE)
if (!is.null(fluorescence) && length(fluorescence) > 1 && !all(is.na(fluorescence))){
flTable <- data.frame(apply(flFitRes[["flFit"]][["flTable"]],2,as.character))
res.table.fl <- cbind(flTable[,1:3], Filter(function(x) !all(is.na(x)),flTable[,-(1:3)]))
if(export.res){
export_Table(table = res.table.fl, out.dir = wd, out.nm = "results.fl1")
message(paste0("\nResults of fluorescence analysis saved as tab-delimited text file in:\n",
"...", gsub(".+/", "", wd), "/results.fl1.txt\n"))
}
# Export grouped results table
if(("l" %in% control$fit.opt) || ("a" %in% control$fit.opt) ){
table_linear_group <- table_group_fluorescence_linear(res.table.fl)
names <- gsub("<sub>", "_", gsub("</sub>|<sup>|</sup>", "", gsub("<br>", " ", colnames(table_linear_group))))
table_linear_group <- as.data.frame(lapply(1:ncol(table_linear_group), function(x) gsub("<strong>", "", gsub("</strong>", "", table_linear_group[,x]))))
colnames(table_linear_group) <- names
if(export.res)
export_Table(table = table_linear_group, out.dir = wd, out.nm = "grouped_results_fluorescence_linear")
}
if(("s" %in% control$fit.opt) || ("a" %in% control$fit.opt) ){
table_spline_group <- table_group_fluorescence_spline(res.table.fl)
names <- gsub("<sub>", "_", gsub("</sub>|<sup>|</sup>", "", gsub("<br>", " ", colnames(table_spline_group))))
table_spline_group <- as.data.frame(lapply(1:ncol(table_spline_group), function(x) gsub("<strong>", "", gsub("</strong>", "", table_spline_group[,x]))))
colnames(table_spline_group) <- names
if(export.res)
export_Table(table = table_spline_group, out.dir = wd, out.nm = "grouped_results_fluorescence_spline")
}
}
# if (!is.null(fluorescence2) && length(fluorescence2) > 1 && !all(is.na(fluorescence2))){
# flTable2 <- data.frame(apply(flFitRes[["flFit2"]][["flTable"]],2,as.character))
# res.table.fl2 <- cbind(flTable2[,1:3], Filter(function(x) !all(is.na(x)),flTable2[,-(1:3)]))
# if(export.res)
# export_Table(table = res.table.fl2, out.dir = wd, out.nm = "results.fl2")
# cat(paste0("Results of fluorescence 2 analysis saved as tab-delimited text file in:\n",
# wd, "/results.fl2.txt\n"))
# # Export grouped results table
# if(("l" %in% control$fit.opt) || ("a" %in% control$fit.opt) ){
# table_linear_group <- table_group_fluorescence_linear(res.table.fl2)
# names <- gsub("<sub>", "_", gsub("</sub>|<sup>|</sup>", "", gsub("<br>", " ", colnames(table_linear_group))))
# table_linear_group <- as.data.frame(lapply(1:ncol(table_linear_group), function(x) gsub("<strong>", "", gsub("</strong>", "", table_linear_group[,x]))))
# colnames(table_linear_group) <- names
# if(export.res)
# export_Table(table = table_linear_group, out.dir = wd, out.nm = "grouped_results_fluorescence2_linear")
# }
#
# if(("s" %in% control$fit.opt) || ("a" %in% control$fit.opt) ){
# table_spline_group <- table_group_fluorescence_spline(res.table.fl2)
# names <- gsub("<sub>", "_", gsub("</sub>|<sup>|</sup>", "", gsub("<br>", " ", colnames(table_spline_group))))
# table_spline_group <- as.data.frame(lapply(1:ncol(table_spline_group), function(x) gsub("<strong>", "", gsub("</strong>", "", table_spline_group[,x]))))
# colnames(table_spline_group) <- names
# if(export.res)
# export_Table(table = table_spline_group, out.dir = wd, out.nm = "grouped_results_fluorescence2_spline")
# }
# }
if (ec50 == TRUE &&
!((dr.parameter.fit.method == "spline" && !any(fit.opt %in% c("s"))) ||
(dr.parameter.fit.method == "linfit" && !any(fit.opt %in% c("l")))
)
) {
if (!is.null(fluorescence) && length(fluorescence) > 1 && !all(is.na(fluorescence))){
if(!is.null(EC50.table1) && length(EC50.table1) > 1) {
res.table.dr_fl1 <- Filter(function(x) !all(is.na(x)),EC50.table1)
if(export.res){
export_Table(table = res.table.dr_fl1, out.dir = wd, out.nm = "results.fl_dr1")
message(paste0("\nResults of EC50 analysis for fluorescence saved as tab-delimited in:\n",
"...", gsub(".+/", "", wd), "/results.fl_dr1.txt\n"))
}
}
}
# if (!is.null(fluorescence2) && length(fluorescence2) > 1 && !all(is.na(fluorescence2))){
# if(!is.null(EC50.table2) && length(EC50.table2) > 1) {
# res.table.dr_fl2 <- Filter(function(x) !all(is.na(x)),EC50.table2)
# if(export.res)
# export_Table(table = res.table.dr_fl2, out.dir = wd, out.nm = "results.fl_dr2")
#
# cat(paste0("Results of EC50 analysis for fluorescence 2 saved as tab-delimited in:\n",
# wd, "/results.fl_dr2.txt\n"))
# }
# }
} else {
res.table.dr_fl1 <- NULL
res.table.dr_fl2 <- NULL
}
# Export RData object
if(export.res)
export_RData(flFitRes, out.dir = wd)
if(any(report %in% c('pdf', 'html'))){
try(fl.report(flFitRes, out.dir = gsub(paste0(getwd(), "/"), "", wd), mean.grp = mean.grp, mean.conc = mean.conc, ec50 = ec50,
export = export.fig, format = report, out.nm = out.nm, parallelize = parallelize))
}
}
invisible(flFitRes)
}
Try the QurvE package in your browser
Any scripts or data that you put into this service are public.
QurvE documentation built on May 29, 2024, 3 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions fl.workflow flFit
Documented in flFit fl.workflow
#' Perform a fluorescence curve analysis on all samples in the provided dataset.
#'
#' \code{flFit} performs all computational fluorescence fitting operations based on the user input.
#'
#' @param fl_data Either... \enumerate{ \item a \code{grodata} object created with \code{\link{read_data}} or \code{\link{parse_data}},
#' \item a list containing a \code{'time'} matrix (for x_type == "time") or \code{'growth'} dataframe (for x_type == "growth") and a \code{'fluorescence'} dataframes,
#' or \item a dataframe containing (normalized) fluorescence values (if a \code{time} matrix or \code{growth} dataframe is provided as separate argument).}
#' @param time (optional) A matrix containing time values for each sample.
#' @param growth (optional) A dataframe containing growth values for each sample and sample identifiers in the first three columns.
#' @param control A \code{fl.control} object created with \code{\link{fl.control}}, defining relevant fitting options.
#' @param parallelize Run linear fits and bootstrapping operations in parallel using all but one available processor cores
#' @param ... Further arguments passed to the shiny app.
#'
#' @return An \code{flFit} object that contains all fluorescence fitting results, compatible with
#' various plotting functions of the QurvE package.
#' \item{raw.x}{Raw x matrix passed to the function as \code{time} (for x_type = 'time') or \code{growth} (for x_type = 'growth').}
#' \item{raw.fl}{Raw growth dataframe passed to the function as \code{data}.}
#' \item{flTable}{Table with fluorescence parameters and related statistics for each fluorescence curve evaluation performed by the function. This table, which is also returned by the generic \code{summary.flFit} method applied to a \code{flFit} object, is used as an input for \code{\link{fl.drFit}}.}
#' \item{flFittedLinear}{List of all \code{flFitLinear} objects, generated by the call of \code{\link{flFitLinear}}. Note: access to each object in the list via double brace: flFittedLinear\[\[#n\]\].}
#' \item{flFittedSplines}{List of all \code{flFitSpline} objects, generated by the call of \code{\link{flFitSpline}}. Note: access to each object via double brace: flFittedSplines\[\[#n\]\].}
#' \item{flBootSplines}{List of all \code{flBootSpline} objects, generated by the call of \code{\link{flBootSpline}}. Note: access to each object via double brace: flFittedSplines\[\[#n\]\].}
#' \item{control}{Object of class \code{fl.control} containing list of options passed to the function as \code{control}.}
#'
#' @details
#' Common response parameters used in dose-response analysis:<br><br><b>Linear fit:</b><br>- max_slope.linfit: Fluorescence increase rate<br>- lambda.linfit: Lag time<br>- dY.linfit: Maximum Fluorescence - Minimum Fluorescence<br>- A.linfit: Maximum fluorescence<br><br><b>Spline fit:</b><br>- max_slope.spline: Fluorescence increase rate<br>- lambda.spline: Lag time<br>- dY.spline: Maximum Fluorescence - Minimum Fluorescence<br>- A.spline: Maximum fluorescence<br>- integral.spline: Integral<br><br><b>Parametric fit:</b><br>- max_slope.model: Fluorescence increase rate<br>- lambda.model: Lag time<br>- dY.model: Maximum Fluorescence - Minimum Fluorescence<br>- A.model: Maximum fluorescence<br>- integral.model: Integral'
#'
#' @family workflows
#' @family fluorescence fitting functions
#' @family dose-response analysis functions
#'
#' @export
#'
#' @importFrom foreach %dopar%
#'
#' @examples
#' # load example dataset
#' input <- read_data(data.growth = system.file("lac_promoters_growth.txt", package = "QurvE"),
#' data.fl = system.file("lac_promoters_fluorescence.txt", package = "QurvE"),
#' csvsep = "\t",
#' csvsep.fl = "\t" )
#'
#' # Define fit controls
#' control <- fl.control(fit.opt = "s",
#' x_type = "time", norm_fl = TRUE,
#' dr.parameter = "max_slope.spline",
#' dr.method = "model",
#' suppress.messages = TRUE)
#'
#' # Run curve fitting workflow
#' res <- flFit(fl_data = input$norm.fluorescence,
#' time = input$time,
#' control = control,
#' parallelize = FALSE)
#'
#' summary(res)
#'
flFit <- function(fl_data, time = NULL, growth = NULL, control= fl.control(), parallelize = TRUE, ...)
{
old.options <- options()
on.exit(options(old.options))
# Define objects based on additional function calls
call <- match.call()
## remove strictly defined arguments
call$time <- call$growth <- call$fl_data <- call$control<- call$parallelize <- NULL
arglist <- sapply(call, function(x) x)
arglist <- unlist(arglist)[-1]
## Assign additional arguments (...) as R objects
if(length(arglist) > 0){
for(i in 1:length(arglist)){
assign(names(arglist)[i], arglist[[i]])
}
}
x_type <- control$x_type
if(!(class(fl_data) %in% c("list", "grodata"))){
if (x_type == "time" && is.numeric(as.matrix(time)) == FALSE)
stop("Need a numeric matrix for 'time' (for x_type = 'time') or a grodata object created with read_data() or parse_data() in the 'fl_data' argument.")
if (x_type == "growth" && is.numeric(as.matrix(growth[-1:-3])) == FALSE)
stop("Need a dataframe for 'growth' (for x_type = 'growth') or a grodata object created with read_data() or parse_data() in the 'fl_data' argument.")
if (is.numeric(as.matrix(fl_data[-1:-3])) == FALSE)
stop("Need a dataframe for 'fl_data' or a grodata object created with read_data() or parse_data() in the 'fl_data' argument.")
} else if(!is.null(fl_data)){
time <- fl_data$time
growth <- fl_data$growth
fl_data <- fl_data$fluorescence
}
# /// check if start growth values are above min.growth in all samples
if(!is.null(growth) && length(growth) > 1){
max.growth <- unlist(lapply(1:nrow(growth), function (x) max(as.numeric(as.matrix(growth[x,-1:-3]))[!is.na(as.numeric(as.matrix(growth[x,-1:-3])))])))
if(is.numeric(control$min.growth) && control$min.growth != 0){
if(!is.na(control$min.growth) && all(as.numeric(max.growth) < control$min.growth)){
stop(paste0("The chosen global start growth value (min.growth) is larger than every value in your dataset.\nThe maximum value in your dataset is: ",
max(as.numeric(max.growth))))
}
}
}
# /// check input parameters
if (is(control)!="fl.control") stop("control must be of class fl.control!")
# Check presence of data for chosen fits
if(x_type == "growth" && is.null(growth))
stop("To perform a fits on fluorescence vs. growth data, please provide a 'growth' data matrix of the same dimensions as 'fl_data'.")
if(x_type == "time" && is.null(time))
stop("To perform a fits on fluorescence vs. time data, please provide a 'time' data matrix of the same dimensions as 'fl_data'.")
# /// check number of datasets
if(control$x_type == "growth"){
if ( (dim(growth)[1])!=(dim(fl_data)[1]) ) stop("flFit: Different number of datasets in fl_data and growth")
x <- growth[,-(1:3)]
}
if(control$x_type == "time"){
if ( (dim(time)[1])!=(dim(fl_data)[1]) ) stop("flFit: Different number of datasets in fl_data and time")
x <- time
}
# /// check fitting options
if (!all(control$fit.opt %in% c("s", "l"))){
options(warn=1)
if(control$suppress.messages==F) message("fit.opt must contain 's', and/or 'l'. Changed to c('s', 'l') (both fit methods)!")
fit.opt=c('s', 'l')
options(warn=0)
}
# /// Initialize some parameters
out.table <- NULL
fitnonpara.all <- list()
fitlinear.all <- list()
boot.all <- list()
fitted.param <- NULL
fitted.nonparam <- NULL
bootstrap.param <- NULL
reliability_tag_linear <- NA
reliability_tag_nonpara <- NA
if(control$interactive == FALSE && parallelize == TRUE &&
dim(fl_data)[1] > 30 &&
(
("l" %in% control$fit.opt) || ("a" %in% control$fit.opt) ||
("s" %in% control$fit.opt && control$nboot.fl > 0)
)
){
x.ls <- lapply(1:nrow(x), function(j) x[j, ][!is.na(x[j, ])][!is.na(fl_data[j, -1:-3])])
wells.ls <- lapply(1:nrow(fl_data), function(j) as.numeric(fl_data[j, -1:-3][!is.na(x[j, ])][!is.na(fl_data[j, -1:-3])]))
IDs.ls <- lapply(1:nrow(fl_data), function(j) as.matrix(fl_data[j, 1:3]))
wellnames.ls <- lapply(1:nrow(fl_data), function(j) paste(as.character(fl_data[j,1]), as.character(fl_data[j,2]),as.character(fl_data[j,3]), sep=" | "))
# Set up computing clusters (all available processor cores - 1)
cl <- parallel::makeCluster(parallel::detectCores(all.tests = FALSE, logical = TRUE)-1)
doParallel::registerDoParallel(cl)
# Perform linear fits in parallel
if (("l" %in% control$fit.opt) || ("a" %in% control$fit.opt)){
fitlinear.all <- foreach::foreach(i = 1:dim(fl_data)[1]
) %dopar% {
if(control$x_type == "growth"){
QurvE::flFitLinear(growth = x.ls[[i]], fl_data = wells.ls[[i]], ID = IDs.ls[[i]], control = control)
} else {
QurvE::flFitLinear(time = x.ls[[i]], fl_data = wells.ls[[i]], ID = IDs.ls[[i]], control = control)
}
}
} else {
# /// generate list with empty objects
fitlinear.all <- lapply(1:nrow(fl_data), function(j) list(raw.x = x.ls[[j]],
raw.fl_data = wells.ls[[j]],
filt.x = NA,
filt.fl_data = NA,
log.fl_data = NA,
ID = IDs.ls[[j]],
FUN = NA,
fit = NA,
par = c(y0 = NA, y0_lm = NA, mumax = 0, mu.se = NA, lag = NA, tmax_start = NA, tmax_end = NA,
t_turn = NA, mumax2 = NA, y0_lm2 = NA, lag2 = NA, tmax2_start = NA,
tmax2_end = NA),
ndx = NA, ndx2 = NA,
quota = NA,
rsquared = NA, rsquared2 = NA,
control = control,
fitFlag = FALSE, fitFlag2 = FALSE)
)
}
# Perform spline bootstrappings in parallel
if ((("s" %in% control$fit.opt) || ("a" %in% control$fit.opt) ) &&
(control$nboot.fl > 10) ){
boot.all <- foreach::foreach(i = 1:dim(fl_data)[1]
) %dopar% {
if(control$x_type == "growth"){
QurvE::flBootSpline(growth = x.ls[[i]], fl_data = wells.ls[[i]], ID = IDs.ls[[i]], control = control)
} else {
QurvE::flBootSpline(time = x.ls[[i]], fl_data = wells.ls[[i]], ID = IDs.ls[[i]], control = control)
}
}
}
else{
# /// create empty gcBootSpline object
boot.all <- lapply(1:nrow(fl_data), function(j) list(raw.x=x.ls[[j]],
raw.fl_data=wells.ls[[j]],
ID =IDs.ls[[j]],
boot.x=NA,
boot.y=NA,
boot.gcSpline=NA,
lambda=NA, mu=NA, A=NA, integral=NA,
bootFlag=FALSE, control=control
)
)
}
parallel::stopCluster(cl = cl)
# Assign classes to list elements
for(i in 1:length(fitlinear.all)){
class(fitlinear.all[[i]]) <- "flFitLinear"
}
# for(i in 1:length(fitpara.all)){
# class(fitpara.all[[i]]) <- "gcFitModel"
# }
# for(i in 1:length(fitnonpara.all)){
# class(fitnonpara.all[[i]]) <- "gcFitSpline"
# }
for(i in 1:length(boot.all)){
class(boot.all[[i]]) <- "flBootSpline"
}
}
reliability_tag <- c()
# /// loop over all wells
for (i in 1:dim(fl_data)[1]){
# Progress indicator for shiny app
if(exists("shiny") && shiny == TRUE){
shiny::incProgress(
amount = 1/(dim(fl_data)[1]),
message = "Computations completed")
}
# /// conversion, to handle even data.frame inputs
actx <-
as.numeric(as.matrix(x[i, ]))[!is.na(as.numeric(as.matrix(x[i, ])))][!is.na(as.numeric(as.matrix((fl_data[i, -1:-3]))))]
actwell <-
as.numeric(as.matrix((fl_data[i, -1:-3])))[!is.na(as.numeric(as.matrix(x[i, ])))][!is.na(as.numeric(as.matrix((fl_data[i, -1:-3]))))]
ID <- as.matrix(fl_data[i,1:3])
wellname <- paste(as.character(fl_data[i,1]), as.character(fl_data[i,2]),as.character(fl_data[i,3]), sep=" | ")
if(control$suppress.messages==FALSE){
cat("\n\n")
cat(paste("=== ", as.character(i), ". [", wellname, "] fluorescence curve =================================\n", sep=""))
cat("----------------------------------------------------\n")
}
if(parallelize == FALSE || control$interactive == TRUE ||
dim(fl_data)[1] <= 30 ||
!("l" %in% control$fit.opt || "a" %in% control$fit.opt || ("s" %in% control$fit.opt && control$nboot.fl > 10))
){
# /// Linear regression fl_data
if ("l" %in% control$fit.opt){
if(control$x_type == "growth"){
fitlinear <- flFitLinear(growth = actx, fl_data = actwell, ID = ID, control = control)
} else {
fitlinear <- flFitLinear(time = actx, fl_data = actwell, ID = ID, control = control)
}
fitlinear.all[[i]] <- fitlinear
}
else{
# /// generate empty object
fitlinear <- list(x.in = actx, fl.in = actwell,
raw.x = actx, raw.fl = actwell,
filt.x = actx, filt.fl = actwell,
ID = ID, FUN = grow_exponential, fit = NA, par = c(
y0 = NA, dY= NA, A = NA, y0_lm = NA, max_slope = 0, tD = NA, slope.se = NA, lag = NA, x.max_start = NA, x.max_end = NA,
x.turn = NA, max_slope2 = NA, tD2 = NA, y0_lm2 = NA, lag2 = NA, x.max2_start = NA,
x.max2_end = NA), ndx = NA, ndx.in = NA, ndx2 = NA, ndx2.in = NA, quota = 0.95, rsquared = NA, rsquared2 = NA, control = control, fitFlag = FALSE, fitFlag2 = FALSE)
class(fitlinear) <- "flFitLinear"
fitlinear.all[[i]] <- fitlinear
}
# /// plot linear fit
if ((control$interactive == TRUE)) {
if (("l" %in% control$fit.opt) || ("a" %in% control$fit.opt)) {
answer_satisfied <- "n"
reliability_tag_linear <- NA
while ("n" %in% answer_satisfied) {
try(plot(fitlinear, log = ""))
mtext(side = 3, line = 0, adj = 0,
outer = FALSE,
cex = 1,
wellname)
answer_satisfied <- readline("Are you satisfied with the linear fit (y/n)?\n\n")
if ("n" %in% answer_satisfied) {
test_answer <- readline("Enter: t0, h, quota, min.growth, R2, RSD, tmax, growth.max >>>>\n\n [Skip (enter 'n'), or adjust fit parameters (see ?flFitLinear).\n Leave {blank} at a given position if standard parameters are desired.]\n\n")
if ("n" %in% test_answer) {
if(control$suppress.messages==FALSE){
cat("\n Tagged the linear fit of this sample as unreliable !\n\n")
}
reliability_tag_linear <- FALSE
fitlinear$reliable <- FALSE
fitlinear.all[[i]]$reliable <- FALSE
answer_satisfied <- "y"
} # end if ("n" %in% test_answer)
else {
new_params <- unlist(strsplit(test_answer, split = ","))
t0_new <- ifelse(!is.na(as.numeric(new_params[1])), as.numeric(new_params[1]), control$t0)
h_new <- if(!is.na(as.numeric(new_params[2]))){
as.numeric(new_params[2])
} else {
control$lin.h
}
quota_new <- ifelse(!is.na(as.numeric(new_params[3])), as.numeric(new_params[3]), 0.95)
min.growth_new <- ifelse(!is.na(as.numeric(new_params[4])), as.numeric(new_params[4]), control$min.growth)
R2_new <- ifelse(!is.na(as.numeric(new_params[5])), as.numeric(new_params[5]), control$lin.R2)
RSD_new <- ifelse(!is.na(as.numeric(new_params[6])), as.numeric(new_params[6]), control$lin.RSD)
tmax_new <- ifelse(!is.na(as.numeric(new_params[7])), as.numeric(new_params[7]), control$tmax)
max.growth_new <- ifelse(!is.na(as.numeric(new_params[8])), as.numeric(new_params[8]), control$max.growth)
control_new <- control
control_new$t0 <- t0_new
control_new$lin.h <- h_new
control_new$lin.R2 <- R2_new
control_new$lin.RSD <- RSD_new
control_new$tmax <- tmax_new
control_new$max.growth <- max.growth_new
if(is.numeric(min.growth_new)){
if(!is.na(min.growth_new) && all(as.vector(actwell) < min.growth_new)){
message(paste0("Start growth values need to be greater than 'min.growth'.\nThe minimum start value in your dataset is: ",
min(as.vector(actwell)),". 'min.growth' was not adjusted."), call. = FALSE)
} else if(!is.na(min.growth_new)){
control_new$min.growth <- min.growth_new
}
}
if ("l" %in% control$fit.opt){
if(control$x_type == "growth"){
fitlinear <- flFitLinear(growth = actx, fl_data = actwell, ID = ID, control = control_new, quota = quota_new)
} else {
fitlinear <- flFitLinear(time = actx, fl_data = actwell, ID = ID, control = control_new, quota = quota_new)
}
fitlinear.all[[i]] <- fitlinear
}
fitlinear.all[[i]] <- fitlinear
} #end else
} # end if ("n" %in% test_answer)
else{
reliability_tag_linear <- TRUE
fitlinear$reliable <- TRUE
fitlinear.all[[i]]$reliable <- TRUE
if(control$suppress.messages==FALSE){
cat("Sample was (more or less) o.k.\n")
}
} # end else
} # end while ("n" %in% answer_satisfied)
} # end if (("l" %in% control$fit.opt) || ("a" %in% control$fit.opt))
} # end if ((control$interactive == TRUE))
else {
reliability_tag_linear <- TRUE
fitlinear$reliable <- TRUE
fitlinear.all[[i]]$reliable <- TRUE
}
} # # control$interactive == TRUE || dim(fl_data)[1] <= 30
# /// Non parametric fit
if ("s" %in% control$fit.opt){
if(control$x_type == "growth"){
nonpara <- flFitSpline(growth = actx, fl_data = actwell, ID = ID, control = control)
} else {
nonpara <- flFitSpline(time = actx, fl_data = actwell, ID = ID, control = control)
}
fitnonpara.all[[i]] <- nonpara
}
else{
# /// generate empty object
nonpara <- list(raw.x = actx, raw.fl = actwell,
fit.x = rep(NA, length(actx)), fit.fl = rep(NA, length(actwell)),
parameters = list(A = NA, dY = NA, max_slope = NA, x.max = NA, lambda = NA, b.tangent = NA, max_slope2 = NA, x.max2 = NA,
lambda2 = NA, b.tangent2 = NA, integral = NA),
spline = NA, reliable = NULL, fitFlag = FALSE, fitFlag2 = FALSE,
control = control)
class(nonpara) <- "flFitSpline"
fitnonpara.all[[i]] <- nonpara
}
# /// plotting parametric fit
if ((control$interactive == TRUE)) {
# /// plotting nonparametric fit
if (("s" %in% control$fit.opt) || ("a" %in% control$fit.opt)) {
if (nonpara$fitFlag == TRUE) {
answer_satisfied <- "n"
reliability_tag_nonpara <- NA
while ("n" %in% answer_satisfied) {
plot(nonpara, add=FALSE, raw=TRUE,slope = TRUE, colData=1, cex=1, plot=T, export=F)
answer_satisfied <- readline("Are you satisfied with the spline fit (y/n)?\n\n")
if ("n" %in% answer_satisfied) {
test_answer <- readline("Enter: smooth.fl, t0, min.growth, tmax, max.growth >>>> \n\n [Skip (enter 'n'), or smooth.fl, t0, and min.growth (see ?fl.control).\n Leave {blank} at a given position if standard parameters are desired.]\n\n ")
if ("n" %in% test_answer) {
if(control$suppress.messages==FALSE){
cat("\n Tagged the linear fit of this sample as unreliable !\n\n")
}
reliability_tag_nonpara <- FALSE
nonpara$reliable <- FALSE
fitnonpara.all[[i]]$reliable <- FALSE
fitnonpara.all[[i]]$FitFlag <- FALSE
answer_satisfied <- "y"
} # end if ("n" %in% test_answer)
else{
new_params <- unlist(strsplit(test_answer, split = ","))
if(!is.na(as.numeric(new_params[2])) && as.numeric(new_params[2]) != ""){
t0_new <- as.numeric(new_params[2])
} else {
t0_new <- control$t0
}
smooth.fl_new <- as.numeric(new_params[1])
control_new <- control
if(!is.na(smooth.fl_new) && smooth.fl_new != ""){
control_new$smooth.fl <- smooth.fl_new
}
control_new$t0 <- t0_new
min.growth_new <- as.numeric(new_params[3])
if(!is.na(min.growth_new)){
if(is.numeric(min.growth_new) && min.growth_new != 0 && all(as.vector(actx) < min.growth_new)){
message(paste0("Start growth values need to be below 'min.growth'.\nThe minimum start value in your dataset is: ",
min(as.vector(growth[,4])),". 'min.growth' was not adjusted."), call. = FALSE)
} else if(!is.na(min.growth_new)){
control_new$min.growth <- min.growth_new
}
}
tmax_new <- as.numeric(new_params[4])
if(!is.na(tmax_new) && tmax_new != ""){
control_new$tmax <- tmax_new
}
max.growth_new <- as.numeric(new_params[5])
if(!is.na(max.growth_new) && max.growth_new != ""){
control_new$max.growth <- max.growth_new
}
if(control$x_type == "growth"){
nonpara <- flFitSpline(growth = actx, fl_data = actwell, ID = ID, control = control_new)
} else {
nonpara <- flFitSpline(time = actx, fl_data = actwell, ID = ID, control = control_new)
}
fitnonpara.all[[i]] <- nonpara
} #end else
} # end if ("n" %in% answer_satisfied)
else{
reliability_tag_nonpara <- TRUE
nonpara$reliable <- TRUE
fitnonpara.all[[i]]$reliable <- TRUE
fitnonpara.all[[i]]$FitFlag <- TRUE
if(control$suppress.messages==FALSE){
cat("Sample was (more or less) o.k.\n")
}
} # end else
} # end while ("n" %in% answer_satisfied)
} # end if (nonpara$fitFlag == TRUE)
} # end if (("s" %in% control$fit.opt) || ("a" %in% control$fit.opt) )
} # end of if((control$interactive == TRUE))
else{
reliability_tag_nonpara <- TRUE
}
if(parallelize == FALSE || control$interactive == TRUE ||
dim(fl_data)[1] <= 30 ||
!("l" %in% control$fit.opt || "a" %in% control$fit.opt || ("s" %in% control$fit.opt && control$nboot.fl > 10))
){
# /// Beginn Bootstrap
if ((("s" %in% control$fit.opt) ) &&
(control$nboot.fl > 0) && (reliability_tag_nonpara ==TRUE) && nonpara$fitFlag == TRUE){
if(control$x_type == "growth") bt <- flBootSpline(growth = actx, fl_data = actwell, ID = ID, control = control)
if(control$x_type == "time") bt <- flBootSpline(time = actx, fl_data = actwell, ID = ID, control = control)
boot.all[[i]] <- bt
} # /// end of if (control$nboot.fl ...)
else{
# /// create empty flBootSpline object
bt <- list(raw.x=actx, raw.fl=actwell, ID =ID, boot.x=NA, boot.y=NA, boot.flSpline=NA,
lambda=NA, mu=NA, A=NA, integral=NA, bootFlag=FALSE, control=control)
class(bt) <- "flBootSpline"
boot.all[[i]] <- bt
}
} # if(interactive == TRUE || dim(fl_data)[1] <= 30 ||
reliability_tag <- c(reliability_tag, any(reliability_tag_linear, reliability_tag_nonpara))
# create output table
# description <- data.frame(TestId=fl_data[i,1], AddId=fl_data[i,2],concentration=fl_data[i,3],
# reliability_tag=reliability_tag,
# log.x.spline=control$log.x.spline, log.y.spline=control$log.y.spline,
# log.x.lin=control$log.x.lin, log.y.spline=control$log.y.lin, nboot.fl=control$nboot.fl)
#
# fitted <- cbind(description, summary.flFitLinear(fitlinear), summary.flFitSpline(nonpara), summary.flBootSpline(bt))
#
# out.table <- rbind(out.table, fitted)
# class(out.table) <- c("data.frame", "flTable")
} # /// end of for (i in 1:dim(fl_data)[1])
# Assign names to list elements
names(fitlinear.all) <- names(fitnonpara.all) <- names(boot.all) <- paste0(as.character(fl_data[,1]), " | ", as.character(fl_data[,2]), " | ", as.character(fl_data[,3]))
# create output table
description <- lapply(1:nrow(fl_data), function(x) data.frame(TestId = fl_data[x,1], AddId = fl_data[x,2],concentration = fl_data[x,3],
reliability_tag = reliability_tag[x],
log.x.spline = control$log.x.spline,
log.y.spline = control$log.y.spline,
log.x.lin = control$log.x.lin,
log.y.lin =control$log.y.lin,
nboot.fl = control$nboot.fl
)
)
fitted <- lapply(1:length(fitlinear.all), function(x) cbind(description[[x]],
summary.flFitLinear(fitlinear.all[[x]]),
summary.flFitSpline(fitnonpara.all[[x]]),
summary.flBootSpline(boot.all[[x]])
)
)
df <- data.frame()
out.table <- do.call(rbind, fitted)
class(out.table) <- c("data.frame", "flTable")
flFit <- list(raw.x = x, raw.fl = fl_data, flTable = out.table, flFittedLinear = fitlinear.all, flFittedSplines = fitnonpara.all, flBootSplines = boot.all, control=control)
class(flFit) <- "flFit"
invisible(flFit)
}
#' Run a complete fluorescence curve analysis and dose-reponse analysis workflow.
#'
#' \code{fl.workflow} runs \code{\link{fl.control}} to create a \code{fl.control} object and then performs all computational fitting operations based on the user input. Finally, if desired, a final report is created in PDF or HTML format that summarizes all results obtained.
#'
#' @param grodata A \code{grodata} object created with \code{\link{read_data}} or \code{\link{parse_data}}, containing fluorescence data and data for the independent variable (i.e., time or growth).
#' @param time (optional) A matrix containing time values for each sample (if a \code{fl_data} dataframe is provided as separate argument).
#' @param growth (optional) A dataframe containing growth data (if a \code{fl_data} matrix is provided as separate argument).
#' @param fl_data (optional) A dataframe containing fluorescence data (if a \code{time} matrix or \code{growth} dataframe is provided as separate argument).
#' @param ec50 (Logical) Perform dose-response analysis (\code{TRUE}) or not (\code{FALSE}).
#' @param mean.grp (\code{"all"}, a string vector, or a list of string vectors) Define groups to combine into common plots in the final report based on sample identifiers (if \code{report == TRUE}). Partial matches with sample/group names are accepted. Note: The maximum number of sample groups (with unique condition/concentration indicators) is 50. If you have more than 50 groups, option \code{"all"} will produce the error \code{! Insufficient values in manual scale. [Number] needed but only 50 provided}.
#' @param mean.conc (A numeric vector, or a list of numeric vectors) Define concentrations to combine into common plots in the final report (if \code{report == TRUE}).
#' @param fit.opt (Character or character vector) Indicates whether the program should perform a linear regression (\code{"l"}), model fit (\code{"m"}), spline fit (\code{"s"}), or all (\code{"a"}). Combinations can be freely chosen by providing a character vector, e.g. \code{fit.opt = c("l", "s")} Default: \code{fit.opt = c("l", "s")}.
#' @param x_type (Character) Which data type shall be used as independent variable? Options are \code{'growth'} and \code{'time'}.
#' @param norm_fl (Logical) use normalized (to growth) fluorescence data in fits. Has an effect only when \code{x_type = 'time'}
#' @param t0 (Numeric) Minimum time value considered for linear and spline fits (if \code{x_type = 'time'}).
#' @param tmax (Numeric) Maximum time value considered for linear and spline fits (if \code{x_type = 'time'})..
#' @param min.growth (Numeric) Indicate whether only values above a certain threshold should be considered for linear regressions or spline fits (if \code{x_type = 'growth'}).
#' @param max.growth (Numeric) Indicate whether only growth values below a certain threshold should be considered for linear regressions or spline fits (if \code{x_type = 'growth'}).
#' @param log.x.lin (Logical) Indicates whether _ln(x+1)_ should be applied to the independent variable for _linear_ fits. Default: \code{FALSE}.
#' @param log.x.spline (Logical) Indicates whether _ln(x+1)_ should be applied to the independent variable for _spline_ fits. Default: \code{FALSE}.
#' @param log.y.lin (Logical) Indicates whether _ln(y/y0)_ should be applied to the fluorescence data for _linear_ fits. Default: \code{FALSE}
#' @param log.y.spline (Logical) Indicates whether _ln(y/y0)_ should be applied to the fluorescence data for _spline_ fits. Default: \code{FALSE}
#' @param lin.h (Numeric) Manually define the size of the sliding window used in \code{\link{flFitLinear}}. If \code{NULL}, h is calculated for each samples based on the number of measurements in the fluorescence increase phase of the plot.
#' @param lin.R2 (Numeric) \ifelse{html}{\out{R<sup>2</sup>}}{\eqn{R^2}} threshold for \code{\link{flFitLinear}}.
#' @param lin.RSD (Numeric) Relative standard deviation (RSD) threshold for the calculated slope in \code{\link{flFitLinear}}.
#' @param lin.dY (Numeric) Threshold for the minimum fraction of growth increase a linear regression window should cover. Default: 0.05 (5%).
#' @param biphasic (Logical) Shall \code{\link{flFitLinear}} and \code{\link{flFitSpline}} try to extract fluorescence parameters for two different phases (as observed with, e.g., regulator-promoter systems with varying response in different growth stages) (\code{TRUE}) or not (\code{FALSE})?
#' @param interactive (Logical) Controls whether the fit for each sample and method is controlled manually by the user. If \code{TRUE}, each fit is visualized in the _Plots_ pane and the user can adjust fitting parameters and confirm the reliability of each fit per sample. Default: \code{TRUE}.
#' @param dr.parameter (Character or numeric) The response parameter in the output table to be used for creating a dose response curve. See \code{\link{fl.drFit}} for further details. Default: \code{"max_slope.spline"}, which represents the maximum slope of the spline fit Typical options include: \code{"max_slope.linfit"}, \code{"dY.linfit"}, \code{"max_slope.spline"}, and \code{"dY.spline"}.
#' @param dr.method (Character) Perform either a smooth spline fit on response parameter vs. concentration data (\code{"spline"}) or fit a biosensor response model (proposed by Meyer et al., 2019).
#' @param dr.have.atleast (Numeric) Minimum number of different values for the response parameter one should have for estimating a dose response curve. Note: All fit procedures require at least six unique values. Default: \code{6}.
#' @param smooth.dr (Numeric) Smoothing parameter used in the spline fit by smooth.spline during dose response curve estimation. Usually (not necessesary) in (0; 1]. See \code{\link{smooth.spline}} for further details. Default: \code{NULL}.
#' @param log.x.dr (Logical) Indicates whether \code{ln(x+1)} should be applied to the concentration data of the dose response curves. Default: \code{FALSE}.
#' @param log.y.dr (Logical) Indicates whether \code{ln(y+1)} should be applied to the response data of the dose response curves. Default: \code{FALSE}.
#' @param nboot.dr (Numeric) Defines the number of bootstrap samples for EC50 estimation. Use \code{nboot.dr = 0} to disable bootstrapping. Default: \code{0}.
#' @param nboot.fl (Numeric) Number of bootstrap samples used for nonparametric curve fitting with \code{\link{flBootSpline}}. Use \code{nboot.fl = 0} to disable the bootstrap. Default: \code{0}
#' @param smooth.fl (Numeric) Parameter describing the smoothness of the spline fit; usually (not necessary) within (0;1]. \code{smooth.gc=NULL} causes the program to query an optimal value via cross validation techniques. Especially for datasets with few data points the option \code{NULL} might cause a too small smoothing parameter. This can result a too tight fit that is susceptible to measurement errors (thus overestimating slopes) or produce an error in \code{\link{smooth.spline}} or lead to overfitting. The usage of a fixed value is recommended for reproducible results across samples. See \code{\link{smooth.spline}} for further details. Default: \code{0.55}
#' @param growth.thresh (Numeric) Define a threshold for growth. Only if any growth value in a sample is greater than \code{growth.thresh} (default: 1.5) times the start growth, further computations are performed. Else, a message is returned.
#' @param suppress.messages (Logical) Indicates whether messages (information about current fluorescence curve, EC50 values etc.) should be displayed (\code{FALSE}) or not (\code{TRUE}). This option is meant to speed up the high-throughput processing data. Note: warnings are still displayed. Default: \code{FALSE}.
#' @param neg.nan.act (Logical) Indicates whether the program should stop when negative fluorescence values or NA values appear (\code{TRUE}). Otherwise, the program removes these values silently (\code{FALSE}). Improper values may be caused by incorrect data or input errors. Default: \code{FALSE}.
#' @param clean.bootstrap (Logical) Determines if negative values which occur during bootstrap should be removed (\code{TRUE}) or kept (\code{FALSE}). Note: Infinite values are always removed. Default: \code{TRUE}.
#' @param report (Character or NULL) Create a PDF (\code{'pdf'}) and/or HTML (\code{'html'}) report after running all computations. Define \code{NULL} if no report should be created. Default: (\code{c('pdf', 'html')})
#' @param out.dir {Character or \code{NULL}} Define the name of a folder in which all result files (tables and reports) are stored. If \code{NULL}, the folder will be named with a combination of "FluorescenceResults_" and the current date and time.
#' @param out.nm {Character or \code{NULL}} Define the name of the report files. If \code{NULL}, the files will be named with a combination of "FluorescenceReport_" and the current date and time.
#' @param export.fig (Logical) Export all figures created in the report as separate PNG and PDF files (\code{TRUE}) or not (\code{FALSE}). Only effective if \code{report = TRUE}.
#' @param export.res (Logical) Create tab-separated TXT files containing calculated parameters and dose-response analysis results as well as an .RData file for the resulting `flFitRes` object.
#' @param parallelize Run linear fits and bootstrapping operations in parallel using all but one available processor cores
#' @param ... Further arguments passed to the shiny app.
#'
#' @return A \code{flFitRes} object that contains all computation results, compatible with various plotting functions of the QurvE package and with \code{\link{fl.report}}.
#' \item{time}{Raw time matrix passed to the function as \code{time} (if no \code{grofit} object is provided. Else, extracted from \code{grofit}).}
#' \item{data}{Raw data dataframe passed to the function as \code{grodata}.}
#' \item{flFit}{\code{flFit} object created with the call of \code{\link{flFit}} on fluorescence data.}
#' \item{drFit}{\code{drFit} or \code{drFitfl} object created with the call of \code{\link{growth.drFit}} or \code{\link{fl.drFit}} for fluorescence data (based on the \code{dr.method} argument in \code{control}; see \code{\link{fl.control}}).}
#' \item{expdesign}{Experimental design table inherited from \code{grodata} or created from the identifier columns (columns 1-3) in \code{data}.}
#' \item{control}{Object of class \code{fl.control} created with the call of \code{\link{fl.control}}.}
#'
#' @export
#'
#' @examples
#' # load example dataset
#' input <- read_data(data.growth = system.file("lac_promoters_growth.txt", package = "QurvE"),
#' data.fl = system.file("lac_promoters_fluorescence.txt", package = "QurvE"),
#' csvsep = "\t",
#' csvsep.fl = "\t")
#'
#' # Run workflow
#' res <- fl.workflow(grodata = input, ec50 = FALSE, fit.opt = "s",
#' x_type = "time", norm_fl = TRUE,
#' dr.parameter = "max_slope.spline",
#' suppress.messages = TRUE,
#' parallelize = FALSE)
#'
#' plot(res, data.type = "raw", legend.ncol = 3, basesize = 15)
#'
fl.workflow <- function(grodata = NULL,
time = NULL,
growth = NULL,
fl_data = NULL,
ec50 = TRUE,
mean.grp = NA,
mean.conc = NA,
fit.opt = c("l", "s"),
x_type = c("growth", "time"),
norm_fl = TRUE,
t0 = 0,
tmax = NA,
min.growth = 0,
max.growth = NA,
log.x.lin = FALSE,
log.x.spline = FALSE,
log.y.lin = FALSE,
log.y.spline = FALSE,
lin.h = NULL,
lin.R2 = 0.97,
lin.RSD = 0.07,
lin.dY = 0.05,
biphasic = FALSE,
interactive = FALSE,
dr.parameter = "max_slope.spline",
dr.method = c("model", "spline"),
dr.have.atleast = 5,
smooth.dr = NULL,
log.x.dr = FALSE,
log.y.dr = FALSE,
nboot.dr = 0,
nboot.fl = 0,
smooth.fl = 0.75,
growth.thresh = 1.5,
suppress.messages = FALSE,
neg.nan.act = FALSE,
clean.bootstrap = TRUE,
report = NULL,
out.dir = NULL,
out.nm = NULL,
export.fig = FALSE,
export.res = FALSE,
parallelize = TRUE,
...)
{
if(ec50 == TRUE){
dr.parameter.fit.method <- gsub(".+\\.", "", dr.parameter)
if((dr.parameter.fit.method == "spline" && !any(fit.opt %in% c("s"))) ||
(dr.parameter.fit.method == "linfit" && !any(fit.opt %in% c("l")))
)
message("The chosen 'dr.parameter' is not compatible with the selected fitting options ('fit.opt'). Dose-response analysis will not be performed.")
}
if(exists("lin.h") && !is.null(lin.h) && (is.na(lin.h) || lin.h == "")) lin.h <- NULL
if ( isTRUE(export.fig) && is.null(report) ){
message(
"The export of plots as separate files (`export.fig = TRUE`) is only valid if `report != NULL`."
)
}
# Define objects based on additional function calls
call <- match.call()
## remove strictly defined arguments
call$grodata <- call$time <- call$growth <- call$fl_data <- call$ec50 <- call$mean.grp <- call$mean.conc <- call$neg.nan.act <- call$clean.bootstrap <- call$suppress.messages <- call$export.res <-
call$fit.opt <- call$t0 <- call$min.growth <- call$log.x.lin <- call$log.x.spline <- call$log.y.spline <- call$log.y.lin <- call$biphasic <- call$norm_fl <- call$x_type <- call$tmax <- call$max.growth <-
call$lin.h <- call$lin.R2 <- call$lin.RSD <- call$lin.dY <- call$interactive <- call$nboot.fl <- call$smooth.fl <- call$dr.method <- call$growth.thresh <- call$parallelize <-
call$dr.have.atleast <- call$dr.parameter <- call$smooth.dr <- call$log.x.dr <- call$log.y.dr <- call$nboot.dr <- call$report <- call$out.dir <- call$out.nm <- call$export.fig <- NULL
arglist <- sapply(call, function(x) x)
arglist <- unlist(arglist)[-1]
## Assign additional arguments (...) as R objects
if(length(arglist) > 0){
for(i in 1:length(arglist)){
assign(names(arglist)[i], arglist[[i]])
}
}
x_type <- match.arg(x_type)
dr.method <- match.arg(dr.method)
if(!is.null(grodata) && !(is(grodata)=="list") && !(is(grodata)=="grodata")){
if (is.numeric(as.matrix(time)) == FALSE)
stop("Need a numeric matrix for 'time' or a grodata object created with read_data() or parse_data().")
if (is.numeric(as.matrix(fl_data[-1:-3])) == FALSE)
stop("Need a numeric matrix for 'fl_data' or a grodata object created with read_data() or parse_data().")
if (is.logical(ec50) == FALSE)
stop("Need a logical value for 'ec50'")
} else {
if(!is.null(grodata$time)) time <- grodata$time
if(!is.null(grodata$growth)) growth <- grodata$growth
if(!is.null(grodata$expdesign)) expdesign <- grodata$expdesign
if(!is.null(grodata$fluorescence)) fluorescence <- grodata$fluorescence
# if(!is.null(grodata$fluorescence2)) fluorescence2 <- grodata$fluorescence2
if(!is.null(grodata$norm.fluorescence)) norm.fluorescence <- grodata$norm.fluorescence
# if(!is.null(grodata$norm.fluorescence2)) norm.fluorescence2 <- grodata$norm.fluorescence2
if(!is.null(time)) time <- time
if(!is.null(growth)) growth <- growth
if(!is.null(fl_data)) fluorescence <- fl_data
}
control <- fl.control(fit.opt = fit.opt, norm_fl = norm_fl, x_type = x_type, t0 = t0, min.growth = min.growth, log.x.lin = log.x.lin,
log.x.spline = log.x.spline, log.y.lin = log.y.lin, log.y.spline = log.y.spline, tmax = tmax, max.growth = max.growth,
lin.h = lin.h, lin.R2 = lin.R2, lin.RSD = lin.RSD, lin.dY = lin.dY, dr.have.atleast = dr.have.atleast,
smooth.dr = smooth.dr, log.x.dr = log.x.dr, log.y.dr = log.y.dr, nboot.dr = nboot.dr,
biphasic = biphasic, interactive = interactive, nboot.fl = nboot.fl, dr.parameter = dr.parameter, dr.method = dr.method, clean.bootstrap = clean.bootstrap,
smooth.fl = smooth.fl, growth.thresh = growth.thresh, suppress.messages = suppress.messages, neg.nan.act = neg.nan.act)
nboot.fl <- control$nboot.fl
nboot.dr <- control$nboot.dr
out.flFit <- NA
out.drFit <- NA
# /// fit of fluorescence curves -----------------------------------
if(norm_fl == TRUE && x_type == "time" && (!is.null(norm.fluorescence) && length(norm.fluorescence) > 1 && !all(is.na(norm.fluorescence)))){
if ((control$suppress.messages==FALSE)){
cat("\n\n")
cat(paste("=== Performing Fits for fluorescence =================================\n"))
cat("----------------------------------------------------\n")
}
if(exists("shiny") && shiny == TRUE){
out.flFit <- flFit(time = time, growth = growth, fl_data = norm.fluorescence, control = control, shiny = TRUE, parallelize = parallelize)
} else {
out.flFit <- flFit(time = time, growth = growth, fl_data = norm.fluorescence, control = control, shiny = FALSE, parallelize = parallelize)
}
} else if (!is.null(fluorescence) && length(fluorescence) > 1 && !all(is.na(fluorescence))){
if(exists("shiny") && shiny == TRUE){
out.flFit <- flFit(time = time, growth = growth, fl_data = fluorescence, control = control, shiny = TRUE, parallelize = parallelize)
} else {
out.flFit <- flFit(time = time, growth = growth, fl_data = fluorescence, control = control, shiny = FALSE, parallelize = parallelize)
}
}
# if(norm_fl == TRUE && x_type == "time" && (!is.null(norm.fluorescence2) && length(norm.fluorescence2) > 1 && !all(is.na(norm.fluorescence2)))){
# if ((control$suppress.messages==FALSE)){
# cat("\n\n")
# cat(paste("=== Performing Fits for Fluorescence 2 =================================\n"))
# cat("----------------------------------------------------\n")
# }
# if(exists("shiny") && shiny == TRUE){
# out.flFit2 <- flFit(time = time, growth = growth, fl_data = norm.fluorescence2, control = control, shiny = TRUE)
# } else {
# out.flFit2 <- flFit(time = time, growth = growth, fl_data = norm.fluorescence2, control = control, shiny = FALSE)
# }
# } else if (!is.null(fluorescence2) && length(fluorescence2) > 1 && !all(is.na(fluorescence2))){
# if ((control$suppress.messages==FALSE)){
# cat("\n\n")
# cat(paste("=== Performing Fits for Fluorescence 2 =================================\n"))
# cat("----------------------------------------------------\n")
# }
# if(exists("shiny") && shiny == TRUE){
# out.flFit2 <- flFit(time = time, growth = growth, fl_data = fluorescence2, control = control, shiny = TRUE)
# } else {
# out.flFit2 <- flFit(time = time, growth = growth, fl_data = fluorescence2, control = control, shiny = FALSE)
# }
# }
# /// Estimate EC50 values
if (ec50 == TRUE &&
!((dr.parameter.fit.method == "spline" && !any(fit.opt %in% c("s"))) ||
(dr.parameter.fit.method == "linfit" && !any(fit.opt %in% c("l")))
)
) {
if (!is.null(fluorescence) && length(fluorescence) > 1 && !all(is.na(fluorescence))){
if ((control$suppress.messages==FALSE)){
cat("\n\n")
cat(paste("=== Performing Dose-Response Analysis for fluorescence =================================\n"))
cat("----------------------------------------------------\n")
}
if(control$dr.method == "spline"){
out.drFit <- growth.drFit(summary.flFit(out.flFit), control)
boot.ec1 <- out.drFit$boot.ec
} else {
out.drFit <- fl.drFit(summary.flFit(out.flFit), control)
boot.ec1 <- NA
}
EC50.table1 <- out.drFit$drTable
}
# if (!is.null(fluorescence2) && length(fluorescence2) > 1 && !all(is.na(fluorescence2))){
# if ((control$suppress.messages==FALSE)){
# cat("\n\n")
# cat(paste("=== Performing Dose-Response Analysis for Fluorescence 2 =================================\n"))
# cat("----------------------------------------------------\n")
# }
# if(control$dr.method == "spline"){
# out.drFit2 <- growth.drFit(summary.flFit(out.flFit2), control)
# boot.ec2 <- out.drFit2$boot.ec
# } else {
# out.drFit2 <- fl.drFit(summary.flFit(out.flFit2), control)
# boot.ec2 <- NA
# }
# EC50.table2 <- out.drFit2$drTable
# }
}
# ///
na.obj <- NA
flFitRes <- list(time = time, data = grodata, flFit = get(ifelse(exists("out.flFit"), "out.flFit", "na.obj")),
# flFit2 = get(ifelse(exists("out.flFit2"), "out.flFit2", "na.obj")),
drFit = get(ifelse(exists("out.drFit"), "out.drFit", "na.obj")),
# drFit2 = get(ifelse(exists("out.drFit2"), "out.drFit2", "na.obj")),
expdesign = expdesign, control = control)
class(flFitRes) <- "flFitRes"
if(!exists("shiny") || shiny != TRUE){
if(!is.null(out.dir)){
wd <- paste0(out.dir)
} else {
wd <- paste(getwd(), "/FluorescenceResults_", format(Sys.time(),
"%Y%m%d_%H%M%S"), sep = "")
}
if(export.res)
dir.create(wd, showWarnings = FALSE)
if (!is.null(fluorescence) && length(fluorescence) > 1 && !all(is.na(fluorescence))){
flTable <- data.frame(apply(flFitRes[["flFit"]][["flTable"]],2,as.character))
res.table.fl <- cbind(flTable[,1:3], Filter(function(x) !all(is.na(x)),flTable[,-(1:3)]))
if(export.res){
export_Table(table = res.table.fl, out.dir = wd, out.nm = "results.fl1")
message(paste0("\nResults of fluorescence analysis saved as tab-delimited text file in:\n",
"...", gsub(".+/", "", wd), "/results.fl1.txt\n"))
}
# Export grouped results table
if(("l" %in% control$fit.opt) || ("a" %in% control$fit.opt) ){
table_linear_group <- table_group_fluorescence_linear(res.table.fl)
names <- gsub("<sub>", "_", gsub("</sub>|<sup>|</sup>", "", gsub("<br>", " ", colnames(table_linear_group))))
table_linear_group <- as.data.frame(lapply(1:ncol(table_linear_group), function(x) gsub("<strong>", "", gsub("</strong>", "", table_linear_group[,x]))))
colnames(table_linear_group) <- names
if(export.res)
export_Table(table = table_linear_group, out.dir = wd, out.nm = "grouped_results_fluorescence_linear")
}
if(("s" %in% control$fit.opt) || ("a" %in% control$fit.opt) ){
table_spline_group <- table_group_fluorescence_spline(res.table.fl)
names <- gsub("<sub>", "_", gsub("</sub>|<sup>|</sup>", "", gsub("<br>", " ", colnames(table_spline_group))))
table_spline_group <- as.data.frame(lapply(1:ncol(table_spline_group), function(x) gsub("<strong>", "", gsub("</strong>", "", table_spline_group[,x]))))
colnames(table_spline_group) <- names
if(export.res)
export_Table(table = table_spline_group, out.dir = wd, out.nm = "grouped_results_fluorescence_spline")
}
}
# if (!is.null(fluorescence2) && length(fluorescence2) > 1 && !all(is.na(fluorescence2))){
# flTable2 <- data.frame(apply(flFitRes[["flFit2"]][["flTable"]],2,as.character))
# res.table.fl2 <- cbind(flTable2[,1:3], Filter(function(x) !all(is.na(x)),flTable2[,-(1:3)]))
# if(export.res)
# export_Table(table = res.table.fl2, out.dir = wd, out.nm = "results.fl2")
# cat(paste0("Results of fluorescence 2 analysis saved as tab-delimited text file in:\n",
# wd, "/results.fl2.txt\n"))
# # Export grouped results table
# if(("l" %in% control$fit.opt) || ("a" %in% control$fit.opt) ){
# table_linear_group <- table_group_fluorescence_linear(res.table.fl2)
# names <- gsub("<sub>", "_", gsub("</sub>|<sup>|</sup>", "", gsub("<br>", " ", colnames(table_linear_group))))
# table_linear_group <- as.data.frame(lapply(1:ncol(table_linear_group), function(x) gsub("<strong>", "", gsub("</strong>", "", table_linear_group[,x]))))
# colnames(table_linear_group) <- names
# if(export.res)
# export_Table(table = table_linear_group, out.dir = wd, out.nm = "grouped_results_fluorescence2_linear")
# }
#
# if(("s" %in% control$fit.opt) || ("a" %in% control$fit.opt) ){
# table_spline_group <- table_group_fluorescence_spline(res.table.fl2)
# names <- gsub("<sub>", "_", gsub("</sub>|<sup>|</sup>", "", gsub("<br>", " ", colnames(table_spline_group))))
# table_spline_group <- as.data.frame(lapply(1:ncol(table_spline_group), function(x) gsub("<strong>", "", gsub("</strong>", "", table_spline_group[,x]))))
# colnames(table_spline_group) <- names
# if(export.res)
# export_Table(table = table_spline_group, out.dir = wd, out.nm = "grouped_results_fluorescence2_spline")
# }
# }
if (ec50 == TRUE &&
!((dr.parameter.fit.method == "spline" && !any(fit.opt %in% c("s"))) ||
(dr.parameter.fit.method == "linfit" && !any(fit.opt %in% c("l")))
)
) {
if (!is.null(fluorescence) && length(fluorescence) > 1 && !all(is.na(fluorescence))){
if(!is.null(EC50.table1) && length(EC50.table1) > 1) {
res.table.dr_fl1 <- Filter(function(x) !all(is.na(x)),EC50.table1)
if(export.res){
export_Table(table = res.table.dr_fl1, out.dir = wd, out.nm = "results.fl_dr1")
message(paste0("\nResults of EC50 analysis for fluorescence saved as tab-delimited in:\n",
"...", gsub(".+/", "", wd), "/results.fl_dr1.txt\n"))
}
}
}
# if (!is.null(fluorescence2) && length(fluorescence2) > 1 && !all(is.na(fluorescence2))){
# if(!is.null(EC50.table2) && length(EC50.table2) > 1) {
# res.table.dr_fl2 <- Filter(function(x) !all(is.na(x)),EC50.table2)
# if(export.res)
# export_Table(table = res.table.dr_fl2, out.dir = wd, out.nm = "results.fl_dr2")
#
# cat(paste0("Results of EC50 analysis for fluorescence 2 saved as tab-delimited in:\n",
# wd, "/results.fl_dr2.txt\n"))
# }
# }
} else {
res.table.dr_fl1 <- NULL
res.table.dr_fl2 <- NULL
}
# Export RData object
if(export.res)
export_RData(flFitRes, out.dir = wd)
if(any(report %in% c('pdf', 'html'))){
try(fl.report(flFitRes, out.dir = gsub(paste0(getwd(), "/"), "", wd), mean.grp = mean.grp, mean.conc = mean.conc, ec50 = ec50,
export = export.fig, format = report, out.nm = out.nm, parallelize = parallelize))
}
}
invisible(flFitRes)
}
Try the QurvE package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.