Merge pull request #325 from USEPA/311-update-tcplPlot-for-loec

Refactored tcplPlot and auxiliary functions to support any number of curves on a comparison plot

Author: cthunes

Diff for: DESCRIPTION

@@ -42,7 +42,9 @@ Imports:
     gridExtra,
     stringr,
     rlang,
-    ctxR (>= 1.1.0)
+    ctxR (>= 1.1.0),
+    viridis,
+    gt
 Suggests:
     roxygen2,
     knitr,
@@ -51,16 +53,16 @@ Suggests:
     htmlTable,
     testthat (>= 3.0.0),
     reshape2,
-    viridis,
     kableExtra,
     colorspace,
     magrittr,
     vdiffr,
     httptest,
-    rmdformats
+    rmdformats,
+    gtable
 License: MIT + file LICENSE
 LazyData: true
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 VignetteBuilder: knitr
 Encoding: UTF-8
 Config/testthat/edition: 3

Diff for: NAMESPACE

@@ -61,6 +61,7 @@ export(tcpldbStats)
 import(DBI)
 import(data.table)
 import(gridExtra)
+import(gt)
 import(stringr)
 importFrom(RColorBrewer,brewer.pal)
 importFrom(RMariaDB,MariaDB)
@@ -85,10 +86,12 @@ importFrom(dplyr,left_join)
 importFrom(dplyr,matches)
 importFrom(dplyr,mutate)
 importFrom(dplyr,mutate_if)
+importFrom(dplyr,n)
 importFrom(dplyr,pull)
 importFrom(dplyr,rowwise)
 importFrom(dplyr,select)
 importFrom(dplyr,summarise)
+importFrom(dplyr,summarize)
 importFrom(dplyr,tibble)
 importFrom(dplyr,ungroup)
 importFrom(ggplot2,aes)
@@ -103,6 +106,7 @@ importFrom(ggplot2,geom_vline)
 importFrom(ggplot2,ggplot)
 importFrom(ggplot2,ggsave)
 importFrom(ggplot2,guide_legend)
+importFrom(ggplot2,is.ggplot)
 importFrom(ggplot2,labs)
 importFrom(ggplot2,margin)
 importFrom(ggplot2,scale_color_manual)
@@ -157,6 +161,8 @@ importFrom(stats,median)
 importFrom(stats,optim)
 importFrom(stats,quantile)
 importFrom(stats,setNames)
+importFrom(stringr,str_count)
+importFrom(stringr,str_trunc)
 importFrom(tcplfit2,tcplfit2_core)
 importFrom(tcplfit2,tcplhit2_core)
 importFrom(tidyr,pivot_longer)
@@ -169,3 +175,4 @@ importFrom(utils,read.csv)
 importFrom(utils,read.table)
 importFrom(utils,tail)
 importFrom(utils,write.table)
+importFrom(viridis,viridis)

Diff for: R/mc5.R

