diff options
Diffstat (limited to 'R/checkexperiment.R')
| -rw-r--r-- | R/checkexperiment.R | 371 |
1 files changed, 185 insertions, 186 deletions
diff --git a/R/checkexperiment.R b/R/checkexperiment.R index 2806317..7fa1a60 100644 --- a/R/checkexperiment.R +++ b/R/checkexperiment.R @@ -1,190 +1,189 @@ -if(getRversion() >= '2.15.1') utils::globalVariables(c("type", "conc", "substance")) -checkexperiment <- function(id, db = "ecotox", endpoint = "%") -{ - databases <- data.frame( - responsename=c("viability","activity","raw_response"), - testtype=c("celltype","enzyme","organism"), - exptype=c("plate","plate","experiment")) - rownames(databases) <- c("cytotox","enzymes","ecotox") - - if (!(db %in% rownames(databases))) stop("Database is not supported") - - if (requireNamespace("RODBC")) { - channel <- RODBC::odbcConnect(db, uid="cytotox", pwd="cytotox", case="tolower") - } else { - stop("For this function, the RODBC package has to be installed and configured.") - } - - responsename = as.character(databases[db,1]) - testtype = as.character(databases[db,2]) - exptype = as.character(databases[db,3]) - - exptable <- paste(exptype, "s", sep="") - commentquery <- paste("SELECT comment FROM ", exptable , - " WHERE ", exptype, " = ", id) - commentdata <- RODBC::sqlQuery(channel,commentquery) - comment <- as.character(commentdata[[1]]) - - expquery <- paste("SELECT experimentator, substance, ", - testtype, ", conc, unit,", responsename, ", type, raw_0, duration, performed, ok", - " FROM ",db," WHERE ",exptype,"=", id, - sep = "") - - if (db == "ecotox") { - expquery <- paste0(expquery, " AND type LIKE '", endpoint, "'") - } - - expdata <- RODBC::sqlQuery(channel,expquery) - - if (db %in% c("cytotox","enzymes")) { - controlquery <- paste("SELECT type,response FROM controls - WHERE plate=",id) - controldata <- RODBC::sqlQuery(channel,controlquery) - } - - RODBC::odbcClose(channel) - - op <- par(ask=TRUE) - on.exit(par(op)) - - if (db %in% c("cytotox","enzymes")) { - blinds <- subset(controldata, type == "blind") - controls <- subset(controldata, type == "control") - QA <- matrix(nrow = 2, ncol = 4, - dimnames = list(c("Blind", "Control (conc = 0)"), - c("Number", "Mean", "Std. Dev.", "% Std. Dev"))) - - QA[1, 1] <- length(blinds$response) - QA[1, 2] <- signif(mean(blinds$response), 2) - QA[1, 3] <- signif(sd(blinds$response), 2) - QA[1, 4] <-signif(QA[1, 3] * 100 / QA[1, 2],2) - } else { - # Use raw response for ecotox - expdata$response <- expdata$raw_response - - if (nlevels(expdata$type) > 1) { - message("There are data for more than one type of raw response in your data.\n", - "The types are ", paste(levels(expdata$type), collapse = " and "), ".\n", - "You should choose one of these types using 'endpoint = \"$type\"'", - "in your call to checkexperiment\n", - "For now, we are continuing with the data for ", levels(expdata$type)[1]) - } - endpoint <- expdata$type[1] - expdata <- subset(expdata, type == endpoint) - - controls <- subset(expdata, conc == 0) - expdata <- subset(expdata, conc != 0) - - QA <- matrix(nrow = 1, ncol = 4, - dimnames = list(c("Control (conc = 0)"), - c("Number", "Mean", "Std. Dev.", "% Std. Dev"))) - - } - - numberOfControls <- length(controls$response) - QA["Control (conc = 0)", 1] <- numberOfControls - if (numberOfControls > 0) { - QA["Control (conc = 0)", 2] <- signif(mean(controls$response),2) - QA["Control (conc = 0)", 3] <- signif(sd(controls$response),2) - QA["Control (conc = 0)", 4] <- signif(QA["Control (conc = 0)", 3] * 100 / - QA["Control (conc = 0)", 2],2) - } +utils::globalVariables(c("type", "conc", "substance")) - if (db == "ecotox") { - if (identical(as.character(levels(expdata$organism)), "Vibrio fischeri")) { - positive <- subset(expdata, substance == "Na Cl") - if (nrow(positive) > 0) { - QA <- rbind(QA, - c(nrow(positive), - signif(mean(positive$raw_response), 2), - signif(sd(positive$raw_response), 2), - signif(100 * sd(positive$raw_response) / - mean(positive$raw_response), 2))) - - rownames(QA) <- c("Control (conc = 0)", - "Positive control (Na Cl)") - } - expdata <- subset(expdata, substance != "Na Cl", drop = TRUE) - } - } - - if (length(expdata$experimentator) < 1) { - stop("There is no response data for ",exptype," ", - id," in database ",db,"\n") - } - exptypestring <- paste(toupper(substring(exptype,1,1)), - substring(exptype,2),sep="") - expdata$experimentator <- factor(expdata$experimentator) - expdata$type <- factor(expdata[[testtype]]) - expdata$performed <- factor(as.character(expdata$performed)) - expdata$substance <- factor(expdata$substance) - expdata$unit <- factor(expdata$unit) - expdata$ok <- factor(expdata$ok) - - # Info on the experiment - cat("\n",exptypestring,id,"from database",db,":\n\n", - "\tExperimentator(s):\t",levels(expdata$experimentator),"\n", - "\tType(s):\t\t",levels(expdata$type),"\n", - "\tPerformed on:\t\t",levels(expdata$performed),"\n", - "\tSubstance(s):\t\t",levels(expdata$substance),"\n", - "\tConcentration unit(s):\t",levels(expdata$unit),"\n", - "\tComment:\t\t",comment,"\n", - "\tOK Levels:\t\t",levels(expdata$ok),"\n\n") - - print(QA) - - # Control growth rate for Lemna and algae - if (endpoint %in% c("cell count", "frond area", "frond number")) { - duration <- unique(expdata$duration) # in hours - if (length(duration) > 1) stop("More than one duration in the data") - response_0 <- unique(expdata$raw_0) - if (length(response_0) > 1) stop("More than one mean response at time 0 in the data") - t_days <- duration / 24 - control_growth_rates <- (log(controls$response) - log(response_0)) / t_days - cat("\nMean growth rate in controls:\t", round(mean(control_growth_rates), 3), "per day\n") - } - - - # Box plot of control data - if (db == "ecotox") { - boxplot(controls$response, - names="controls", - ylab=endpoint, - ylim=range(controls$response, na.rm = TRUE), - boxwex=0.4, - main=paste("Plate ",id)) - } else { - boxplot(blinds$response,controls$response, - names=c("blinds","controls"), - ylab="Response", - boxwex=0.4, - main=paste("Plate ",id)) +checkexperiment <- function(id, + db = c("ecotox", "cytotox", "enzymes"), + endpoint = "%") +{ + db = match.arg(db) + + databases <- data.frame( + responsename = c("viability", "activity", "raw_response"), + testtype = c("celltype", "enzyme", "organism"), + exptype = c("plate", "plate", "experiment"), + row.names = c("cytotox", "enzymes", "ecotox"), + stringsAsFactors = FALSE) + + con <- dbConnect(odbc(), "cytotox", database = db) + + responsename <- databases[db, 1] + testtype <- databases[db, 2] + exptype <- databases[db, 3] + + exptable <- paste(exptype, "s", sep = "") + + commentquery <- paste0("SELECT comment FROM ", exptable, " ", + "WHERE ", exptype, " = ", id) + commentdata <- dbGetQuery(con, commentquery) + comment <- as.character(commentdata[[1]]) + + expquery <- paste0("SELECT ", + "experimentator, substance, ", testtype, ", conc, unit,", responsename, ", + type, raw_0, duration, performed, ok ", + "FROM ", db, " ", + "WHERE ", exptype, "=", id) + + if (db == "ecotox") { + expquery <- paste0(expquery, " AND type LIKE '", endpoint, "'") + } + + expdata <- dbGetQuery(con, expquery) + + if (db %in% c("cytotox", "enzymes")) { + controlquery <- paste0("SELECT type, response FROM controls ", + " WHERE plate=", id) + controldata <- dbGetQuery(con, controlquery) + } + + op <- par(ask=TRUE) + on.exit(par(op)) + + if (db %in% c("cytotox","enzymes")) { + blinds <- subset(controldata, type == "blind") + controls <- subset(controldata, type == "control") + QA <- matrix(nrow = 2, ncol = 4, + dimnames = list(c("Blind", "Control (conc = 0)"), + c("Number", "Mean", "Std. Dev.", "% Std. Dev"))) + + QA[1, 1] <- length(blinds$response) + QA[1, 2] <- signif(mean(blinds$response), 2) + QA[1, 3] <- signif(sd(blinds$response), 2) + QA[1, 4] <- signif(QA[1, 3] * 100 / QA[1, 2], 2) + } else { + # Use raw response for ecotox + expdata$response <- expdata$raw_response + + if (nlevels(expdata$type) > 1) { + message("There are data for more than one type of raw response in your data.\n", + "The types are ", paste(levels(expdata$type), collapse = " and "), ".\n", + "You should choose one of these types using 'endpoint = \"$type\"'", + "in your call to checkexperiment\n", + "For now, we are continuing with the data for ", levels(expdata$type)[1]) } - - # Plot of dose response data - drdata <- expdata[c(2,4,6)] - drdata$substance <- factor(drdata$substance) - substances <- levels(drdata$substance) - - lld <- log10(min(subset(drdata,conc!=0)$conc)) - lhd <- log10(max(drdata$conc)) - - ylab <- if (db == "ecotox") endpoint - else responsename - - plot(1,type="n", - xlim = c(lld - 0.5, lhd + 2), - ylim = range(expdata[responsename], na.rm = TRUE), - xlab = paste("decadic logarithm of the concentration in ",levels(expdata$unit)), - ylab = ylab) - - drdatalist <- split(drdata,drdata$substance) - - for (i in 1:length(drdatalist)) { - points(log10(drdatalist[[i]]$conc),drdatalist[[i]][[responsename]],col=i); + endpoint <- expdata$type[1] + expdata <- subset(expdata, type == endpoint) + + controls <- subset(expdata, conc == 0) + expdata <- subset(expdata, conc != 0) + + QA <- matrix(nrow = 1, ncol = 4, + dimnames = list(c("Control (conc = 0)"), + c("Number", "Mean", "Std. Dev.", "% Std. Dev"))) + } + + numberOfControls <- length(controls$response) + QA["Control (conc = 0)", 1] <- numberOfControls + if (numberOfControls > 0) { + QA["Control (conc = 0)", 2] <- signif(mean(controls$response),2) + QA["Control (conc = 0)", 3] <- signif(sd(controls$response),2) + QA["Control (conc = 0)", 4] <- signif(QA["Control (conc = 0)", 3] * 100 / + QA["Control (conc = 0)", 2],2) + } + + if (db == "ecotox") { + if (identical(as.character(levels(expdata$organism)), "Vibrio fischeri")) { + positive <- subset(expdata, substance == "Na Cl") + if (nrow(positive) > 0) { + QA <- rbind(QA, + c(nrow(positive), + signif(mean(positive$raw_response), 2), + signif(sd(positive$raw_response), 2), + signif(100 * sd(positive$raw_response) / + mean(positive$raw_response), 2))) + + rownames(QA) <- c("Control (conc = 0)", + "Positive control (Na Cl)") + } + expdata <- subset(expdata, substance != "Na Cl", drop = TRUE) } - - legend("topright",substances, pch=1, col=1:length(substances), inset=0.05) - title(main=paste(levels(expdata$experimentator), - " - ",levels(expdata$type))) + } + + if (length(expdata$experimentator) < 1) { + stop("There is no response data for ", exptype, " ", + id, " in database ", db, "\n") + } + exptypestring <- paste0(toupper(substring(exptype, 1, 1)), + substring(exptype, 2)) + expdata$experimentator <- factor(expdata$experimentator) + expdata$type <- factor(expdata[[testtype]]) + expdata$performed <- factor(as.character(expdata$performed)) + expdata$substance <- factor(expdata$substance) + expdata$unit <- factor(expdata$unit) + expdata$ok <- factor(expdata$ok) + + # Info on the experiment + cat("\n", + exptypestring, id, "from database", db, ":\n\n", + "\tExperimentator(s):\t",levels(expdata$experimentator),"\n", + "\tType(s):\t\t",levels(expdata$type),"\n", + "\tPerformed on:\t\t",levels(expdata$performed),"\n", + "\tSubstance(s):\t\t",levels(expdata$substance),"\n", + "\tConcentration unit(s):\t",levels(expdata$unit),"\n", + "\tComment:\t\t",comment,"\n", + "\tOK Levels:\t\t",levels(expdata$ok),"\n\n") + + print(QA) + + # Control growth rate for Lemna and algae + if (endpoint %in% c("cell count", "frond area", "frond number")) { + duration <- as.numeric(unique(expdata$duration)) # in hours + if (length(duration) > 1) stop("More than one duration in the data") + response_0 <- unique(expdata$raw_0) + if (length(response_0) > 1) stop("More than one mean response at time 0 in the data") + t_days <- duration / 24 + control_growth_rates <- (log(controls$response) - log(response_0)) / t_days + cat("\nMean growth rate in controls:\t", round(mean(control_growth_rates), 3), "per day\n") + } + + # Box plot of control data + if (db == "ecotox") { + boxplot(controls$response, + names="controls", + ylab=endpoint, + ylim=range(controls$response, na.rm = TRUE), + boxwex=0.4, + main=paste("Plate ",id)) + } else { + boxplot(blinds$response,controls$response, + names=c("blinds","controls"), + ylab="Response", + boxwex=0.4, + main=paste("Plate ",id)) + } + + # Plot of dose response data + drdata <- expdata[c(2,4,6)] + drdata$substance <- factor(drdata$substance) + substances <- levels(drdata$substance) + + lld <- log10(min(subset(drdata,conc!=0)$conc)) + lhd <- log10(max(drdata$conc)) + + ylab <- if (db == "ecotox") endpoint + else responsename + + plot(1,type="n", + xlim = c(lld - 0.5, lhd + 2), + ylim = range(expdata[responsename], na.rm = TRUE), + xlab = paste("decadic logarithm of the concentration in ",levels(expdata$unit)), + ylab = ylab) + + drdatalist <- split(drdata,drdata$substance) + + for (i in 1:length(drdatalist)) { + points(log10(drdatalist[[i]]$conc),drdatalist[[i]][[responsename]],col=i); + } + + legend("topright",substances, pch=1, col=1:length(substances), inset=0.05) + title(main=paste(levels(expdata$experimentator), + " - ",levels(expdata$type))) } |