R/Plot_synergy.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.
#

library(SpatialExtremes)
library(gplots)
library(lattice)

#' Drug interaction landscape
#'
#' A function to visualize the synergy scores for drug combinations as 2D or 3D
#' interaction landscape over the dose-response matrix.
#'
#' The interaction landscapes are displayed on the current graphics device no
#' matter what.  If the save.file parameter is TRUE, they are also saved to PDF
#' files.
#'
#' @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 len a var that influences statistical smoothing of generated plots.
#' @param save.file when TRUE, saves the interaction landscape a pdf file in the
#'   current working directory.  When FALSE, displays the landscape in a local
#'   graphics window. Note that local display will always use a window with
#'   specified height and width, and will never use the RStudio plot viewer.  By
#'   default, it is FALSE.
#' @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 default, 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
#'   default, it is NULL so the legend ending point is fixed by the data
#'   automatically.
#' @param x.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 y.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 NULL
#' @author Liye He \email{liye.he@helsinki.fi}
#' @author David L. Roxe \email{dlroxe@google.com}
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' scores <- CalculateSynergy(data)
#' PlotSynergy(scores, "2D")
PlotSynergy <- function(data,
                        type = "2D",
                        len = 3,
                        save.file = FALSE,
                        pair.index = NULL,
                        legend.start = NULL,
                        legend.end = NULL,
                        x.range = NULL,
                        y.range = NULL) {
  if(!is.list(data)) {
    stop("Input data is not of type 'list'!")
  }

  plotDataList <- lapply(
    if(is.null(pair.index)) 1:length(data$scores) else pair.index,
    function(i) {
      .GetPlotData(data, i, legend.start, legend.end, y.range, x.range, len)
    })

  lapply(
    plotDataList,
    function(plotData) {
      # R doesn't seem to allow generic method dispatch based on string value of
      # 'type'.  So, we just assign a specialized function with if/else-if/else.
      PlotFunc <-
        if(type=="2D") .Plot2D else if(type=="3D") .Plot3D else .PlotAll

      if(save.file) {
        pdf(
          file=paste(plotData$metadata$file.name.prefix, "synergy",
                     plotData$metadata$drug.blockID, data$method, type, "pdf",
                     sep="."),
          width=if(type == "all") 12 else 10,
          height=if(type == "all") 6 else 10)
        print(do.call(PlotFunc, plotData$data))
        dev.off()
      } else {
        dev.new(noRStudioGD = TRUE)
        print(do.call(PlotFunc, plotData$data))
      }
    })
  return(NULL)
}

# Returns nested associative lists with all data needed for plots of any type.
.GetPlotData <- function(
  data, index, legend.start, legend.end, row.range, col.range, len) {
  scores.dose <- data$scores[[index]]
  color.range <- .RoundUpToNextTen(scores.dose)
  unit.text <- paste("(", data$drug.pairs$concUnit[index], ")", sep = "")
  return(list(
    metadata = list(  # useful for naming output files
      drug.blockID = data$drug.pairs$blockIDs[index],
      file.name.prefix = paste(data$drug.pairs$drug.row[index],
                               data$drug.pairs$drug.col[index], sep = '.')
    ),
    data = list(  # useful for actual plot content
      plot.title =
        .ComputePlotTitle(scores.dose, row.range, col.range, data$method),
      # Note, when considering application of ranges, that with len=3 mat.tmp will
      # have 9 times as many elements as scores.dose.  It is best to invoke
      # .ExtendedScores only after ranges have been applied.
      mat.tmp =
        .ExtendedScores(.ApplyRanges(scores.dose, row.range, col.range), len),
      len = len,
      row.conc = .NumericApplyRange(rownames(scores.dose), row.range),
      col.conc = .NumericApplyRange(colnames(scores.dose), col.range),
      drug.row = paste(data$drug.pairs$drug.row[index], unit.text, sep = " "),
      drug.col = paste(data$drug.pairs$drug.col[index], unit.text, sep = " "),
      start.point = if(is.null(legend.start)) -color.range else legend.start,
      end.point = if(is.null(legend.end)) color.range else legend.end)))
}

# Returns the highest absolute value in 'scores', rounded up to next ten.
.RoundUpToNextTen <- function(scores) {
  return(round(5 + max(abs(max(scores)),
                       abs(min(scores))),
               -1))
}

