R/plot_two_drugs.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 in this page:
# Plot2DrugHeatmap: Heatmap Plot for 2-drug Combination Dose-Response/Synergy
#                   Scores
# Plot2DrugSurface: 3D surface Plot for 2-drug Combination Dose-Response/Synergy
#                   Scores
# Auxiliary functions:
# .Extract2DrugPlotData: Extract Data for 2 Drug combination plots
# .RoundValues: Round the Numbers for Plotting
# .Pt2mm: Round the Numbers for Plotting
# .ExtendedScores: Convert Font Size from pt to mm

#' Heatmap Plot for 2-drug Combination Dose-Response/Synergy Scores
#' 
#' This function will generate a plot for 2-drug combinations. The axes are the
#' dosage at which drugs are tested. The values could be observed response,
#' synergy scores or the reference effects calculated from different models.
#'
#' @param data A list object generated by function \code{\link{ReshapeData}}.
#' @param plot_block A character/integer. It indicates the block ID for the
#'   block to visualize.
#' @param drugs A vector of characters or integers with length of 2. It contains
#'   the index for two drugs to plot. For example, \code{c(1, 2)} indicates to
#'   plot "drug1" and "drug2" in the input \code{data}.
#' @param plot_value A character value. It indicates the value to be visualized.
#'   If the \code{data} is the direct output from \link{ReshapeData}, the values
#'   for this parameter are:
#'   \itemize{
#'     \item \strong{response_origin} The original response value in input data.
#'     It might be \% inhibition or \% viability.
#'     \item \strong{response} The \% inhibition after preprocess by function 
#'     \link{ReshapeData}
#'   }
#'   If the \code{data} is the output from \link{CalculateSynergy}, following
#'   values are also available:
#'   \itemize{
#'     \item \strong{ZIP_ref, Bliss_ref, HSA_ref, Loewe_ref} The reference
#'     additive effects calculated by ZIP, Bliss, HSA or Loewe model,
#'     respectively.
#'     \item \strong{ZIP_synergy, Bliss_synergy, HSA_synergy, Loewe_synergy}
#'     The synergy score calculated by ZIP, Bliss, HSA or Loewe model,
#'     respectively.
#'     \item \strong{ZIP_fit} The response fitted by ZIP model.
#'   }
#' @param statistic A character or NULL. It indicates the statistics printed
#'   in the plot while there are replicates in input data. Available values are:
#'   \itemize{
#'     \item \strong{sd} Standard deviation;
#'     \item \strong{sem} Standard error of mean;
#'     \item \strong{ci} 95\% confidence interval.
#'   }
#'   If it is \code{NULL}, no statistics will be printed.
#' @param summary_statistic A vector of characters or NULL. It indicates the
#'   summary statistics for all the \code{plot_value} in whole combination
#'   matrix. Available values are:
#'   \itemize{
#'     \item \strong{mean} Median value for all the responses or synergy
#'     scores in the matrix and the p-value if it is valid;
#'     \item \strong{median} Median value for all the responses or synergy
#'     scores in the matrix;
#'     \item \strong{quantile_90} 90\% quantile. User could change the number to
#'     print different sample quantile. For example quantile_50 equal to median. 
#'   }
#'   If it is \code{NULL}, no statistics will be printed.
#' @param plot_title A character value. It specifies the plot title. If it is
#'   \code{NULL}, the function will automatically generate a title.
#' @param dynamic A logical value. If it is \code{TRUE}, this function will
#'   use \link[plotly]{plot_ly} to generate an interactive plot. If it is
#'   \code{FALSE}, this function will use \link[lattice]{wireframe} to generate
#'   a static plot.
#' @param col_range A vector of two integers. They 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 \code{NULL} so all the concentrations are used.
#' @param row_range A vector of two integers. They 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 \code{NULL} so all the concentrations are used.
#' @param color_range A vector of two numeric values. They specify the range
#'   of the color bars. The first item (lower bounder) must be less than the
#'   second one (upper bounder). The plotted values larger than defined upper
#'   bounder will be filled in color \code{high_value_color}. The plotted values
#'   less than defined lower bounder will be filled in color
#'   \code{low_value_color}. If the defined range includes 0, value 0 will be
#'   filled in color "white". By default, it is set as \code{NULL} which
#'   means the function will automatically set the color range according to
#'   the plotted values.
#' @param high_value_color An R color value. It indicates the color for the
#'   high values.
#' @param low_value_color An R color value. It indicates the color for low
#'   values.
#' @param text_label_size_scale A numeric value. It is used to control the size
#'   text labels on the heatmap. The text size will multiply by this scale
#'   factor.
#' @param text_label_color NULL or an R color value. It indicates the color for
#'   text labels on the heatmap. If it is \code{NULL}, no text will be shown.
#' @param title_text_size_scale A numeric value. It is used to control the size
#'   of legend title, legend text, plot title, axis title, axis tick, subtitle.
#'   All the text size will multiply by this scale factor.
#'   
#' @return A ggplot plot object.
#'
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#' 
#' @export
#'
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' Plot2DrugHeatmap(data)
Plot2DrugHeatmap <- function(data,
                             plot_block = 1,
                             drugs = c(1, 2),
                             plot_value = "response",
                             statistic = NULL,
                             summary_statistic = NULL,
                             dynamic = FALSE,
                             plot_title = NULL,
                             col_range = NULL,
                             row_range = NULL,
                             color_range = NULL,
                             high_value_color = "#FF0000",
                             low_value_color = "#00FF00",
                             text_label_size_scale = 1,
                             text_label_color = "#000000",
                             title_text_size_scale = 1
                             ) {
  # Extract plot data
  plot_data <- .Extract2DrugPlotData(
    data = data,
    plot_block = plot_block,
    drugs = drugs,
    plot_value = plot_value,
    statistic = statistic
  )
  plot_table <- plot_data$plot_table
  drug_pair <- plot_data$drug_pair
  
  # Subset plot matrix
  if (!is.null(row_range)) {
    selected_rows <- levels(plot_table$conc2)[row_range[1]:row_range[2]]
    plot_table <- plot_table[plot_table$conc2 %in% selected_rows, ]
    plot_table$conc2 <- factor(plot_table$conc2)
  }
  
  if (!is.null(col_range)) {
    selected_cols <- levels(plot_table$conc1)[col_range[1]:col_range[2]]
    plot_table <- plot_table[plot_table$conc1 %in% selected_cols, ]
    plot_table$conc1 <- factor(plot_table$conc1)
  }
  
  # Generate plot title and legend title
  if (plot_value == "response") {
    if (is.null(plot_title)){
      plot_title <- paste(
        "Dose Response Matrix",
        sep = " "
      )
    }
    legend_title <- "Inhibition (%)"
  } else if (plot_value == "response_origin") {
    if (is.null(plot_title)){
      plot_title <- paste(
        "Dose Response Matrix",
        sep = " "
      )
    }
    legend_title <- paste(stringr::str_to_title(drug_pair$input_type), "%")
  } else {
    if (is.null(plot_title)){
      plot_title <- switch(
        sub(".*_", "", plot_value),
        "ref" = sub("_ref", " Reference Additive Effect", plot_value),
        "fit" = sub("_fit", " Fitted Effect", plot_value),
        "synergy" = sub("_synergy", " Synergy Score", plot_value)
      )
    }
    legend_title <- switch(
      sub(".*_", "", plot_value),
      "ref" = "Inhibition (%)",
      "fit" = "Inhibition (%)",
      "synergy" = "Synergy Score"
    )
  }
  # plot subtitle (summary statistics)
  plot_subtitle <- c()
  if (!is.null(summary_statistic)) {
    concs <- plot_table[, grepl("conc\\d+", colnames(plot_table))]
    if (endsWith(plot_value, "_synergy")) {
      concs_zero <- apply(concs, 2, function(x){x == 0})
      index <- rowSums(concs_zero) < 1
      summary_value_table <- plot_table[index, ]
    } else {
      summary_value_table <- plot_table
    }
    avail_value <- grepl("mean|median|quantile_\\d+", summary_statistic)
    if (!any(avail_value)){
      warning("Input value for parameter summary_statistic is not available.")
      plot_subtitle <- ""
    } else {
      if ("mean" %in% summary_statistic) {
        value <- .RoundValues(mean(summary_value_table$value))
        if (length(concs == 2) & 
            (drug_pair$replicate | 
             !plot_value %in% c("response", "response_origin"))) {# & drug_pair$replicate) {
          p_value <- data$drug_pairs[data$drug_pairs$block_id == plot_block,
                                     paste0(plot_value, "_p_value")]
          if (p_value != "< 2e-324") {
            p_value <- paste0("= ", p_value)
          }
          plot_subtitle <- c(
            plot_subtitle, 
            paste0(
              "Mean: ",
              value,
              " (p ",
              p_value,
              ")"
            )
          )
        } else {
          plot_subtitle <- c(plot_subtitle, paste0("Mean: ", value))
        }
      }
      if ("median" %in% summary_statistic) {
        value <- .RoundValues(stats::median(summary_value_table$value))
        plot_subtitle <-  c(plot_subtitle, paste0("Median: ", value))
      }
      qua <- grep("quantile_\\d+", summary_statistic, value = TRUE)
      if (length(qua) > 0) {
        for (q in qua) {
          pro <- as.numeric(sub("quantile_", "", q))
          value <- .RoundValues(stats::quantile(summary_value_table$value, 
                                                probs = pro / 100))
          plot_subtitle <-  c(plot_subtitle, paste0(pro, "% Quantile: ", value))
        }
      }
    }
  }
  plot_subtitle <- paste(plot_subtitle, collapse = " | ")
  
  # Color range
  if (is.null(color_range)){
    color_range <- round(max(abs(plot_table$value)), -1) + 10
    start_point <- -color_range
    end_point <- color_range
  } else {
    if (length(color_range) != 2 | class(color_range) != "numeric"){
      stop(
        "The variable 'color_range' should be a vector with exact 2 numeric ",
        "values."
      )
    } else if (color_range[1] >= color_range[2]){
      stop(
        "The first item in 'color_range' vector should be less than the ",
        "second item."
      )
    } else {
      if (color_range[1] > max(plot_table$value) |
          color_range[2] < min(plot_table$value)){
        stop(
          "There is no overlap between 'color_range' (",
          paste(color_range, collapse = ", "),
          ") and the range of 'plot_value' (",
          paste(range(plot_table$value), collapse = ", "), ")")
      }
      start_point <- color_range[1]
      end_point <- color_range[2]
    }
  }
  
  # plot heatmap for dose-response matrix
  if (dynamic){
    conc1 <- unique(plot_table$conc1)
    conc2 <- unique(plot_table$conc2)
    mat <- reshape2::acast(plot_table, conc1~conc2, value.var = "value")
    colnames(mat) <- seq(1, ncol(mat))
    rownames(mat) <- seq(1, nrow(mat))
    x_ticks_text <- as.character(sort(conc1))
    y_ticks_text <- as.character(sort(conc2))
    x <- data.frame(x = seq(1, length(conc1)),
                    ticks = as.character(sort(conc1)))
    y <- data.frame(y = seq(1, length(conc2)),
                    ticks = as.character(sort(conc2)))
    plot_table <- plot_table %>% 
      dplyr::left_join(x, by = c("conc1" = "ticks")) %>% 
      dplyr::left_join(y, by = c("conc2" = "ticks")) 
    x_axis_title <- paste0(drug_pair$drug1, " (", drug_pair$conc_unit1, ")")
    y_axis_title <- paste0(drug_pair$drug2, " (", drug_pair$conc_unit2, ")")
    
    # Hover text
    concs <- expand.grid(conc1, conc2)
    hover_text <- NULL
    for (i in 1:nrow(concs)) {
      hover_text <- c(
        hover_text,
        paste0(
          drug_pair[, c("drug1", "drug2")],
          ": ",
          sapply(concs[i, ], as.character),
          " ",
          drug_pair[, c("conc_unit1", "conc_unit1")],
          "<br>",
          collapse = ""
        )
      )
    }
    hover_text <- paste0(
      hover_text,
      "Value: ",
      .RoundValues(plot_table$value),
      sep = ""
    )
    hover_text <- matrix(hover_text, nrow = length(conc1))
    
    # Color scale
    if (start_point < 0 & end_point > 0){
      zero_pos <- -start_point/(end_point - start_point)
      color_scale <- list(
        c(0, low_value_color),
        c(zero_pos, "white"),
        c(1, high_value_color)
      )
    } else {
      color_scale <- list(
        c(0, low_value_color),
        c(1, high_value_color)
      )
    }
    
    p <- plotly::plot_ly(
      x = ~plot_table$x,
      y = ~plot_table$y,
      z = ~plot_table$value,
      zmin = start_point,
      zmid = 0,
      zmax = end_point,
      type = "heatmap",
      # text = hover_text,
      hoverinfo = "",
      autocolorscale = FALSE,
      colorscale = color_scale,
      colorbar = list(
        x = 1,
        y = 0.75,
        align = "center",
        outlinecolor = "#FFFFFF",
        tickcolor = "#FFFFFF",
        title = legend_title,
        titlefont = list(size = 12 * title_text_size_scale, family = "arial"),
        tickfont = list(size = 12 * title_text_size_scale, family = "arial")
      )
      ) %>% 
      plotly::layout(
        title = list(
          text = paste0("<b>", plot_title, "</b>"),
          font = list(size = 18 * title_text_size_scale, family = "arial"),
          y = 1.3
        ),
        xaxis = list(
          title = paste0(x_axis_title),
          tickfont = list(size = 12 * title_text_size_scale, family = "arial"),
          titlefont = list(size = 12 * title_text_size_scale, family = "arial"),
          ticks = "none",
          showspikes = FALSE,
          showgrid = FALSE,
          tickmode = "array", 
          tickvals = x$x,
          ticktext = x$ticks
        ),
        yaxis = list(
          title = paste0(y_axis_title),
          tickfont = list(size = 12 * title_text_size_scale, family = "arial"),
          titlefont = list(size = 12 * title_text_size_scale, family = "arial"),
          ticks = "none",
          showspikes = FALSE,
          showgrid = FALSE,
          tickmode = "array", 
          tickvals = y$y,
          ticktext = y$ticks
        ),
        margin = list(
          l = 50,
          r = 50,
          b = 50,
          t = 90,
          pad = 4
        )
      ) %>%
      plotly::add_annotations(
        text = plot_subtitle,
        x = 0.3,
        y = 1.05,
        yref = "paper",
        xref = "paper",
        xanchor = "middle",
        yanchor = "top",
        showarrow = FALSE,
        font = list(size = 15 * title_text_size_scale, family = "arial")
      )
      
    if (!is.null(text_label_color)) {
      p <- p %>% 
        plotly::add_annotations(
          x = ~plot_table$x,
          y = ~plot_table$y,
          text = ~plot_table$text,
          showarrow = FALSE,
          font = list(
            color = text_label_color,
            size = 12 * text_label_size_scale
          ),
          ax = 20,
          ay = -20
        ) 
    }
    p <- p %>%
      plotly::config(
        toImageButtonOptions = list(
          format = "svg",
          filename = plot_title,
          width = NULL,# 500,
          height = NULL,# 500,
          scale = 1
        )
      ) 
  } else{
    warn = getOption("warn")
    options(warn=-1)
    p <- ggplot2::ggplot(
      data = plot_table,
      aes(x = conc1, y = conc2, fill = value)
    ) +
      ggplot2::geom_tile() +
      ggplot2::geom_text(
        ggplot2::aes(label = text),
        size = .Pt2mm(7) * text_label_size_scale
      )
    
    # Manage color bar
    if (start_point < 0 & end_point > 0){
      colour_breaks <- c(start_point, 0, end_point)
      colours <- c(low_value_color, "#FFFFFF", high_value_color)
      p <- p +
        scale_fill_gradientn(
          limits  = c(start_point, end_point),#range(plot_table$value),
          colours = colours[c(1, seq_along(colours), length(colours))],
          values  = c(
            0,
            scales::rescale(colour_breaks, from = c(start_point, end_point)), #range(plot_table$value)),
            1
          ),
          oob = scales::oob_squish_any
        )
      
    } else {
      p <- p +
        ggplot2::scale_fill_gradient(
          high= high_value_color,
          low = low_value_color,
          name = legend_title,
          limits = c(start_point, end_point),
          oob = scales::oob_squish_any
        )
    }
    
    p <- p +
      ggplot2::guides(
        fill = ggplot2::guide_colorbar(
          barheight = 10,
          barwidth = 1.5,
          ticks = FALSE
        )
      ) +
      # Add the title for heatmap
      ggplot2::labs(
        title = plot_title,
        subtitle = plot_subtitle,
        x = paste0(drug_pair$drug1, " (", drug_pair$conc_unit1, ")"),
        y = paste0(drug_pair$drug2, " (", drug_pair$conc_unit2, ")"),
        fill = legend_title
      ) +
      ggplot2::theme(
        plot.title = ggplot2::element_text(
          size = 13.5 * title_text_size_scale,
          face = "bold",
          hjust = 0.5
        ),
        plot.subtitle = ggplot2::element_text(
          size = 12 * title_text_size_scale,
          hjust = 0.5
        ),
        panel.background = ggplot2::element_blank(),
        # Set label's style of heatmap
        axis.text = ggplot2::element_text(
          size = 10 * title_text_size_scale
        ),
        axis.title = ggplot2::element_text(
          size = 10 * title_text_size_scale
        ),
        legend.title = ggplot2::element_text(
          size = 10 * title_text_size_scale
        ),
        legend.text = ggplot2::element_text(
          size = 10 * title_text_size_scale
        ),
        legend.background = element_rect(color = NA)
      )
  }
  return(p)
  options(warn=warn)
}

