R/checkSS.R

In shirewoman2/Consultancy: Standardized tools for using R within the Simcyp Consultancy Team

#' Create a graph of simulated concentrations to check for whether a perpetrator
#' drug is at steady-state when the victim drug is dosed
#'
#' \code{checkSS} creates a graph with time on the x axis and then some useful
#' concentration point -- C0, Cmax, Cmin (see notes on usage), or Clast for each
#' dose -- on the y axis. You can optionally overlay this graph with the full
#' concentration-time data for that same compound or a different one. For
#' example, you could create a graph with points for Clast for an inhibitor and
#' then lines showing the substrate concentration-time data overlaid on top of
#' that. This way, you can check whether the substrate is dosed after the
#' inhibitor reaches steady state, and you could also check whether the
#' inhibitor was present the entire time the substrate was getting eliminated.
#' Additionally, if you supply a data.frame of enzyme-abundance data, that will
#' be graphed above the concentration-time data, giving another check for
#' whether everything important is at steady state. 
#'
#' @param t0 start time for compound being plotted
#' @param conc_point Concentration point to plot. Options are: \describe{
#'
#'   \item{"C0"}{concentration for the first time point available for that dose}
#'
#'   \item{"Clast"}{(default) concentration for the last time point available
#'   for that dose}
#'
#'   \item{"Cmin"}{the minimum concentration for that dosing interval, which may
#'   not be the last time point; this is somewhat prone to sampling artifacts
#'   and we generally recommend using "Clast" instead}
#'
#'   \item{"Cmax"}{the maximum concentration for that dose number}}
#' @param mark_dosing_substrate optionally mark substrate dosing intervals on
#'   the graph as "none" (default) to have no marks for the dosing intervals or
#'   a combination of a color in R and a named linetype, e.g., "red dotted" or
#'   "blue dashed" or even "#FFBE33 longdash".
#' @param mark_dosing_inhibitor1 optionally mark inhibitor 1 dosing intervals on
#'   the graph as "none" (default) to have no marks for the dosing intervals or
#'   a combination of a color in R and a named linetype, e.g., "red dotted" or
#'   "blue dashed" or even "#FFBE33 longdash".
#' @param mark_dosing_inhibitor2 optionally mark inhibitor 2 dosing intervals on
#'   the graph as "none" (default) to have no marks for the dosing intervals or
#'   a combination of a color in R and a named linetype, e.g., "red dotted" or
#'   "blue dashed" or even "#FFBE33 longdash".
#' @param diff_cutoff what percent difference cutoff would you like to use to
#'   color the points? The default is for points with less than a 5\% difference
#'   from the previous point to be blue and points with a larger percent
#'   difference to be red.
#' @param ct_dataframe the input concentration-time data generated by running
#'   the function \code{\link{extractConcTime}} or, if you'd also like to see an
#'   overlay of the substrate with or without any inhibitors, by running
#'   \code{\link{extractConcTime_mult}}. This function assumes there's data for
#'   only one tissue and that it's appropriate to compare all the data included,
#'   so you'll get graph artifacts if those are not the case.
#' @param mean_type the mean type to use since this function only displays
#'   summary data; defaults to the arithmetic mean
#' @param accum_compoundID the compound ID to monitor for accumulation. Defaults
#'   to "inhibitor 1". The time point requested will be shown as points and will
#'   be colored by percent difference from the previous point.
#' @param overlay_compoundID the compound ID to overlay the complete
#'   concentration-time data for. Defaults to "substrate"; use "none" for no
#'   overlaid plot. This was designed for the following scenario: Monitor
#'   Ctrough for an inhibitor to make sure that it's at steady state, plot that
#'   as points, and then plot the concentration-time profile of a substrate over
#'   that to make sure that the inhibitor is present the whole time the
#'   substrate is being eliminated.
#' @param sim_enz_dataframe optionally include a data.frame of enzyme abundances
#'   to plot below the main accumulation plot. If provided, all enzymes and
#'   tissues in the data will be included in the plot.
#' @param x_axis_interval optionally set the x-axis major tick-mark interval.
#'   Acceptable input: any number or leave as NA to accept default values, which
#'   are generally reasonable guesses as to aesthetically pleasing and
#'   PK-relevant intervals.
#' @param graph_title optionally specify a title that will be centered across
#'   your graph or set of graphs
#' @param graph_title_size the font size for the graph title if it's included;
#'   default is 14. This also determines the font size of the graph labels.
#' @param save_graph optionally save the output graph by supplying a file name
#'   in quotes here, e.g., "My conc time graph.png". If you do not designate a
#'   file extension, it will be saved as a png file, but if you specify a
#'   different file extension, it will be saved as that file format. Acceptable
#'   extensions are "eps", "ps", "jpeg", "jpg", "tiff", "png", "bmp", or "svg".
#'   Do not include any slashes, dollar signs, or periods in the file name.
#'   Leaving this as NA means the file will not be automatically saved to disk.
#' @param fig_height figure height in inches; default is 6
#' @param fig_width figure width in inches; default is 5
#'
#' @return a ggplot2 graph
#' @export
#'
#' @examples
#'
#' checkSS(ct_dataframe = MDZ_Keto)
#'
#' checkSS(ct_dataframe = MDZ_Keto, conc_point = "Cmax")
#'
#' checkSS(ct_dataframe = MDZ_Keto, conc_point = "Cmax",
#'    accum_compoundID = "inhibitor 1", overlay_compoundID = "none",
#'    mark_dosing_substrate = "pink dotted", diff_cutoff = 0.01,
#'    save_graph = "MDZ keto accumulation check.png")
#'
#' 
checkSS <- function(ct_dataframe,
                    accum_compoundID = "inhibitor 1", 
                    overlay_compoundID = "substrate",
                    sim_enz_dataframe = NA,
                    t0 = 0,
                    conc_point = "Clast",
                    mark_dosing_substrate = "none",
                    mark_dosing_inhibitor1 = "none", 
                    mark_dosing_inhibitor2 = "none", 
                    diff_cutoff = 0.05,
                    mean_type = "arithmetic", 
                    x_axis_interval = NA,
                    graph_title = NA,
                    graph_title_size = 14, 
                    save_graph = NA,
                    fig_height = 4,
                    fig_width = 8){
   
   # Error catching ----------------------------------------------------------
   
   # Check whether tidyverse is loaded
   if("package:tidyverse" %in% search() == FALSE){
      stop("The SimcypConsultancy R package requires the package tidyverse to be loaded, and it doesn't appear to be loaded yet. Please run\nlibrary(tidyverse)\n    ...and then try again.", 
           call. = FALSE)
   }
   
   if(nrow(ct_dataframe) == 0){
      stop("Please check your input. The data.frame you supplied for ct_dataframe doesn't have any rows.", 
           call. = FALSE)
   }
   
   if(length(sort(unique(ct_dataframe$File))) > 1){
      stop(paste0("The checkSS function is for graphing only one simulator file at a time, but you have ",
                  length(sort(unique(ct_dataframe$File))), 
                  " simulator files. An example of how to rectify this without re-extracting any data, where `CT` is the data.frame of concentration-time data you've already got:
checkSS(ct_dataframe = CT %>% filter(File == `mysim-01.xlsx`))"),
call. = FALSE)
   }
   
   
   # Main body of function --------------------------------------------------
   
   MyMeanType <- ct_dataframe %>%
      filter(CompoundID == accum_compoundID & 
                Trial %in% c("geomean", "mean", "median")) %>% 
      pull(Trial) %>% unique() %>% 
      factor(levels = c("mean", "geomean", "median")) %>% 
      sort()
   
   if(switch(mean_type, "arithmetic" = "mean", "geometric" = "geomean",
             "median" = "median") %in% ct_dataframe$Trial == FALSE){
      
      warning(wrapn(paste0("You requested the ", 
                           switch(mean_type, "arithmetic" = "arithmetic means",
                                  "geometric" = "geometric means", 
                                  "median" = "medians"), 
                           ", but those are not included in your data. Instead, the ",
                           ifelse(MyMeanType[1] == "mean", 
                                  "arithmetic mean", MyMeanType[1]),
                           "s will be used.")), 
              call. = FALSE)
      MyMeanType <- MyMeanType[1] %>% as.character()
      
   } else {
      MyMeanType <- switch(mean_type, "arithmetic" = "mean", "geometric" = "geomean",
                           "median" = "median")
   }
   
   # Adjusting conc units to ng/mL. At some point, we could make this detect
   # the units and then adjust the graph titles to whatever is in the data, but
   # for now, just making everything be ng/mL for convenience.
   ct_dataframe <- convert_units(ct_dataframe, conc_units = "ng/mL")
   
   # suppressMessages(
   SScheck <- ct_dataframe %>% 
      filter(Trial == MyMeanType & 
                DoseNum > 0 & 
                CompoundID == accum_compoundID) %>%  
      group_by(CompoundID, DoseNum, Inhibitor) %>% 
      # switch doesn't seem to work with summarize. Calculating each value.
      summarize(t0 = min(Time),
                tlast = max(Time),
                tmin = Time[which.min(Conc)],
                tmax = Time[which.max(Conc)], 
                Cmin = min(Conc), 
                Cmax = max(Conc), 
                C0 = Conc[which.min(Time)], 
                Clast = Conc[which.max(Time)]) %>% 
      ungroup() %>% 
      mutate(Conc = switch(conc_point, 
                           "Cmin" = Cmin,
                           "Cmax" = Cmax, 
                           "C0" = C0, 
                           "Clast" = Clast),
             Time = switch(conc_point, 
                           "Cmin" = tmin, 
                           "Cmax" = tmax,
                           "C0" = t0, 
                           "Clast" = tlast))
   
   SScheck <- split(SScheck, f = list(SScheck$CompoundID, 
                                      SScheck$Inhibitor))
   
   for(i in 1:length(SScheck)){
      SScheck[[i]] <- SScheck[[i]] %>% 
         mutate(PercDiff = c(NA, diff(Conc, lag = 1))/Conc, 
                DiffCriterion = abs(PercDiff) < diff_cutoff,
                DiffCriterion = ifelse(is.na(DiffCriterion), FALSE, DiffCriterion), 
                DiffCriterion_lab = ifelse(DiffCriterion, 
                                           paste0("<", diff_cutoff*100, "%"),
                                           paste0("\u2265", diff_cutoff*100, "%")), 
                DiffCriterion_lab = factor(DiffCriterion_lab, 
                                           levels = c(paste0("\u2265", diff_cutoff*100, "%"),
                                                      paste0("<", diff_cutoff*100, "%"))))
   }
   
   SScheck <- bind_rows(SScheck) %>% 
      mutate(Perp = case_when(Inhibitor != "none" ~ "DDI", 
                              CompoundID == "substrate" & Inhibitor == "none" ~ "baseline"))
   
   # Noting whether perpetrator present and what it is
   MyPerpetrator <- unique(ct_dataframe$Inhibitor[ct_dataframe$Inhibitor != "none"])
   
   if(length(MyPerpetrator) > 0){
      ct_dataframe <- ct_dataframe %>%
         mutate(Inhibitor = factor(Inhibitor, levels = c("none", MyPerpetrator)))
   }
   
   # Setting up graphs
   LineAES_substrate <- str_split(mark_dosing_substrate, pattern = " ")[[1]]
   LineAES_inhibitor1 <- str_split(mark_dosing_inhibitor1, pattern = " ")[[1]]
   LineAES_inhibitor2 <- str_split(mark_dosing_inhibitor2, pattern = " ")[[1]]
   
   G <- ggplot(SScheck %>% filter(CompoundID == accum_compoundID),
               aes(x = Time, y = Conc, color = DiffCriterion_lab, 
                   shape = Perp))
   
   if(mark_dosing_substrate != "none"){
      G <- G + 
         geom_vline(xintercept = SScheck$t0[SScheck$CompoundID == "substrate" &
                                               SScheck$DoseNum != 0], 
                    color = LineAES_substrate[1], linetype = LineAES_substrate[2])
   }
   
   if(mark_dosing_inhibitor1 != "none"){
      G <- G + 
         geom_vline(xintercept = SScheck$t0[SScheck$CompoundID == "inhibitor 1" &
                                               SScheck$DoseNum != 0], 
                    color = LineAES_inhibitor1[1], linetype = LineAES_inhibitor1[2])
   }
   
   if(mark_dosing_inhibitor2 != "none"){
      G <- G + 
         geom_vline(xintercept = SScheck$t0[SScheck$CompoundID == "inhibitor 2" &
                                               SScheck$DoseNum != 0], 
                    color = LineAES_inhibitor2[1], linetype = LineAES_inhibitor2[2])
   }
   
   G <- G +
      geom_point(size = 2) + 
      labs(color = paste(str_to_title(accum_compoundID),
                         "difference\nfrom previous point"), 
           linetype = paste(str_to_title(overlay_compoundID), "concentration"), 
           shape = NULL) +
      xlab("Time (h)") +
      scale_x_time(x_axis_interval = x_axis_interval, 
                   time_range = c(0, 
                                  max(ct_dataframe$Time[
                                     ct_dataframe$CompoundID %in% c(accum_compoundID, 
                                                                    overlay_compoundID)]))) +
      scale_color_brewer(palette = "Set1") +
      theme_consultancy()
   
   if(length(unique(SScheck$Perp[
      SScheck$CompoundID == accum_compoundID])) == 1){
      G <- G + guides(shape = "none")
   }
   
   if(overlay_compoundID != "none"){
      
      OverlayCT <- ct_dataframe %>% 
         filter(Trial == MyMeanType & CompoundID == overlay_compoundID) %>% 
         mutate(Inhibitor = ifelse(Inhibitor == "none", 
                                   "without perpetrator",
                                   paste("with", prettify_compound_name(MyPerpetrator))), 
                Inhibitor = factor(Inhibitor, 
                                   levels = c("without perpetrator",
                                              paste("with", 
                                                    prettify_compound_name(MyPerpetrator)))))
      
      # Checking for differences in scale between accum_compoundID and
      # overlay_compoundID
      ScaleChallenges <- abs(max(SScheck$Conc, na.rm = T) - 
                                max(OverlayCT$Conc, na.rm = T))
      
      if(complete.cases(ScaleChallenges) && ScaleChallenges > 100){
         ScaleFactor <- round_up(max(SScheck$Conc, na.rm = T) /
                                    max(OverlayCT$Conc, na.rm = T))
         OverlayCT <- OverlayCT %>% 
            mutate(Conc = Conc * ScaleFactor)
      } 
      
      if(length(MyPerpetrator) > 0){
         G <- G + 
            geom_line(data = OverlayCT, 
                      aes(x = Time, y = Conc, linetype = Inhibitor), 
                      inherit.aes = F, color = "gray50") +
            scale_linetype_manual(values = c("solid", "dashed"))
      } else {
         
         G <- G + 
            geom_line(data = OverlayCT, 
                      aes(x = Time, y = Conc),
                      inherit.aes = F, color = "gray50")
      }
      
      # Adding 2nd y axis if there were big differences in scale.
      if(complete.cases(ScaleChallenges) && ScaleChallenges> 100){
         G <- G +
            scale_y_continuous(
               paste(str_to_title(accum_compoundID), 
                     switch(conc_point, 
                            "Cmin" = "Cmin (ng/mL)", 
                            "Cmax" = "Cmax (ng/mL)", 
                            "C0" = "C0 (ng/mL)",
                            "Clast" = "Clast (ng/mL)")), 
               # FIXME - I've got something screwed up here. I don't remember
               # why I set this up to have a label of "ng/mL" here w/out noting
               # what the concentration was adjusted by, but the numbers on the
               # 2ndary axis are complete fiction.
               sec.axis = sec_axis(~ . * ScaleFactor,
                                   name = paste(str_to_title(overlay_compoundID),
                                                "concentration (ng/mL)"))) +
            theme(axis.line.y.right = element_line(color = "black"))
      } else {
         
         G <- G +
            scale_y_continuous(limits = c(0, max(c(OverlayCT$Conc, 
                                                   SScheck$Conc)))) +
            ylab(switch(conc_point, 
                        "Cmin" = expression(C[min]~"(ng/mL)"), 
                        "Cmax" = expression(C[max]~"(ng/mL)"), 
                        "C0" = expression(C[0]~"(ng/mL)"),
                        "Clast" = expression(C[last]~"(ng/mL)")))
      }
      
      # Adding legend item saying that points are for accum_compoundID and
      # lines are for overlay_compoundID. 
      Empty <- data.frame(Time = mean(OverlayCT$Time),
                          Conc = as.numeric(NA), 
                          Type = c(accum_compoundID, overlay_compoundID))
      G <- G +
         geom_point(data = Empty, aes(x = Time, y = Conc, alpha = Type),
                    inherit.aes = F) +
         geom_line(data = Empty, aes(x = Time, y = Conc, alpha = Type),
                   inherit.aes = F)  +
         scale_alpha_manual(
            name = "Compound", values = c(1,1),
            breaks = c(accum_compoundID, overlay_compoundID))
      
   } else {
      G <- G + 
         scale_y_continuous(limits = c(0, max(SScheck$Conc[
            SScheck$CompoundID == accum_compoundID]))) +
         ylab(switch(conc_point, 
                     "Cmin" = expression(C[min]~"(ng/mL)"), 
                     "Cmax" = expression(C[max]~"(ng/mL)"), 
                     "C0" = expression(C[0]~"(ng/mL)"),
                     "Clast" = expression(C[last]~"(ng/mL)")))
      
   }
   
   # Setting order of legend items 
   G <- G + 
      guides(alpha = guide_legend(order = 1, 
                                  override.aes =
                                     list(linetype = c("blank", "solid"),
                                          shape = c(16, NA),
                                          size = c(2, 0.5),
                                          color = c("#7030A0", "gray50"))), 
             color = guide_legend(order = 2),
             linetype = guide_legend(order = 3)) +
      theme(axis.title = element_text(face = "plain"))
   
   if(complete.cases(graph_title) & 
      "data.frame" %in% class(sim_enz_dataframe) == FALSE){
      G <- G + ggtitle(graph_title) +
         theme(plot.title = element_text(hjust = 0.5, size = graph_title_size), 
               plot.title.position = "panel")
   }
   
   
   if("data.frame" %in% class(sim_enz_dataframe)){
      sim_enz_dataframe <- sim_enz_dataframe %>% filter(File == unique(ct_dataframe$File))
      
      suppressWarnings(suppressMessages(
         EnzPlot <- enz_plot_overlay(
            sim_enz_dataframe = sim_enz_dataframe %>% filter(Inhibitor != "none"), 
            colorBy_column = Enzyme, 
            linetype_column = Tissue, 
            y_axis_label = ifelse(length(unique(sim_enz_dataframe$Enzyme)) == 1, 
                                  paste("Relative", unique(sim_enz_dataframe$Enzyme), 
                                        "abundance"), 
                                  "Relative enzyme abundance"))
      ))
      
      # Removing box around graph and making axis titles not be bold for
      # consistency w/G
      EnzPlot <- EnzPlot +
         theme_consultancy(border = FALSE) +
         theme(axis.title = element_text(face = "plain"))
      
      if(mark_dosing_substrate != "none"){
         EnzPlot <- EnzPlot +
            geom_vline(xintercept = SScheck$t0[SScheck$CompoundID == "substrate" &
                                                  SScheck$DoseNum != 0], 
                       color = LineAES_substrate[1], linetype = LineAES_substrate[2])
      }
      
      if(mark_dosing_inhibitor1 != "none"){
         EnzPlot <- EnzPlot +
            geom_vline(xintercept = SScheck$t0[SScheck$CompoundID == "inhibitor 1" &
                                                  SScheck$DoseNum != 0], 
                       color = LineAES_inhibitor1[1], linetype = LineAES_inhibitor1[2])
      }
      
      if(mark_dosing_inhibitor2 != "none"){
         EnzPlot <- EnzPlot +
            geom_vline(xintercept = SScheck$t0[SScheck$CompoundID == "inhibitor 2" &
                                                  SScheck$DoseNum != 0], 
                       color = LineAES_inhibitor2[1], linetype = LineAES_inhibitor2[2])
      }
      
      G <- ggpubr::ggarrange(EnzPlot, G, ncol = 1, align = "hv")
      
      if(complete.cases(graph_title)){
         G <- ggpubr::annotate_figure(G, 
                                      top = ggpubr::text_grob(graph_title, 
                                                              size = graph_title_size))
      }
      
   }
   
   
   # Saving ------------------------------------------------------------------
   
   if(complete.cases(save_graph)){
      FileName <- save_graph
      if(str_detect(FileName, "\\.")){
         # Making sure they've got a good extension
         Ext <- sub("\\.", "", str_extract(FileName, "\\..*"))
         FileName <- sub(paste0(".", Ext), "", FileName)
         Ext <- ifelse(Ext %in% c("eps", "ps", "jpeg", "tiff",
                                  "png", "bmp", "svg", "jpg"), 
                       Ext, "png")
         FileName <- paste0(FileName, ".", Ext)
      } else {
         FileName <- paste0(FileName, ".png")
      }
      
      ggsave(FileName, height = fig_height, width = fig_width, dpi = 600,
             plot = G)
   }
   
   return(G)
   
}