@@ -33,7 +33,7 @@ mc5 <- function(ae, wr = FALSE) {
   modl_prob <- cnst_prob <- hill_er <- modl_tp <- modl_ga <- modl_gw <- NULL
   hill_rmse <- hill_prob <- modl_acb <- modl_acc <- gnls_er <- modl_la <- NULL
   gnls_la <- modl_lw <- gnls_lw <- gnls_rmse <- gnls_prob <- actp <- NULL
-  modl_ac10 <- model_type <- NULL 
+  modl_ac10 <- model_type <- conc <- loec_hitc <- NULL 
 
   overwrite_osd <- FALSE
 

Diff for: R/tcplLoadData.R

@@ -218,14 +218,11 @@ tcplLoadData <- function(lvl, fld = NULL, val = NULL, type = "mc", add.fld = TRU
 
     if (length(colnames(dat))) {
       if (lvl == 3 | lvl == "agg") {
-        dat$resp <- lapply(dat$resp, unlist)
-        dat$logc <- lapply(dat$logc, unlist)
         if (lvl == 3) dat <- unnest_longer(dat, c(conc, logc, resp)) %>% as.data.table() 
         else dat <- unnest_longer(dat, c(logc, resp)) %>% as.data.table() 
       }
 
       if (lvl == 6) {
-        dat$flag <- lapply(dat$flag, unlist)
         dat <- unnest_longer(dat, flag) %>% filter(flag != "NULL") %>% as.data.table()
       }
     }

Diff for: R/tcplPlot.R

@@ -59,10 +59,14 @@
 #' @seealso \code{\link{tcplPlot}}
 #'
 #' @import data.table
+#' @importFrom dplyr mutate group_by summarize rowwise
+#' @importFrom stringr str_count
 #' @export
 tcplPlotLoadData <- function(type = "mc", fld = "m4id", val, flags = FALSE){
-  #variable binding
-  lvl <- m4id <- conc <- resp <- conc_unit <- NULL
+  
+  # variable binding for R CMD check
+  lvl <- m4id <- s2id <- conc <- resp <- conc_unit <- flag_count <- NULL
+  bmd <- top <- bmr <- model_type <- coff <- max_med <- hitc <- NULL
 
   # Validate vars based on some assumed properties
   validated_vars <- tcplPlotValidate(type = type,flags = flags)
@@ -73,13 +77,6 @@ tcplPlotLoadData <- function(type = "mc", fld = "m4id", val, flags = FALSE){
   dat <- tcplLoadData(lvl = lvl, fld = fld, val = val, type = type)
   if (nrow(dat) == 0) stop("No data for fld/val provided")
 
-  # set order to given order
-  dat <- dat[order(match(get(fld[1]), if(is.list(val)) val[[1]] else val))]
-  if (getOption("TCPL_DRVR") == "API" && tolower(fld) == "aeid") {
-    dat <- dat %>% arrange(m4id)
-  }
-  dat$order <- 1:nrow(dat)
-  
   mcLoadDat <- function(m4id = NULL,flags) {
     l4 <- tcplLoadData(lvl = 4, fld = "m4id", val = m4id, add.fld = T)
     dat <- l4[dat, on = "m4id"]
@@ -103,8 +100,10 @@ tcplPlotLoadData <- function(type = "mc", fld = "m4id", val, flags = FALSE){
   if (getOption("TCPL_DRVR") != "API") {
     if (type == "mc") {
       dat <- mcLoadDat(dat$m4id,flags = flags)
+      setorder(dat, m4id)
       agg <- tcplLoadData(lvl = "agg", fld = "m4id", val = dat$m4id)
     } else { # type == 'sc'
+      setorder(dat, s2id)
       agg <- tcplLoadData(lvl = "agg", fld = "s2id", val = dat$s2id, type = "sc")
     }
 
@@ -113,11 +112,15 @@ tcplPlotLoadData <- function(type = "mc", fld = "m4id", val, flags = FALSE){
 
     #determine if we're single conc or multiconc based on dat
     join_condition <- c("m4id","s2id")[c("m4id","s2id") %in% colnames(dat)]
-    conc_resp_table <- agg %>% group_by(.data[[join_condition]]) %>% summarise(conc = list(conc), resp = list(resp)) %>% as.data.table()
+    conc_resp_table <- agg %>% group_by(.data[[join_condition]]) %>% summarize(conc = list(conc), resp = list(resp)) %>% as.data.table()
     dat <- dat[conc_resp_table, on = join_condition]
 
+    # get chemical and sample information
     dat <- tcplPrepOtpt(dat)
 
+    # determine missing chem info and replace with string description of well type(s)
+    dat <- tcplPlotLoadWllt(dat, type)
+    
   } else {
     # fix flags from API for plotting
     if (flags == TRUE) {
@@ -129,6 +132,12 @@ tcplPlotLoadData <- function(type = "mc", fld = "m4id", val, flags = FALSE){
     dat$conc_unit <- dat$tested_conc_unit
   }
 
+  # add flag_count
+  if (flags == TRUE) {
+    dat$flag_count <- 0
+    dat[flag != "None", flag_count := stringr::str_count(flag, "\n") + 1]
+  }
+  
   # add normalized data type for y axis
   ndt <- tcplLoadAeid(fld = "aeid", val = dat$aeid, add.fld = "normalized_data_type")
   dat <- dat[ndt, on = "aeid"]
@@ -138,5 +147,31 @@ tcplPlotLoadData <- function(type = "mc", fld = "m4id", val, flags = FALSE){
   dat <- dat[conc_unit=="uM", conc_unit:="\u03BCM"]
   dat <- dat[conc_unit=="mg/l", conc_unit:="mg/L"]
 
+  # remove zero width space if exists
+  dat$chnm <- gsub("\u200b", "", dat$chnm)
+  
+  # replace null bmd in dat table
+  dat <- dat[is.null(dat$bmd), bmd:=NA]
+  
+  #replace null top with 0
+  dat[is.null(dat$top), top := 0]
+  dat[is.na(top), top := 0]
+  
+  # correct bmr and coff direction
+  if (type == "mc") {
+    # if top if less than 0, flip bmr no matter what
+    dat[top < 0, bmr := bmr * -1]
+    # if model type is loss, flip cut off
+    dat[model_type == 4, coff := coff * -1]
+    # if model type is bidirectional, flip cut off if top is less than 0
+    dat[model_type == 2 & top < 0, coff := coff * -1]
+  } else { # sc
+    # if max median is less than 0, flip cut off to align with it
+    dat[max_med < 0, coff := coff * -1]
+    # if hitc is less than 0, max median is in the opposite of intended direction, 
+    # so flip cut off (possibly again)
+    dat[hitc < 0, coff := coff * -1]
+  }
+  
   dat
 }

Diff for: R/tcplPlotLoadData.R

@@ -1,3 +1,46 @@
+#' tcplPlotLoadWllt
+#' Replaces NA dtxsid and chnm with a string description of the sample's well type(s)
+#'
+#' @param dat dataset
+#' @param type mc or sc
+#'
+#' @return dat with updated dtxsid/chnm if they are NA
+#' @import data.table
+#' @importFrom dplyr group_by summarize select left_join case_when 
+tcplPlotLoadWllt <- function(dat = NULL, type = "mc") {
+  
+  # variable binding for R CMD check
+  dsstox_substance_id <- chnm <- aeid <- acid <- wllt <- wllt_desc <- NULL
+  
+  # determine missing chemical info
+  missing <- dat[is.na(dsstox_substance_id) | is.na(chnm), .(spid, aeid)]
+  if (nrow(missing) > 0) {
+    acid_map <- tcplLoadAcid(fld = "aeid", val = missing$aeid)[,.(aeid,acid)]
+    missing <- missing[acid_map, on="aeid"]
+    l0_dat <- tcplLoadData(type = type, lvl = 0, fld = list("spid","acid"), 
+                           list(missing$spid, missing$acid))
+    
+    # replace with string describing the well type(s)
+    wllts <- l0_dat |> select(acid, spid, wllt) |> group_by(acid, spid) |> 
+      summarize(wllt_desc = dplyr::case_when(
+        all(unique(wllt) %in% c('c', 'p')) ~ "Gain-of-signal control",
+        all(unique(wllt) %in% c('m', 'o')) ~ "Loss-of-signal control",
+        all(unique(wllt) == 'n') ~ "Neutral/negative control",
+        all(unique(wllt) == 'b') ~ "Blank",
+        all(unique(wllt) == 'v') ~ "Viability control",
+        all(!is.na(unique(wllt))) && all(!is.null(unique(wllt))) ~ paste0("Well type: ", paste(unique(wllt), collapse = ", ")),
+        .default = NA
+      ), .groups="drop") |> left_join(acid_map, by = "acid") |> as.data.table()
+    no_wllt <- wllts[is.na(wllt_desc), spid]
+    if (length(no_wllt)) warning(paste0("wllt for SPID(s): ", paste(no_wllt, collapse = ", "), " missing. Leaving dsstox_substance_id and chnm as NA."))
+    dat <- left_join(dat, wllts, by = c("aeid", "spid"))
+    dat[is.na(dsstox_substance_id) | is.na(chnm), c("dsstox_substance_id", "chnm") := wllt_desc]
+  }
+  dat
+  
+}
+
+
 #' tcplPlotSetYRange
 #'
 #' @param dat dataset
@@ -47,17 +90,69 @@ tcplPlotSetYRange <- function(dat, yuniform, yrange, type) {
   yrange
 }
 
+#' tcplPlotCalcAspectRatio
+#'
+#' @param type string of mc or sc indicating if it is single or multi conc
+#' @param verbose should the plot return a table with parameters
+#' @param multi Boolean, by default TRUE for "pdf" if the number of plots exceeds
+#' one. Prints variable number of plots per page depending on 'verbose' and 'type' 
+#' settings.
+#' @param nrows Integer, number of rows each compare plot uses
+#' @param output How should the plot be presented. To work with the plot in 
+#' environment, use "ggplot"; to interact with the plot in application, use 
+#' "console"; or to save as a file type, use "pdf", "jpg", "png", "svg", or "tiff".
+#' @param group.threshold Integer of length 1, minimum number of samples in a 
+#' given plot where comparison plots, instead of coloring models by sample, should
+#' delineate curve color by a given group.fld. By default 9.
+#' @param nrow Integer, number of rows in multiplot. By default 2.
+#' @param ncol Integer, number of columns in multiplot. By default 3, 2 if verbose, 
+#' 1 for compare plots.
+#' @param flags Boolean, by default FALSE. If TRUE, level 6 flags are displayed
+#' within output.
+#'
+#' @return a list of validated parameters for plotting
+tcplPlotCalcAspectRatio <- function(type = "mc", verbose = FALSE, multi = FALSE, nrows = 0, 
+                                    output = c("ggplot", "console", "pdf", "png", "jpg", "svg", "tiff"), 
+                                    group.threshold = 9, nrow = NULL, ncol = NULL, flags = FALSE) {
+  
+  # assign nrow = ncol = 1 for output="pdf" and multi=FALSE to plot one plot per page
+  if(output == "pdf" && multi == FALSE)
+    nrow <- ncol <- 1
+  w <- 7
+  h <- 5
+  if (verbose && any(nrows > 1 & nrows < group.threshold)) 
+    w <- 12
+  if (type == "sc" && any(!verbose, all(nrows == 1 | nrows >= group.threshold))) 
+    w <- 5
+  if (type == "sc")
+    h <- 6
+  if (verbose && any(nrows > 1 & nrows < group.threshold)) 
+    h <- max(2 + max(nrows[nrows > 1 & nrows < group.threshold]), 6)
+  if (any(nrows >= group.threshold) && flags)
+    h <- max(h, 8)
+  if(is.null(nrow))
+    nrow <- max(round(10/h), 1)
+  if(is.null(ncol))
+    ncol <- max(round(14/w), 1)
+  list(w=w,h=h,nrow=nrow,ncol=ncol)
+}
 
 #' tcplPlotValidate
 #'
+#' @param dat data.table containing plot-prepared data
 #' @param type string of mc or sc indicating if it is single or multi conc
+#' @param compare Character vector, the field(s) to join samples on to create comparison
+#' plots
+#' @param by Parameter to divide files into e.g. "aeid".
 #' @param flags bool - should we return flags
 #' @param output how should the plot be formatted
 #' @param multi are there multiple plots
 #' @param verbose should the plot return a table with parameters
 #'
 #' @return a list of validated parameters for plotting
-tcplPlotValidate <- function(type = "mc", flags = NULL, output = "none", multi = NULL, verbose = FALSE) {
+tcplPlotValidate <- function(dat = NULL, type = "mc", compare = "m4id", by = NULL, 
+                             flags = NULL, output = c("ggplot", "console", "pdf", "png", "jpg", "svg", "tiff"), 
+                             multi = NULL, verbose = FALSE) {
   # set lvl based on type
   lvl <- 5
   if (type == "sc") {
@@ -66,28 +161,39 @@ tcplPlotValidate <- function(type = "mc", flags = NULL, output = "none", multi =
       warning("'flags' was set to TRUE - no flags exist for plotting single concentration")
       flags <- FALSE
     }
+    if (compare == "m4id") compare <- "s2id"
   }
-
-  # default assign multi=TRUE for output="pdf"
-  if (output == "pdf" && is.null(multi)) {
-    multi <- TRUE
+  
+  if (!is.null(by) && length(by) > 1) stop("'by' must be of length 1.")
+  if (length(output) > 1) output <- output[1] else if (length(output) == 0) stop("'output' cannot be NULL")
+  
+  if (!is.null(dat) && !is.data.table(dat)) {
+    if (!is.list(dat) || !is.data.table(dat[[1]])) {
+      stop("'dat' must be a data.table or a list of data.tables.")
+    }
+    if (!compare %in% c("m4id", "s2id")) {
+      warning("'dat' provided as list of list of data tables, meaning compare plots are already subdivided. 'compare' field will be ignored and have no effect.")
+    }
+    if (!is.null(by)) {
+      warning("Using 'by' can have unintended consequences when 'dat' is provided as a list of data.tables. Instead, consider adding a custom field to group comparison plots, and specify using the 'compare' parameter. Then, use 'by' to split plots into files.")
+    }
   }
+  
   # forced assign multi=FALSE for output = c("console","png","jpg","svg","tiff"), verbose=FALSE for output="console"
   if (output != "pdf") {
     multi <- FALSE
-    if (output == "console") {
-      verbose <- FALSE
-    }
+    if (output == "console") verbose <- FALSE
   }
 
-  list(lvl = lvl, type = type, flags = flags, output = output, multi = multi, verbose = verbose)
+  list(lvl = lvl, compare = compare, flags = flags, output = output, multi = multi, verbose = verbose)
 }
 
 
 #' tcplLegacyPlot
 #'
 #' @return a ggplot based on old plotting methodology
 tcplLegacyPlot <- function() {
+  stop("Legacy plotting is currently unsupported.")
   # VARIABLE BINDING
   fld <- val <- lvl <- multi <- fileprefix <- NULL
   if (length(output) > 1) output <- output[1]

Diff for: R/tcplPlotUtils.R

@@ -44,8 +44,14 @@ tcplQueryAPI <- function(resource = "data", fld = NULL, val = NULL, return_flds
 
     dat$dsstox_substance_id <- dat$dtxsid
 
-    # unlist logc to conc
-    dat <- dat %>% rowwise() %>% mutate(conc = list(10^unlist(logc))) %>% as.data.table()
+    # unlist logc, resp, flag, mc6mthdid
+    unlist_cols <- c("logc", "resp", "flag", "mc6MthdId")
+    for (c in unlist_cols) {
+      dat[[c]] <- lapply(dat[[c]], unlist)
+    }
+    
+    # unlog logc to conc
+    dat <- dat %>% rowwise() %>% mutate(conc = list(10^logc)) %>% as.data.table()
 
   } else if (resource == "assay") {
 

Diff for: R/tcplQueryAPI.R

@@ -58,7 +58,7 @@
 tcplSubsetChid <- function(dat, flag = TRUE, type = "mc", export_ready = FALSE) {
   ## Variable-binding to pass R CMD Check
   chit <- hitc <- aeid <- casn <- fitc <- fitc.ordr <- m4id <- nflg <- NULL
-  chid <- conc <- minc <- NULL
+  chid <- conc <- minc <- actc <- NULL
 
   if (!type %in% c("mc", "sc")) {
     stop("type must be sc (single concentration) or mc (multi-concentration)")