#' 2D Contour Plot for 2-drug Combination Dose-Response/Synergy Scores
#' 
#' This function will generate a contour level plot for 2-drug combinations.
#' The axes are the dosage for each drug. The values could be observed response,
#' synergy scores or the reference effects calculated from different models.
#'
#' @param data A list object generated by function \code{\link{ReshapeData}}.
#' @param plot_block A character/integer. It indicates the block ID for the
#'   block to visualize.
#' @param drugs A vector of characters or integers with length of 2. It contains
#'   the index for two drugs to plot. For example, \code{c(1, 2)} indicates to
#'   plot "drug1" and "drug2" in the input \code{data}.
#' @param plot_value A character value. It indicates the score or response value
#'   to be visualized. If the \code{data} is the direct output from
#'   \link{ReshapeData}, the available values for this parameter are:
#'   \itemize{
#'     \item \strong{response_origin} The original response value in input data.
#'     It might be \% inhibition or \% viability.
#'     \item \strong{response} The \% inhibition after preprocess by function 
#'     \link{ReshapeData}
#'   }
#'   If the \code{data} is the output from \link{CalculateSynergy}, following
#'   values are also available:
#'   \itemize{
#'     \item \strong{ZIP_ref, Bliss_ref, HSA_ref, Loewe_ref} The reference
#'     additive effects calculated by ZIP, Bliss, HSA or Loewe model,
#'     respectively.
#'     \item \strong{ZIP_synergy, Bliss_synergy, HSA_synergy, Loewe_synergy}
#'     The synergy score calculated by ZIP, Bliss, HSA or Loewe model,
#'     respectively.
#'     \item \strong{ZIP_fit} The response fitted by ZIP model.
#'   }
#' @param summary_statistic A vector of characters or NULL. It indicates the
#'   summary statistics for all the \code{plot_value} in whole combination
#'   matrix. Available values are:
#'   \itemize{
#'     \item \strong{mean} Median value for all the responses or synergy
#'     scores in the matrix and the p-value if it is valid;
#'     \item \strong{median} Median value for all the responses or synergy
#'     scores in the matrix;
#'     \item \strong{quantile_90} 90\% quantile. User could change the number to
#'     print different sample quantile. For example quantile_50 equal to median. 
#'   }
#'   If it is \code{NULL}, no statistics will be printed.
#' @param plot_title A character value. It specifies the plot title. If it is
#'   \code{NULL}, the function will automatically generate a title.
#' @param axis_line A logical value. Whether to show the axis lines and ticks.
#' @param interpolate_len An integer. It specifies how many values need to be
#'    interpolated between two concentrations. It is used to control the 
#'    smoothness of the synergy surface.
#' @param col_range A vector of two integers. They 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 \code{NULL} so all the concentrations are used.
#' @param row_range A vector of two integers. They 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 \code{NULL} so all the concentrations are used.
#' @param color_range A vector of two numeric values. They specify the range
#'   of the color bars. The first item (lower bounder) must be less than the
#'   second one (upper bounder). The plotted values larger than defined upper
#'   bounder will be filled in color \code{high_value_color}. The plotted values
#'   less than defined lower bounder will be filled in color
#'   \code{low_value_color}. If the defined range includes 0, value 0 will be
#'   filled in color "white". By default, it is set as \code{NULL} which
#'   means the function will automatically set the color range according to
#'   the plotted values.
#' @param dynamic A logical value. If it is \code{TRUE}, this function will
#'   use \link[plotly]{plot_ly} to generate an interactive plot. If it is
#'   \code{FALSE}, this function will use \link[lattice]{wireframe} to generate
#'   a static plot.
#' @param grid A logical value. It indicates whether to add grids.
#' @param high_value_color An R color value. It indicates the color for the
#'   high values.
#' @param low_value_color An R color value. It indicates the color for low
#'   values.
#' @param text_size_scale A numeric value. It is used to control the size
#'   of text in the plot. All the text size will multiply by this scale factor.
#'   
#' @return If \code{dynamic = FALSE}, this function will return a plot project 
#'   recorded by \link[grDevices]{recordPlot}. If \code{dynamic = FALSE}, this
#'   function will return a plotly plot object.
#'
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#' 
#' @export
#'
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' Plot2DrugContour(data)
Plot2DrugContour <- function(data,
                             plot_block = 1,
                             drugs = c(1, 2),
                             plot_value = "response",
                             interpolate_len = 2,
                             summary_statistic = NULL,
                             dynamic = FALSE,
                             grid = TRUE,
                             plot_title = NULL,
                             axis_line = FALSE,
                             col_range = NULL,
                             row_range = NULL,
                             color_range = NULL,
                             # width = 600,
                             # height = 400,
                             high_value_color = "#FF0000",
                             low_value_color = "#00FF00",
                             text_size_scale = 1) {
  # Extract plot data
  plot_data <- .Extract2DrugPlotData(
    data = data,
    plot_block = plot_block,
    drugs = drugs,
    plot_value = plot_value,
    statistic = NULL
  )
  plot_table <- plot_data$plot_table
  drug_pair <- plot_data$drug_pair
  
  # Subset plot matrix
  if (!is.null(row_range)) {
    selected_rows <- levels(plot_table$conc2)[row_range[1]:row_range[2]]
    plot_table <- plot_table[plot_table$conc2 %in% selected_rows, ]
    plot_table$conc2 <- factor(plot_table$conc2)
  }
  
  if (!is.null(col_range)) {
    selected_cols <- levels(plot_table$conc1)[col_range[1]:col_range[2]]
    plot_table <- plot_table[plot_table$conc1 %in% selected_cols, ]
    plot_table$conc1 <- factor(plot_table$conc1)
  }
  # Generate plot title and legend title
  if (plot_value == "response") {
    if (is.null(plot_title)){
      plot_title <- paste(
        "Dose Response Matrix",
        sep = " "
      )
    }
    legend_title <- "Inhibition (%)"
    z_axis_title <- "Response (% inhibition)"
  } else if (plot_value == "response_origin") {
    if (is.null(plot_title)){
      plot_title <- paste(
        "Dose Response Matrix",
        sep = " "
      )
    }
    legend_title <- paste0(stringr::str_to_title(drug_pair$input_type), " (%)")
    z_axis_title <- paste0("Response (% ",drug_pair$input_type,")")
  } else {
    if (is.null(plot_title)){
      plot_title <- switch(
        sub(".*_", "", plot_value),
        "ref" = sub("_ref", " Reference Additive Effect", plot_value),
        "fit" = sub("_fit", " Fitted Effect", plot_value),
        "synergy" = sub("_synergy", " Synergy Score", plot_value)
      )
    }
    legend_title <- switch(
      sub(".*_", "", plot_value),
      "ref" = "Inhibition (%)",
      "fit" = "Inhibition (%)",
      "synergy" = "Synergy Score"
    )
    z_axis_title <- switch(
      sub(".*_", "", plot_value),
      "ref" = "Response (% inhibition)",
      "fit" = "Response (% inhibition)",
      "synergy" = "Synergy Score"
    )
  }
  # plot subtitle (summary statistics)
  plot_subtitle <- c()
  if (!is.null(summary_statistic)) {
    concs <- plot_table[, grepl("conc\\d+", colnames(plot_table))]
    if (endsWith(plot_value, "_synergy")) {
      concs_zero <- apply(concs, 2, function(x){x == 0})
      index <- rowSums(concs_zero) < 1
      summary_value_table <- plot_table[index, ]
    } else {
      summary_value_table <- plot_table
    }
    avail_value <- grepl("mean|median|quantile_\\d+", summary_statistic)
    if ("mean" %in% summary_statistic) {
      value <- .RoundValues(mean(summary_value_table$value))
      if (length(concs == 2) & 
          (drug_pair$replicate | 
           !plot_value %in% c("response", "response_origin"))) {# & drug_pair$replicate) {
        p_value <- data$drug_pairs[data$drug_pairs$block_id == plot_block,
                                   paste0(plot_value, "_p_value")]
        if (p_value != "< 2e-324") {
          p_value <- paste0("= ", p_value)
        }
        plot_subtitle <- c(
          plot_subtitle, 
          paste0(
            "Mean: ",
            value,
            " (p ",
            p_value,
            ")"
          )
        )
      } else {
        plot_subtitle <- c(plot_subtitle, paste0("Mean: ", value))
      }
    }
    if ("median" %in% summary_statistic) {
      value <- .RoundValues(stats::median(summary_value_table$value))
      plot_subtitle <-  c(plot_subtitle, paste0("Median: ", value))
    }
    qua <- grep("quantile_\\d+", summary_statistic, value = TRUE)
    if (length(qua) > 0) {
      for (q in qua) {
        pro <- as.numeric(sub("quantile_", "", q))
        value <- .RoundValues(stats::quantile(summary_value_table$value, 
                                              probs = pro / 100))
        plot_subtitle <-  c(plot_subtitle, paste0(pro, "% Quantile: ", value))
      }
    }
  }
  plot_subtitle <- paste(plot_subtitle, collapse = " | ")
  
  conc1 <- unique(plot_table$conc1)
  conc2 <- unique(plot_table$conc2)
  
  # Interpolation by kriging
  mat <- reshape2::acast(plot_table, conc1~conc2, value.var = "value")
  extended_mat <- .ExtendedScores(mat, len = interpolate_len)
  colnames(extended_mat) <- seq(1, ncol(extended_mat))
  rownames(extended_mat) <- seq(1, nrow(extended_mat))
  x_ticks <- seq(1, nrow(extended_mat), by = interpolate_len + 1)
  y_ticks <- seq(1, ncol(extended_mat), by = interpolate_len + 1)
  x_ticks_text <- as.character(sort(conc1))
  y_ticks_text <- as.character(sort(conc2))
  x_axis_title <- paste0(drug_pair$drug1, " (", drug_pair$conc_unit1, ")")
  y_axis_title <- paste0(drug_pair$drug2, " (", drug_pair$conc_unit2, ")")
  
  # Color range
  if (is.null(color_range)){
    color_range <- round(max(abs(plot_table$value)), -1) + 10
    start_point <- -color_range
    end_point <- color_range
  } else {
    if (length(color_range) != 2 | class(color_range) != "numeric"){
      stop(
        "The variable 'color_range' should be a vector with exact 2 numeric ",
        "values."
      )
    } else if (color_range[1] >= color_range[2]){
      stop(
        "The first item in 'color_range' vector should be less than the ",
        "second item."
      )
    } else {
      if (color_range[1] > max(plot_table$value) |
          color_range[2] < min(plot_table$value)){
        stop(
          "There is no overlap between 'color_range' (",
          paste(color_range, collapse = ", "),
          ") and the range of 'plot_value' (",
          paste(range(plot_table$value), collapse = ", "), ")")
      }
      start_point <- color_range[1]
      end_point <- color_range[2]
    }
  }
  
  if (dynamic){
    y <- seq(1, ncol(extended_mat))
    x <- seq(1, nrow(extended_mat))

    # Hover text
    concs <- expand.grid(conc1, conc2)
    hover_text <- NULL
    for (i in 1:nrow(concs)) {
      hover_text <- c(
        hover_text,
        paste0(
          drug_pair[, c("drug1", "drug2")],
          ": ",
          sapply(concs[i, ], as.character),
          " ",
          drug_pair[, c("conc_unit1", "conc_unit1")],
          "<br>",
          collapse = ""
        )
      )
    }
    hover_text <- paste0(
      hover_text,
      "Value: ",
      .RoundValues(plot_table$value),
      sep = ""
    )
    hover_text <- matrix(hover_text, nrow = length(conc1))
    hover_text_mat <- matrix(rep(NA, length(extended_mat)),
                             nrow = nrow(extended_mat))
    for (i in 1:length(x_ticks)) {
      for (j in 1:length(y_ticks)) {
        hover_text_mat[x_ticks[i], y_ticks[j]] <- hover_text[i, j]
      }
    }
    
    # Color scale
    if (start_point < 0 & end_point > 0){
      zero_pos <- -start_point/(end_point - start_point)
      color_scale <- list(
        c(0, low_value_color),
        c(zero_pos, "white"),
        c(1, high_value_color)
      )
    } else {
      color_scale <- list(
        c(0, low_value_color),
        c(1, high_value_color)
      )
    }
    
    p <- plotly::plot_ly(
      x = x,
      y = y,
      z = t(extended_mat),
      zmin = start_point,
      zmax = end_point,
      type = "contour",
      line = list(color = "white", smoothing = 0),
      text = t(hover_text_mat),
      hoverinfo = "text",
      # autocolorscale = FALSE,
      colorscale = color_scale,
      autocontour = F,
      # ncontours = 21,
      colorbar = list(
        x = 1,
        y = 0.75,
        title = legend_title,
        titlefont = list(size = 16 * text_size_scale, family = "arial"),
        showticklabels = TRUE,
        align = "center",
        outlinewidth = 0,
        ticks = "inside",
        tickcolor = "#FFFFFFF",
        tickfont = list(size = 12 * text_size_scale, family = "arial")
      ),
      contours = list(
        coloring = "heatmap",
        start = start_point,
        end = end_point,
        size = 10,
        showlabels = FALSE,
        x = list(
          # highlight = FALSE,
          highlightwidth = 2,
          show = FALSE,
          color = 'black',
          width = 5,
          start = 1,
          end = max(x),
          size = 1 * text_size_scale
        ),
        y = list(
          # highlight = FALSE,
          highlightwidth = 2,
          show = FALSE,
          color = 'black',
          width = 5,
          start = 1,
          end = max(y),
          size = 1 * text_size_scale
        )
      )
      ) %>% 
      plotly::add_annotations(
        text = plot_subtitle,
        x = 0.3,
        y = 1.05,
        yref = "paper",
        xref = "paper",
        xanchor = "middle",
        yanchor = "top",
        showarrow = FALSE,
        font = list(size = 15 * text_size_scale)
      ) %>% 
      plotly::layout(
        title = list(
          text = paste0("<b>", plot_title, "</b>"),
          font = list(size = 18 * text_size_scale, family = "arial"),
          y = 1.3
        ),
        xaxis = list(
          title = paste0(x_axis_title),
          tickfont = list(size = 12 * text_size_scale, family = "arial"),
          titlefont = list(size = 16 * text_size_scale, family = "arial"),
          ticks = ifelse(axis_line, "outside", "none"),
          showline = axis_line,
          showspikes = FALSE,
          tickmode = "array", 
          tickvals = x_ticks,
          ticktext = x_ticks_text,
          range = c(min(x_ticks), max(x_ticks))
        ),
        yaxis = list(
          title = paste0(y_axis_title),
          tickfont = list(size = 12 * text_size_scale, family = "arial"),
          titlefont = list(size = 16 * text_size_scale, family = "arial"),
          ticks = ifelse(axis_line, "outside", "none"),
          showline = axis_line,
          showspikes = FALSE,
          tickmode = "array", 
          tickvals = y_ticks,
          ticktext = y_ticks_text,
          range = c(min(y_ticks), max(y_ticks))
        ),
        margin = list(
          l = 50,
          r = 50,
          b = 50,
          t = 90,
          pad = 0
        )
      ) %>% 
      plotly::config(
        toImageButtonOptions = list(
          format = "svg",
          filename = plot_title,
          width = NULL,# 500,
          height = NULL,# 500,,
          scale = 1
        )
      )
    if (grid) {
      for(x in x_ticks){
        p <- plotly::add_segments(
          p,
          x = x,
          xend = x,
          y = min(y_ticks),
          yend = max(y_ticks),
          line = list(color = "black", dash = "dot", width = 1),
          inherit = FALSE,
          showlegend = FALSE
        )
      }
      for(y in y_ticks){
        p <- plotly::add_segments(
          p,
          y = y,
          yend = y,
          x = min(x_ticks),
          xend = max(x_ticks),
          line = list(color = "black", dash = "dot", width = 1),
          inherit = FALSE,
          showlegend = FALSE
        )
      }
    }
  } else{
    warn = getOption("warn")
    options(warn=-1)
    extended_df <- reshape2::melt(extended_mat)
    colnames(extended_df) <- c("x", "y", "value")
    p <- ggplot2::ggplot(
        extended_df,
        ggplot2::aes(x = x, y = y, z = value),
        na.rm = TRUE
      ) +
      metR::geom_contour_fill(
        ggplot2::aes(
          x = x, y = y, z = value),
        color = "#FFFFFF50",
        binwidth = 10
      ) +
      ggplot2::scale_x_continuous(
        expand = c(0, 0),
        breaks = x_ticks,
        labels = x_ticks_text,
        limits = c(min(x_ticks), max(x_ticks))
      ) +
      ggplot2::scale_y_continuous(
        expand = c(0, 0),
        breaks = y_ticks,
        labels = y_ticks_text,
        limits = c(min(y_ticks), max(y_ticks))
      )
      
    # Manage color bar
    if (start_point < 0 & end_point > 0){
      colour_breaks <- c(start_point, 0, end_point)
      colours <- c(low_value_color, "#FFFFFF", high_value_color)
      p <- p +
        scale_fill_gradientn(
          limits  = c(start_point, end_point), # range(plot_table$value),
          colours = colours[c(1, seq_along(colours), length(colours))],
          values  = c(
            0,
            scales::rescale(colour_breaks, from = c(start_point, end_point)), #range(plot_table$value)),
            1
          ),
          oob = scales::oob_squish_any
        )
      
    } else {
      p <- p +
        ggplot2::scale_fill_gradient(
          high= high_value_color,
          low = low_value_color,
          name = legend_title,
          limits = c(start_point, end_point),
          oob = scales::oob_squish_any
        )
    }
    
    p <- p +
    ggplot2::guides(
      fill = ggplot2::guide_colorbar(
        title = legend_title,
        barheight = 10,
        barwidth = 1.5,
        ticks = FALSE
      )
    ) +
    ggplot2::xlab(x_axis_title) +
    ggplot2::ylab(y_axis_title) +
    ggplot2::ggtitle(
      label = paste0(
        plot_title
      ),
      subtitle = plot_subtitle
    ) +
    ggplot2::theme(
      plot.title = ggplot2::element_text(
        size = 13.5 * text_size_scale,
        face = "bold",
        hjust = 0.5
      ),
      plot.subtitle = ggplot2::element_text(
        size = 12 * text_size_scale,
        hjust = 0.5
      ),
      panel.background = ggplot2::element_blank(),
      # Set label's style of heatmap
      axis.text = ggplot2::element_text(
        size = 10 * text_size_scale
      ),
      axis.title = ggplot2::element_text(
        size = 10 * text_size_scale
      ),
      legend.title = ggplot2::element_text(
        size = 10 * text_size_scale
      ),
      legend.text = ggplot2::element_text(
        size = 10 * text_size_scale
      ),
      legend.background = element_rect(color = NA)
    )
    
    if (axis_line){
      p <- p +
        ggplot2::theme(
          axis.ticks = ggplot2::element_line(),
          axis.line = ggplot2::element_line()
        )
    } else {
      p <- p +
        ggplot2::theme(
          axis.ticks = ggplot2::element_blank(),
          axis.line = ggplot2::element_blank()
        )
    }
    if (grid) {
      p <- p +
        geom_vline(
          xintercept = x_ticks,
          linetype = "dotted"
        ) +
        geom_hline(
          yintercept = y_ticks,
          linetype = "dotted"
        )
    }
  }
  return(p)
  options(warn=warn)
}