# Returns a string that may be used as a plot title.  The string includes
# the mean value of all scores in the selected ranges, and the name of the
# method used to compute them.
.ComputePlotTitle <- function(scores.dose, row.range, col.range, data.method) {
  # The first row and column are removed no matter what,
  # presumably (? TODO/REVIEW) on the assumptions that:
  #  1. For each drug, the first recorded effect is for dose=0.
  #  2. When dose=0 for either drug, the synergy score is also 0.
  return(paste("Average synergy: ",
               round(mean(scores.dose[.SliceAtLeastFirst(row.range),
                                      .SliceAtLeastFirst(col.range)],
                          na.rm = TRUE),
                     3),
               " (",data.method,")",
               sep = ""))
}

# Convert 'val' to a slice.  However, guarantee that the slice starts at
# 2 or higher.  So, both for val=c(1,5) and val=c(2,5), the function should
# return 2:5.  The first two elements of list 'val' are taken to be the slice
# endpoints; any other elements are ignored.  If 'val' is NULL, then the
# returned slice is '-1' -- 'remove the first element'.
.SliceAtLeastFirst <- function(val) {
  return(if(is.null(val)) -1 else max(val[1],2):val[2])
}

# Returns a slice of 'val' specified by 'range', as numeric.
# If 'range' is null, 'val' is wholly returned.  The first two elements of
# 'range' are taken to be slice endpoints; any other values are ignored.
# The desired portion of 'val' is converted to numeric form, then returned.  An
# error occurs if any element of the selected slice of 'val' cannot be converted
# to a number.
.NumericApplyRange <- function(val, range) {
  return(as.numeric(if(is.null(range)) val else val[range[1]:range[2]]))
}

# Returns a 2D slice of 'data' specified by 'row.range' and 'col.range'.  If
# either range is NULL, the data for the corresponding axis of 'data' is
# returned in full.  The row.range and col.range inputs are each treated as
# lists whose first two elements are range endpoints, and whose remaining
# elements (if any) are ignored.
.ApplyRanges <- function(data, row.range, col.range) {
  .RangeWithDefault <- function(range, default) {
    return(if(is.null(range)) default else range[1]:range[2])
  }
  # Default value is 'remove the element just after the end of the list',
  # for each dimension of the array (a no-op that permits cleaner syntax).
  return(data[.RangeWithDefault(row.range, default=-(nrow(data)+1)),
              .RangeWithDefault(col.range, default=-(ncol(data)+1))])
}

.GetKrigingObservationCoords <- function(data) {
  nr <- nrow(data)
  nc <- ncol(data)
  if (nr > 0 && nc > 0) {
    xtics <- 1:nr
    ytics <- 1:nc
    return(cbind(rep(xtics, nc),
                 rep(ytics, each = nr)))
  }
}

.GetKrigingInterpolationParams <- function(data, len) {
  nr <- nrow(data)
  nc <- ncol(data)
  if (nr > 0 && nc > 0) {
    krig.row.len <- nr + len * (nr - 1)
    krig.col.len <- nc + len * (nc - 1)
    krig.xtics <- seq(1, nr, length=krig.row.len)
    krig.ytics <- seq(1, nc, length=krig.col.len)
    return(list(rows = krig.row.len,
                cols = krig.col.len,
                coord = (cbind(rep(krig.xtics, krig.col.len),
                              rep(krig.ytics, each = krig.row.len)))))
  }
}

# Perform statistical interpolation to smooth plotted values.
.ExtendedScores <- function(scores.dose, len) {
  # len: how many values need to be predicted between two adjacent elements
  #      of scores.dose
  kriging.params = .GetKrigingInterpolationParams(scores.dose, len)
  # Note, not all kriging implementations are suitable for all purposes.  In
  # particular, SpatialExtremes offers an implementation that works with small
  # data sets.  Earlier versions of this implementation that used other kriging
  # packages only worked well with large ones.
  return(
    matrix(
      SpatialExtremes::kriging(
        data = c(scores.dose),
        data.coord = .GetKrigingObservationCoords(scores.dose),
        krig.coord = kriging.params$coord,
        cov.mod = "whitmat",
        grid = FALSE,
        sill = 1,
        range = 10,
        smooth = 0.8
        )$krig.est,
      nrow=kriging.params$rows,
      ncol=kriging.params$cols))
}

.Plot2D <- function(plot.title, mat.tmp, len, row.conc, col.conc,
                    drug.row, drug.col, start.point, end.point) {
  # Note a different layout() here than in .PlotAll()
  layout(matrix(c(1, 2), nrow = 2L, ncol = 1L), heights = c(0.1, 1))
  .Plot2DTitle(plot.title, start.point, end.point)
  .Plot2DContourMap(
    mat.tmp, row.conc, col.conc, drug.row, drug.col, start.point, end.point)
  return(recordPlot())
}

