params <- list(my_css = "css/rmdformats.css") ## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.align = 'center' ) ## ----setup,class.source="fold-show",warning = FALSE, message = FALSE---------- # Primary Packages # library(tcplfit2) library(tcpl) # Data Formatting Packages # library(data.table) library(DT) library(htmlTable) library(dplyr) library(stringr) # Plotting Packages # library(ggplot2) library(gridExtra) ## ----ex1_concRespCore,class.source="fold-show",warning=FALSE------------------ # tested concentrations conc <- list(.03,.1,.3,1,3,10,30,100) # observed responses at respective concentrations resp <- list(0,.2,.1,.4,.7,.9,.6, 1.2) # row object with relevant parameters row = list(conc = conc,resp = resp,bmed = 0,cutoff = 1,onesd = 0.5,name="some chemical") # execute concentration-response modeling through potency estimation res <- concRespCore(row, fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", "pow", "exp2", "exp3", "exp4", "exp5"), conthits = T) ## ----echo=FALSE--------------------------------------------------------------- htmlTable::htmlTable(head(res), align = 'l', align.header = 'l', rnames = FALSE , css.cell = ' padding-bottom: 5px; vertical-align:top; padding-right: 10px;min-width: 5em ') ## ----ex1_concRespPlot2, fig.height = 4.55, fig.width = 8---------------------- # plot the winning curve from example 1, add a title concRespPlot2(res, log_conc = TRUE) + ggtitle("Example 1: Chemical A") ## ----example2,class.source="fold-show",warning=FALSE-------------------------- # read in the data # Loading in the level 3 example data set from invitrodb stored in tcplfit2 data("mc3") # view the first 6 rows of the mc3 data # dtxsid = unique chemical identifier from EPA's DSSTox Database # casrn = unique chemical identifier from Chemical Abstracts Service # name = chemical name # spid = sample id # logc = log_10 concentration value # resp = response # assay = assay name head(mc3) # estimate the background variability # assume the two lowest concentrations (logc <= -2) for baseline in this example # Note: The baseline may be assay/application specific temp <- mc3[mc3$logc<= -2,"resp"] # obtain response in the two lowest concentrations bmad <- mad(temp) # obtain the baseline median absolute deviation onesd <- sd(temp) # obtain the baseline standard deviation cutoff <- 3*bmad # estimate the cutoff, use the typical cutoff=3*BMAD # select six chemical samples # Note: there may be more than one sample processed for a given chemical spid.list <- unique(mc3$spid) spid.list <- spid.list[1:6] # create empty objects to store fitting results and plots model_fits <- NULL result_table <- NULL plt_lst <- NULL # loop over the samples to perform concentration-response modeling & hitcalling for(spid in spid.list) { # select the data for just this sample temp <- mc3[is.element(mc3$spid,spid),] # The data file stores concentrations in log10 units, so back-transform to "raw scale" conc <- 10^temp$logc # Save the response values resp <- temp$resp # pull out all of the chemical identifiers and the assay name dtxsid <- temp[1,"dtxsid"] casrn <- temp[1,"casrn"] name <- temp[1,"name"] assay <- temp[1,"assay"] # Execute curve fitting # Input concentrations, responses, cutoff, a list of models to fit, and other model fitting requirements # force.fit is set to true so that all models will be fit regardless of cutoff # bidirectional = FALSE indicates only fit models in the positive direction. # if using bidirectional = TRUE the coff only needs to be specified in the positive direction. model_fits[[spid]] <- tcplfit2_core(conc, resp, cutoff, force.fit = TRUE, fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", "pow", "exp2","exp3", "exp4", "exp5"), bidirectional = FALSE) # Get a plot of all curve fits plt_lst[[spid]] <- plot_allcurves(model_fits[[spid]], conc = conc, resp = resp, log_conc = TRUE) # Pass the output from 'tcplfit2_core' to 'tcplhit2_core' along with # cutoff, onesd, and any identifiers out <- tcplhit2_core(model_fits[[spid]], conc, resp, bmed = 0, cutoff = cutoff, onesd = onesd, identifiers = c(dtxsid = dtxsid, casrn = casrn, name = name, assay = assay)) # store all results in one table result_table <- rbind(result_table,out) } ## ----example 2 fit results---------------------------------------------------- # shows the structure of the output object from tcplfit2_core (only top level) str(model_fits[[1]],max.lev = 1) ## ----hill_model_fit_str------------------------------------------------------- # structure of the model fit list - hill model results str(model_fits[[1]][["hill"]]) ## ----example2 plot1, fig.height = 9, fig.width = 7---------------------------- grid.arrange(grobs=plt_lst,ncol=2) ## ----echo=FALSE--------------------------------------------------------------- htmlTable::htmlTable(result_table, align = 'l', align.header = 'l', rnames = FALSE , css.cell = ' padding-bottom: 5px; vertical-align:top; padding-right: 10px;min-width: 5em ') ## ----example2 plot2----------------------------------------------------------- # plot the first row concRespPlot2(result_table[1,],log_conc = TRUE) + # add a descriptive title to the plot ggtitle(paste(result_table[1,"dtxsid"], result_table[1,"name"])) ## ----example3_init, fig.height = 6, fig.width = 7, message=FALSE, warning = FALSE---- # Loading in the Level 0 example data set from invitrodb data("mc0") data.table::setDTthreads(2) dat <- mc0 ## ----echo=FALSE--------------------------------------------------------------- # only show the top 6 rows for the treatment samples htmlTable::htmlTable(head(dat[wllt=='t',]), align = 'l', align.header = 'l', rnames = FALSE , css.cell = ' padding-bottom: 5px; vertical-align:top; padding-right: 10px;min-width: 5em ') ## ----example3_cndx, class.source="scroll-100", fig.height = 6, fig.width = 7, warning=FALSE---- # Order by the following columns setkeyv(dat, c('acid', 'srcf', 'apid', 'coli', 'rowi', 'spid', 'conc')) # Define a temporary replicate ID (rpid) column for test compound wells # rpid consists of the sample ID, well type (wllt), source file, assay plate ID, and # concentration. # the := operator is a data.table function to add/update rows nconc <- dat[wllt == "t" , ## denotes test well as the well type (wllt) list(n = lu(conc)), # total number of unique concentrations by = list(acid, apid, spid)][ , list(nconc = min(n)), by = acid] dat[wllt == "t" & acid %in% nconc[nconc > 1, acid], rpid := paste(acid, spid, wllt, srcf, apid, "rep1", conc, sep = "_")] dat[wllt == "t" & acid %in% nconc[nconc == 1, acid], rpid := paste(acid, spid, wllt, srcf, "rep1", conc, sep = "_")] # Define rpid column for non-test compound wells dat[wllt != "t", rpid := paste(acid, spid, wllt, srcf, apid, "rep1", conc, sep = "_")] # set the replicate index (repi) based on rowid # increment repi every time a replicate ID is duplicated dat[, dat_rpid := rowid(rpid)] dat[, rpid := sub("_rep[0-9]+.*", "",rpid, useBytes = TRUE)] dat[, rpid := paste0(rpid,"_rep",dat_rpid)] # For each replicate, define concentration index # by ranking the unique concentrations indexfunc <- function(x) as.integer(rank(unique(x))[match(x, unique(x))]) dat[ , cndx := indexfunc(conc), by = list(rpid)] ## ----example3_mc2, fig.height = 6, fig.width = 7------------------------------ # If no adjustments are required for the data, the corrected value (cval) should be set as original rval dat[,cval := rval] # Poor well quality (wllq) wells should be removed dat <- dat[!wllq == 0,] ##Fitting generally cannot occur if response values are NA therefore values need to be removed dat <- dat[!is.na(cval),] ## ----echo=FALSE--------------------------------------------------------------- htmlTable::htmlTable(head(tcpl::tcplMthdList(3)), align = 'l', align.header = 'l', rnames = FALSE , css.cell = ' padding-bottom: 5px; vertical-align:top; padding-right: 10px;min-width: 5em ') ## ----example3_normalize------------------------------------------------------- # calculate bval of the median of all the wells that have a type of n dat[, bval := median(cval[wllt == "n"]), by = list(apid)] # calculate pval based on the wells that have type of m or o excluding any NA wells dat[, pval := median(cval[wllt %in% c("m","o")], na.rm = TRUE), by = list(apid, wllt, conc)] # take pval as the minimum per assay plate (apid) dat[, pval := min(pval, na.rm = TRUE), by = list(apid)] # Calculate normalized responses dat[, resp := ((cval - bval)/(pval - bval) * 100)] ## ----echo=FALSE--------------------------------------------------------------- htmlTable::htmlTable(head(tcpl::tcplMthdList(4)), align = 'l', align.header = 'l', rnames = FALSE , css.cell = ' padding-bottom: 5px; vertical-align:top; padding-right: 10px;min-width: 5em ') ## ----example3_get_bmad.and.onesd---------------------------------------------- bmad <- mad(dat[cndx %in% c(1, 2) & wllt == "t", resp]) onesd <- sd(dat[cndx %in% c(1, 2) & wllt == "t", resp]) ## ----example3_fitting,class.source="fold-show"-------------------------------- #do tcplfit2 fitting myfun <- function(y) { res <- tcplfit2::tcplfit2_core(y$conc, y$resp, cutoff = 3*bmad, bidirectional = TRUE, verbose = FALSE, force.fit = TRUE, fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", "pow", "exp2", "exp3", "exp4", "exp5") ) list(list(res)) #use list twice because data.table uses list(.) to look for values to assign to columns } ## ----example3_fitting_full, eval=FALSE---------------------------------------- # # only want to run tcplfit2 for test wells in this case # # this chunk doesn't run, fit the curves on the subset below # dat[wllt == 't',params:= myfun(.SD), by = .(spid)] ## ----example3_fitting_subset,class.source="fold-show"------------------------- # create a subset that contains 6 samples and run curve fitting subdat <- dat[spid %in% unique(spid)[10:15],] subdat[wllt == 't',params:= myfun(.SD), by = .(spid)] ## ----hitc_plots,fig.height = 6, fig.width = 7,warning=FALSE,message=FALSE----- #### Data Set-Up #### # obtain the base example data DATA_CASE <- tcplfit2::signatures[1,] conc <- strsplit(DATA_CASE[,"conc"],split = "[|]") %>% unlist() %>% as.numeric() resp <- strsplit(DATA_CASE[,"resp"],split = "[|]") %>% unlist() %>% as.numeric() OG_data <- data.frame(xval = conc,yval = resp) %>% # obtain the concentrations that are outside the cutoff band dplyr::mutate(type = ifelse(abs(resp)>=abs(DATA_CASE[,"cutoff"]),"Extreme Responses",NA)) %>% mutate(.,df = "OG_data") # obtain the fit and best fitting/hitcalling information fit <- tcplfit2::tcplfit2_core(conc = conc,resp = resp, cutoff = DATA_CASE[,"cutoff"]) hit <- tcplfit2::tcplhit2_core(params = fit, conc = conc,resp = resp, cutoff = DATA_CASE[,"cutoff"], onesd = DATA_CASE[,"onesd"]) # obtain the continuous curve from fit information XC <- seq(from = min(conc),to = max(conc),length.out = 100) YC <- tcplfit2::exp4(x = XC,ps = unlist(fit$exp4[fit$exp4$pars])) # set up a continuous curve dataset cont_fit <- # best fit data.frame(xval = XC,yval = YC,type = "Best Fit") %>% # constant (flat) fit rbind.data.frame(data.frame(xval = XC,yval = rep(0,length(XC)),type = "Constant Fit")) ## prop weight 3 - continuous curve dataset addition ## # set up temporary data needed for re-scaling plot tmp_cutoff <- DATA_CASE[,"cutoff"] # cutoff value tmp_top <- fit$exp4$top # maximal predicted response from baseline tmp_ps <- unlist(fit$exp4[fit$exp4$pars]) # model parameters # code from toplikelihood.R lines 51-56 for the "exp4" model if (tmp_top == tmp_ps[1]) { # check if the top and tp are the same tmp_ps[1] = tmp_cutoff } else { x_top = acy(y = tmp_top, modpars = list(tp=tmp_ps[1],ga=tmp_ps[2],er=tmp_ps[3]),type="exp4") tmp_ps[1] = tmp_cutoff/( 1 - 2^(-x_top/tmp_ps[2])) } # obtain the rescaled predicted response YC_rescale <- tcplfit2::exp4(x = XC,ps = tmp_ps) # add the continuous rescaled curve to the continuous curve dataset cont_fit <- rbind.data.frame( cont_fit, data.frame(xval = XC,yval = YC_rescale,type = "Rescaled Best Fit") ) %>% mutate(.,df = "cont_fit") # dataset with reference lines (e.g. cutoff, bmr, top, etc.) ref_df <- data.frame( xval = rep(0,6), yval = c(hit$cutoff*c(-1,1), hit$bmr*c(-1,1), fit$exp4$top, hit$cutoff), type = c(rep("Cutoff",2),rep("BMR",2),"Top","Top at Cutoff") ) %>% mutate(.,df = "ref_df") ## plotting dataframe combined plot_highlight_df <- rbind.data.frame(OG_data,cont_fit,ref_df) #### Generate Plots #### ## Generate a Base Plot for the Concentration-Response ## base_plot <- ggplot2::ggplot()+ geom_point(data = dplyr::filter(plot_highlight_df,df == "OG_data"), aes(x = log10(xval),y = yval))+ geom_line(data = dplyr::filter(plot_highlight_df,df == "cont_fit" & type == "Best Fit"), aes(x = log10(xval),y = yval))+ geom_hline(data = dplyr::filter(plot_highlight_df,df == "ref_df" & type %in% c("Cutoff","BMR")), aes(yintercept = yval,linetype = type,colour = type))+ ggplot2::ylim(c(-1,1))+ scale_colour_manual(breaks = c("Cutoff","BMR"),values = rep("black",2))+ scale_linetype_manual(breaks = c("Cutoff","BMR"),values = c("dashed","dotted"))+ theme_bw()+ theme(axis.title.x = element_blank(),axis.title.y = element_blank()) ## Proportional Weight 1 Plot ## p1_plot <- base_plot+ # add a title for the subplot ggplot2::ggtitle("p1",subtitle = "AIC Weight")+ # add the constant (reference) and winning model (comparison) - highlighted geom_line(data = dplyr::filter(plot_highlight_df,df == "cont_fit" & type != "Rescaled Best Fit"), aes(x = log10(xval),y = yval,colour = type,linetype = type))+ scale_colour_manual(name = "", breaks = c("Constant Fit","Best Fit","Cutoff","BMR"), values = c("blue","red",rep("black",2)))+ scale_linetype_manual(name = "", breaks = c("Constant Fit","Best Fit","Cutoff","BMR"), values = c("solid","solid","dashed","dotted"))+ theme(legend.position = "inside", legend.position.inside = c(0.5,0.15), legend.key.size = unit(0.5,"cm"), legend.text = element_text(size = 7), legend.title = element_blank(), legend.background = element_rect(fill = alpha("lemonchiffon",0.5))) ## Proportional Weight 2 Plot ## p2_plot <- base_plot+ # add a title for the subplot ggplot2::ggtitle("p2",subtitle = "Responses Outside Cutoff")+ # add the concentrations with median responses outside the cutoff band - highlighted geom_point(data = dplyr::filter(plot_highlight_df,df == "OG_data" & type == "Extreme Responses"), aes(x = log10(xval),y = yval,shape = type),col = "red")+ # add the cutoff band - highlighted geom_hline(data = dplyr::filter(plot_highlight_df,df == "ref_df" & type %in% c("Cutoff","BMR")), aes(yintercept = yval,linetype = type,colour = type))+ scale_colour_manual(name = "", breaks = c("Cutoff","BMR"), values = c("blue","black"))+ scale_linetype_manual(name = "", breaks = c("Cutoff","BMR"), values = c("dashed","dotted"))+ scale_shape(name = "")+ theme(legend.position = "inside", legend.position.inside = c(0.5,0.15), legend.key.size = unit(0.5,"cm"), legend.spacing.y = unit(-4,"lines"), legend.text = element_text(size = 7), legend.title = element_blank(), legend.background = element_rect(fill = alpha("lemonchiffon",0.5))) ## Proportional Weight 3 Plot ## p3_plot <- base_plot+ # add a title for the subplot ggplot2::ggtitle("p3",subtitle = "Top Likelihood Ratio")+ # add the original predicted curve & the re-scaled predicted curve - highlighted ggplot2::geom_line(data = dplyr::filter(plot_highlight_df,df == "cont_fit" & type != "Constant Fit"), aes(x = log10(xval),y = yval,colour = type,linetype = type))+ # add the 'top' (maximal predicted change in response from baseline) & the cutoff band - highlighted ggplot2::geom_hline(data = dplyr::filter(plot_highlight_df,df == "ref_df"), aes(yintercept = yval,colour = type,linetype = type))+ scale_linetype_manual(name = "", breaks = c("Best Fit","Rescaled Best Fit","Cutoff","BMR","Top","Top at Cutoff"), values = c(rep("solid",2),"dashed","dotted",rep("dashed",2)))+ scale_colour_manual(name = "", breaks = c("Best Fit","Rescaled Best Fit","Cutoff","BMR","Top","Top at Cutoff"), values = c("blue","red",rep("black",2),"skyblue","hotpink"))+ theme(legend.position = "inside", legend.position.inside = c(0.5,0.175), legend.key.size = unit(0.5,"cm"), legend.text = element_text(size = 7), legend.title = element_blank(), legend.background = element_rect(fill = alpha("lemonchiffon",0.5))) ## All Plots ## grid.arrange(p1_plot,p2_plot,p3_plot, ncol = 3, top = paste(DATA_CASE[,"signature"],DATA_CASE[,"dtxsid"],sep = "\n"), left = "response", bottom = paste("log10(conc)", paste(paste("hitc:",signif(hit[,"hitcall"],3)), paste("log10(bmd):",signif(log10(hit[,"bmd"]),3)),sep = ", "), sep = "\n") ) ## ----example3_hitcalling------------------------------------------------------ #do tcplfit2 hitcalling myfun2 <- function(y) { res <- tcplfit2::tcplhit2_core(params = y$params[[1]], conc = y$conc, resp = y$resp, cutoff = 3*bmad, onesd = onesd ) list(list(res)) } # continue with hitcalling res <- subdat[wllt == 't', myfun2(.SD), by = .(spid)] # pivot wider res_wide <- rbindlist(Map(cbind, spid = res$spid, res$V1)) # add a binary hitcall column to the data res_wide[,hitb := ifelse(hitcall >= 0.9,1,0)] ## ----echo=FALSE--------------------------------------------------------------- htmlTable::htmlTable(head(res_wide), align = 'l', align.header = 'l', rnames = FALSE , css.cell = ' padding-bottom: 5px; vertical-align:top; padding-right: 10px;min-width: 5em ') ## ----example3_plot, fig.height = 8, fig.width = 7----------------------------- # allocate a place-holder object plt_list <- NULL # plot results using `concRespPlot` for(i in 1:nrow(res_wide)){ plt_list[[i]] <- concRespPlot2(res_wide[i,]) } # compile and display winning model plots for concentration-response series grid.arrange(grobs=plt_list,ncol=2) ## ----ex4_lower,warning=FALSE-------------------------------------------------- # We'll use data from mc3 in this section data("mc3") # determine the background variation # background is defined per the assay. In this case we use logc <= -2 # However, background should be defined in a way that makes sense for your application temp <- mc3[mc3$logc<= -2,"resp"] bmad <- mad(temp) onesd <- sd(temp) cutoff <- 3*bmad # load example data spid <- unique(mc3$spid)[94] ex_df <- mc3[is.element(mc3$spid,spid),] # The data file has stored concentration in log10 form, fix it conc <- 10^ex_df$logc # back-transforming concentrations on log10 scale resp <- ex_df$resp # modify the data for demonstration purposes conc2 <- conc[conc>0.41] resp2 <- resp[which(conc>0.41)] # pull out all of the chemical identifiers and the name of the assay dtxsid <- ex_df[1,"dtxsid"] casrn <- ex_df[1,"casrn"] name <- ex_df[1,"name"] assay <- ex_df[1,"assay"] # create the row object row_low <- list(conc = conc2, resp = resp2, bmed = 0, cutoff = cutoff, onesd = onesd, assay=assay, dtxsid=dtxsid,casrn=casrn,name=name) # run the concentration-response modeling for a single sample res_low <- concRespCore(row_low,fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", "pow", "exp2", "exp3", "exp4", "exp5"), bidirectional=F) # plotting the results min_conc <- min(conc2) concRespPlot2(res_low, log_conc = T) + geom_vline(aes(xintercept = log10(min_conc)),lty = "dashed")+ geom_rect(aes(xmin = log10(res_low[1, "bmdl"]), xmax = log10(res_low[1, "bmdu"]),ymin = 0,ymax = 30), alpha = 0.05,fill = "skyblue") + geom_segment(aes(x = log10(res_low[, "bmd"]), xend = log10(res_low[, "bmd"]), y = 0, yend = 30),col = "blue")+ ggtitle(label = paste(name,"-",assay),subtitle = dtxsid) ## ----ex4_lower-res------------------------------------------------------------ # function results res_low['Min. Conc.'] <- min(conc2) res_low['Name'] <- name res_low[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")] <- round(res_low[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")], 3) ## ----ex4_table, echo=FALSE---------------------------------------------------- DT::datatable(res_low[1, c("Name","Min. Conc.", "bmd", "bmdl", "bmdu")],rownames = FALSE) ## ----ex4_lower-demo,class.source="fold-show"---------------------------------- # using the argument to set a lower bound for BMD res_low2 <- concRespCore(row_low,fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", "pow", "exp2", "exp3", "exp4", "exp5"), bidirectional=F, bmd_low_bnd = 0.8) ## ----ex4_lower-res-bnd-------------------------------------------------------- # print out the new results # include previous results side by side for comparison res_low2['Min. Conc.'] <- min(conc2) res_low2['Name'] <- paste(name, "after `bounding`", sep = "-") res_low['Name'] <- paste(name, "before `bounding`", sep = "-") res_low2[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")] <- round(res_low2[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")], 3) output_low <- rbind(res_low[1, c('Name', "Min. Conc.", "bmd", "bmdl", "bmdu")], res_low2[1, c('Name', "Min. Conc.", "bmd", "bmdl", "bmdu")]) ## ----example_4_lower_res_table, echo = FALSE---------------------------------- DT::datatable(output_low,rownames = FALSE) ## ----ex4_lower-plot-bnd, class.source="scroll-100"---------------------------- # generate some concentrations for the fitted curve logc_plot <- seq(from=-3,to=2,by=0.05) conc_plot <- 10^logc_plot # initiate the plot plot(conc2,resp2,xlab="conc (uM)",ylab="Response",xlim=c(0.001,100),ylim=c(-5,60), log="x",main=paste(name,"\n",assay),cex.main=0.9) # add vertical lines to mark the minimum concentration in the data and the lower threshold set by bmd_low_bnd abline(v=min(conc2), lty = 1, col = "brown", lwd = 2) abline(v=res_low2$bmd, lty = 2, col = "darkviolet", lwd = 2) # add markers for BMD and its boundaries before `bounding` lines(c(res_low$bmd,res_low$bmd),c(0,50),col="green",lwd=2) rect(xleft=res_low$bmdl,ybottom=0,xright=res_low$bmdu,ytop=50,col=rgb(0,1,0, alpha = .5), border = NA) points(res_low$bmd, -0.5, pch = "x", col = "green") # add markers for BMD and its boundaries after `bounding` lines(c(res_low2$bmd,res_low2$bmd),c(0,50),col="blue",lwd=2) rect(xleft=res_low2$bmdl,ybottom=0,xright=res_low2$bmdu,ytop=50,col=rgb(0,0,1, alpha = .5), border = NA) points(res_low2$bmd, -0.5, pch = "x", col = "blue") # add the fitted curve lines(conc_plot, exp4(ps = c(res_low$tp, res_low$ga), conc_plot)) legend(1e-3, 60, legend=c("Lowest Dose Tested", "Boundary", "BMD-before", "BMD-after"), col=c("brown", "darkviolet", "green", "blue"), lty=c(1,2,1,1)) ## ----ex5_upper,warning=FALSE-------------------------------------------------- # load example data spid <- unique(mc3$spid)[26] ex_df <- mc3[is.element(mc3$spid,spid),] # The data file has stored concentration in log10 form, so fix that conc <- 10^ex_df$logc # back-transforming concentrations on log10 scale resp <- ex_df$resp # pull out all of the chemical identifiers and the name of the assay dtxsid <- ex_df[1,"dtxsid"] casrn <- ex_df[1,"casrn"] name <- ex_df[1,"name"] assay <- ex_df[1,"assay"] # create the row object row_up <- list(conc = conc, resp = resp, bmed = 0, cutoff = cutoff, onesd = onesd,assay=assay, dtxsid=dtxsid,casrn=casrn,name=name) # run the concentration-response modeling for a single sample res_up <- concRespCore(row_up,fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", "pow", "exp2", "exp3", "exp4", "exp5"), bidirectional=F) # plotting the results max_conc <- max(conc) concRespPlot2(res_up, log_conc = T) + # geom_vline(aes(xintercept = max(log10(conc))),lty = "dashed")+ geom_vline(aes(xintercept = log10(max_conc)),lty = "dashed")+ geom_rect(aes(xmin = log10(res_up[1, "bmdl"]), xmax = log10(res_up[1, "bmdu"]),ymin = 0,ymax = 125), alpha = 0.05,fill = "skyblue") + geom_segment(aes(x = log10(res_up[, "bmd"]), xend = log10(res_up[, "bmd"]), y = 0, yend = 125),col = "blue")+ ggtitle(label = paste(name,"-",assay),subtitle = dtxsid) ## ----ex5_upper-res------------------------------------------------------------ # max conc res_up['Max Conc.'] <- max(conc) res_up['Name'] <- name res_up[1, c("Max Conc.", "bmd", "bmdl", "bmdu")] <- round(res_up[1, c("Max Conc.", "bmd", "bmdl", "bmdu")], 3) # function results ## ----example_5_table, echo = FALSE-------------------------------------------- DT::datatable(res_up[1, c('Name','Max Conc.', "bmd", "bmdl", "bmdu")],rownames = FALSE) ## ----ex5_upper-demo,class.source="fold-show"---------------------------------- # using bmd_up_bnd = 2 res_up2 <- concRespCore(row_up,fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", "pow", "exp2", "exp3", "exp4", "exp5"), bidirectional=F, bmd_up_bnd = 2) ## ----ex5_upper-bnd------------------------------------------------------------ # print out the new results # include previous results side by side for comparison res_up2['Max Conc.'] <- max(conc) res_up2['Name'] <- paste(name, "after `bounding`", sep = "-") res_up['Name'] <- paste(name, "before `bounding`", sep = "-") res_up2[1, c("Max Conc.", "bmd", "bmdl", "bmdu")] <- round(res_up2[1, c("Max Conc.", "bmd", "bmdl", "bmdu")], 3) output_up <- rbind(res_up[1, c('Name', "Max Conc.", "bmd", "bmdl", "bmdu")], res_up2[1, c('Name', "Max Conc.", "bmd", "bmdl", "bmdu")]) ## ----example_upper_2_table, echo = FALSE-------------------------------------- DT::datatable(output_up,rownames = FALSE) ## ----ex5_upper-bnd-plot, class.source="scroll-100"---------------------------- # generate some concentration for the fitting curve logc_plot <- seq(from=-3,to=2,by=0.05) conc_plot <- 10^logc_plot # initiate plot plot(conc,resp,xlab="conc (uM)",ylab="Response",xlim=c(0.001,500),ylim=c(-5,150), log="x",main=paste(name,"\n",assay),cex.main=0.9) # add vertical lines to mark the maximum concentration in the data and the upper boundary set by bmd_up_bnd abline(v=max(conc), lty = 1, col = "brown", lwd=2) abline(v=160, lty = 2, col = "darkviolet", lwd=2) # add marker for BMD and its boundaries before `bounding` lines(c(res_up$bmd,res_up$bmd),c(0,125),col="green",lwd=2) rect(xleft=res_up$bmdl,ybottom=0,xright=res_up$bmdu,ytop=125,col=rgb(0,1,0, alpha = .5), border = NA) points(res_up$bmd, -0.5, pch = "x", col = "green") # add marker for BMD and its boundaries after `bounding` lines(c(res_up2$bmd,res_up2$bmd),c(0,125),col="blue",lwd=2) rect(xleft=res_up2$bmdl,ybottom=0,xright=res_up2$bmdu,ytop=125,col=rgb(0,0,1, alpha = .5), border = NA) points(res_up2$bmd, -0.5, pch = "x", col = "blue") # add the fitting curve lines(conc_plot, poly1(ps = c(res_up$a), conc_plot)) legend(1e-3, 150, legend=c("Maximum Dose Tested", "Boundary", "BMD-before", "BMD-after"), col=c("brown", "darkviolet", "green", "blue"), lty=c(1,2,1,1)) ## ----ex6_hitcore,class.source="fold-show"------------------------------------- # using the same data, fit curves param <- tcplfit2_core(conc2, resp2, cutoff = cutoff) hit_res <- tcplhit2_core(param, conc2, resp2, cutoff = cutoff, onesd = onesd, bmd_low_bnd = 0.8) ## ----ex6_hitcore-res---------------------------------------------------------- # adding the result from tcplhit2_core to the output table for comparison hit_res["Name"]<- paste("Chlorothalonil", "tcplhit2_core", sep = "-") hit_res['Min. Conc.'] <- min(conc2) hit_res[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")] <- round(hit_res[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")], 3) output_low <- rbind(output_low, hit_res[1, c('Name', "Min. Conc.", "bmd", "bmdl", "bmdu")]) ## ----ex6_res-hit_table, echo = FALSE------------------------------------------ DT::datatable(output_low,rownames = FALSE) ## ----ex7_lower-bnd,class.source="fold-show"----------------------------------- res_low3 <- concRespCore(row_low,fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", "pow", "exp2", "exp3", "exp4", "exp5"), conthits = T, aicc = F, bidirectional=F, bmd_low_bnd = 0.4) ## ----ex7_lower-bnd-res-------------------------------------------------------- # print out the new results # add to previous results for comparison res_low3['Min. Conc.'] <- min(conc2) res_low3['Name'] <- paste("Chlorothalonil", "after `bounding` (two fifths)", sep = "-") res_low3[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")] <- round(res_low3[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")], 3) output_low <- rbind(output_low[-3, ], res_low3[1, c('Name', "Min. Conc.", "bmd", "bmdl", "bmdu")]) ## ----ex7_lower-bnd-res-table, echo = FALSE------------------------------------ DT::datatable(output_low,rownames = FALSE) ## ----ex7_lower-bnd-plot, class.source="scroll-100"---------------------------- # initiate the plot plot(conc2,resp2,xlab="conc (uM)",ylab="Response",xlim=c(0.001,100),ylim=c(-5,60), log="x",main=paste(name,"\n",assay),cex.main=0.9) # add vertical lines to mark the minimum concentration in the data and the lower boundary set by bmd_low_bnd abline(v=min(conc2), lty = 1, col = "brown", lwd = 2) abline(v=0.4*min(conc2), lty = 2, col = "darkviolet", lwd = 2) # add markers for BMD and its boundaries before `bounding` lines(c(res_low$bmd,res_low$bmd),c(0,50),col="green",lwd=2) rect(xleft=res_low$bmdl,ybottom=0,xright=res_low$bmdu,ytop=50,col=rgb(0,1,0, alpha = .5), border = NA) points(res_low$bmd, 0, pch = "x", col = "green") # add markers for BMD and its boundaries after `bounding` lines(c(res_low3$bmd,res_low3$bmd),c(0,50),col="blue",lwd=2) rect(xleft=res_low3$bmdl,ybottom=0,xright=res_low3$bmdu,ytop=50,col=rgb(0,0,1, alpha = .5), border = NA) points(res_low3$bmd, 0, pch = "x", col = "blue") # add the fitted curve lines(conc_plot, exp4(ps = c(res_low$tp, res_low$ga), conc_plot)) legend(1e-3, 60, legend=c("Lowest Dose Tested", "Boundary Dose", "BMD-before", "BMD-after"), col=c("brown", "darkviolet", "green", "blue"), lty=c(1,2,1,1)) ## ----appendix_plt1, fig.height = 6, fig.width = 7, warning = FALSE------------ # read in the file data("signatures") # set up a 3 x 2 grid for the plots oldpar <- par(no.readonly = TRUE) on.exit(par(oldpar)) par(mfrow=c(3,2),mar=c(4,4,5,2)) # fit 6 observations in signatures for(i in 1:nrow(signatures)){ # set up input data row = list(conc=as.numeric(str_split(signatures[i,"conc"],"\\|")[[1]]), resp=as.numeric(str_split(signatures[i,"resp"],"\\|")[[1]]), bmed=0, cutoff=signatures[i,"cutoff"], onesd=signatures[i,"onesd"], name=signatures[i,"name"], assay=signatures[i,"signature"]) # run concentration-response modeling (1st plotting option) out = concRespCore(row,conthits=F,do.plot=T) if(i==1){ res <- out }else{ res <- rbind.data.frame(res,out) } } ## ----appendix plt2, fig.height = 8, fig.width = 7----------------------------- # set up a 3 x 2 grid for the plots oldpar <- par(no.readonly = TRUE) on.exit(par(oldpar)) par(mfrow=c(3,2),mar=c(4,4,2,2)) # plot results using `concRespPlot` for(i in 1:nrow(res)){ concRespPlot(res[i,],ymin=-1,ymax=1) } ## ----------------------------------------------------------------------------- # Load the example data set data("signatures") # using the first row of signature as an example conc <- as.numeric(str_split(signatures[1,"conc"],"\\|")[[1]]) resp <- as.numeric(str_split(signatures[1,"resp"],"\\|")[[1]]) cutoff <- signatures[1,"cutoff"] # run curve fitting output <- tcplfit2_core(conc, resp, cutoff) # show the structure of the output summary(output) ## ----class.source="fold-show",fig.height=8,fig.width=7------------------------ # get plots in the original and in log-10 concentration scale basic <- plot_allcurves(output, conc, resp) basic_log <- plot_allcurves(output, conc, resp, log_conc = T) # arrange the ggplot2 output into a grid grid.arrange(basic, basic_log) ## ----------------------------------------------------------------------------- # prepare the 'row' object for concRespCore row <- list(conc=conc, resp=resp, bmed=0, cutoff=cutoff, onesd=signatures[1,"onesd"], name=signatures[1,"name"], assay=signatures[1,"signature"]) # run concentration-response modeling out <- concRespCore(row,conthits=F) # show the output out ## ----class.source="fold-show"------------------------------------------------- # pass the output to the plotting function basic_plot <- concRespPlot2(out) basic_log <- concRespPlot2(out, log_conc = TRUE) # arrange the ggplot2 output into a grid grid.arrange(basic_plot, basic_log) ## ----------------------------------------------------------------------------- # Using the fitted result and plot from the example in the last section # get the cutoff from the output cutoff <- out[, "cutoff"] basic_plot + # Cutoff Band - a transparent rectangle geom_rect(aes(xmin = 0,xmax = 30,ymin = -cutoff,ymax = cutoff), alpha = 0.1,fill = "skyblue") + # Titles ggtitle( label = paste("Best Model Fit", out[, "name"], sep = "\n"), subtitle = paste("Assay Endpoint: ", out[, "assay"])) + ## Add BMD and BMR labels geom_hline( aes(yintercept = out[, "bmr"]), col = "blue") + geom_segment( aes(x = out[, "bmd"], xend = out[, "bmd"], y = -0.5, yend = out[, "bmr"]), col = "blue" ) + geom_point(aes(x = out[, "bmd"], y = out[, "bmr"], fill = "BMD"), shape = 21, cex = 2.5) ## ----------------------------------------------------------------------------- # Get all potency estimates and the corresponding y value on the curve estimate_points <- out %>% select(bmd, acc, ac50, ac10, ac5) %>% tidyr::pivot_longer(everything(), names_to = "Potency Estimates") %>% mutate(`Potency Estimates` = toupper(`Potency Estimates`)) y <- c(out[, "bmr"], out[, "cutoff"], rep(out[, "top"], 3)) y <- y * c(1, 1, .5, .1, .05) estimate_points <- cbind(estimate_points, y = y) # add Potency Estimate Points and set colors basic_plot + geom_point( data = estimate_points, aes(x = value, y = y, fill = `Potency Estimates`), shape = 21, cex = 2.5 ) ## ----------------------------------------------------------------------------- # add Potency Estimate Points and set colors - with plot in log-10 concentration basic_log + geom_point( data = estimate_points, aes(x = log10(value), y = y, fill = `Potency Estimates`), shape = 21, cex = 2.5 ) ## ----------------------------------------------------------------------------- # maybe want to extract and use the same x's in the base plot # to calculate predicted responses conc_plot <- basic_plot[["layers"]][[2]][["data"]][["conc_plot"]] basic_plot + # fitted parameter values of another curve you want to add geom_line(data=data.frame(x=conc_plot, y=tcplfit2::exp5(c(0.5, 10, 1.2), conc_plot)), aes(x,y,color = "exp5"))+ # add different colors for comparisons scale_colour_manual(values=c("#CC6666", "#9999CC"), labels = c("Curve 1-exp4", "Curve 2-exp5")) + labs(title = "Curve 1 v.s. Curve 2") ## ----ex1_AUC,class.source="fold-show"----------------------------------------- # some example data conc <- list(.03, .1, .3, 1, 3, 10, 30, 100) resp <- list(0, .2, .1, .4, .7, .9, .6, 1.2) row <- list(conc = conc, resp = resp, bmed = 0, cutoff = 1, onesd = .5) # AUC is included in the output concRespCore(row, conthits = TRUE, AUC = TRUE) ## ----ex2_AUC, fig.height = 4.55, fig.width = 8-------------------------------- # This is taken from the example under tcplfit2_core conc_ex2 <- c(0.03, 0.1, 0.3, 1, 3, 10, 30, 100) resp_ex2 <- c(0, 0.1, 0, 0.2, 0.6, 0.9, 1.1, 1) # fit all available models in the package # show all fitted curves output_ex2 <- tcplfit2_core(conc_ex2, resp_ex2, 0.8) # arrange the ggplot2 output into a grid grid.arrange(plot_allcurves(output_ex2, conc_ex2, resp_ex2), plot_allcurves(output_ex2, conc_ex2, resp_ex2, log_conc = TRUE), ncol = 2) ## ----ex2_AUC_posthit,class.source="fold-show"--------------------------------- # hitcalling results out <- tcplhit2_core(output_ex2, conc_ex2, resp_ex2, 0.8, onesd = 0.4) out # perform AUC estimation post_hit_AUC(out) ## ----ex2_AUC-getAUC,class.source="fold-show"---------------------------------- fit_method <- "hill" # extract the parameters modpars <- output_ex2[[fit_method]][output_ex2[[fit_method]]$pars] # plug into get_AUC function estimated_auc1 <- get_AUC(fit_method, min(conc_ex2), max(conc_ex2), modpars) estimated_auc1 # extract the predicted responses from the model pred_resp <- output_ex2[[fit_method]][["modl"]] ## ----ex2_AUC-getAUC-plot,fig.height = 6, fig.width = 6------------------------ # plot to see if the result make sense # the shaded area is what the function tries to find plot(conc_ex2, pred_resp,ylim = c(0,1), xlab = "Concentration",ylab = "Response",main = "Positive Response AUC") lines(conc_ex2, pred_resp) polygon(c(conc_ex2, max(conc_ex2)), c(pred_resp, min(pred_resp)), col=rgb(1,0,0,0.5)) ## ----ex2_AUC-other-models----------------------------------------------------- # list of models fitmodels <- c("gnls", "poly1", "poly2", "pow", "exp2", "exp3", "exp4", "exp5") mylist <- list() for (model in fitmodels){ fit_method <- model # extract corresponding model parameters modpars <- output_ex2[[fit_method]][output_ex2[[fit_method]]$pars] # get AUC mylist[[fit_method]] <- get_AUC(fit_method, min(conc_ex2), max(conc_ex2), modpars) } # print AUC's for other models data.frame(mylist,row.names = "AUC") ## ----ex3_AUC, fig.height = 4.55, fig.width = 8-------------------------------- # use row 5 in the data conc <- as.numeric(str_split(signatures[5,"conc"],"\\|")[[1]]) resp <- as.numeric(str_split(signatures[5,"resp"],"\\|")[[1]]) cutoff <- signatures[5,"cutoff"] # plot all models, this is an example of negative curves output_negative <- tcplfit2_core(conc, resp, cutoff) grid.arrange(plot_allcurves(output_negative, conc, resp), plot_allcurves(output_negative, conc, resp, log_conc = TRUE), ncol = 2) ## ----ex3_AUC-getAUC,class.source="fold-show"---------------------------------- # choose fit method fit_method <- "exp3" # extract corresponding model parameters and predicted response modpars <- output_negative[[fit_method]][output_negative[[fit_method]]$pars] pred_resp <- output_negative[[fit_method]][["modl"]] estimated_auc2 <- get_AUC(fit_method, min(conc), max(conc), modpars) estimated_auc2 ## ----ex3_AUC-plot,fig.height = 6, fig.width = 6------------------------------- # plot this curve pred_resp <- pred_resp[order(conc)] plot(conc[order(conc)], pred_resp,ylim = c(-1,0), xlab = "Concentration",ylab = "Response",main = "Negative Response AUC") lines(conc[order(conc)], pred_resp) polygon(c(conc[order(conc)], max(conc)), c(pred_resp, max(pred_resp)), col=rgb(1,0,0,0.5)) ## ----ex3_AUC-abs-neg,class.source="fold-show"--------------------------------- get_AUC(fit_method, min(conc), max(conc), modpars, return.abs = TRUE) ## ----ex4_AUC, fig.height = 6, fig.width = 6----------------------------------- # simulate a poly2 curve conc_sim <- seq(0,3, length.out = 100) ## biphasic poly2 parameters b1 <- -1.3 b2 <- 0.7 ## converted to tcplfit2's poly2 parameters a <- b1^2/b2 b <- b1/b2 c(a,b) ## plot the curve resp_sim <- poly2(c(a, b, 0.1), conc_sim) plot(conc_sim, resp_sim, type = "l", xlab = "Concentration",ylab = "Response",main = "Biphasic Response") abline(h = 0) ## ----ex4_AUC-cont,class.source="fold-show"------------------------------------ # get AUC for the simulated Polynomial 2 curve get_AUC("poly2", min(conc_sim), max(conc_sim), ps = c(a, b)) ## ----fig.height = 6, fig.width = 6-------------------------------------------- ## plot the curve for the AUC plot(conc_sim, resp_sim, type = "l", xlab = "Concentration",ylab = "Response",main = "Biphasic Response AUC") abline(h = 0) polygon(c(conc_sim[which(resp_sim <= 0)], max(conc_sim[which(resp_sim <= 0)])), c(resp_sim[which(resp_sim <= 0)], max(resp_sim[which(resp_sim <= 0)])), col="skyblue") polygon(c(conc_sim[c(max(which(resp_sim <= 0)),which(resp_sim > 0))], max(conc_sim[which(resp_sim > 0)])), c(0,resp_sim[which(resp_sim > 0)], 0), col="indianred") ## ----setup-2, warning=FALSE--------------------------------------------------- # prepare concentration data for demonstration ex_conc <- seq(0, 100, length.out = 500) ex2_conc <- seq(0, 3, length.out = 100) ## ----poly-1,fig.width=5,fig.height=5,warning=FALSE---------------------------- poly1_plot <- ggplot(mapping=aes(ex_conc)) + geom_line(aes(y = 55*ex_conc, color = "a=55")) + geom_line(aes(y = 10*ex_conc, color = "a=10")) + geom_line(aes(y = 0.05*ex_conc, color = "a=0.05")) + geom_line(aes(y = -5*ex_conc, color = "a=(-5)")) + labs(x = "Concentration", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.8)) + scale_color_manual(name='a values', breaks=c('a=(-5)', 'a=0.05', 'a=10', 'a=55'), values=c('a=(-5)'='black', 'a=0.05' = 'red', 'a=10'='blue', 'a=55'='darkviolet')) poly1_plot ## ----poly-2, fig.width=8, fig.height=5, warning=FALSE------------------------- fits_poly <- data.frame( # change a y1 = poly2(ps = c(a = 40, b = 2),x = ex_conc), y2 = poly2(ps = c(a = 6, b = 2),x = ex_conc), y3 = poly2(ps = c(a = 0.1, b = 2),x = ex_conc), y4 = poly2(ps = c(a = -2, b = 2),x = ex_conc), y5 = poly2(ps = c(a = -20, b = 2),x = ex_conc), # change b y6 = poly2(ps = c(a = 4,b = 1.8),x = ex_conc), y7 = poly2(ps = c(a = 4,b = 7),x = ex_conc), y8 = poly2(ps = c(a = 4,b = 16),x = ex_conc) ) # shows how changes in parameter 'a' affect the shape of the curve poly2_plot1 <- ggplot(fits_poly, aes(ex_conc)) + geom_line(aes(y = y1, color = "a=40")) + geom_line(aes(y = y2, color = "a=6")) + geom_line(aes(y = y3, color = "a=0.1")) + geom_line(aes(y = y4, color = "a=(-2)")) + geom_line(aes(y = y5, color = "a=(-20)")) + labs(x = "Concentration", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.8)) + scale_color_manual(name='a values', breaks=c('a=(-20)', 'a=(-2)', 'a=0.1', 'a=6', 'a=40'), values=c('a=(-20)'='black', 'a=(-2)'='red', 'a=0.1'='blue', 'a=6'='darkviolet', 'a=40'='darkgoldenrod1')) # shows how changes in parameter 'b' affect the shape of the curve poly2_plot2 <- ggplot(fits_poly, aes(ex_conc)) + geom_line(aes(y = y6, color = "b=1.8")) + geom_line(aes(y = y7, color = "b=7")) + geom_line(aes(y = y8, color = "b=16")) + labs(x = "Concentration", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.8)) + scale_color_manual(name='b values', breaks=c('b=1.8', 'b=7', 'b=16'), values=c('b=1.8'='black', 'b=7'='red', 'b=16'='blue')) grid.arrange(poly2_plot1, poly2_plot2, ncol = 2) ## ----pow, fig.width=8, fig.height=5, warning=FALSE---------------------------- fits_pow <- data.frame( # change a y1 = pow(ps = c(a = 0.48,p = 1.45),x = ex2_conc), y2 = pow(ps = c(a = 7.2,p = 1.45),x = ex2_conc), y3 = pow(ps = c(a = -3.2,p = 1.45),x = ex2_conc), # change p y4 = pow(ps = c(a = 1.2,p = 0.3),x = ex2_conc), y5 = pow(ps = c(a = 1.2,p = 1.6),x = ex2_conc), y6 = pow(ps = c(a = 1.2,p = 3.2),x = ex2_conc) ) # shows how changes in parameter 'a' affect the shape of the curve pow_plot1 <- ggplot(fits_pow, aes(ex2_conc)) + geom_line(aes(y = y1, color = "a=0.48")) + geom_line(aes(y = y2, color = "a=7.2")) + geom_line(aes(y = y3, color = "a=(-3.2)")) + labs(x = "Concentration", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.8)) + scale_color_manual(name='a values', breaks=c('a=(-3.2)', 'a=0.48', 'a=7.2'), values=c('a=(-3.2)'='black', 'a=0.48'='red', 'a=7.2'='blue')) # shows how changes in parameter 'p' affect the shape of the curve pow_plot2 <- ggplot(fits_pow, aes(ex2_conc)) + geom_line(aes(y = y4, color = "p=0.3")) + geom_line(aes(y = y5, color = "p=1.6")) + geom_line(aes(y = y6, color = "p=3.2")) + labs(x = "Concentration", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.8)) + scale_color_manual(name='p values', breaks=c('p=0.3', 'p=1.6', 'p=3.2'), values=c('p=0.3'='black', 'p=1.6'='red', 'p=3.2'='blue')) grid.arrange(pow_plot1, pow_plot2, ncol = 2) ## ----Hill, fig.height=5, fig.width=8, warning=FALSE--------------------------- fits_hill <- data.frame( # change tp y1 = hillfn(ps = c(tp = -200,ga = 5,p = 1.76), x = ex_conc), y2 = hillfn(ps = c(tp = 200,ga = 5,p = 1.76), x = ex_conc), y3 = hillfn(ps = c(tp = 850,ga = 5,p = 1.76), x = ex_conc), # change ga y4 = hillfn(ps = c(tp = 120,ga = 4,p = 1.76), x = ex_conc), y5 = hillfn(ps = c(tp = 120,ga = 12,p = 1.76), x = ex_conc), y6 = hillfn(ps = c(tp = 120,ga = 20,p = 1.76), x = ex_conc), # change p y7 = hillfn(ps = c(tp = 120,ga = 5,p = 0.5), x = ex_conc), y8 = hillfn(ps = c(tp = 120,ga = 5,p = 2), x = ex_conc), y9 = hillfn(ps = c(tp = 120,ga = 5,p = 5), x = ex_conc) ) # shows how changes in parameter 'tp' affect the shape of the curve hill_plot1 <- ggplot(fits_hill, aes(log10(ex_conc))) + geom_line(aes(y = y1, color = "tp=(-200)")) + geom_line(aes(y = y2, color = "tp=200")) + geom_line(aes(y = y3, color = "tp=850")) + labs(x = "Concentration in Log-10 Scale", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.7), legend.key.size = unit(0.5, 'cm')) + scale_color_manual(name='tp values', breaks=c('tp=(-200)', 'tp=200', 'tp=850'), values=c('tp=(-200)'='black', 'tp=200'='red', 'tp=850'='blue')) # shows how changes in parameter 'ga' affect the shape of the curve hill_plot2 <- ggplot(fits_hill, aes(log10(ex_conc))) + geom_line(aes(y = y4, color = "ga=4")) + geom_line(aes(y = y5, color = "ga=12")) + geom_line(aes(y = y6, color = "ga=20")) + labs(x = "Concentration in Log-10 Scale", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.7), legend.key.size = unit(0.4, 'cm')) + scale_color_manual(name='ga values', breaks=c('ga=4', 'ga=12', 'ga=20'), values=c('ga=4'='black', 'ga=12'='red', 'ga=20'='blue')) # shows how changes in parameter 'p' affect the shape of the curve hill_plot3 <- ggplot(fits_hill, aes(log10(ex_conc))) + geom_line(aes(y = y7, color = "p=0.5")) + geom_line(aes(y = y8, color = "p=2")) + geom_line(aes(y = y9, color = "p=5")) + labs(x = "Concentration in Log-10 Scale", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.7), legend.key.size = unit(0.4, 'cm')) + scale_color_manual(name='p values', breaks=c('p=0.5', 'p=2', 'p=5'), values=c('p=0.5'='black', 'p=2'='red', 'p=5'='blue')) grid.arrange(hill_plot1, hill_plot2, hill_plot3, ncol = 2, nrow = 2) ## ----gnls, fig.width=8, fig.height=5, warning=FALSE--------------------------- fits_gnls <- data.frame( # change la y1 = gnls(ps = c(tp = 750,ga = 15,p = 1.45,la = 17,q = 1.34), x = ex_conc), y2 = gnls(ps = c(tp = 750,ga = 15,p = 1.45,la = 50,q = 1.34), x = ex_conc), y3 = gnls(ps = c(tp = 750,ga = 15,p = 1.45,la = 100,q = 1.34), x = ex_conc), # change q y4 = gnls(ps = c(tp = 750,ga = 15,p = 1.45,la = 20,q = 0.3), x = ex_conc), y5 = gnls(ps = c(tp = 750,ga = 15,p = 1.45,la = 20,q = 1.2), x = ex_conc), y6 = gnls(ps = c(tp = 750,ga = 15,p = 1.45,la = 20,q = 8), x = ex_conc) ) # shows how changes in parameter 'la' affect the shape of the curve gnls_plot1 <- ggplot(fits_gnls, aes(log10(ex_conc))) + geom_line(aes(y = y1, color = "la=17")) + geom_line(aes(y = y2, color = "la=50")) + geom_line(aes(y = y3, color = "la=100")) + labs(x = "Concentration in Log-10 Scale", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.8)) + scale_color_manual(name='la values', breaks=c('la=17', 'la=50', 'la=100'), values=c('la=17'='black', 'la=50'='red', 'la=100'='blue')) # shows how changes in parameter 'q' affect the shape of the curve gnls_plot2 <- ggplot(fits_gnls, aes(log10(ex_conc))) + geom_line(aes(y = y4, color = "q=0.3")) + geom_line(aes(y = y5, color = "q=1.2")) + geom_line(aes(y = y6, color = "q=8")) + labs(x = "Concentration in Log-10 Scale", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.8)) + scale_color_manual(name='q values', breaks=c('q=0.3', 'q=1.2', 'q=8'), values=c('q=0.3'='black', 'q=1.2'='red', 'q=8'='blue')) grid.arrange(gnls_plot1, gnls_plot2, ncol = 2) ## ----exp2, fig.width=8, fig.height=5, warning=FALSE--------------------------- fits_exp2 <- data.frame( # change a y1 = exp2(ps = c(a = 20,b = 12), x = ex2_conc), y2 = exp2(ps = c(a = 9,b = 12), x = ex2_conc), y3 = exp2(ps = c(a = 0.1,b = 12), x = ex2_conc), y4 = exp2(ps = c(a = -3,b = 12), x = ex2_conc), # change b y5 = exp2(ps = c(a = 0.45,b = 4), x = ex2_conc), y6 = exp2(ps = c(a = 0.45,b = 9), x = ex2_conc), y7 = exp2(ps = c(a = 0.45,b = 20), x = ex2_conc) ) # shows how changes in parameter 'a' affect the shape of the curve exp2_plot1 <- ggplot(fits_exp2, aes(ex2_conc)) + geom_line(aes(y = y1, color = "a=20")) + geom_line(aes(y = y2, color = "a=9")) + geom_line(aes(y = y3, color = "a=0.1")) + geom_line(aes(y = y4, color = "a=(-3)")) + labs(x = "Concentration", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.8)) + scale_color_manual(name='a values', breaks=c('a=(-3)', 'a=0.1', 'a=9', 'a=20'), values=c('a=(-3)'='black', 'a=0.1'='red', 'a=9'='blue', 'a=20'='darkviolet')) # shows how changes in parameter 'b' affect the shape of the curve exp2_plot2 <- ggplot(fits_exp2, aes(ex2_conc)) + geom_line(aes(y = y5, color = "b=4")) + geom_line(aes(y = y6, color = "b=9")) + geom_line(aes(y = y7, color = "b=20")) + labs(x = "Concentration", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.8)) + scale_color_manual(name='b values', breaks=c('b=4', 'b=9', 'b=20'), values=c('b=4'='black', 'b=9'='red', 'b=20'='blue')) grid.arrange(exp2_plot1, exp2_plot2, ncol = 2) ## ----exp3, fig.width=5, fig.height=5, warning=FALSE--------------------------- fits_exp3 <- data.frame( # change p y1 = exp3(ps = c(a = 1.67,b = 12.5,p = 0.3), x = ex2_conc), y2 = exp3(ps = c(a = 1.67,b = 12.5,p = 0.9), x = ex2_conc), y3 = exp3(ps = c(a = 1.67,b = 12.5,p = 1.2), x = ex2_conc) ) # shows how changes in parameter 'p' affect the shape of the curve exp3_plot <- ggplot(fits_exp3, aes(ex2_conc)) + geom_line(aes(y = y1, color = "p=0.3")) + geom_line(aes(y = y2, color = "p=0.9")) + geom_line(aes(y = y3, color = "p=1.2")) + labs(x = "Concentration", y = "Response") + theme_bw()+ theme(legend.position = c(0.15,0.8)) + scale_color_manual(name='p values', breaks=c('p=0.3', 'p=0.9', 'p=1.2'), values=c('p=0.3'='black', 'p=0.9'='red', 'p=1.2'='blue')) exp3_plot ## ----exp4, fig.width=8, fig.height=5, warning=FALSE--------------------------- fits_exp4 <- data.frame( # change tp y1 = exp4(ps = c(tp = 895,ga = 15),x = ex_conc), y2 = exp4(ps = c(tp = 200,ga = 15),x = ex_conc), y3 = exp4(ps = c(tp = -500,ga = 15),x = ex_conc), # change ga y4 = exp4(ps = c(tp = 500,ga = 0.4),x = ex_conc), y5 = exp4(ps = c(tp = 500,ga = 10),x = ex_conc), y6 = exp4(ps = c(tp = 500,ga = 20),x = ex_conc) ) # shows how changes in parameter 'tp' affect the shape of the curve exp4_plot1 <- ggplot(fits_exp4, aes(ex_conc)) + geom_line(aes(y = y1, color = "tp=895")) + geom_line(aes(y = y2, color = "tp=200")) + geom_line(aes(y = y3, color = "tp=(-500)")) + labs(x = "Concentration", y = "Response") + theme_bw()+ theme(legend.position = c(0.8,0.2)) + scale_color_manual(name='tp values', breaks=c('tp=(-500)', 'tp=200', 'tp=895'), values=c('tp=(-500)'='black', 'tp=200'='red', 'tp=895'='blue')) # shows how changes in parameter 'ga' affect the shape of the curve exp4_plot2 <- ggplot(fits_exp4, aes(ex_conc)) + geom_line(aes(y = y4, color = "ga=0.4")) + geom_line(aes(y = y5, color = "ga=10")) + geom_line(aes(y = y6, color = "ga=20")) + labs(x = "Concentration", y = "Response") + theme_bw()+ theme(legend.position = c(0.8,0.2)) + scale_color_manual(name='ga values', breaks=c('ga=0.4', 'ga=10', 'ga=20'), values=c('ga=0.4'='black', 'ga=10'='red', 'ga=20'='blue')) grid.arrange(exp4_plot1, exp4_plot2, ncol = 2) ## ----exp5, fig.width=5, fig.height=5, warning=FALSE--------------------------- fits_exp5 <- data.frame( # change p y1 = exp5(ps = c(tp = 793,ga = 6.25,p = 0.3), x = ex_conc), y2 = exp5(ps = c(tp = 793,ga = 6.25,p = 3.4), x = ex_conc), y3 = exp5(ps = c(tp = 793,ga = 6.25,p = 8), x = ex_conc) ) # shows how changes in parameter 'p' affect the shape of the curve exp5_plot <- ggplot(fits_exp5, aes(ex_conc)) + geom_line(aes(y = y1, color = "p=0.3")) + geom_line(aes(y = y2, color = "p=3.4")) + geom_line(aes(y = y3, color = "p=8")) + labs(x = "Concentration", y = "Response") + theme_bw()+ theme(legend.position = c(0.8,0.2)) + scale_color_manual(name='p values', breaks=c('p=0.3', 'p=3.4', 'p=8'), values=c('p=0.3'='black', 'p=3.4'='red', 'p=8'='blue')) exp5_plot ## ----appendix-table, echo=FALSE, warning=FALSE-------------------------------- # First column - tcplfit2 available models. Model <- c( "Constant", "Linear", "Quadratic","Power", "Hill", "Gain-Loss", "Exponential 2", "Exponential 3","Exponential 4", "Exponential 5" ) # Second column - model abbreviations used in invitrodb & tcplfit2. Abbreviation <- c( "cnst", "poly1", "poly2","pow", "hill", "gnls", "exp2", "exp3", "exp4", "exp5" ) # Third column - model equations. Equations <- c( "$f(x) = 0$", # constant "$f(x) = ax$", # linear "$f(x) = a(\\frac{x}{b}+(\\frac{x}{b})^{2})$", # quadratic "$f(x) = ax^p$", # power "$f(x) = \\frac{tp}{1 + (\\frac{ga}{x})^{p}}$", # hill "$f(x) = \\frac{tp}{(1 + (\\frac{ga}{x})^{p} )(1 + (\\frac{x}{la})^{q} )}$", # gain-loss "$f(x) = a*(exp(\\frac{x}{b}) - 1)$", # exp 2 "$f(x) = a*(exp((\\frac{x}{b})^{p}) - 1)$", # exp 3 "$f(x) = tp*(1-2^{\\frac{-x}{ga}})$", # exp 4 "$f(x) = tp*(1-2^{-(\\frac{x}{ga})^{p}})$" # exp 5 ) # Fourth column - model parameter descriptions. OutputParameters <- c( "", # constant "a (y-scale)", # linear, "a (y-scale)
b (x-scale)", # quadratic "a (y-scale)
p (power)", # power "tp (top parameter)
ga (gain AC50)
p (gain-power)", # hill "tp (top parameter)
ga (gain AC50)
p (gain power)
la (loss AC50)
q (loss power)", # gain-loss "a (y-scale)
b (x-scale)", # exp2 "a (y-scale)
b (x-scale)
p (power)", # exp3 "tp (top parameter)
ga (AC50)", # exp4 "tp (top parameter)
ga (AC50)
p (power)" # exp5 ) # Fifth column - additional model details. Details <- c( "Parameters always equals 'er'.", # constant "", # linear "", # quadratic "", # power "Concentrations are converted internally to log10 units and optimized with f(x) = tp/(1 + 10^(p*(gax))), then ga and ga_sd are converted back to regular units before returning.", # hill "Concentrations are converted internally to log10 units and optimized with f(x) = tp/[(1 + 10^(p*(gax)))(1 + 10^(q*(x-la)))], then ga, la, ga_sd, and la_sd are converted back to regular units before returning." , # gain-loss "", # exp2 "", # exp3 "", # exp4 "") # exp5 # Consolidate all columns into a table. output <- data.frame(Model, Abbreviation, Equations, OutputParameters, Details) # Export/print the table into an html rendered table. htmlTable(output, align = 'l', align.header = 'l', rnames = FALSE , css.cell = ' padding-bottom: 5px; vertical-align:top; padding-right: 10px;min-width: 5em ', caption="*tcplfit2* model details.", tfoot = "Model descriptions are pulled from tcplFit2 manual at ." )