#' 3D Surface Plot for 2-drug Combination Dose-Response/Synergy Scores
#' 
#' This function will generate a surface plot for 2-drug combinations. The axes
#' are the dosage for each drug. The values could be observed response,
#' synergy scores or the reference effects calculated from different models.
#'
#' @param data A list object generated by function \code{\link{ReshapeData}}.
#' @param plot_block A character/integer. It indicates the block ID for the
#'   block to visualize.
#' @param drugs A vector of characters or integers with length of 2. It contains
#'   the index for two drugs to plot. For example, \code{c(1, 2)} indicates to
#'   plot "drug1" and "drug2" in the input \code{data}.
#' @param plot_value A character value. It indicates the score or response value
#'   to be visualized. If the \code{data} is the direct output from
#'   \link{ReshapeData}, the available values for this parameter are:
#'   \itemize{
#'     \item \strong{response_origin} The original response value in input data.
#'     It might be \% inhibition or \% viability.
#'     \item \strong{response} The \% inhibition after preprocess by function 
#'     \link{ReshapeData}
#'   }
#'   If the \code{data} is the output from \link{CalculateSynergy}, following
#'   values are also available:
#'   \itemize{
#'     \item \strong{ZIP_ref, Bliss_ref, HSA_ref, Loewe_ref} The reference
#'     additive effects calculated by ZIP, Bliss, HSA or Loewe model,
#'     respectively.
#'     \item \strong{ZIP_synergy, Bliss_synergy, HSA_synergy, Loewe_synergy}
#'     The synergy score calculated by ZIP, Bliss, HSA or Loewe model,
#'     respectively.
#'     \item \strong{ZIP_fit} The response fitted by ZIP model.
#'   }
#' @param summary_statistic A vector of characters or NULL. It indicates the
#'   summary statistics for all the \code{plot_value} in whole combination
#'   matrix. Available values are:
#'   \itemize{
#'     \item \strong{mean} Median value for all the responses or synergy
#'     scores in the matrix and the p-value if it is valid;
#'     \item \strong{median} Median value for all the responses or synergy
#'     scores in the matrix;
#'     \item \strong{quantile_90} 90\% quantile. User could change the number to
#'     print different sample quantile. For example quantile_50 equal to median. 
#'   }
#'   If it is \code{NULL}, no statistics will be printed.
#' @param plot_title A character value. It specifies the plot title. If it is
#'   \code{NULL}, the function will automatically generate a title.
#' @param interpolate_len An integer. It specifies how many values need to be
#'   interpolated between two concentrations. It is used to control the 
#'   smoothness of the synergy surface.
#' @param axis_line A logical value. It specifies whether to show the axis lines
#'   and ticks in dynamic plot. Axis lines are always shown in the static plot.
#' @param col_range A vector of two integers. They 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 \code{NULL} so all the concentrations are used.
#' @param row_range A vector of two integers. They 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 \code{NULL} so all the concentrations are used.
#' @param z_range A vector of two numeric values. They specify the range of
#'   z-axis plotted.Default value is \code{NULL}. The function
#'   automatically set the range.
#' @param color_range A vector of two numeric values. They specify the range
#'   of the color bars. The first item (lower bounder) must be less than the
#'   second one (upper bounder). The plotted values larger than defined upper
#'   bounder will be filled in color \code{high_value_color}. The plotted values
#'   less than defined lower bounder will be filled in color
#'   \code{low_value_color}. If the defined range includes 0, value 0 will be
#'   filled in color "white". By default, it is set as \code{NULL} which
#'   means the function will automatically set the color range according to
#'   the plotted values.
#' @param dynamic A logical value. If it is \code{TRUE}, this function will
#'   use \link[plotly]{plot_ly} to generate an interactive plot. If it is
#'   \code{FALSE}, this function will use \link[lattice]{wireframe} to generate
#'   a static plot.
#' @param grid A logical value. It indicates whether to add grids on the 
#'   surface.
#' @param high_value_color An R color value. It indicates the color for the
#'   high values.
#' @param low_value_color An R color value. It indicates the color for low
#'   values.

