R/Plot_dose_response.R

In google/synergyfinderengineered: Calculate and Visualize Synergy Scores for Drug Combinations

# Copyright 2018 Google LLC
#
# Use of this source code is governed by a MIT-style
# license that can be found in the LICENSE file or at
# https://opensource.org/licenses/MIT.

#' Visualize the drug combination dose-response data
#'
#' A function to visualize the drug combination dose-response data
#'
#' @param data a list object generated by function \code{\link{ReshapeData}}.
#' @param save.file a parameter to specify wheter the visualization results are
#'   saved as pdf files in current working directory or not. If it is FALSE, the
#'   results are displayed locally. It is FALSE by default.
#' @param pair.index a parameter to specify which drug combination if there are
#'   many drug combinations in the data. By default, it is NULL so that the
#'   visualization of all the drug combinations in the data is returned.
#' @param Emin 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. Defaults to NA.
#' @param Emax 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. Defaults to NA.
#' @param ... further graphical parameters from \code{\link{plot}} for plotting
#'   the single drug dose-response curve. Use e.g., cex.lab to change the axis
#'   label size and cex.axis to change the tick size of axises.
#' @return there is no return value
#' @author Liye He \email{liye.he@helsinki.fi}
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' PlotDoseResponse(data)
PlotDoseResponse <- function (data, save.file = FALSE, pair.index = NULL, Emin = NA, Emax = NA, ...) {
  if(!is.list(data)) {
    stop("Input data is not a list format!")
  }
  dose.response.mats <- data$dose.response.mats
  drug.pairs <- data$drug.pairs
  num.pairs <- 1:length(dose.response.mats)
  plots <- list()
  if(!is.null(pair.index)) {
      num.pairs <- pair.index
  }
  for (i in num.pairs) {
    if(save.file) {
      pdf(file=paste(drug.pairs$drug.row[i],
                     drug.pairs$drug.col[i], "dose.response",
                     drug.pairs$blockIDs[i],
                     "pdf",
                     sep = "."),
          width=12,
          height=6)
    } else {
      dev.new(noRStudioGD = TRUE)
    }
    response.mat <- dose.response.mats[[i]]
    num.row <- length(response.mat)
    data.plot <- data.frame(x = numeric(num.row), y = numeric(num.row),
                            Inhibition = numeric(num.row))
    data.plot$Inhibition <- round(c(response.mat), 2)
    data.plot$y <- rep(c(1:ncol(response.mat)), nrow(response.mat))
    data.plot$x <- rep(1:nrow(response.mat), each = ncol(response.mat))
    data.plot$x <- as.factor(data.plot$x)
    data.plot$y <- as.factor(data.plot$y)
    conc.unit <- drug.pairs$concUnit[i] ## concentration unit
    unit.text <- paste("(", conc.unit, ")", sep = "")
    drug.row <- drug.pairs$drug.row[i]
    drug.col <- drug.pairs$drug.col[i]
    plot.title <- paste("Dose-response matrix (inhibition)", "\n BlockID:",
                        drug.pairs$blockIDs[i], sep = " ")


    # plot dose-response matrix
    axis.x.text <- round(as.numeric(colnames(response.mat)), 1)
    axis.y.text <- round(as.numeric(rownames(response.mat)), 1)
    dose.response.p <- ggplot2::ggplot(
      data.plot,
      ggplot2::aes_string(x = "x", y = "y")) +
      ggplot2::geom_tile(ggplot2::aes_string(fill = 'Inhibition')) +
      ggplot2::geom_text(ggplot2::aes_string(fill = 'Inhibition', label = 'Inhibition')) +
      ggplot2::scale_fill_gradient2(low = "green", high = "red", midpoint = 0, name = "Inhibiton (%)") +
      ggplot2::scale_x_discrete(labels = axis.x.text) +
      ggplot2::scale_y_discrete(labels = axis.y.text) +
      ggplot2::xlab(paste(drug.col, unit.text, sep = " ")) +
      ggplot2::ylab(paste(drug.row, unit.text, sep = " "))
    dose.response.p <- dose.response.p + ggplot2::theme(axis.text.x = ggplot2::element_text(color = "red", face = "bold", size = 15))
    dose.response.p <- dose.response.p + ggplot2::theme(axis.text.y = ggplot2::element_text(color = "red", face = "bold", size = 15))
    dose.response.p <- dose.response.p + ggplot2::theme(axis.title = ggplot2::element_text(size = 15))
    dose.response.p <- dose.response.p + ggplot2::ggtitle(plot.title) + ggplot2::theme(plot.title =
                                                                                         ggplot2::element_text(size = 25))


    single.fitted <- FittingSingleDrug(response.mat, fixed = c(NA, Emin, Emax, NA))


    layout(matrix(c(1, 3, 2, 3), 2, 2, byrow = TRUE))
    # plot the curve for the row drug
    suppressWarnings(par(mgp=c(3, .5, 0)))
    x.lab <- paste("Concentration", unit.text, sep = " ")
    plot(single.fitted$drug.row.model, xlab = x.lab, ylab = "Inhibition (%)", type = "obs", col = "red",
         cex = 1.5, pch = 16, ...)
    plot(single.fitted$drug.row.model, xlab = x.lab, ylab = "Inhibition (%)", type = "none",
         cex = 1.5, add = TRUE, lwd = 3)
    title(paste("Dose-response curve for drug:", drug.row, "in Block", drug.pairs$blockIDs[i]), cex.main = 1)




    # plot the curve for the col drug

    plot(single.fitted$drug.col.model, xlab = x.lab, ylab = "Inhibition (%)", type = "obs", col = "red",
         cex = 1.5, pch = 16, ...)
    plot(single.fitted$drug.col.model, xlab = x.lab, ylab = "Inhibition (%)", type = "none", cex = 1.5, add = TRUE, lwd = 3)
    title(paste("Dose-response curve for drug:", drug.col, "in Block", drug.pairs$blockIDs[i]), cex.main = 1)

    plot.new()
    #vps <- baseViewports()
    #pushViewport()
    print(dose.response.p, vp = grid::viewport(height = grid::unit(1, "npc"),
                                               width = grid::unit(0.5, "npc"),
                                               just = c("left","top"),
                                               y = 1, x = 0.5))
    #popViewport()
    merge.plot <- recordPlot()
    plots[[i]] <- merge.plot
    if(save.file) {
      dev.off()
    }
  }
}