aboutsummaryrefslogtreecommitdiff
path: root/R/drfit.R
diff options
context:
space:
mode:
Diffstat (limited to 'R/drfit.R')
-rw-r--r--R/drfit.R1051
1 files changed, 530 insertions, 521 deletions
diff --git a/R/drfit.R b/R/drfit.R
index 7ffb1a0..bc0bddd 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -1,521 +1,530 @@
-drdata <- function(substances, experimentator = "%", db = "cytotox",
- celltype="IPC-81",enzymetype="AChE",whereClause="1",
- ok="'ok'")
-{
- library(RODBC)
- channel <- odbcConnect(db,uid="cytotox",pwd="cytotox",case="tolower")
- slist <- paste(substances,collapse="','")
- if (db == "cytotox") {
- responsetype <- "viability"
- testtype <- "celltype"
- type <- celltype
- } else {
- responsetype <- "activity"
- testtype <- "enzyme"
- type <- enzymetype
- }
-
- query <- paste("SELECT conc,",responsetype,",unit,experimentator,substance,",testtype,
- ",plate,ok FROM ", db, " WHERE substance IN ('",
- slist,"') AND experimentator LIKE '",
- experimentator,"' AND ",testtype," LIKE '",
- type,"' AND ",
- whereClause," AND ok in (",
- ok,")",sep="")
- data <- sqlQuery(channel,query)
- odbcClose(channel)
- names(data)[[1]] <- "dose"
- names(data)[[2]] <- "response"
- data$dosefactor <- factor(data$dose)
- data$substance <- factor(data$substance,levels=substances)
- return(data)
-}
-
-linearlogisf <- function(x,k,f,mu,b)
-{
- k*(1 + f*x) / (1 + ((2*f*(10^mu) + 1) * ((x/(10^mu))^b)))
-}
-
-drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
- lognorm = TRUE, logis = FALSE,
- linearlogis = FALSE, linearlogisWrong = NA,
- b0 = 2, f0 = 0)
-{
- substances <- levels(data$substance)
- unit <- levels(as.factor(data$unit))
-
- ri <- rix <- 0 # ri is the index over the result rows
- # rix is used later to check if any
- # model result was appended
- rsubstance <- array() # the substance names in the results
- rn <- vector() # number of dose-response curves
- rlhd <- rlld <- vector() # highest and lowest doses tested
- mtype <- array() # the modeltypes
- sigma <- array() # the standard deviation of the residuals
- logEC50 <- vector()
- stderrlogEC50 <- vector()
- slope <- vector()
- b <- vector()
- f <- vector()
-
- splitted <- split(data,data$substance)
- for (i in substances) {
- tmp <- splitted[[i]]
- fit <- FALSE
- n <- round(length(tmp$response)/9)
- if (length(tmp$response) == 0) {
- nodata = TRUE
- } else {
- nodata = FALSE
- }
- rix <- ri
- if (!nodata) {
- if (is.na(startlogEC50[i])){
- w <- 1/abs(tmp$response - 0.3)
- startlogEC50[[i]] <- sum(w * log10(tmp$dose))/sum(w)
- }
- highestdose <- max(tmp$dose)
- lowestdose <- min(tmp$dose)
- lhd <- log10(highestdose)
- lld <- log10(lowestdose)
- responseathighestdose <- mean(subset(tmp,dose==highestdose)$response)
- if (responseathighestdose < 0.5) {
- inactive <- FALSE
-
- if (linearlogis &&
- length(subset(linearlogisWrong,linearlogisWrong == i))==0) {
- m <- try(nls(response ~ linearlogisf(dose,1,f,logEC50,b),
- data=tmp,
- start=list(f=f0,logEC50=startlogEC50[[i]],b=b0)))
- if (!inherits(m, "try-error")) {
- fit <- TRUE
- ri <- ri + 1
- s <- summary(m)
- sigma[[ri]] <- s$sigma
- rsubstance[[ri]] <- i
- rn[[ri]] <- n
- rlld[[ri]] <- log10(lowestdose)
- rlhd[[ri]] <- log10(highestdose)
- mtype[[ri]] <- "linearlogis"
- logEC50[[ri]] <- coef(m)[["logEC50"]]
- slope[[ri]] <- NA
- if (logEC50[[ri]] > rlhd[[ri]]) {
- logEC50[[ri]] <- NA
- stderrlogEC50[[ri]] <- NA
- b[[ri]] <- NA
- f[[ri]] <- NA
- } else {
- stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
- b[[ri]] <- coef(m)[["b"]]
- f[[ri]] <- coef(m)[["f"]]
- }
- }
- }
-
- if (logis) {
- m <- try(nls(response ~ plogis(-log10(dose),-logEC50,slope),
- data=tmp,
- start=list(logEC50=startlogEC50[[i]],slope=1)))
- if (chooseone==FALSE || fit==FALSE) {
- if (!inherits(m, "try-error")) {
- fit <- TRUE
- ri <- ri + 1
- s <- summary(m)
- sigma[[ri]] <- s$sigma
- rsubstance[[ri]] <- i
- rn[[ri]] <- n
- rlld[[ri]] <- log10(lowestdose)
- rlhd[[ri]] <- log10(highestdose)
- mtype[[ri]] <- "logis"
- logEC50[[ri]] <- coef(m)[["logEC50"]]
- b[[ri]] <- NA
- f[[ri]] <- NA
- if (logEC50[[ri]] > rlhd[[ri]]) {
- logEC50[[ri]] <- NA
- slope[[ri]] <- NA
- stderrlogEC50[[ri]] <- NA
- } else {
- slope[[ri]] <- coef(m)[["slope"]]
- stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
- }
- }
- }
- }
-
- if (lognorm) {
- m <- try(nls(response ~ pnorm(-log10(dose),-logEC50,slope),
- data=tmp,
- start=list(logEC50=startlogEC50[[i]],slope=1)))
- if (chooseone==FALSE || fit==FALSE) {
- if (!inherits(m, "try-error")) {
- fit <- TRUE
- ri <- ri + 1
- s <- summary(m)
- sigma[[ri]] <- s$sigma
- rsubstance[[ri]] <- i
- rn[[ri]] <- n
- rlld[[ri]] <- log10(lowestdose)
- rlhd[[ri]] <- log10(highestdose)
- mtype[[ri]] <- "lognorm"
- logEC50[[ri]] <- coef(m)[["logEC50"]]
- b[[ri]] <- NA
- f[[ri]] <- NA
- if (logEC50[[ri]] > rlhd[[ri]]) {
- logEC50[[ri]] <- NA
- slope[[ri]] <- NA
- stderrlogEC50[[ri]] <- NA
- } else {
- slope[[ri]] <- coef(m)[["slope"]]
- stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
- }
- }
- }
- }
-
- } else {
- inactive <- TRUE
- }
- }
- if (ri == rix) { # if no entry was appended for this substance
- ri <- ri + 1
- rsubstance[[ri]] <- i
- rn[[ri]] <- n
- if (nodata) {
- rlld[[ri]] <- rlhd[[i]] <- NA
- mtype[[ri]] <- "no data"
- } else {
- rlld[[ri]] <- log10(lowestdose)
- rlhd[[i]] <- log10(highestdose)
- if (inactive) {
- mtype[[ri]] <- "inactive"
- } else {
- mtype[[ri]] <- "no fit"
- }
- }
- sigma[[ri]] <- NA
- logEC50[[ri]] <- NA
- stderrlogEC50[[ri]] <- NA
- slope[[ri]] <- NA
- b[[ri]] <- NA
- f[[ri]] <- NA
- }
- }
- results <- data.frame(rsubstance,rn, rlhd, mtype, logEC50, stderrlogEC50, unit, sigma)
- names(results) <- c("Substance","n", "lhd","mtype","logEC50","std","unit","sigma")
- if (lognorm || logis) {
- results$slope <- slope
- }
- if (linearlogis) {
- results$b <- b
- results$f <- f
- }
- return(results)
-}
-
-drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
- path = "./", fileprefix = "drplot", overlay = FALSE,
- postscript = FALSE, png = FALSE, bw = TRUE,
- colors = 1:8,devoff=T)
-{
- unitlevels <- levels(as.factor(drresults$unit))
- if (length(unitlevels) == 1) {
- unit <- unitlevels
- } else {
- unit <- "different units"
- }
-
- # Get the plot limits on the x-axis (log of the dose)
- if(is.data.frame(data)) {
- if (min(data$dose == 0)) {
- cat("At least one of the dose levels is 0 - this is not a valid dose.")
- } else {
- lld <- log10(min(data$dose))
- }
- lhd <- log10(max(data$dose))
- hr <- max(data$response)
- dsubstances <- levels(data$substance)
- } else {
- lld <- min(drresults[["logEC50"]],na.rm=TRUE) - 2
- lhd <- max(drresults[["logEC50"]],na.rm=TRUE) + 2
- if (length(subset(drresults,mtype=="linearlogis")$Substance) != 0) {
- hr <- 1.8
- } else {
- hr <- 1.0
- }
- }
-
- # Prepare overlay plot if requested
- if (overlay)
- {
- if (postscript) {
- filename = paste(path,fileprefix,".eps",sep="")
- postscript(file=filename,
- paper="special",width=7,height=7,horizontal=FALSE,pointsize=12)
- cat("Created File: ",filename,"\n")
- }
- if (png) {
- filename = paste(path,fileprefix,".png",sep="")
- png(filename=filename,
- width=500, height=500,pointsize=12)
- cat("Created File: ",filename,"\n")
- }
- if (!postscript && !png) {
- x11(width=7,height=7,pointsize=12)
- }
-
- plot(0,type="n",
- xlim=c(lld - 0.5, lhd + 1),
- ylim= c(-0.1, hr + 0.2),
- xlab=paste("Decadic Logarithm of the dose in ", unit),
- ylab="Normalized response")
- }
-
- # Plot the data either as raw data or as error bars
- if(is.data.frame(data)) {
- splitted <- split(data,data$substance)
- n <- 0
- if (bw) colors <- rep("black",length(dsubstances))
- for (i in dsubstances) {
- n <- n + 1
- tmp <- splitted[[i]]
- if (length(tmp$response) != 0) {
- color <- colors[[n]]
- # Prepare the single graphs if an overlay is not requested
- if (!overlay)
- {
- if (postscript) {
- filename = paste(path,fileprefix,sub(" ","_",gsub("([\(\) ])", "", i)),".eps",sep="")
- postscript(file=filename,
- paper="special",width=7,height=7,horizontal=FALSE,pointsize=12)
- cat("Created File: ",filename,"\n")
- }
- if (png) {
- filename = paste(path,fileprefix,sub(" ","_",gsub("([\(\) ])", "", i)),".png",sep="")
- png(filename=filename,
- width=500, height=500,pointsize=12)
- cat("Created File: ",filename,"\n")
- }
- if (!postscript && !png) {
- x11(width=7,height=7,pointsize=12)
- }
-
- plot(0,type="n",
- xlim=c(lld - 0.5, lhd + 2),
- ylim= c(-0.1, hr + 0.2),
- xlab=paste("Decadic Logarithm of the dose in ", unit),
- ylab="Normalized response")
- }
- if (!overlay) legend(lhd - 1, hr + 0.1, i,lty = 1, col = color)
- tmp$dosefactor <- factor(tmp$dose) # necessary because the old
- # factor has all levels, not
- # only the ones tested with
- # this substance
- # Plot the data, if requested
- if (dtype != "none") {
- if (dtype == "raw") {
- points(log10(tmp$dose),tmp$response,col=color)
- } else {
- splitresponses <- split(tmp$response,tmp$dosefactor)
- means <- sapply(splitresponses,mean)
- lengths <- sapply(splitresponses,length)
- vars <- sapply(splitresponses,var)
- standarddeviations <- sqrt(vars)
- }
- if (dtype == "std")
- {
- tops <- means + standarddeviations
- bottoms <- means - standarddeviations
- }
- if (dtype == "conf")
- {
- confidencedeltas <- qt((1 + alpha)/2, lengths - 1) * sqrt(vars)
- tops <- means + confidencedeltas
- bottoms <- means - confidencedeltas
- }
- if (dtype != "raw")
- {
- x <- log10(as.numeric(levels(tmp$dosefactor)))
- segments(x,bottoms,x,tops,col=color)
- points(x,means,col=color)
- smidge <- 0.05
- segments(x - smidge,bottoms,x + smidge,bottoms,col=color)
- segments(x - smidge,tops,x + smidge,tops,col=color)
- }
- }
-
- # Plot the fits, if there are any
- fits <- subset(drresults,Substance == i)
- nf <- length(fits$Substance) # number of fits to plot
- if (nf > 0) {
- for (j in 1:nf)
- {
- logEC50 <- fits[j,"logEC50"]
- mtype <- as.character(fits[j, "mtype"])
- if (mtype == "lognorm") {
- slope <- fits[j,"slope"]
- plot(function(x) pnorm(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=color)
- }
- if (mtype == "logis") {
- slope <- fits[j,"slope"]
- plot(function(x) plogis(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=color)
- }
- if (mtype == "linearlogis") {
- plot(function(x) linearlogisf(10^x,1,fits[j,"f"],fits[j,"logEC50"],fits[j,"b"]),
- lld - 0.5, lhd + 2,
- add=TRUE,col=color)
- }
- }
- }
- if (!overlay && (postscript || png)) dev.off()
- } else {
- cat("No data for ",i,"\n")
- }
- }
- }
- if (overlay) legend(lhd - 1, hr + 0.1, dsubstances,lty = 1, col = colors)
- if (overlay && (postscript || png)) {
- if (devoff) {
- dev.off()
- }
- }
-}
-
-checkplate <- function(plate,db="cytotox")
-{
- library(RODBC)
- channel <- odbcConnect(db,uid="cytotox",pwd="cytotox",case="tolower")
-
- if (db == "cytotox") {
- responsetype <- "viability"
- testtype <- "celltype"
- } else {
- responsetype <- "activity"
- testtype <- "enzyme"
- }
-
- platequery <- paste("SELECT experimentator,substance,",testtype,",conc,unit,",responsetype,",performed,ok",
- "FROM ",db," WHERE plate=", plate)
-
- controlquery <- paste("SELECT type,response FROM controls WHERE plate=",plate)
-
- platedata <- sqlQuery(channel,platequery)
- controldata <- sqlQuery(channel,controlquery)
-
- odbcClose(channel)
-
- if (length(platedata$experimentator) < 1) {
- cat("There is no response data for plate ",plate," in database ",db,"\n")
- } else {
- platedata$experimentator <- factor(platedata$experimentator)
- platedata$type <- factor(platedata[[testtype]])
- platedata$substance <- factor(platedata$substance)
- platedata$unit <- factor(platedata$unit)
- platedata$performed <- factor(platedata$performed)
- platedata$ok <- factor(platedata$ok)
-
- blinds <- subset(controldata,type=="blind")
- controls <- subset(controldata,type=="control")
-
- numberOfBlinds <- length(blinds$response)
- numberOfControls <- length(controls$response)
- meanOfBlinds <- mean(blinds$response)
- meanOfControls <- mean(controls$response)
- stdOfBlinds <- sd(blinds$response)
- stdOfControls <- sd(controls$response)
-
- cat("Plate ",plate," from database ",db,"\n",
- "\tExperimentator: ",levels(platedata$experimentator),"\n",
- "\tType(s): ",levels(platedata$type),"\n",
- "\tPerformed on : ",levels(platedata$performed),"\n",
- "\tSubstance(s): ",levels(platedata$substance),"\n",
- "\tConcentration unit: ",levels(platedata$unit),"\n",
- "\tOK: ",levels(platedata$ok),"\n",
- "\t\tNumber \tMean \tStandard Deviation\n",
- "blind\t\t",numberOfBlinds,"\t",meanOfBlinds,"\t",stdOfBlinds,"\n",
- "control\t",numberOfControls,"\t",meanOfControls,"\t",stdOfControls,"\n")
-
- par(ask=TRUE)
-
- boxplot(blinds$response,controls$response,names=c("blinds","controls"),ylab="Response",main=paste("Plate ",plate))
-
- drdata <- platedata[c(2,4,6)]
- drdata$substance <- factor(drdata$substance)
- substances <- levels(drdata$substance)
-
- plot(log10(drdata$conc),drdata$viability,
- xlim=c(-2.5, 4.5),
- ylim= c(-0.1, 2),
- xlab=paste("decadic logarithm of the concentration in ",levels(platedata$unit)),
- ylab=responsetype)
-
- drdatalist <- split(drdata,drdata$substance)
-
- for (i in 1:length(drdatalist)) {
- points(log10(drdatalist[[i]]$conc),drdatalist[[i]][[responsetype]],col=i);
- }
-
- legend(3.0,1.5,substances, pch=1, col=1:length(substances))
- title(main=paste("Plate ",plate," - ",levels(platedata$experimentator)," - ",levels(platedata$type)))
- }
-}
-
-checksubstance <- function(substance,db="cytotox",experimentator="%",celltype="%",enzymetype="%",whereClause="1",ok="%")
-{
- library(RODBC)
- channel <- odbcConnect(db,uid="cytotox",pwd="cytotox",case="tolower")
-
- if (db == "cytotox") {
- responsetype <- "viability"
- testtype <- "celltype"
- type <- celltype
- } else {
- responsetype <- "activity"
- testtype <- "enzyme"
- type <- enzymetype
- }
- query <- paste("SELECT experimentator,substance,",testtype,",plate,conc,unit,",responsetype,",ok",
- " FROM ",db," WHERE substance LIKE '",
- substance,"' AND experimentator LIKE '",
- experimentator,"' AND ",testtype," LIKE '",
- type,"' AND ",
- whereClause," AND ok LIKE '",ok,"'",sep="")
-
- data <- sqlQuery(channel,query)
- odbcClose(channel)
-
- data$experimentator <- factor(data$experimentator)
- data$substance <- factor(data$substance)
- substances <- levels(data$substance)
- data$type <- factor(data[[testtype]])
- data$plate <- factor(data$plate)
- plates <- levels(data$plate)
- concentrations <- split(data$conc,data$conc)
- concentrations <- as.numeric(names(concentrations))
- data$unit <- factor(data$unit)
- data$ok <- factor(data$ok)
-
- if (length(plates)>6) {
- palette(rainbow(length(plates)))
- }
-
- plot(log10(data$conc),data[[responsetype]],
- xlim=c(-2.5, 4.5),
- ylim= c(-0.1, 2),
- xlab=paste("decadic logarithm of the concentration in ",levels(data$unit)),
- ylab=responsetype)
-
- platelist <- split(data,data$plate)
-
- for (i in 1:length(platelist)) {
- points(log10(platelist[[i]]$conc),platelist[[i]][[responsetype]],col=i);
- }
-
- legend(3.5,1.7,plates, pch=1, col=1:length(plates))
- title(main=paste(substance," - ",levels(data$experimentator)," - ",levels(data$type)))
-
- cat("Substanz ",substance,"\n",
- "\tExperimentator(s):",levels(data$experimentator),"\n",
- "\tType(s):\t",levels(data$type),"\n",
- "\tSubstance(s):\t",levels(data$substance),"\n",
- "\tPlate(s):\t",plates,"\n\n")
-}
+drdata <- function(substances, experimentator = "%", db = "cytotox",
+ celltype="IPC-81",enzymetype="AChE",whereClause="1",
+ ok="'ok'")
+{
+ library(RODBC)
+ channel <- odbcConnect(db,uid="cytotox",pwd="cytotox",case="tolower")
+ slist <- paste(substances,collapse="','")
+ if (db == "cytotox") {
+ responsetype <- "viability"
+ testtype <- "celltype"
+ type <- celltype
+ } else {
+ responsetype <- "activity"
+ testtype <- "enzyme"
+ type <- enzymetype
+ }
+
+ query <- paste("SELECT conc,",responsetype,",unit,experimentator,substance,",testtype,
+ ",plate,ok FROM ", db, " WHERE substance IN ('",
+ slist,"') AND experimentator LIKE '",
+ experimentator,"' AND ",testtype," LIKE '",
+ type,"' AND ",
+ whereClause," AND ok in (",
+ ok,")",sep="")
+ data <- sqlQuery(channel,query)
+ odbcClose(channel)
+ names(data)[[1]] <- "dose"
+ names(data)[[2]] <- "response"
+ data$dosefactor <- factor(data$dose)
+ data$substance <- factor(data$substance,levels=substances)
+ return(data)
+}
+
+linearlogisf <- function(x,k,f,mu,b)
+{
+ k*(1 + f*x) / (1 + ((2*f*(10^mu) + 1) * ((x/(10^mu))^b)))
+}
+
+drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
+ lognorm = TRUE, logis = FALSE,
+ linearlogis = FALSE, linearlogisWrong = NA,
+ b0 = 2, f0 = 0)
+{
+ substances <- levels(data$substance)
+ unit <- levels(as.factor(data$unit))
+
+ ri <- rix <- 0 # ri is the index over the result rows
+ # rix is used later to check if any
+ # model result was appended
+ rsubstance <- array() # the substance names in the results
+ rn <- vector() # number of dose-response curves
+ rlhd <- rlld <- vector() # highest and lowest doses tested
+ mtype <- array() # the modeltypes
+ sigma <- array() # the standard deviation of the residuals
+ logEC50 <- vector()
+ stderrlogEC50 <- vector()
+ slope <- vector()
+ b <- vector()
+ f <- vector()
+
+ splitted <- split(data,data$substance)
+ for (i in substances) {
+ tmp <- splitted[[i]]
+ fit <- FALSE
+ n <- round(length(tmp$response)/9)
+ if (length(tmp$response) == 0) {
+ nodata = TRUE
+ } else {
+ nodata = FALSE
+ }
+ rix <- ri
+ if (!nodata) {
+ if (is.na(startlogEC50[i])){
+ w <- 1/abs(tmp$response - 0.3)
+ startlogEC50[[i]] <- sum(w * log10(tmp$dose))/sum(w)
+ }
+ highestdose <- max(tmp$dose)
+ lowestdose <- min(tmp$dose)
+ lhd <- log10(highestdose)
+ lld <- log10(lowestdose)
+ responseathighestdose <- mean(subset(tmp,dose==highestdose)$response)
+ if (responseathighestdose < 0.5) {
+ inactive <- FALSE
+
+ if (linearlogis &&
+ length(subset(linearlogisWrong,linearlogisWrong == i))==0) {
+ m <- try(nls(response ~ linearlogisf(dose,1,f,logEC50,b),
+ data=tmp,
+ start=list(f=f0,logEC50=startlogEC50[[i]],b=b0)))
+ if (!inherits(m, "try-error")) {
+ fit <- TRUE
+ ri <- ri + 1
+ s <- summary(m)
+ sigma[[ri]] <- s$sigma
+ rsubstance[[ri]] <- i
+ rn[[ri]] <- n
+ rlld[[ri]] <- log10(lowestdose)
+ rlhd[[ri]] <- log10(highestdose)
+ mtype[[ri]] <- "linearlogis"
+ logEC50[[ri]] <- coef(m)[["logEC50"]]
+ slope[[ri]] <- NA
+ if (logEC50[[ri]] > rlhd[[ri]]) {
+ logEC50[[ri]] <- NA
+ stderrlogEC50[[ri]] <- NA
+ b[[ri]] <- NA
+ f[[ri]] <- NA
+ } else {
+ stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
+ b[[ri]] <- coef(m)[["b"]]
+ f[[ri]] <- coef(m)[["f"]]
+ }
+ }
+ }
+
+ if (logis) {
+ m <- try(nls(response ~ plogis(-log10(dose),-logEC50,slope),
+ data=tmp,
+ start=list(logEC50=startlogEC50[[i]],slope=1)))
+ if (chooseone==FALSE || fit==FALSE) {
+ if (!inherits(m, "try-error")) {
+ fit <- TRUE
+ ri <- ri + 1
+ s <- summary(m)
+ sigma[[ri]] <- s$sigma
+ rsubstance[[ri]] <- i
+ rn[[ri]] <- n
+ rlld[[ri]] <- log10(lowestdose)
+ rlhd[[ri]] <- log10(highestdose)
+ mtype[[ri]] <- "logis"
+ logEC50[[ri]] <- coef(m)[["logEC50"]]
+ b[[ri]] <- NA
+ f[[ri]] <- NA
+ if (logEC50[[ri]] > rlhd[[ri]]) {
+ logEC50[[ri]] <- NA
+ slope[[ri]] <- NA
+ stderrlogEC50[[ri]] <- NA
+ } else {
+ slope[[ri]] <- coef(m)[["slope"]]
+ stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
+ }
+ }
+ }
+ }
+
+ if (lognorm) {
+ m <- try(nls(response ~ pnorm(-log10(dose),-logEC50,slope),
+ data=tmp,
+ start=list(logEC50=startlogEC50[[i]],slope=1)))
+ if (chooseone==FALSE || fit==FALSE) {
+ if (!inherits(m, "try-error")) {
+ fit <- TRUE
+ ri <- ri + 1
+ s <- summary(m)
+ sigma[[ri]] <- s$sigma
+ rsubstance[[ri]] <- i
+ rn[[ri]] <- n
+ rlld[[ri]] <- log10(lowestdose)
+ rlhd[[ri]] <- log10(highestdose)
+ mtype[[ri]] <- "lognorm"
+ logEC50[[ri]] <- coef(m)[["logEC50"]]
+ b[[ri]] <- NA
+ f[[ri]] <- NA
+ if (logEC50[[ri]] > rlhd[[ri]]) {
+ logEC50[[ri]] <- NA
+ slope[[ri]] <- NA
+ stderrlogEC50[[ri]] <- NA
+ } else {
+ slope[[ri]] <- coef(m)[["slope"]]
+ stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
+ }
+ }
+ }
+ }
+
+ } else {
+ inactive <- TRUE
+ }
+ }
+ if (ri == rix) { # if no entry was appended for this substance
+ ri <- ri + 1
+ rsubstance[[ri]] <- i
+ rn[[ri]] <- n
+ if (nodata) {
+ rlld[[ri]] <- rlhd[[i]] <- NA
+ mtype[[ri]] <- "no data"
+ } else {
+ rlld[[ri]] <- log10(lowestdose)
+ rlhd[[i]] <- log10(highestdose)
+ if (inactive) {
+ mtype[[ri]] <- "inactive"
+ } else {
+ mtype[[ri]] <- "no fit"
+ }
+ }
+ sigma[[ri]] <- NA
+ logEC50[[ri]] <- NA
+ stderrlogEC50[[ri]] <- NA
+ slope[[ri]] <- NA
+ b[[ri]] <- NA
+ f[[ri]] <- NA
+ }
+ }
+ results <- data.frame(rsubstance,rn, rlhd, mtype, logEC50, stderrlogEC50, unit, sigma)
+ names(results) <- c("Substance","n", "lhd","mtype","logEC50","std","unit","sigma")
+ if (lognorm || logis) {
+ results$slope <- slope
+ }
+ if (linearlogis) {
+ results$b <- b
+ results$f <- f
+ }
+ return(results)
+}
+
+drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
+ path = "./", fileprefix = "drplot", overlay = FALSE,
+ postscript = FALSE, png = FALSE, bw = TRUE,
+ colors = 1:8,devoff=T,lpos=FALSE)
+{
+ unitlevels <- levels(as.factor(drresults$unit))
+ if (length(unitlevels) == 1) {
+ unit <- unitlevels
+ } else {
+ unit <- "different units"
+ }
+
+ # Get the plot limits on the x-axis (log of the dose)
+ if(is.data.frame(data)) {
+ if (min(data$dose == 0)) {
+ cat("At least one of the dose levels is 0 - this is not a valid dose.")
+ } else {
+ lld <- log10(min(data$dose))
+ }
+ lhd <- log10(max(data$dose))
+ hr <- max(data$response)
+ dsubstances <- levels(data$substance)
+ } else {
+ lld <- min(drresults[["logEC50"]],na.rm=TRUE) - 2
+ lhd <- max(drresults[["logEC50"]],na.rm=TRUE) + 2
+ if (length(subset(drresults,mtype=="linearlogis")$Substance) != 0) {
+ hr <- 1.8
+ } else {
+ hr <- 1.0
+ }
+ }
+
+ # Legend position
+ if (!lpos[[1]]) {
+ lx <- lhd - 1
+ ly <- hr + 0.1
+ } else {
+ lx <- lpos[[1]]
+ ly <- lpos[[2]]
+ }
+
+ # Prepare overlay plot if requested
+ if (overlay)
+ {
+ if (postscript) {
+ filename = paste(path,fileprefix,".eps",sep="")
+ postscript(file=filename,
+ paper="special",width=7,height=7,horizontal=FALSE,pointsize=12)
+ cat("Created File: ",filename,"\n")
+ }
+ if (png) {
+ filename = paste(path,fileprefix,".png",sep="")
+ png(filename=filename,
+ width=500, height=500,pointsize=12)
+ cat("Created File: ",filename,"\n")
+ }
+ if (!postscript && !png) {
+ x11(width=7,height=7,pointsize=12)
+ }
+
+ plot(0,type="n",
+ xlim=c(lld - 0.5, lhd + 1),
+ ylim= c(-0.1, hr + 0.2),
+ xlab=paste("Decadic Logarithm of the dose in ", unit),
+ ylab="Normalized response")
+ }
+
+ # Plot the data either as raw data or as error bars
+ if(is.data.frame(data)) {
+ splitted <- split(data,data$substance)
+ n <- 0
+ if (bw) colors <- rep("black",length(dsubstances))
+ for (i in dsubstances) {
+ n <- n + 1
+ tmp <- splitted[[i]]
+ if (length(tmp$response) != 0) {
+ color <- colors[[n]]
+ # Prepare the single graphs if an overlay is not requested
+ if (!overlay)
+ {
+ if (postscript) {
+ filename = paste(path,fileprefix,sub(" ","_",gsub("([\(\) ])", "", i)),".eps",sep="")
+ postscript(file=filename,
+ paper="special",width=7,height=7,horizontal=FALSE,pointsize=12)
+ cat("Created File: ",filename,"\n")
+ }
+ if (png) {
+ filename = paste(path,fileprefix,sub(" ","_",gsub("([\(\) ])", "", i)),".png",sep="")
+ png(filename=filename,
+ width=500, height=500,pointsize=12)
+ cat("Created File: ",filename,"\n")
+ }
+ if (!postscript && !png) {
+ x11(width=7,height=7,pointsize=12)
+ }
+
+ plot(0,type="n",
+ xlim=c(lld - 0.5, lhd + 2),
+ ylim= c(-0.1, hr + 0.2),
+ xlab=paste("Decadic Logarithm of the dose in ", unit),
+ ylab="Normalized response")
+ }
+ if (!overlay) legend(lx, ly, i,lty = 1, col = color)
+ tmp$dosefactor <- factor(tmp$dose) # necessary because the old
+ # factor has all levels, not
+ # only the ones tested with
+ # this substance
+ # Plot the data, if requested
+ if (dtype != "none") {
+ if (dtype == "raw") {
+ points(log10(tmp$dose),tmp$response,col=color)
+ } else {
+ splitresponses <- split(tmp$response,tmp$dosefactor)
+ means <- sapply(splitresponses,mean)
+ lengths <- sapply(splitresponses,length)
+ vars <- sapply(splitresponses,var)
+ standarddeviations <- sqrt(vars)
+ }
+ if (dtype == "std")
+ {
+ tops <- means + standarddeviations
+ bottoms <- means - standarddeviations
+ }
+ if (dtype == "conf")
+ {
+ confidencedeltas <- qt((1 + alpha)/2, lengths - 1) * sqrt(vars)
+ tops <- means + confidencedeltas
+ bottoms <- means - confidencedeltas
+ }
+ if (dtype != "raw")
+ {
+ x <- log10(as.numeric(levels(tmp$dosefactor)))
+ segments(x,bottoms,x,tops,col=color)
+ points(x,means,col=color)
+ smidge <- 0.05
+ segments(x - smidge,bottoms,x + smidge,bottoms,col=color)
+ segments(x - smidge,tops,x + smidge,tops,col=color)
+ }
+ }
+
+ # Plot the fits, if there are any
+ fits <- subset(drresults,Substance == i)
+ nf <- length(fits$Substance) # number of fits to plot
+ if (nf > 0) {
+ for (j in 1:nf)
+ {
+ logEC50 <- fits[j,"logEC50"]
+ mtype <- as.character(fits[j, "mtype"])
+ if (mtype == "lognorm") {
+ slope <- fits[j,"slope"]
+ plot(function(x) pnorm(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=color)
+ }
+ if (mtype == "logis") {
+ slope <- fits[j,"slope"]
+ plot(function(x) plogis(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=color)
+ }
+ if (mtype == "linearlogis") {
+ plot(function(x) linearlogisf(10^x,1,fits[j,"f"],fits[j,"logEC50"],fits[j,"b"]),
+ lld - 0.5, lhd + 2,
+ add=TRUE,col=color)
+ }
+ }
+ }
+ if (!overlay && (postscript || png)) dev.off()
+ } else {
+ cat("No data for ",i,"\n")
+ }
+ }
+ }
+ if (overlay) legend(lx, ly, dsubstances,lty = 1, col = colors)
+ if (overlay && (postscript || png)) {
+ if (devoff) {
+ dev.off()
+ }
+ }
+}
+
+checkplate <- function(plate,db="cytotox")
+{
+ library(RODBC)
+ channel <- odbcConnect(db,uid="cytotox",pwd="cytotox",case="tolower")
+
+ if (db == "cytotox") {
+ responsetype <- "viability"
+ testtype <- "celltype"
+ } else {
+ responsetype <- "activity"
+ testtype <- "enzyme"
+ }
+
+ platequery <- paste("SELECT experimentator,substance,",testtype,",conc,unit,",responsetype,",performed,ok",
+ "FROM ",db," WHERE plate=", plate)
+
+ controlquery <- paste("SELECT type,response FROM controls WHERE plate=",plate)
+
+ platedata <- sqlQuery(channel,platequery)
+ controldata <- sqlQuery(channel,controlquery)
+
+ odbcClose(channel)
+
+ if (length(platedata$experimentator) < 1) {
+ cat("There is no response data for plate ",plate," in database ",db,"\n")
+ } else {
+ platedata$experimentator <- factor(platedata$experimentator)
+ platedata$type <- factor(platedata[[testtype]])
+ platedata$substance <- factor(platedata$substance)
+ platedata$unit <- factor(platedata$unit)
+ platedata$performed <- factor(platedata$performed)
+ platedata$ok <- factor(platedata$ok)
+
+ blinds <- subset(controldata,type=="blind")
+ controls <- subset(controldata,type=="control")
+
+ numberOfBlinds <- length(blinds$response)
+ numberOfControls <- length(controls$response)
+ meanOfBlinds <- mean(blinds$response)
+ meanOfControls <- mean(controls$response)
+ stdOfBlinds <- sd(blinds$response)
+ stdOfControls <- sd(controls$response)
+
+ cat("Plate ",plate," from database ",db,"\n",
+ "\tExperimentator: ",levels(platedata$experimentator),"\n",
+ "\tType(s): ",levels(platedata$type),"\n",
+ "\tPerformed on : ",levels(platedata$performed),"\n",
+ "\tSubstance(s): ",levels(platedata$substance),"\n",
+ "\tConcentration unit: ",levels(platedata$unit),"\n",
+ "\tOK: ",levels(platedata$ok),"\n",
+ "\t\tNumber \tMean \tStandard Deviation\n",
+ "blind\t\t",numberOfBlinds,"\t",meanOfBlinds,"\t",stdOfBlinds,"\n",
+ "control\t",numberOfControls,"\t",meanOfControls,"\t",stdOfControls,"\n")
+
+ par(ask=TRUE)
+
+ boxplot(blinds$response,controls$response,names=c("blinds","controls"),ylab="Response",main=paste("Plate ",plate))
+
+ drdata <- platedata[c(2,4,6)]
+ drdata$substance <- factor(drdata$substance)
+ substances <- levels(drdata$substance)
+
+ plot(log10(drdata$conc),drdata$viability,
+ xlim=c(-2.5, 4.5),
+ ylim= c(-0.1, 2),
+ xlab=paste("decadic logarithm of the concentration in ",levels(platedata$unit)),
+ ylab=responsetype)
+
+ drdatalist <- split(drdata,drdata$substance)
+
+ for (i in 1:length(drdatalist)) {
+ points(log10(drdatalist[[i]]$conc),drdatalist[[i]][[responsetype]],col=i);
+ }
+
+ legend(3.0,1.5,substances, pch=1, col=1:length(substances))
+ title(main=paste("Plate ",plate," - ",levels(platedata$experimentator)," - ",levels(platedata$type)))
+ }
+}
+
+checksubstance <- function(substance,db="cytotox",experimentator="%",celltype="%",enzymetype="%",whereClause="1",ok="%")
+{
+ library(RODBC)
+ channel <- odbcConnect(db,uid="cytotox",pwd="cytotox",case="tolower")
+
+ if (db == "cytotox") {
+ responsetype <- "viability"
+ testtype <- "celltype"
+ type <- celltype
+ } else {
+ responsetype <- "activity"
+ testtype <- "enzyme"
+ type <- enzymetype
+ }
+ query <- paste("SELECT experimentator,substance,",testtype,",plate,conc,unit,",responsetype,",ok",
+ " FROM ",db," WHERE substance LIKE '",
+ substance,"' AND experimentator LIKE '",
+ experimentator,"' AND ",testtype," LIKE '",
+ type,"' AND ",
+ whereClause," AND ok LIKE '",ok,"'",sep="")
+
+ data <- sqlQuery(channel,query)
+ odbcClose(channel)
+
+ data$experimentator <- factor(data$experimentator)
+ data$substance <- factor(data$substance)
+ substances <- levels(data$substance)
+ data$type <- factor(data[[testtype]])
+ data$plate <- factor(data$plate)
+ plates <- levels(data$plate)
+ concentrations <- split(data$conc,data$conc)
+ concentrations <- as.numeric(names(concentrations))
+ data$unit <- factor(data$unit)
+ data$ok <- factor(data$ok)
+
+ if (length(plates)>6) {
+ palette(rainbow(length(plates)))
+ }
+
+ plot(log10(data$conc),data[[responsetype]],
+ xlim=c(-2.5, 4.5),
+ ylim= c(-0.1, 2),
+ xlab=paste("decadic logarithm of the concentration in ",levels(data$unit)),
+ ylab=responsetype)
+
+ platelist <- split(data,data$plate)
+
+ for (i in 1:length(platelist)) {
+ points(log10(platelist[[i]]$conc),platelist[[i]][[responsetype]],col=i);
+ }
+
+ legend(3.5,1.7,plates, pch=1, col=1:length(plates))
+ title(main=paste(substance," - ",levels(data$experimentator)," - ",levels(data$type)))
+
+ cat("Substanz ",substance,"\n",
+ "\tExperimentator(s):",levels(data$experimentator),"\n",
+ "\tType(s):\t",levels(data$type),"\n",
+ "\tSubstance(s):\t",levels(data$substance),"\n",
+ "\tPlate(s):\t",plates,"\n\n")
+}

Contact - Imprint