R/plot_dose_response_curve.R

In shuyuzheng/synergyfinder: Calculate and Visualize Synergy Scores for Drug Combinations

# Copyright Shuyu Zheng and Jing Tang - All Rights Reserved
# Unauthorized copying of this file, via any medium is strictly prohibited
# Proprietary and confidential
# Written by Shuyu Zheng <shuyu.zheng@helsinki.fi>, March 2021
#
# Functions on this page:
#
# PlotDoseResponseCurve: Plot Dose Response Curve for Single Drug

#' Plot Dose Response Curve for Single Drug
#' 
#' This function will pot the dose response curve fitted by 4 parameters 
#' log-logistic curve.
#'
#' @param data A list object generated by function \code{\link{ReshapeData}}.
#' @param plot_block A character/integer. It indicates the block ID for the
#'   block to visualize.
#' @param drug_index A character/integer. It indicates the index of the drug to 
#'   plot. For example, 1 or "1" indicates to plot "drug1" in the input 
#' \code{data}.
#' @param adjusted A logical value. If it is \code{FALSE}, original response
#'   value will be used to fit the curve. If it is \code{TRUE}, the response
#'   values adjusted by (adding noise and/or imputation) will be plotted.
#' @param Emin A numeric or \code{NA}. the minimal effect of the drug used in
#'   the 4-parameter log-logistic function to fit the dose-response curve. If
#'   it is not NA, it is fixed the value assigned by the user. Default setting
#'   is \code{NA}.
#' @param Emax A numeric or \code{NA}. the maximal effect of the drug used in
#'   the 4-parameter log-logistic function to fit the dose-response curve. If
#'   it is not NA, it is fixed the value assigned by the user. Default setting
#'   is \code{NA}.
#' @param grid A logical value. It indicates whether to show the grids in the
#'   plots. The default value is \code{TRUE}. If \code{grid = NULL} or
#'   \code{grid = FALSE}, there is no grid.
#' @param point_color An R color value. It indicates the color for points in
#'   dose response curve plots.
#' @param curve_color An R color value. It indicates the color for curves in 
#'   dose response curve plots.
#' @param plot_title A character value or NULL. It specifies the plot title.
#'   If it is \code{NULL}, the function will automatically generate a title.
#' @param plot_subtitle A character or NULL. It indicates the subtitle for the
#'   plot. If it is \code{NULL}, the function will automatically generate a 
#'   subtitle.
#' @param plot_setting A list of graphical arguments. The arguments are passed 
#'   to \link[graphics]{par} function to modify the appearance of plots.
#' @param plot_new A logic value. If it is \code{TRUE}, a new device will be
#'   initiate with \link[graphics]{plot.new}. You might want to set it as
#'   \code{FALSE} while combining with other plots by using
#'    \link[graphics]{layout} function.
#' @param record_plot A logic value. If it is \code{TRUE}, a plot object 
#'   recorded by \link[grDevices]{recordPlot} will be returned. If it is
#'   \code{FALSE}, this function will return \code{NULL}.
#' @param text_size_scale A numeric value. It is used to control the size
#'   of text in the plot. All the text size will multiply by this scale factor.
#' @param ... Additional graphical arguments that are inherited from
#'   link[drc]{plot.drc} function. For example, use xlim = c(0.5, 500) or
#'   ylim = (0, 100) to control the ranges of x-axis or y-axis, respectively. 
#' @return A plot object recorded by \link[grDevices]{recordPlot} or NULL.
#' 
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#' 
#' @export
#'
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' p <- PlotDoseResponseCurve(data, plot_block = 1, drug_index = 2, grid = NULL)
#' replayPlot(p)
PlotDoseResponseCurve <- function(data,
                                  plot_block = 1,
                                  drug_index = 1,
                                  adjusted = TRUE,
                                  Emin = NA,
                                  Emax = NA,
                                  grid = TRUE,
                                  point_color = "#C24B40",
                                  curve_color = "black",
                                  text_size_scale = 1,
                                  plot_title = NULL,
                                  plot_subtitle = NULL,
                                  plot_setting = list(
                                    cex.lab = 1 * text_size_scale,
                                    mgp = c(2, 0.5, 0),
                                    font.main = 2,
                                    font.lab = 1,
                                    cex.main = 14 / 12 * text_size_scale,
                                    bty = "l",
                                    lwd = 1.5
                                  ),
                                  plot_new = TRUE,
                                  record_plot = TRUE,
                                  ...) {
  
  # 1. Check the input data
  # Data structure of 'data'
  if (!is.list(data)) {
    stop("Input data is not in list format!")
  }
  if (!all(c("drug_pairs", "response") %in% names(data))) {
    stop("Input data should contain at least tow elements: 'drug_pairs' and 
         'response'. Please prepare your data with 'ReshapeData' function.")
  }
  # Parameter 'plot_block'
  if (length(drug_index) != 1) {
    stop("The length of 'plot_block' parameter is not 1. Please chosed only one
         block for visualization.")
  }
  # Parameter 'drug'
  concs <- grep("conc\\d", colnames(data$response), value = TRUE)
  if (length(drug_index) != 1) {
    stop("The length of 'drug' parameter is not 1. Please chosed only one
         drug in the block for visualization.")
  } else if (!drug_index %in% sub("conc", "", concs)) {
    stop("The input drug_index '", drug_index, "' is not found in input data. ",
         "Available indexes are '", 
         paste(sub("conc", "", concs), collapse = "', '"),
         "'.")
  }
  
  # Annotation data
  drug_anno <- data$drug_pairs[data$drug_pairs$block_id == plot_block, ] %>% 
    dplyr::select(drug = paste0("drug", drug_index),
                  unit = paste0("conc_unit", drug_index))
  
  # Extract single drug dose response
  if (!adjusted) {
    response <- data$response %>% 
      dplyr::select(-response) %>% 
      dplyr::rename(response = response_origin) %>% 
      dplyr::filter(block_id == plot_block)
  } else {
    response <- data$response %>% 
      dplyr::filter(block_id == plot_block)
  }
  single_drug_data <- ExtractSingleDrug(response)
  single_drug_data <- single_drug_data[[paste0("conc", drug_index)]]
  # Fit model for the row drug
  drug_model <- suppressWarnings(
    FitDoseResponse(
      single_drug_data,
      Emin = Emin,
      Emax = Emax
    )
  )
  
  if (is.null(plot_subtitle)) {
    plot_subtitle <- paste(
      drug_anno$drug,
      "in Block",
      plot_block
    )
  }
  if (is.null(plot_title)) {
    plot_title <- "Dose-Response Curve"
  }
  
  # plot the curve for the drug
  # For all of R's graphical devices, the default text size is 12 points but it
  # can be reset by including a pointsize argument to the function that opens
  #the graphical device. From ?pdf:
  # curve_pred <- data.frame(
  #   dose = seq(0, max(single_drug_data$dose), length.out = 700),
  #   response = PredictModelSpecify(
  #     drug_model, 
  #     seq(0, max(single_drug_data$dose), length.out = 700)
  #   ),
  #   stringsAsFactors = FALSE
  # )
  # p1 <- ggplot(single_drug_data) +
  #   geom_point(aes(log10(dose), response)) +
  #   geom_line(aes(log10(dose), response), data = curve_pred) + 
  #   scale_y_continuous(trans = log10_trans())
  # p1
  
  model_type <- FindModelType(drug_model)
  if (model_type == "LL.4"){
    # Set break point of x-axis (numbers smaller than it are set as 0)
    bp <- round(min(log10(drug_model$origData$dose[drug_model$origData$dose > 1e-10])))-1
    max <- round(max(log10(drug_model$origData$dose)))
    step <- (max - bp)%/%4
    if (step < 1){
      step <- 1
    }
    xt <- 10 ^ seq(bp, max, by = step)
    xtlab <- xt
    xtlab[1] <- 0
  } else { # model type is L.4
    # Set break point of x-axis (numbers smaller than it are set as 0)
    bp <- round(min(log10(drug_model$origData$dose[drug_model$origData$dose > 1e-10])))-1
    # if (bp == 0) {
    #   bp <- -1
    # }
    max <- round(max(log10(drug_model$origData$dose)))
    step <- (max - bp)%/%4
    if (step < 1){
      step <- 1
    }
    # Set x-axis tick markders
    xt <- 10 ^ seq(bp, max, by = step)
    xtlab <- xt
    xtlab[1] <- 0
    drug_model$dataList$dose <- exp(drug_model$dataList$dose)
    drug_model$dataList$dose[drug_model$dataList$dose < 10^bp] <- 10^bp
  }
  
  if (is.null(grid)){
    grid_exp <- NULL
  } else if (grid){
    grid_exp <- expression(
      {
        graphics::grid(nx = NA, ny = NULL, col = "#DFDFDF", lty = 1)
        graphics::abline(col = "#DFDFDF", v = xt)
      }
    )
  } else {
    grid_exp <- NULL
  }
  
  if (plot_new) {
    graphics::plot.new()
    grDevices::dev.control("enable")
  }

  suppressWarnings(graphics::par(plot_setting))
  
  # Plot dots
  graphics::plot(
    x = drug_model,
    xlab = paste0("Concentration (", drug_anno$unit, ")"),
    ylab = "Inhibition (%)",
    type = "obs",
    log = "x",
    col = point_color,
    pch = 16,
    cex.axis = 1 * text_size_scale,
    panel.first = eval(grid_exp),
    xttrim = FALSE,
    conName = NULL,
    bp = 10 ^ bp,
    xt = xt,
    xtlab = xtlab,
    ...
  )
  
  # Plot curve
  graphics::plot(
    x = drug_model,
    type = "none",
    log = "x",
    col = curve_color,
    cex.axis = 1 * text_size_scale,
    add = TRUE,
    lwd = 3,
    xttrim = FALSE,
    conName = NULL,
    bp = 10 ^ bp,
    xt = xt,
    xtlab = xtlab,
    ...
  )

  # Plot title
  graphics::title(plot_title)
  graphics::mtext(
    plot_subtitle,
    cex = 7/9 * graphics::par()$cex.main * text_size_scale
  )
  
  if (record_plot) {
    p <- grDevices::recordPlot()
    grDevices::dev.off()
    print(p)
    return(p)
  } else {
    return(NULL)
  }
}