R/Plot_synergy.R

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

#'Drug interaction landscape
#'
#'A function to visualize the synergy scores for drug combinations as 2D or 3D
#'interaction landscape over the dose-response matrix.
#'
#' @param data a list object generated by function \code{\link{CalculateSynergy}}.
#' @param type a parameter to specify the type of the interaction landscape, 2D, 3D
#' or both. By default, 2D interaction landscape is returned.
#' @param save.file a logical parameter to specify if the interaction landscape is saved as a pdf file
#' in the current working directory or returned as an R object. By default, it is FALSE.
#' @param len a parameter to specify how many values need to be predicted between two concentrations
#' @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 synergy score visualization of all the drug combinations in the data is returned.
#' @param legend.start a parameter to specify the starting point of the legend. By defualt, it is NULL so the legend 
#' starting point is fixed by the data automatically.
#' @param legend.end a parameter to specify the ending point of the legend. By defualt, it is NULL so the legend 
#' ending point is fixed by the data automatically.
#' @param col.range a parameter to specify the starting and ending concentration of the drug on x-axis. Use e.g., c(1, 3) to 
#' specify that only from 1st to 3rd concentrations of the drug on x-axis are used. By default, it is NULl so all the 
#' concentrations are used.
#' @param row.range a parameter to specify the starting and ending concentration of the drug on y-axis. Use e.g., c(1, 3) to 
#' specify that only from 1st to 3rd concentrations of the drug on y-axis are used. By default, it is NULl so all the 
#' concentrations are used.
#' @return a pdf file or the interaction landscapes are only displayed depending on
#' the save.file parameter.
#' @author Liye He \email{liye.he@helsinki.fi}
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' scores <- CalculateSynergy(data)
#' PlotSynergy(scores, "2D")