.Plot2DTitle <- function (plot.title, start.point, end.point) {
  par(mar = c(0, 6.1, 2.1, 4.1))
  suppressWarnings(par(mgp = c(3, -1, 0)))
  levels <- seq(start.point, end.point, by = 2)
  plot.new()
  plot.window(ylim = c(0, 1), xlim = range(levels), xaxs = "i", yaxs = "i")
  rect(xleft = levels[-length(levels)],
       ybottom = 0,
       xright = levels[-1L],
       ytop = 0.3,
       col = .GetColors(start.point, end.point),
       border = NA)
  title(plot.title)
  axis(3,
       tick = FALSE,
       at = lattice::do.breaks(c(start.point, end.point), end.point/10))
  par(mar = c(5.1,4.1,2.1,2.1))
  suppressWarnings(par(mgp = c(2,1,0)))
}

.Plot2DContourMap <- function(
  mat.tmp, row.conc, col.conc, drug.row, drug.col, start.point, end.point) {
  plot.new()
  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),
              log = "",  # no axis uses log scale
              xaxs = "i",
              yaxs = "i")
  # TODO: Find a way to stop depending on private method .filled.contour().
  .filled.contour(x = x.2D,
                  y = y.2D,
                  z = mat.tmp,
                  levels = seq(start.point, end.point, by = 2),
                  col = .GetColors(start.point, end.point))
  box()
  mtext(drug.col, 1, cex = 1, padj = 3)
  mtext(drug.row, 2, cex = 1, las = 3, padj = -3)

  axis(side = 1,
       at = seq(0, 1, by = 1/(length(col.conc) - 1)),
       labels = round(col.conc, 1))
  axis(side = 2,
       at = seq(0, 1, by = 1/(length(row.conc) - 1)),
       labels = round(row.conc, 1))
}

# Default values for 'extra' args position and new.page let us take calls via
# the PlotFunc indirection in .PlotScores(), yet still re-use this function from
# .PlotAll(), which must specify a layout that places this graph beside the 2D
# contour plot.
.Plot3D <- function(
  plot.title, mat.tmp, len, row.conc, col.conc, drug.row, drug.col,
  start.point, end.point, position = c(0,0,1,1), new.page = TRUE) {
  print(
    lattice::wireframe(
      x = mat.tmp,
      scales = list(arrows=FALSE,
                    distance=c(0.8,0.8,0.8),
                    col = 1,
                    cex = 0.8,
                    z = list(tick.number = 6),
                    x = list(at = seq(1, nrow(mat.tmp), by = len + 1),
                             labels = round(col.conc, 3)),
                    y = list(at = seq(1, nrow(mat.tmp), by = len + 1),
                             labels = round(row.conc, 3))),
      drape = TRUE,
      colorkey = list(space="top",width=0.5),
      screen = list(z = 30, x = -55),
      zlab=list(expression("Synergy score"),rot=90,cex=1,axis.key.padding=0),
      xlab=list(as.character(drug.col),cex=1, rot=20),
      ylab=list(as.character(drug.row),cex=1,rot=-50),
      zlim=c(start.point, end.point),
      # .GetColors(-100,100) =~ colorRampPalette(c("green","white","red"))(101)
      col.regions=.GetColors(-100,100),
      main = plot.title,
      at=lattice::do.breaks(c(start.point, end.point),100),
      par.settings = list(axis.line=list(col="transparent")),
      zoom = 1,
      aspect = 1
    ),
    position = position,
    newpage = new.page)
  return(recordPlot())
}

.PlotAll <- function(
  plot.title, mat.tmp, len, row.conc, col.conc, drug.row, drug.col,
  start.point, end.point) {
  # Note a different layout() here than in .Plot2D().
  layout(matrix(c(1, 2, 3, 3), nrow = 2L, ncol = 2L), heights = c(0.1, 1))
  .Plot2DTitle(plot.title, start.point, end.point)
  .Plot2DContourMap(
    mat.tmp, row.conc, col.conc, drug.row, drug.col, start.point, end.point)
  .Plot3D(plot.title, mat.tmp, len, row.conc, col.conc, drug.row, drug.col,
          start.point, end.point, position = c(0.5,0, 1, 1), new.page = FALSE)
  return(recordPlot())
}

.GetColors <- function(start.point, end.point) {
  levels <- seq(start.point, end.point, by = 2)
  col1 <- colorRampPalette(c('green', 'white'))(length(which(levels<=0)))
  col2 <- colorRampPalette(c('white', 'red'))(length(which(levels>=0)))
  return(c(col1, col2[-1]))
}