#' @param text_size_scale A numeric value. It is used to control the size
#'   of text in the plot. All the text size will multiply by this scale factor.
#'
#' @return If \code{dynamic = FALSE}, this function will return a plot project 
#'   recorded by \link[grDevices]{recordPlot}. If \code{dynamic = FALSE}, this
#'   function will return a plotly plot object.
#'
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#' 
#' @export
#'
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' Plot2DrugSurface(data)
#' Plot2DrugSurface(data, dynamic = TRUE)
Plot2DrugSurface <- function(data,
                             plot_block = 1,
                             drugs = c(1, 2),
                             plot_value = "response",
                             summary_statistic = NULL,
                             plot_title = NULL,
                             interpolate_len = 2,
                             axis_line = FALSE,
                             col_range = NULL,
                             row_range = NULL,
                             z_range = NULL,
                             color_range = NULL,
                             dynamic = FALSE,
                             grid = TRUE,
                             # width = 600,
                             # height = 400,
                             high_value_color = "#FF0000",
                             low_value_color = "#00FF00",
                             text_size_scale = 1) {
  # Extract plot data
  plot_data <- .Extract2DrugPlotData(
    data = data,
    plot_block = plot_block,
    drugs = drugs,
    plot_value = plot_value,
    statistic = NULL
  )
  plot_table <- plot_data$plot_table
  drug_pair <- plot_data$drug_pair
  
  # Subset plot matrix
  if (!is.null(row_range)) {
    selected_rows <- levels(plot_table$conc2)[row_range[1]:row_range[2]]
    plot_table <- plot_table[plot_table$conc2 %in% selected_rows, ]
    plot_table$conc2 <- factor(plot_table$conc2)
  }
  
  if (!is.null(col_range)) {
    selected_cols <- levels(plot_table$conc1)[col_range[1]:col_range[2]]
    plot_table <- plot_table[plot_table$conc1 %in% selected_cols, ]
    plot_table$conc1 <- factor(plot_table$conc1)
  }
  # Generate plot title and legend title
  if (plot_value == "response") {
    if (is.null(plot_title)){
      plot_title <- paste(
        "Dose Response Matrix",
        sep = " "
      )
    }
    legend_title <- "Inhibition (%)"
    z_axis_title <- "Response (% inhibition)"
  } else if (plot_value == "response_origin") {
    if (is.null(plot_title)){
      plot_title <- paste(
        "Dose Response Matrix",
        sep = " "
      )
    }
    legend_title <- paste0(stringr::str_to_title(drug_pair$input_type), " (%)")
    z_axis_title <- paste0("Response (% ",drug_pair$input_type,")")
  } else {
    if (is.null(plot_title)){
      plot_title <- switch(
        sub(".*_", "", plot_value),
        "ref" = sub("_ref", " Reference Additive Effect", plot_value),
        "fit" = sub("_fit", " Fitted Effect", plot_value),
        "synergy" = sub("_synergy", " Synergy Score", plot_value)
      )
    }
    legend_title <- switch(
      sub(".*_", "", plot_value),
      "ref" = "Inhibition (%)",
      "fit" = "Inhibition (%)",
      "synergy" = "Synergy Score"
    )
    z_axis_title <- switch(
      sub(".*_", "", plot_value),
      "ref" = "Response (% inhibition)",
      "fit" = "Response (% inhibition)",
      "synergy" = "Synergy Score"
    )
  }
  
  # plot subtitle (summary statistics)
  plot_subtitle <- c()
  if (!is.null(summary_statistic)) {
    concs <- plot_table[, grepl("conc\\d+", colnames(plot_table))]
    if (endsWith(plot_value, "_synergy")) {
      concs_zero <- apply(concs, 2, function(x){x == 0})
      index <- rowSums(concs_zero) < 1
      summary_value_table <- plot_table[index, ]
    } else {
      summary_value_table <- plot_table
    }
    avail_value <- grepl("mean|median|quantile_\\d+", summary_statistic)
    if ("mean" %in% summary_statistic) {
      value <- .RoundValues(mean(summary_value_table$value))
      if (length(concs == 2) & 
          (drug_pair$replicate | 
           !plot_value %in% c("response", "response_origin"))) {# & drug_pair$replicate) {
        p_value <- data$drug_pairs[data$drug_pairs$block_id == plot_block,
                                   paste0(plot_value, "_p_value")]
        if (p_value != "< 2e-324") {
          p_value <- paste0("= ", p_value)
        }
        plot_subtitle <- c(
          plot_subtitle, 
          paste0(
            "Mean: ",
            value,
            " (p ",
            p_value,
            ")"
          )
        )
      } else {
        plot_subtitle <- c(plot_subtitle, paste0("Mean: ", value))
      }
    }
    if ("median" %in% summary_statistic) {
      value <- .RoundValues(stats::median(summary_value_table$value))
      plot_subtitle <-  c(plot_subtitle, paste0("Median: ", value))
    }
    qua <- grep("quantile_\\d+", summary_statistic, value = TRUE)
    if (length(qua) > 0) {
      for (q in qua) {
        pro <- as.numeric(sub("quantile_", "", q))
        value <- .RoundValues(
          stats::quantile(summary_value_table$value, probs = pro / 100)
        )
        plot_subtitle <-  c(plot_subtitle, paste0(pro, "% Quantile: ", value))
      }
    }
  }
  plot_subtitle <- paste(plot_subtitle, collapse = " | ")
  
  conc1 <- unique(plot_table$conc1)
  conc2 <- unique(plot_table$conc2)
  
  # kriging with kriging from  (which )
  mat <- reshape2::acast(plot_table, conc1~conc2, value.var = "value")
  extended_mat <- .ExtendedScores(mat, len = interpolate_len)
  colnames(extended_mat) <- seq(1, ncol(extended_mat))
  rownames(extended_mat) <- seq(1, nrow(extended_mat))
  x_ticks <- seq(1, nrow(extended_mat), by = interpolate_len + 1)
  y_ticks <- seq(1, ncol(extended_mat), by = interpolate_len + 1)
  x_ticks_text <- as.character(sort(conc1))
  y_ticks_text <- as.character(sort(conc2))
  x_axis_title <- paste0(drug_pair$drug1, " (", drug_pair$conc_unit1, ")")
  y_axis_title <- paste0(drug_pair$drug2, " (", drug_pair$conc_unit2, ")")
  
  # Color range
  if (is.null(color_range)){
    color_range <- round(max(abs(plot_table$value)), -1) + 10
    start_point <- -color_range
    end_point <- color_range
  } else {
    if (length(color_range) != 2 | class(color_range) != "numeric"){
      stop(
        "The variable 'color_range' should be a vector with exact 2 numeric ",
        "values."
      )
    } else if (color_range[1] >= color_range[2]){
      stop(
        "The first item in 'color_range' vector should be less than the ",
        "second item."
      )
    } else {
      if (color_range[1] > max(plot_table$value) |
          color_range[2] < min(plot_table$value)){
        stop(
          "There is no overlap between 'color_range' (",
          paste(color_range, collapse = ", "),
          ") and the range of 'plot_value' (",
          paste(range(plot_table$value), collapse = ", "), ")")
      }
      start_point <- color_range[1]
      end_point <- color_range[2]
    }
  }
  
  # z axis range
  if (is.null(z_range)){
    z_range <- round(max(abs(plot_table$value)), -1) + 10
    z_range <- c(-z_range, z_range)
  } else {
    if (length(z_range) != 2 | class(z_range) != "numeric"){
      stop(
        "The variable 'z_range' should be a vector with exact 2 numeric ",
        "values."
      )
    } else if (z_range[1] >= z_range[2]){
      stop(
        "The first item in 'z_range' vector should be less than the ",
        "second item."
      )
    } else {
      if (z_range[1] > max(plot_table$value) |
          z_range[2] < min(plot_table$value)){
        stop(
          "There is no overlap between 'color_range' (",
          paste(z_range, collapse = ", "),
          ") and the range of 'plot_value' (",
          paste(range(plot_table$value), collapse = ", "), ")")
      }
    }
  }
  
  if (dynamic) {
    y <- seq(1, ncol(extended_mat))
    x <- seq(1, nrow(extended_mat))
    
    # Hover text
    concs <- expand.grid(conc1, conc2)
    hover_text <- NULL
    for (i in 1:nrow(concs)) {
      hover_text <- c(
        hover_text,
        paste0(
          drug_pair[, c("drug1", "drug2")],
          ": ",
          sapply(concs[i, ], as.character),
          " ",
          drug_pair[, c("conc_unit1", "conc_unit1")],
          "<br>",
          collapse = ""
        )
      )
    }
    hover_text <- paste0(
      hover_text,
      "Value: ",
      .RoundValues(plot_table$value),
      sep = ""
    )
    hover_text <- matrix(hover_text, nrow = length(conc1))
    hover_text_mat <- matrix(rep(NA, length(extended_mat)),
                             nrow = nrow(extended_mat))
    for (i in 1:length(x_ticks)) {
      for (j in 1:length(y_ticks)) {
        hover_text_mat[x_ticks[i], y_ticks[j]] <- hover_text[i, j]
      }
    }
    
    # Color scale
    if (start_point < 0 & end_point > 0){
      zero_pos <- -start_point/(end_point - start_point)
      color_scale <- list(
        c(0, low_value_color),
        c(zero_pos, "white"),
        c(1, high_value_color)
      )
    } else {
      color_scale <- list(
        c(0, low_value_color),
        c(1, high_value_color)
      )
    }
    
    # z-axis setting
    zaxis_setting <- list(
      title = paste0(z_axis_title),
      tickfont = list(size = 12 * text_size_scale, family = "arial"),
      titlefont = list(size = 12 * text_size_scale, family = "arial"),
      ticks = ifelse(axis_line, "outside", "none"),
      showline = axis_line,
      ticks = "none",
      tickmode = "array",
      showspikes = FALSE
    )
    
    if (!is.null(z_range)){
      zaxis_setting$range <- z_range
    }
    
    p <- plotly::plot_ly() %>% 
      plotly::add_surface(
        name = "surface",
        x = ~x,
        y = ~y,
        z = t(extended_mat),
        text = t(hover_text_mat),
        hoverinfo = "text",
        colorscale = color_scale,
        cauto = FALSE,
        colorbar = list(
          x = 1,
          y = 0.75,
          align = "center",
          outlinecolor = "#FFFFFF",
          tickcolor = "#FFFFFF",
          title = legend_title,
          titlefont = list(size = 12 * text_size_scale, family = "arial"),
          tickfont = list(size = 12 * text_size_scale, family = "arial")
        ),
        cmin = start_point,
        cmax = end_point,
        contours = list(
          x = list(
            # highlight = FALSE,
            show = grid,
            color = 'black',
            width = 1,
            start = 1,
            end = max(x),
            size = 1
          ),
          y = list(
            # highlight = FALSE,
            show = grid,
            color = 'black',
            width = 1,
            start = 1,
            end = max(y),
            size = 1
          )#,
          # z = list(highlight = FALSE)
        )
      ) %>% 
      plotly::add_annotations(
        text = plot_subtitle,
        x = 0.3,
        y = 1.05,
        yref = "paper",
        xref = "paper",
        xanchor = "middle",
        yanchor = "top",
        showarrow = FALSE,
        font = list(size = 15 * text_size_scale)
      ) %>% 
      plotly::layout(
        title = list(
          text = paste0("<b>", plot_title, "</b>"),
          font = list(size = 18 * text_size_scale, family = "arial"),
          y = 1.3
        ),
        scene = list(
          aspectratio = list(x=1, y=1, z=1),
          xaxis = list(
            title = paste0(x_axis_title),
            tickfont = list(size = 12 * text_size_scale, family = "arial"),
            titlefont = list(size = 12 * text_size_scale, family = "arial"),
            ticks = ifelse(axis_line, "outside", "none"),
            showline = axis_line,
            showspikes = FALSE,
            tickmode = "array", 
            tickvals = x_ticks,
            ticktext = x_ticks_text
          ),
          yaxis = list(
            title = paste0(y_axis_title),
            tickfont = list(size = 12 * text_size_scale, family = "arial"),
            titlefont = list(size = 12 * text_size_scale, family = "arial"),
            ticks = ifelse(axis_line, "outside", "none"),
            showline = axis_line,
            showspikes = FALSE,
            tickmode = "array", 
            tickvals = y_ticks,
            ticktext = y_ticks_text
          ),
          zaxis = zaxis_setting,
          camera = list(eye = list(x = -1.25, y = -1.25, z = 1.25))
        ),
        margin = list(
          l = 50,
          r = 50,
          b = 50,
          t = 90,
          pad = 4
        )
      ) %>% 
      plotly::config(
        toImageButtonOptions = list(
          format = "svg",
          filename = plot_title,
          width = NULL,# 500,
          height = NULL,# 500,
          scale = 1
        )
      ) 
  } else { # static plot
    
    # Color scale
    if (start_point < 0 & end_point > 0){
      mid_point <- -start_point/(end_point - start_point)
      color_level <- round(seq(start_point, end_point, by = 2), 0)
      col1 <- grDevices::colorRampPalette(c(
        "white",
        low_value_color
      ))(length(which(color_level <= 0)))
      col2 <- grDevices::colorRampPalette(c(
        "white",
        high_value_color
      ))(length(which(color_level >= 0)))
      col <- c(
        rev(col1), 
        col2[-1]
      )
    } else {
      color_level <- (end_point - start_point + 1)%/%2 + 1
      col <- grDevices::colorRampPalette(c(
        low_value_color,
        high_value_color
      ))(color_level)
    }
    col_break <- c(start_point, end_point)
    
    additional_low_tick <- (start_point - min(plot_table$value)) %/% 2
    additional_high_tick <- (max(plot_table$value) - end_point) %/% 2
    col_expand <- col
    col_break_expand <- col_break
    if (additional_low_tick > 0) {
      col_expand <- c(
        rep(low_value_color, additional_low_tick),
        col_expand
      )
      col_break_expand[1] <- min(plot_table$value)
    }
    if (additional_high_tick > 0) {
      col_expand <- c(
        col_expand,
        rep(high_value_color, additional_high_tick)
      )
      col_break_expand[2] <- max(plot_table$value)
    }
    
    scale_par <- list(
      arrows = FALSE,
      distance = c(0.8, 0.8, 0.8),
      col = 1,
      z = list(
        tick.number = 6,
        cex = 0.75 * text_size_scale),
      x =  list(
        at = x_ticks,
        labels = x_ticks_text,
        cex = 0.75 * text_size_scale
      ),
      y = list(
        at = y_ticks,
        labels = y_ticks_text,
        cex = 0.75 * text_size_scale
      )
    )
  
    # Let the text able to shown outside the plot panel
    clip <- lattice::trellis.par.get("clip")
    clip$panel <- "off"
    lattice::trellis.par.set("clip", clip)
    p <- lattice::wireframe(
      extended_mat,
      scales = scale_par,
      drape = TRUE,
      colorkey = list(
        space = "right",
        width = 2,
        height = 0.4,
        labels = list(cex = 0.75 * text_size_scale),
        col = col,
        at = lattice::do.breaks(
          c(start_point, end_point),
          length(col)
        )
      ),
      screen = list(z = 30, x = -55),
      zlab = list(
        z_axis_title,
        axis.key.padding = 0,
        cex = 0.75 * text_size_scale,
        rot = 93
      ),
      xlab = list(
        x_axis_title,
        cex = 0.75 * text_size_scale,
        rot = 23
      ),
      ylab = list(
        y_axis_title,
        cex = 0.75 * text_size_scale,
        rot = -53
      ),
      zlim = z_range,
      col.regions = col_expand,
      main = list(
        label = plot_title,
        fontsize = 13.5 * text_size_scale
      ),
      at = lattice::do.breaks(
        col_break_expand,
        length(col_expand)
      ),
      par.settings = list(
        axis.line = list(col = "transparent")
      ),
      zoom = 1,
      aspect = 1,
      panel = function(...) {
        lattice::panel.wireframe(...)
        # subtitle
        grid::grid.text(
          label = plot_subtitle,
          x=unit(0.55, "npc"),
          y=unit(1, "npc"),
          gp = grid::gpar(
            col = "black",
            fontsize = 10.5 * text_size_scale
          ) # pt
        )
        # legend title
        grid::grid.text(
          label = legend_title,
          x = unit(1.05, "npc"),
          y = unit(0.75, "npc"),
          gp = grid::gpar(
            col = "black",
            fontsize = 8.5 * text_size_scale
          ) # pt
        )
        },
      pretty = TRUE
    )
    print(p)
    p <- grDevices::recordPlot()
  }
  return(p)
}

