R/Plot_dose_response.R

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

#' 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 if the visualization results are saved as pdf files in current working directory or not. If it is FALSE, the results are
#' returned as a list of the plots. 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 if save.file parameter is TRUE, pdf files are returned. Otherwise, the plots are only displayed.
#' @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) {
    #dev.new()
    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.runit <- drug.pairs$concRUnit[i] ## concentration row unit
    conc.cunit <- drug.pairs$concCUnit[i] ## concentration col unit
    runit.text <- paste("(", conc.runit, ")", sep = "")
    cunit.text <- paste("(", conc.cunit, ")", 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 <- ggplot(data.plot, aes_string(x = "x", y = "y")) + geom_tile(aes_string(fill = 'Inhibition')) +
      geom_text(aes_string(fill = 'Inhibition', label = 'Inhibition')) +
      scale_fill_gradient2(low = "green", high = "red", midpoint = 0, name = "Inhibition (%)") +
      scale_x_discrete(labels = axis.x.text) + scale_y_discrete(labels = axis.y.text) +
      xlab(paste(drug.col, runit.text, sep = " ")) + ylab(paste(drug.row, cunit.text, sep = " "))
    dose.response.p <- dose.response.p + theme(axis.text.x = element_text(color = "red", face = "bold", size = 15))
    dose.response.p <- dose.response.p + theme(axis.text.y = element_text(color = "red", face = "bold", size = 15))
    dose.response.p <- dose.response.p + theme(axis.title = element_text(size = 15))
    dose.response.p <- dose.response.p + ggtitle(plot.title) + theme(plot.title = 
                                                                    element_text(size = 20))
    

    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", runit.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
    x.lab <- paste("Concentration", cunit.text, sep = " ")
    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 = viewport(height = unit(1, "npc"), width = unit(0.5, "npc"), just = c("left","top"),
                                         y = 1, x = 0.5))
    #popViewport()
    merge.plot <- recordPlot()
    plots[[i]] <- merge.plot
    if(save.file) {
      file.name <- paste(drug.row, drug.col, "dose.response", drug.pairs$blockIDs[i], "pdf", sep = ".")
      pdf(file.name, width = 12, height = 6)
      replayPlot(merge.plot)
      dev.off()
    }
  }
  if(!save.file) {
    for(i in num.pairs) {
      dev.new(noRStudioGD = TRUE)
      replayPlot(plots[[i]])
    }
  }
}