PlotSynergy <- function(data, type = "2D", save.file = FALSE, len = 3, pair.index = NULL, legend.start = NULL, legend.end = NULL, row.range = NULL, col.range = NULL){
  if(!is.list(data)) {
    stop("Input data is not a list format!")
  }
  scores <- data$scores
  drug.pairs <- data$drug.pairs
  num.pairs <- 1:length(scores)
  plots <- list()
  if(!is.null(pair.index)) {
      num.pairs <- pair.index
  }
  for (i in num.pairs) {
    scores.dose <- scores[[i]]
    drug.col <- drug.pairs$drug.col[i]
    drug.row <- drug.pairs$drug.row[i]
    
    scores.tmp <- scores.dose
    if(!is.null(col.range)) {
        if(col.range[1] == 1) {
            # the first column is not included while calculating the average score
            scores.tmp <- scores.tmp[, (col.range[1] + 1):col.range[2]]
        } else {
            scores.tmp <- scores.tmp[,col.range[1]:col.range[2]]
        }
        
    } else {
      # delete the first column
      scores.tmp <- scores.tmp[, -1]
    }
    if(!is.null(row.range)) {
        if(row.range[1] == 1) {
            # the first row is not included while calculating the average score
            scores.tmp <- scores.tmp[(row.range[1] + 1):row.range[2],]
        } else {
            scores.tmp <- scores.tmp[row.range[1]:row.range[2], ]
        }
        
    } else {
      # delete the first row
      scores.tmp <- scores.tmp[-1, ]
    }

    
    
    summary.score <- round(mean(scores.tmp, na.rm = TRUE), 3)
    # drug.col: the col-wise drug
    # drug.row: the row-wise drug
    
    # kriging with kriging from SpatialExtremes package!
    
    row.conc <- as.numeric(rownames(scores.dose))  ## concentrations on the row
    col.conc <- as.numeric(colnames(scores.dose))  ## concentrations on the column
    
    nr <- nrow(scores.dose)
    nc <- ncol(scores.dose)
    
    # coordinates for the predicted values of the matrix
    
    # len: how many values need to be predicted between two concentrations
    scores.extended <- .ExtendedScores(scores.dose, len)
    
    # get the subset of the extended matrix with predicted values
    mat.tmp <- scores.extended
    if(!is.null(col.range)){
      # selecting cols
      mat.tmp[, which(pred.cory >= col.range[1] & pred.cory <= col.range[2])]
    }
    if(!is.null(row.range)){
      # selecting rows
      mat.tmp[which(pred.corx >= row.range[1] & pred.corx <= row.range[2]), ]
    }
    
    # the matrix for plotting
    pp <- matrix(NA, nrow(mat.tmp)*ncol(mat.tmp), 3)
    pp <- data.frame(pp)
    colnames(pp) <- c("x", "y", "z")
    pp$x <- rep(colnames(mat.tmp), each = nrow(mat.tmp))
    pp$y <- rep(rownames(mat.tmp), ncol(mat.tmp))
    pp$z <- c(mat.tmp)
    
    plot.title <- paste("Average synergy: ", summary.score, " (",data$method,")", sep = "")

    conc.runit <- drug.pairs$concRUnit[i] ## concentration unit row
    conc.cunit <- drug.pairs$concCUnit[i] ## concentration unit col
    unit.rtext <- paste("(", conc.runit, ")", sep = "")
    unit.ctext <- paste("(", conc.cunit, ")", sep = "")
    file.name <- paste(drug.row, drug.col, "synergy", drug.pairs$blockIDs[i], data$method, "pdf", sep = ".")
    drug.row <- paste(drug.row, unit.rtext, sep = " ")
    drug.col <- paste(drug.col, unit.ctext, sep = " ")
    max.dose <- max(abs(max(scores.dose)), abs(min(scores.dose)))
    color.range <- round(max.dose + 5, -1)
    if(is.null(legend.start)){
        start.point <- -color.range
    } else {
        start.point <- legend.start
    }
    
    if(is.null(legend.end)) {
        end.point <- color.range
    } else {
        end.point <- legend.end
    }
    
    # colors
    levels <- seq(start.point, end.point, by = 2)
    col1 <- colorRampPalette(c('green', "#FFFFFF"))(length(which(levels<=0)))
    col2 <- colorRampPalette(c("#FFFFFF", 'red'))(length(which(levels>=0)))
    col <- c(col1, col2[-1])
    
    
    

    if (type == "3D") {
      # x-axis ticks settings
      if(!is.null(col.range)) {
        xaxis <- list(at = seq(1, ncol(mat.tmp), by = len + 1),
                      labels = round(col.conc[col.range[1]:col.range[2]], 3))
      } else {
        xaxis <- list(at = seq(1, ncol(mat.tmp), by = len + 1),
                      labels = round(col.conc, 3))
      }
      
      # y-axis ticks settings
      if(!is.null(row.range)) {
        yaxis <- list(at = seq(1, nrow(mat.tmp), by = len + 1),
                   labels = round(row.conc[row.range[1]:row.range[2]], 3))
      } else {
        yaxis <- list(at = seq(1, nrow(mat.tmp), by = len + 1),
                   labels = round(row.conc, 3))
      }
      
      
      par1 <- list(arrows=FALSE, distance=c(0.8,0.8,0.8), col = 1,
                   cex = 0.8, z = list(tick.number = 6),
                   x=xaxis, y = yaxis)
      zlabs <- list(expression("Synergy score"), rot = 90, 
                    cex = 1, axis.key.padding = 0)
      xpar <- list(as.character(drug.col), cex = 1, rot = 20)
      ypar <- list(as.character(drug.row), cex = 1, rot = -50)
      
      fig <- wireframe(t(mat.tmp),scales = par1,
                       drape = TRUE, colorkey = list(space = "top",width = 0.5),
                       screen = list(z = 30, x = -55),
                       zlab = zlabs, xlab = xpar, ylab = ypar,
                       zlim = c(start.point, end.point),
                       col.regions = col,
                       main = plot.title,
                       at = do.breaks(c(start.point, end.point), length(col)),
                       par.settings = list(axis.line=list(col="transparent")),
                       zoom = 1, aspect = 1)
      print(fig)
      fig <- recordPlot()
    } else if (type == "2D") {
      dev.new(noRStudioGD = TRUE)
      layout(matrix(c(1, 2), nrow = 2L, ncol = 1L), heights = c(0.1, 1))
      par(mar = c(0, 6.1, 2.1, 4.1))
      suppressWarnings(par(mgp = c(3, -1.5, 0)))
      #levels <- seq(start.point, end.point, by = 2)
      #col <- colorpanel(end.point, low = "green", mid = "white", high = "red")
      plot.new()
      plot.window(ylim = c(0, 1), xlim = range(levels), xaxs = "i", yaxs = "i")
      rect(levels[-length(levels)],0, levels[-1L],0.3, col = col, border = NA)
      axis(3,tick = FALSE, at = do.breaks(c(start.point, end.point), end.point/10))
      title(plot.title)
      par(mar = c(5.1,4.1,1.1,2.1))
      suppressWarnings(par(mgp = c(2,1,0)))
      plot.new()
      mat.tmp <- t(mat.tmp)
      x.2D <- (1:dim(mat.tmp)[1] - 1) / (dim(mat.tmp)[1] - 1)
      y.2D <- (1:dim(mat.tmp)[2] - 1) / (dim(mat.tmp)[2] - 1)
      plot.window(asp = NA, xlim = range(x.2D), ylim = range(y.2D), "", xaxs = "i", yaxs = "i")
      .filled.contour(x.2D, y.2D, z = mat.tmp, levels, col = col)
      # #grid(dim(c)[1] - 1, dim(c)[2] - 1, col = "gray", lty = "solid", lwd = 0.3)
      box()
      mtext(drug.col, 1, cex = 1, padj = 3)
      mtext(drug.row, 2, cex = 1, las = 3, padj = -3)
      if(!is.null(row.range)) {
        yconc <- round(row.conc[row.range[1]:row.range[2]], 3)
      } else {
        yconc <- round(row.conc, 3)
      }
      
      if(!is.null(col.range)) {
        xconc <- round(col.conc[col.range[1]:col.range[2]], 3)
      } else {
        xconc <- round(col.conc, 3)
      }

      axis(side = 1, at = seq(0, 1, by = 1/(length(xconc) - 1)), labels = xconc)
      axis(side = 2, at = seq(0, 1, by = 1/(length(yconc) - 1)), labels = yconc)
      fig <- recordPlot()
      dev.off()
    } else {
      # 3d plot
      # x-axis ticks settings
      if(!is.null(col.range)) {
        xaxis <- list(at = seq(1, ncol(mat.tmp), by = len + 1),
                      labels = round(col.conc[col.range[1]:col.range[2]], 3))
      } else {
        xaxis <- list(at = seq(1, ncol(mat.tmp), by = len + 1),
                      labels = round(col.conc, 3))
      }
      
      # y-axis ticks settings
      if(!is.null(row.range)) {
        yaxis <- list(at = seq(1, nrow(mat.tmp), by = len + 1),
                      labels = round(row.conc[row.range[1]:row.range[2]], 3))
      } else {
        yaxis <- list(at = seq(1, nrow(mat.tmp), by = len + 1),
                      labels = round(row.conc, 3))
      }
      
      
      par1 <- list(arrows=FALSE, distance=c(0.8,0.8,0.8), col = 1,
                   cex = 0.8, z = list(tick.number = 6),
                   x=xaxis, y = yaxis)
      zlabs <- list(expression("Synergy score"), rot = 90, 
                    cex = 1, axis.key.padding = 0)
      xpar <- list(as.character(drug.col), cex = 1, rot = 20)
      ypar <- list(as.character(drug.row), cex = 1, rot = -50)
      
      syn.3d.plot <- wireframe(t(mat.tmp),scales = par1,
                       drape = TRUE, colorkey = list(space = "top",width = 0.5),
                       screen = list(z = 30, x = -55),
                       zlab = zlabs, xlab = xpar, ylab = ypar,
                       zlim = c(start.point, end.point),
                       col.regions = col,
                       main = plot.title,
                       at = do.breaks(c(start.point, end.point), length(col)),
                       par.settings = list(axis.line=list(col="transparent")),
                       zoom = 1, aspect = 1)
      # 2d plot
      layout(matrix(c(1, 2, 3, 3), nrow = 2L, ncol = 2L), heights = c(0.1, 1))
      par(mar = c(0, 6.1, 2.1, 4.1))
      suppressWarnings(par(mgp = c(3, -0.8, 0)))
      plot.new()
      plot.window(ylim = c(0, 1), xlim = range(levels), xaxs = "i", yaxs = "i")
      rect(levels[-length(levels)],0, levels[-1L],0.3, col = col, border = NA)
      axis(3,tick = FALSE, at = do.breaks(c(start.point, end.point), end.point/10))
      title(plot.title)
      par(mar = c(5.1,4.1,1.1,2.1))
      suppressWarnings(par(mgp = c(2,1,0)))
      plot.new()
      mat.tmp <- t(mat.tmp)
      x.2D <- (1:dim(mat.tmp)[1] - 1) / (dim(mat.tmp)[1] - 1)
      y.2D <- (1:dim(mat.tmp)[2] - 1) / (dim(mat.tmp)[2] - 1)
      plot.window(asp = NA, xlim = range(x.2D), ylim = range(y.2D), "", xaxs = "i", yaxs = "i")
      .filled.contour(x.2D, y.2D, z = mat.tmp, levels, col = col)
      box()
      mtext(drug.col, 1, cex = 1, padj = 3)
      mtext(drug.row, 2, cex = 1, las = 3, padj = -3)
      if(!is.null(row.range)) {
        yconc <- round(row.conc[row.range[1]:row.range[2]], 3)
      } else {
        yconc <- round(row.conc, 3)
      }
      
      if(!is.null(col.range)) {
        xconc <- round(col.conc[col.range[1]:col.range[2]], 3)
      } else {
        xconc <- round(col.conc, 3)
      }
      
      axis(side = 1, at = seq(0, 1, by = 1/(length(xconc) - 1)), labels = xconc)
      axis(side = 2, at = seq(0, 1, by = 1/(length(yconc) - 1)), labels = yconc)
      print(syn.3d.plot, position = c(0.5,0, 1, 1), newpage = FALSE)
      fig <- recordPlot()
    }
    plots[[i]] <- fig
    if(save.file) {
      if (type %in% c("2D", "3D")) {
        pdf(file.name, width = 10, height = 10)
      } else {
        pdf(file.name, width = 12, height = 6)
      }
      print(fig)
      dev.off()
    }

  }

  if(!save.file) {
    for(i in num.pairs) {
      dev.new(noRStudioGD = TRUE)
      replayPlot(plots[[i]])
    }
  }


}

.ExtendedScores <- function(scores.dose, len) {
  # len: how many values need to be predicted between two adjacent elements
  #      of scores.dose
  nr <- nrow(scores.dose)
  nc <- ncol(scores.dose)
  extended.row.idx <- seq(1, nr, length = (nr - 1)*(len + 2) - (nr - 2))
  extended.col.idx <- seq(1, nc, length = (nc - 1)*(len + 2) - (nc - 2))
  krig.coord <- cbind(rep(extended.row.idx, each = length(extended.col.idx)), 
                      rep(extended.col.idx, length(extended.row.idx)))
  extended.scores <- SpatialExtremes::kriging(
    data = c(scores.dose),
    data.coord = cbind(rep(1:nr, nc), rep(1:nc, each = nr)),
    krig.coord = krig.coord,
    cov.mod = "whitmat",
    grid = FALSE,
    sill = 1,
    range = 10,
    smooth = 0.8
  )$krig.est
  extended.scores <- matrix(extended.scores, 
                            nrow = length(extended.row.idx),
                            ncol = length(extended.col.idx),
                            byrow = TRUE)	
  extended.scores
}