# Auxiliary functions -----------------------------------------------------

#' Extract Data for 2 Drug Combination Plots
#'
#' @param data A list object generated by function \code{\link{ReshapeData}}.
#' @param plot_block A character/integer. It indicates the block ID for the
#'   block to visualize.
#' @param drugs A vector of characters or integers with length of 2. It contains
#'   the index for two drugs to plot. For example, \code{c(1, 2)} indicates to
#'   plot "drug1" and "drug2" in the input \code{data}.
#' @param plot_value A character value. It indicates the value to be visualized.
#'   If the \code{data} is the direct output from \link{ReshapeData}, the values
#'   for this parameter are:
#'   \itemize{
#'     \item \strong{response_origin} The original response value in input data.
#'     It might be \% inhibition or \% viability.
#'     \item \strong{response} The \% inhibition after preprocess by function 
#'     \link{ReshapeData}
#'   }
#'   If the \code{data} is the output from \link{CalculateSynergy}, following
#'   values are also available:
#'   \itemize{
#'     \item \strong{ZIP_ref, Bliss_ref, HSA_ref, Loewe_ref} The reference
#'     additive effects calculated by ZIP, Bliss, HSA or Loewe model,
#'     respectively.
#'     \item \strong{ZIP_synergy, Bliss_synergy, HSA_synergy, Loewe_synergy}
#'     The synergy score calculated by ZIP, Bliss, HSA or Loewe model,
#'     respectively.
#'     \item \strong{ZIP_fit} The response fitted by ZIP model.
#'   }
#' @param statistic A character or NULL. It indicates the statistics printed
#'   in the plot while there are replicates in input data. Available values are:
#'   \itemize{
#'     \item \strong{sem} Standard error of mean;
#'     \item \strong{ci} 95\% confidence interval.
#'   }
#'   If it is \code{NULL}, no statistics will be printed.
#'
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#' 
#' @return A data frame. It contains the concentrations for selected drugs, the
#'   selected values for plotting, and the text for printing on the heatmap.
.Extract2DrugPlotData <- function(data,
                                  plot_block = 1,
                                  drugs = c(1, 2),
                                  plot_value = "response",
                                  statistic = NULL){
  # 1. Check the input data
  # Data structure of 'data'
  if (!is.list(data)) {
    stop("Input data is not in list format!")
  }
  if (!all(c("drug_pairs", "response") %in% names(data))) {
    stop("Input data should contain at least tow elements: 'drug_pairs' and 
         'response'. Please prepare your data with 'ReshapeData' function.")
  }
  # Parameter 'drugs'
  if (length(drugs) != 2) {
    stop("The length of 'drugs' parameter is not 2. Please chosed exactly 2
         drugs for heatmap.")
  }
  # Parameter 'plot_value'
  avail_value <- c("response", "response_origin", "ZIP_ref", "ZIP_fit",
                   "ZIP_synergy", "HSA_ref", "HSA_synergy", "Bliss_ref",
                   "Bliss_synergy", "Loewe_ref", "Loewe_synergy")
  if (!plot_value %in% avail_value) {
    stop("The parameter 'plot_value = ", plot_value, "' is not available.",
         "Avaliable values are '", paste(avail_value, collapse = ", "), "'.")
  }
  
  # Annotation data
  drug_pair <- data$drug_pairs[data$drug_pairs$block_id == plot_block, ] %>% 
    dplyr::select(
      drug1 = paste0("drug", drugs[1]),
      drug2 = paste0("drug", drugs[2]),
      conc_unit1 = paste0("conc_unit", drugs[1]),
      conc_unit2 = paste0("conc_unit", drugs[2]),
      replicate,
      input_type
    )
  
  # Parameter 'statistic'
  if (is.null(statistic)){
    statistic_table <- drug_pair$replicate
  } else {
    avail_statistic <- c("sd", "sem", "ci")
    if (!drug_pair$replicate) {
      warning("The selected block ", plot_block,
              " doesn't have the replicate data. Statistics is not available.")
      statistic_table <- FALSE
    } else if(!statistic %in% avail_statistic) {
      warning("The parameter 'statistic = ", statistic, "' is not available.",
              "Avaliable values are ", paste(avail_statistic, sep = ", "), ".")
      statistic_table <- FALSE
    } else {
      statistic_table <- TRUE
    }
  }
  
  # 1. Extract tables for plotting
  
  # Data table
  concs <- grep("conc\\d", colnames(data$response), value = TRUE)
  selected_concs <- paste0("conc", drugs)
  
  if (statistic_table){
    if (startsWith(plot_value, "response")){
      plot_table <- data$response_statistics
    } else {
      if (!"synergy_scores" %in% names(data)){
        stop("The synergy scores are not calculated. Please run function ",
             "'CalculateSynergy' first.")
      }
      plot_table <- data$synergy_scores_statistics
    }
    plot_table <- plot_table %>% 
      dplyr::filter(block_id == plot_block) %>% 
      dplyr::ungroup()
    if (is.null(statistic)) {
      plot_table <- plot_table %>% 
        dplyr::select(
          dplyr::starts_with("conc"), 
          value = !!paste0(plot_value, "_mean")
        ) %>%
        dplyr::mutate(
          text = as.character(.RoundValues(value))
        )
    } else if (statistic == "sd") {
      plot_table <- plot_table %>% 
        dplyr::select(
          dplyr::starts_with("conc"), 
          value = !!paste0(plot_value, "_mean"),
          statistic = !!paste0(plot_value, "_sd")
        ) %>%
        dplyr::mutate(
          text = paste(
            .RoundValues(value), "\n" ,
            "\u00B1",
            .RoundValues(statistic)
          )
        )
    } else if (statistic == "sem") {
      plot_table <- plot_table %>% 
        dplyr::select(
          dplyr::starts_with("conc"), 
          value = !!paste0(plot_value, "_mean"),
          statistic = !!paste0(plot_value, "_sem")
        ) %>%
        dplyr::mutate(
          text = paste(
            .RoundValues(value), "\n" ,
            "\u00B1",
            .RoundValues(statistic)
          )
        )
    } else if (statistic == "ci") {
      plot_table <- plot_table %>% 
        dplyr::select(
          dplyr::starts_with("conc"), 
          value = !!paste0(plot_value, "_mean"),
          left = !!paste0(plot_value, "_ci_left"),
          right = !!paste0(plot_value, "_ci_right")
        ) %>%
        dplyr::mutate(
          text = paste0(
            .RoundValues(value), "\n",
            "[",
            .RoundValues(left),
            ", ",
            .RoundValues(right),
            "]"
          )
        )
    }
  } else {
    if (startsWith(plot_value, "response")){
      plot_table <- data$response
    } else {
      if (!"synergy_scores" %in% names(data)){
        stop("The synergy scores are not calculated. Please run function ",
             "'CalculateSynergy' first.")
      }
      plot_table <- data$synergy_scores
    }
    plot_table <- plot_table %>% 
      dplyr::filter(block_id == plot_block) %>% 
      dplyr::select(
        dplyr::starts_with("conc"),
        value = !!plot_value
      ) %>%
      dplyr::mutate(
        text = as.character(.RoundValues(value))
      )
  }
  
  # Extract data for selected two drugs from multi-drug combination
  other_concs <- setdiff(concs, selected_concs)
  if (length(other_concs) > 0) {
    conc_zero <- apply(
      dplyr::select(plot_table, dplyr::all_of(concs)), 
      1,
      function(x) {
        sum(x != 0) <= 2
      })
    plot_table <- plot_table[conc_zero, ]
    other_concs_sum <- plot_table %>% 
      dplyr::ungroup() %>% 
      dplyr::select(dplyr::all_of(other_concs)) %>% 
      rowSums()
    plot_table <- plot_table[other_concs_sum == 0, ] %>% 
      dplyr::select(-dplyr::all_of(other_concs))
    colnames(plot_table) <- sapply(colnames(plot_table), function(x){
      if (x == selected_concs[1]){
        return("conc1")
      } else if (x == selected_concs[2]) {
        return("conc2")
      } else {
        return(x)
      }
    })
  }
  
  # Transform conc into factor
  plot_table[, c("conc1", "conc2")] <- lapply(
    plot_table[, c("conc1", "conc2")],
    function(x) {
      factor(.RoundValues(x))
    }
  )
  plot_table <- plot_table %>% 
    dplyr::select(conc1, conc2, value, text)
  return(list(plot_table = plot_table, drug_pair = drug_pair))
}

#' Round the Numbers for Plotting
#' 
#' This function will round the input numbers by 2 digits, if the absolute of 
#' number is larger than or equal to 1. It will take 2 significant digits, if
#' the absolute of number is less than 1. 
#'
#' @param numbers A vector of numeric values. It contains the numbers need to
#' be rounded.
#'
#' @return A vector of rounded numbers.
#'
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#' 
.RoundValues <- function(numbers) {
  numbers[abs(numbers) >= 1 & !is.na(numbers)] <- round(
    numbers[abs(numbers) >= 1 & !is.na(numbers)],
    2
  )
  numbers[abs(numbers) < 1 & !is.na(numbers)] <- signif(
    numbers[abs(numbers) < 1 & !is.na(numbers)],
    2
  )
  return(numbers)
}

#' Convert Font Size from pt to mm
#'
#' This function converts font sizes from "pt" unite to "mm" unite.
#'  
#' @param x A numerical value. It is the font size in "pt" unite.
#'
#' @return A numerical value in "mm" unite
#' 
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#' 
.Pt2mm <- function(x) {
  5 * x / 14
}

#' Make a Smooth Surface for Scores
#'
#' @param scores_mat a matrix contains scores which will be visualized
#' @param len length of the interval between plotted data points.
#'
#' @return a matrix which which contains interpolated points for input
#'   \code{scores_mat}.
#' 
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#' 
.ExtendedScores <-  function (scores_mat, len) {
  # len: how many values need to be predicted between two adjacent elements
  #      of scores.dose
  options(scipen = 999)
  nr <- nrow(scores_mat)
  nc <- ncol(scores_mat)
  
  ext.row.len <- (nr - 1) * (len + 2) - (nr - 2)
  ext.col.len <- (nc - 1) * (len + 2) - (nc - 2)
  
  extended.row.idx <- seq(1, nr, length = ext.row.len)
  extended.col.idx <- seq(1, nc, length = ext.col.len)
  
  krig.coord <- cbind(rep(extended.row.idx, each = ext.col.len),
                      rep(extended.col.idx, times = ext.row.len))
  extended.scores <- SpatialExtremes::kriging(data = c(scores_mat),
                             data.coord = cbind(rep(seq_len(nr), nc),
                                                rep(seq_len(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 = ext.row.len,
                            ncol = ext.col.len, byrow = TRUE)
  
  # extended.scores = data.frame(extended.scores)
  extended.scores <- round(extended.scores, 3)
  
  row.dose <- as.numeric(rownames(scores_mat))
  col.dose <- as.numeric(colnames(scores_mat))
  
  extend.row.dose <- mapply(function(x, y){seq(from = x, to = y,
                                               length.out = len + 2)},
                            row.dose[-nr], row.dose[-1])
  #extend.row.dose <- unique(round(c(extend.row.dose), 8))
  extend.row.dose <- unique(signif(c(extend.row.dose), 8))
  
  extend.col.dose <- mapply(function(x, y){seq(from = x, to = y,
                                               length.out = len + 2)},
                            col.dose[-nc], col.dose[-1])
  #extend.col.dose <- unique(round(c(extend.col.dose), 8))
  extend.col.dose <- unique(signif(c(extend.col.dose), 8))
  
  rownames(extended.scores) <- extend.row.dose
  colnames(extended.scores) <- extend.col.dose
  
  return(extended.scores)
}