R/plot_barometer.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
#
# SynergyFinder
#
# Functions on this page:
#
# PlotBarometer: Plot Barometer for Responses at One Data Point

#' Plot Barometer for Responses at One Data Point
#' 
#' This function will plot a barometer. The needle will point at the response
#' (\% inhibition) at the data point specified by \code{plot_concs}. The
#' reference additive effects calculated by different models will be marked on
#' the bar if they are included in input \code{data}.
#' 
#' @param data A list object generated by function 
#'   \code{\link{CalculateSynergy}}.
#' @param plot_block An integer or character. It indicates the block id for the
#'   combination matrix to visualize.
#' @param plot_concs A vector of numeric values with the length same as the
#'   number of drugs in selected block. It contains the concentrations  for
#'   "drug1", "drug2", ... The data point selected by these concentrations will
#'   be highlighted in the plot.
#' @param needle_color An R color value. It indicates the color of the needle.
#' @param needle_text_size A numeric value. It indicates the size of the text
#'   near the center of barometer which showing the response value. The unit is
#'   "mm".
#' @param needle_text_offset A numeric value. It is used to set the position of
#' the response values text. Smaller value means the text is closer to the
#' center.
#' @param graduation_color An R color value. It indicates the color of the
#'   graduation texts and ticks.
#' @param graduation_label_size A numeric value. It indicates the size of the
#'   graduation texts. The unit is "mm".
#' @param graduation_label_offset A numeric value. It is used to set the
#'   position of graduation texts. Smaller values means the graduation texts is
#'   closer to the ticks. It ranges from 0 to 1.
#' @param annotation_label_size A numeric value. It indicates the size of the
#'   labels for the additive effects at the out-most layer. The unit is "mm".
#' @param annotation_label_offset A numeric value. It is used to set the
#'   position of additive effect labels. Smaller values means the labels is
#'   closer to the color bar. It ranges from 0 to 1.
#' @param annotation_label_color An R color value. It indicates the color of the
#'   additive effects at the out-most layer.
#' @param font_family The font family for all the texts in the plot.
#' @param color_bar_color An R color value. It indicates the color of the
#'   largest value in the color bar.
#' @param color_bar_outer A numerical value. It indicates the proportion of the
#'   radius for the outer side of color bar comparing to the outermost edge for
#'   plotting area. It ranges from 0 to 1.
#' @param color_bar_inner A numerical value. It indicates the proportion of the
#'   radius for the inner side of color bar comparing to the outermost edge for
#'   plotting area. It ranges from 0 to 1.
#' @param major_graduation_outer A numerical value. It indicates the proportion
#'   of the radius for the outer side of graduation comparing to the outermost
#'   edge for plotting area. It ranges from 0 to 1.
#' @param minor_graduation_inner A numerical value. It indicates the proportion
#'   of the radius for the outer side of graduation comparing to the outermost
#'   edge for plotting area. It ranges from 0 to 1.
#' @param major_graduation_inner A numerical value. It indicates the proportion
#'   of the radius for the outer side of graduation comparing to the outermost
#'   edge for plotting area. It ranges from 0 to 1.
#' @param show_concs A logical value. If it is \code{TRUE}, the concentration of
#'   drugs will be shown on the plot.
#'   
#' @return A ggplot object.
#' 
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#' 
#' @references Tang J, Wennerberg K and Aittokallio T (2015) 
#' href{https://www.frontiersin.org/articles/10.3389/fphar.2015.00181/full}{What
#'  is synergy? The Saariselkä agreement revisited}. Front. Pharmacol. 6:181. 
#'  doi: 10.3389/fphar.2015.00181
#' 
#' @export
#'
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' data <- CalculateSynergy(data, method = c("ZIP", "HSA", "Bliss", "Loewe"))
#' p <- PlotBarometer(data, plot_block = 1, c(625, 50))
#' p
PlotBarometer <- function(data,
                          plot_block = 1,
                          plot_concs,
                          graduation_color = "#6C6C6C",
                          needle_color = "#6C6C6C",
                          needle_text_size = 5,
                          needle_text_offset = 2,
                          graduation_label_size = 4,
                          graduation_label_offset = 0.7,
                          annotation_label_size = 4,
                          annotation_label_offset = 0.6,
                          annotation_label_color = "#6C6C6C",
                          font_family = "",
                          color_bar_color = "#8f1b01",
                          color_bar_outer = 9,
                          color_bar_inner = 8,
                          major_graduation_outer = 7.8,
                          minor_graduation_inner = 7.5,
                          major_graduation_inner = 7,
                          show_concs = TRUE) {
  # Check plot_block
  if (!plot_block %in% data$drug_pairs$block_id) {
    stop("The input block id '", plot_block, "' could not be found in the input
         data.")
  }
  # Prepare data tables for plot
  drug_pair <- data$drug_pairs[which(data$drug_pairs$block_id == plot_block), ]
  concs <- grep("conc", colnames(data$response), value = TRUE)
  concs <- sort(concs)
  if (drug_pair$replicate) {
    plot_table <- data$response_statistics %>% 
      dplyr::filter(block_id == plot_block) %>% 
      dplyr::select(block_id, dplyr::all_of(concs), response = response_mean)
  } else {
    plot_table <- data$response %>% 
      dplyr::filter(block_id == plot_block) %>% 
      dplyr::select(block_id, dplyr::all_of(concs), response)
  }
  if ("synergy_scores" %in% names(data)) {
    plot_table <- plot_table %>% 
      dplyr::left_join(
        dplyr::filter(data$synergy_scores, block_id == plot_block),
        by = c("block_id", concs)
      ) %>% 
      dplyr::ungroup() %>% 
      dplyr::select(-block_id) %>% 
      dplyr::relocate(dplyr::any_of(concs))
  } else {
    plot_table <- plot_table %>% 
      dplyr::ungroup() %>% 
      dplyr::select(-block_id) %>% 
      dplyr::relocate(dplyr::any_of(concs))
  }
  
  # Check input "plot_concs"
  if (length(plot_concs) != length(concs)) {
    stop("The length of input 'plot_concs' is not equal to the number of ",
         "drugs in data. Please specify ", length(concs),
         " concentrations in 'plot_concs'.")
  }
  conc_exist <- NULL
  for (i in 1:length(concs)) {
    conc_exist[i] <- plot_concs[i] %in%
      plot_table[paste0("conc", i)][[1]]
  }
  if (!all(conc_exist)) {
    stop("The concentrations for drug ", paste(which(!conc_exist), collapse = ", "), 
         " specified by 'highlight_row' are not in data.")
  }
  selected_data <- plot_table[
    apply(
      plot_table[, concs],
      1,
      function(x) {
        all(x == plot_concs)
      }
    ),
  ]
  
  # Data table for color bar
  start_angle <- - pi * 1 / 4
  end_angle <- pi * 5 / 4
  angle_slice = (end_angle - start_angle)/100
  
  color_bar_data <- data.frame(
    start = seq(start_angle, end_angle, length.out = 101),
    end = seq(start_angle, end_angle, length.out = 101) + angle_slice,
    label = seq(0, 100, length.out = 101),
    stringsAsFactors = FALSE)
  
  # Data table for needle (The coordinate for vertex of triangle)
  needle_value <- selected_data$response
  
  # Shrink the values out of the range [0, 100]
  if (needle_value < 0) {
    needle_label_value <- needle_value
    needle_value <- needle_value / 10
  } else if (needle_value > 100) {
    needle_label_value <- needle_value
    needle_value <- 100 + (needle_value - 100) / 10
  } else {
    needle_label_value <- needle_value
  }
  
  theta_radius = angle_slice * needle_value + start_angle
  

  needle_length = (color_bar_outer + color_bar_inner) / 2
  needle <- data.frame(
    theta = c(theta_radius, theta_radius - pi /2, theta_radius + pi /2),
    r = c(needle_length, 0.15, 0.15),
    stringsAsFactors = FALSE
  ) %>% 
    dplyr::mutate(x = r * cos(theta), y = r * sin(theta))
  
  # Data table for reference effects
  ref <- grep(".*_ref", colnames(selected_data), value = TRUE)
  
  if (length(ref) == 0) {
    reference <- NULL
  } else {
    reference <- selected_data %>% 
      dplyr::select(dplyr::all_of(ref)) %>% 
      tidyr::gather(key = "label", value = "value")
    reference$value[which(reference$value < 0)] <- reference$value[reference$value < 0] /10
    reference$value[reference$value > 100] <- 100 + (reference$value[reference$value > 100] - 100)/10
    reference <- reference %>% 
      dplyr::mutate(
        adjust = 0,
        angle = value * angle_slice + start_angle
      ) %>% 
      dplyr::arrange(value)
    # Shrink the values out of the range [0, 100]

    reference$label <- sub("_ref", "", reference$label)
    # Separate overlapped labels
    if (nrow(reference) > 1) {
      for (i in 2:nrow(reference)){
        dif <- abs(reference$value[1:(i-1)] - reference$value[i]) < 
          0.7 * annotation_label_size
        if (any(dif)) {
          if (length(dif) == 1) {
            reference$adjust[i] <- 0.15 * annotation_label_size
          } else {
            reference$adjust[i] <- min(
              setdiff(
                seq(
                  0,
                  max(reference$adjust[1:(i-1)][which(dif)]) +
                    0.15 * annotation_label_size,
                  0.15 * annotation_label_size
                ),
                reference$adjust[1:(i-1)][which(dif)]
              )
            )
          }
        }
      }
    }
  }
  if (show_concs){
    # Generate text for concentrations
    conc_text <- sapply(1:length(plot_concs), function(i){
      paste0(
        drug_pair[, paste0("drug", i)], ": ",
        .RoundValues(plot_concs[i]),
        " (", drug_pair[, paste0("conc_unit", i)], ")"
      )
    })
    conc_text <- paste(conc_text, collapse = "\n")
    center_texts <- paste0(
      "[ ", 
      .RoundValues(needle_label_value), 
      "% ]\n",
      conc_text)
  } else {
    center_texts <- paste0(
      "[ ", 
      .RoundValues(needle_label_value), 
      "% ]\n")
  }
  p <- ggplot(color_bar_data) +
    ggforce::geom_arc_bar(
      data = color_bar_data[-101,],
      aes(
        x0 = 0,
        y0 = 0,
        r0 = color_bar_inner,
        r = color_bar_outer,
        start = start - pi/2,
        end = end - pi/2,
        fill = 100 - label,
        color = 100 - label
      )
    ) +
    ggplot2::scale_fill_gradient(
      high= "grey90",
      low = color_bar_color
    ) +
    ggplot2::scale_color_gradient(
      high= "grey90",
      low = color_bar_color
    ) +
    # minor graduations
    geom_segment(
      aes(
        x = - major_graduation_outer * cos(start),
        xend = - minor_graduation_inner * cos(start),
        y = major_graduation_outer * sin(start),
        yend = minor_graduation_inner * sin(start)
      ),
      color = graduation_color
    ) + 
    # major graduations
    geom_segment(
      data = subset(color_bar_data, label %% 10 == 0),
      aes(
        x = - major_graduation_outer * cos(start),
        xend = - major_graduation_inner * cos(start),
        y = major_graduation_outer * sin(start),
        yend = major_graduation_inner * sin(start)
      ),
      color = graduation_color
    ) +
    # graduation text
    geom_text(
      data = subset(color_bar_data, label %% 10 == 0),
      family = font_family,
      size = graduation_label_size, color = graduation_color,
      aes(
        x = -(major_graduation_inner - graduation_label_offset) * cos(start),
        y = (major_graduation_inner - graduation_label_offset) * sin(start),
        label = paste0(label, "%")
      )
    ) +
    theme(
      axis.title = element_blank(),
      axis.text = element_blank(),
      axis.ticks = element_blank(),
      panel.background = element_blank(),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      legend.position = "none"
    ) +
    coord_fixed() 
  # Add needle
  p <- p + 
    # needle
    geom_polygon(
      data = needle,
      aes(x = -x, y = y),
      fill = needle_color
    ) +
    geom_point(
      aes(x = 0, y = 0),
      colour = needle_color,
      size = 4
    ) + 
    # response text
    annotate(
      "text",
      x = 0,
      y = -needle_text_offset,
      size = needle_text_size,
      family = font_family,
      label = center_texts,
      color = needle_color
    )
  # Mark reference
  if (!is.null(reference)) {
    p <- p + 
      # mark reference
      geom_text(
        data = reference,
        size = annotation_label_size,
        family = font_family,
        color = annotation_label_color,
        aes(
          x = -(color_bar_outer + annotation_label_offset + adjust) * cos(angle),
          y = (color_bar_outer + annotation_label_offset + adjust) * sin(angle),
          label = label, angle = (pi /2 - angle) * 180 / pi
        )
      ) +
      geom_segment(
        data = reference,
        color = "grey90",
        aes(
          x = -color_bar_inner * cos(value * angle_slice + start_angle),
          xend = -color_bar_outer * cos(value * angle_slice + start_angle),
          y = color_bar_inner * sin(value * angle_slice + start_angle),
          yend = color_bar_outer * sin(value * angle_slice + start_angle)
        )
      )
  }
  return(p)
}