R/calculate_synergy_score.R

In 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:
#
# CalculateSynergy: Calculate the Synergy Scores for Drug Combinations
# ZIP/Bliss/HSA/Loewe: 4 functions to calculate synergy scores
# CorrectBaseLine: Base Line Correction for Drug Combination Matrix
#
# Auxiliary functions:
# .Distance: Calculate Distance from a Point to a Plane
# .SolveExpDose/L4/LL4: 3 functions to solve the expected dose of drug at which
#                      it could achieve the given % inhibition
# .Bootstrapping: Bootstraping Sample from Replicates in Response Data

#' Calculate the Synergy Scores for Drug Combinations
#'
#' \code{CalculateSynergy} is the main function for calculating synergy scores
#'   based on model(ZIP, Bliss, Loewe, and HSA) from one dose-response
#' \strong{matrix}.
#'
#' @param data A list object generated by function \code{\link{ReshapeData}}.
#' @param method A parameter to specify which models to use to calculate the
#'   synergy scores. Choices are "ZIP", "Bliss", "HSA" and "Loewe". Defaults to
#'   "ZIP".
#' @param Emax The expected maximum response value in the 4 parameter 
#'   log-logistic model. It is used while calling \code{\link{ZIP}} and
#'   \code{\link{Loewe}}.
#' @param Emin The expected minimum response value in the 4 parameter 
#'   log-logistic model. It is used while calling \code{\link{ZIP}} and
#'   \code{\link{Loewe}}.
#' @param adjusted A logical value. If it is \code{TRUE}, the
#'   'adjusted.response.mats' will be used to calculate synergy scores. If it is
#'   \code{FALSE}, the raw data ('dose.response.mats') will be used to calculate
#'   synergy scores.
#' @param correct_baseline  A character value. It indicates the method used for
#'   baseline correction. Available values are:
#'   \itemize{
#'     \item \strong{non} No baseline correction.
#'     \item \strong{part} Adjust only the negative values in the matrix.
#'     \item \strong{all} Adjust all values in the matrix.
#'   }
#' @param iteration An integer value. It indicates the number of iterations for
#'   synergy scores calculation on data with replicates.
#' @param seed An integer or NULL. It is used to set the random seed in synergy
#'   scores calculation on data with replicates. 
#'
#' @return This function will add 1 or 2 elements into inputted \code{data}
#'   list:
#'   \itemize{
#'     \item \strong{scores} A data frame. It contains synergy scores, 
#'     reference effects and fitted response values (only for "ZIP" model) 
#'     calculated by selected \code{method}. If there are replicates in the
#'     block, the mean values across all iterations for all the metrics
#'     mentioned above will be output.
#'     \item \strong{scores_statistics} A data frame. It will be output if
#'     there is block have replicated response values. It contains the
#'     the statistics (including mean, standard deviation, standard error of 
#'     mean and 95% confidence interval) for synergy scores, reference effects
#'     and fitted response values (only for "ZIP" model) across results from
#'     iterations.
#'  }
#'  This function also add mean of synergy scores across whole combination
#'  matrix to the \code{data$drug_pair} table.
#'
#' @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)
#' scores <- CalculateSynergy(data)
CalculateSynergy <- function(data,
                             method = c("ZIP", "HSA", "Bliss", "Loewe"),
                             Emin = NA,
                             Emax = NA,
                             adjusted = TRUE,
                             correct_baseline = "non",
                             iteration = 10,
                             seed = 123) {
  options(scipen = 999)
  # 1. Check the input 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.")
  }

  # 2. Select the dose response table for plotting.
  if (adjusted) {
    response <- data$response %>%
      dplyr::select(-response_origin)
  } else {
    response <- data$response %>%
      dplyr::select(-response) %>% 
      dplyr::rename(response = response_origin)
  }
  
  # 3. Calculate synergy scores
  if (!all(method %in% c("ZIP", "HSA", "Bliss", "Loewe"))) {
    stop("The method parameter can only be one of the following: ZIP, HSA, Bliss
       and Loewe.")
  }

  blocks <- unique(response$block_id)
  scores <- NULL
  scores_statistics <- NULL
  for (b in blocks) {      
    response_one_block <- response %>% 
      dplyr::filter(block_id == b) %>% 
      dplyr::select(-block_id) %>% 
      dplyr::ungroup()
    # Correct baseline
    concs <- grep("conc\\d", colnames(response_one_block), value = TRUE)
    tmp <- dplyr::select(response_one_block, dplyr::all_of(concs)) %>% 
      unique()
    tmp_scores_statistic <- tmp
    # Calculate synergy score. Different work flow for replicated data
    replicate <- data$drug_pairs$replicate[data$drug_pairs$block_id == b]
    if (replicate) { # Block with replicate
      for (m in method) {
        set.seed(seed)
        message("Calculating ", m, " score for block ", b, "...")
        iter <- pbapply::pblapply(seq(1, iteration), function(x){
          response_boot <- .Bootstrapping(response_one_block)
          response_boot <- CorrectBaseLine(
            response_boot,
            method = correct_baseline
          )
          s <- eval(call(m, response_boot))
          return(s)
          }) %>% 
          purrr::reduce(dplyr::left_join, by = concs) %>% 
          dplyr::ungroup()
        # column names for scores
        s <- grep("conc\\d", colnames(iter), value = TRUE, invert = TRUE)
        s <- unique(gsub("\\..*$", "", s))
        tmp_m <- dplyr::select(iter, concs)
        tmp_score_statistic_m <- tmp_m
        # combo cell index
        concs_iter <- iter[, grepl("conc\\d", colnames(iter))]
        concs_iter <- apply(concs_iter, 2, function(x){x == 0})
        index <- rowSums(concs_iter) < 1
        for (i in s){
          tmp_m[[i]] <- rowMeans(dplyr::select(iter, dplyr::starts_with(i)))
          tmp_score_statistic_m[[paste0(i, "_mean")]] <- tmp_m[[i]]
          tmp_score_statistic_m[[paste0(i, "_sem")]] <- 
            apply(dplyr::select(iter, dplyr::starts_with(i)), 1, stats::sd)
          tmp_score_statistic_m[[paste0(i, "_ci_left")]] <- 
            apply(dplyr::select(iter, dplyr::starts_with(i)), 1, 
                  function(x) stats::quantile(x, probs = 0.025))
          tmp_score_statistic_m[[paste0(i, "_ci_right")]] <- 
            apply(dplyr::select(iter, dplyr::starts_with(i)), 1, 
                  function(x) stats::quantile(x, probs = 0.975))
          # calculate P value for synergy scores
          if (endsWith(i, "synergy")) {
            matrix_mean <- colMeans(iter[index, startsWith(colnames(iter), i)])
            z <- abs(mean(matrix_mean)) / stats::sd(matrix_mean)
            p <- exp(-0.717 * z - 0.416 * z ^2)
            p <- formatC(p, format = "e", digits = 2, zero.print = "< 2e-324")
            data$drug_pairs[data$drug_pairs$block_id == b, 
                            paste0(i, "_p_value")] <- p
          }
        }
        tmp_scores_statistic <- dplyr::left_join(tmp_scores_statistic, 
                                                 tmp_score_statistic_m,
                                                 by = concs)
        tmp <- dplyr::left_join(tmp, tmp_m, by = concs)
      }
      tmp$block_id <- b
      tmp_scores_statistic$block_id <- b
      scores <- rbind.data.frame(scores, tmp)
      scores_statistics <- rbind.data.frame(scores_statistics, 
                                           tmp_scores_statistic)
    } else { # Blocks without replicates
      # Correct base line
      response_one_block <- CorrectBaseLine(
        response_one_block,
        method = correct_baseline
      )
      tmp <- dplyr::select(response_one_block, dplyr::all_of(concs)) %>% 
        unique()
      for (m in method) {
        message("Calculating ", m, " score for block ", b, "...")
        if (m %in% c("Bliss", "HSA")) {
          fun <- call(
            m,
            response_one_block
          )
        } else {
          fun <- call(
            m,
            response_one_block,
            Emin = Emin,
            Emax = Emax
          )
        }
        tmp_m <- eval(fun)
        tmp <- tmp %>% 
          dplyr::left_join(tmp_m, by = concs)
      }
      tmp$block_id <- b
      tmp <- dplyr::select(tmp, block_id, dplyr::everything())
      scores <- rbind.data.frame(scores, tmp)
    }
  }
  
  ## 4. Save data into the list
  data$synergy_scores <- dplyr::select(scores, block_id, dplyr::everything())
  if (length(scores_statistics) != 0){
    data$synergy_scores_statistics <- dplyr::select(scores_statistics, block_id,
                                            dplyr::everything())
  }
  # combo cell index
  concs <- data$synergy_scores[, grepl("conc\\d", colnames(data$synergy_scores))]
  concs <- apply(data$synergy_scores, 2, function(x){x == 0})
  index <- rowSums(concs) < 1
  summarized_score <- data$synergy_scores %>% 
    dplyr::select(block_id, dplyr::ends_with("_synergy")) %>% 
    dplyr::filter(index) %>% 
    dplyr::group_by(block_id) %>%
    dplyr::summarise_all(mean) %>% 
    dplyr::ungroup()
  data$drug_pairs <- data$drug_pairs %>% 
    dplyr::select(-dplyr::ends_with("_synergy")) %>% 
    dplyr::left_join(summarized_score, by = "block_id")
  return(data)
}

#' Calculate Delta Synergy Score Based on ZIP Model
#'
#' \code{ZIP} calculates the Delta score matrix from a dose-response
#' matrix by using Zero Interaction Potency (ZIP) method.
#'
#' @details Zero Interaction Potency (ZIP) is a reference model for evaluating
#' the conbimation effect of two drugs. It captures the effect of drug
#' combination by comparing the change in the potency of the dose-response
#' curves between individual drugs and their combinations. \cr
#' \cr
#' The optional arguments \code{drug.col.model}, \code{drug.row.model} are
#' designed for reuse the single drug dose response model fitting results, if
#' it has been down before. Functions \code{\link{FitDoseResponse}} and
#' \code{\link{ExtractSingleDrug}} could be used to calculate these arguments.
#'
#' @param response A data frame. It must contain the columns: "conc1", "conc2",
#'   ..., for the concentration of the combined drugs and "response" for the
#'   observed \%inhibition at certain combination.
#' @param Emax The expected maximum response value in the 4 parameter 
#'   log-logistic model.
#' @param Emin The expected minimum response value in the 4 parameter 
#'   log-logistic model.
#' @param quiet A logical value. If it is \code{TRUE} then the warning message
#'   will not show during calculation.
#'  
#' @return A data frame containing the concentrations for drugs, reference
#'   effect, fitted response and synergy score estimated by ZIP model.
#'
#' @author
#'   \itemize{
#'     \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'     \item Jing Tang \email{jing.tang@helsinki.fi}
#'  }
#'
#' @references 
#'   \itemize{
#'     \item{Yadav B, Wennerberg K, Aittokallio T, Tang J. (2015).
#'     \href{https://doi.org/10.1016/j.csbj.2015.09.001}{Searching for Drug
#'     Synergy in Complex Dose-Response Landscape Using an Interaction Potency
#'     Model.} Comput Struct Biotechnol J, 13:504– 513.}
#'   }
#' @export
#'
#' @examples
#' # No single drug fitted modle before
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' response <- data$response[data$response$block_id == 1, 
#'                           c("conc1", "conc2", "response")]
#' ZIP_score <- ZIP(response)
#' 
#' \dontrun{
#' # Parallel processing:
#' if (future::supportsMulticore()) {
#'   future::plan(future::multicore)
#' } else {
#'   future::plan(future::multisession)
#' }
#' ZIP(response)
#' # future::plan(future::sequential) # Turn off the multicore setting
#' }
ZIP <- function(response,
                Emin = NA,
                Emax = NA,
                quiet = TRUE) {
  if (quiet) {
    options(warn = -1)
  }
  # 1. Calculate ZIP reference effect
  
  # Get predicted response for all single drugs
  single_drug_data <- ExtractSingleDrug(response)
  single_drug_model <- lapply(
    single_drug_data,
    function(x) FitDoseResponse(x, Emin = Emin, Emax = Emax)
  )
  
  single_drug_pred <- lapply(single_drug_model, 
                         function(x) {
                           data.frame(dose = x$origData$dose,
                                      pred = stats::predict(x),
                                      stringsAsFactors = FALSE)
                           })
  
  ref_zip <- expand.grid(lapply(single_drug_pred, function(x) x$dose))
  ref_zip$ZIP_ref <- apply(
    expand.grid(lapply(single_drug_pred, function(x) x$pred)),
    1, function(x) {
      (1 - prod(1 - x / 100)) * 100
    }
  )
  
  # 2. Calculate ZIP fitted effect
  concs <- grep("conc\\d", colnames(response), value = TRUE)
  response <- purrr::map(dplyr::all_of(concs), function(conc) {
    response %>% 
      dplyr::ungroup() %>% 
      dplyr::rename(dose = conc) %>% 
      # Group table by conditions, and wrap as a input data for model fitting
      tidyr::nest(dplyr::any_of(c("dose", "response"))) %>% 
      dplyr::mutate(pred = furrr::future_map(data, function(x, conc) {
        condition_baseline <- x$response[which(x$dose == 0)]
        model <- FitDoseResponse(x, Emin = condition_baseline,
                                 Emax = NA) # Fit dose response curve
        pred <- stats::predict(model) # Predict response on corresponding dosage
        return(pred)
      })) %>% 
      tidyr::unnest() %>% 
      dplyr::rename(!!conc := "dose")
  })
  
  # Take the average of fitted response as the fit_zip
  response <- response %>% 
    purrr::reduce(dplyr::left_join, by = c(concs, "response")) %>% 
    dplyr::mutate(ZIP_fit = rowMeans(dplyr::select(., dplyr::starts_with("pred"))))
  
  # 3. Calculate synergy score
  response <- response %>% 
    dplyr::left_join(ref_zip, by = concs)
  
  # Assign response value to fit_zip in non-combination wells (single drug or DMSO)
  no_comb_rows <- apply(
    response[, grep("conc\\d", colnames(response), value = TRUE)], 1,
    function(x) {
      (length(x) - sum(x == 0)) <= 1
    }
  )
  response$ZIP_fit[which(no_comb_rows)] <- response$response[which(no_comb_rows)]
  response$ZIP_ref[which(no_comb_rows)] <- response$response[which(no_comb_rows)]
  
  response <- response %>% 
    dplyr::mutate(ZIP_synergy = ZIP_fit - ZIP_ref) %>% 
    dplyr::select(!!concs, ZIP_fit, ZIP_ref, ZIP_synergy)
  
  return(response)
  options(warn = 0)
}

#' Calculate Bliss Synergy Score
#'
#' \code{Bliss} calculates the synergy score matrix for a block of
#' drug combination by using a drug interaction reference model introduced by
#' C. I. Bliss in 1939.
#'
#' This model is a reference model for evaluating the combination effect of two
#' drugs. The basic assumption of this model is "The expected effect of two
#' drugs acting independently". The Bliss reference effect 
#' y = 1 - product_all_drug(1-\%Inhibition) * 100.
#'
#' @param response A data frame. It must contain the columns: "conc1", "conc2",
#'   ..., for the concentration of the combined drugs and "response" for the
#'   observed \%inhibition at certain combination.
#' @param single_drug_data A list or NULL. It contains the monotherapy dose
#'   response data for all the tested drugs in the inputted block. If it is
#'   \code{NULL}, the data will extract the dose response table from inputted
#'   \code{response} table.
#'
#' @return  A data frame containing the concentrations for drugs, reference
#'   effect and synergy score estimated by Bliss model.
#'
#' @author
#'   \itemize{
#'     \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'     \item Jing Tang \email{jing.tang@helsinki.fi}
#'  }
#'
#' @references 
#'   \itemize{
#'      \item{Yadav B, Wennerberg K, Aittokallio T, Tang J. (2015).
#'      \href{https://doi.org/10.1016/j.csbj.2015.09.001}{Searching for Drug
#'      Synergy in Complex Dose-Response Landscape Using an Interaction Potency
#'      Model.} Comput Struct Biotechnol J, 13:504– 513.}
#'      \item{Bliss, C. I. (1939).
#'      \href{https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1744-7348.1939.tb06990.x}{The toxicity of poisons applied jointly.}
#'      Annals of Applied Biology, 26(3):585–615.}
#' }
#'
#' @export
#'
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' response <- data$response[data$response$block_id == 1,
#'                           c("conc1", "conc2", "response")]
#' Bliss.score <- Bliss(response)
Bliss <- function(response, single_drug_data) {
  # Get all possible combinations of drug dosages
  single_drug_data <- ExtractSingleDrug(response = response)
  bliss <- expand.grid(lapply(single_drug_data, function(x) x$dose))

  # Calculate the Bliss reference effect
  ref <- expand.grid(lapply(single_drug_data, function(x) x$response))
  bliss$Bliss_ref <- apply(ref, 1, function(x) (1 - prod(1 - x / 100)) * 100)
  concs <- grep("conc\\d", colnames(response), value = TRUE)
  bliss <- response %>%
    dplyr::left_join(bliss, by = concs)
  # Assign response value to reference additive effects in non-combination wells
  # (single drug or DMSO)
  no_comb_rows <- apply(
    bliss[, concs], 1,
    function(x) {
      (length(x) - sum(x == 0)) <= 1
    }
  )
  bliss$Bliss_ref[which(no_comb_rows)] <- bliss$response[which(no_comb_rows)]
  
  # Calculate Bliss synergy score
  bliss <- bliss %>%
    dplyr::mutate(Bliss_synergy = response - Bliss_ref) %>% 
    dplyr::select(-response)
  return(bliss)
}

#' Calculate HSA Synergy Score
#'
#' \code{HSA} calculates the synergy score matrix for a block of
#' drug combination by using Highest Single Agent (HSA) reference model.
#'
#' This model is a reference model for evaluating the combination effect of two
#' drugs. The basic assumption of this model is "The reference effect of drug
#' combination is the maximal single drug effect".
#'
#' @param response A data frame. It must contain the columns: "conc1", "conc2",
#'   ..., for the concentration of the combined drugs and "response" for the
#'   observed \%inhibition at certain combination.
#'
#' @return  A data frame containing the concentrations for drugs, reference
#'   effect and synergy score estimated by HSA model.
#' 
#' @author
#'   \itemize{
#'     \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'     \item Jing Tang \email{jing.tang@helsinki.fi}
#'  }
#'
#' @references \itemize{
#'    \item{Yadav B, Wennerberg K, Aittokallio T, Tang J.(2015).
#'    \href{https://doi.org/10.1016/j.csbj.2015.09.001}{Searching for Drug
#'    Synergy in Complex Dose-Response Landscape Using an Interaction Potency
#'    Model.}Comput Struct Biotechnol J, 13:504– 513.}
#'    \item{Berenbaum MC. (1989).
#'    \href{https://www.ncbi.nlm.nih.gov/pubmed/2692037}{What is synergy?}
#'    Pharmacol Rev 1990 Sep;41(3):422.
#'    }
#' }
#'
#' @export
#'
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' response <- data$response[data$response$block_id == 1, 
#'                           c("conc1", "conc2", "response")]
#' HSA.score <- HSA(response)
HSA <- function(response) {
  
  # Get all possible combinations of drug dosages
  single_drug_data <- ExtractSingleDrug(response = response)
  hsa <- expand.grid(lapply(single_drug_data, function(x) x$dose))

  # Calculate the HSA reference effect (max effect in all single drugs)
  ref <- expand.grid(lapply(single_drug_data, function(x) x$response))
  hsa$HSA_ref <- apply(ref, 1, max)
  concs <- grep("conc\\d", colnames(response), value = TRUE)

  hsa <- response %>%
    dplyr::left_join(hsa, by = concs)
  
  # Assign response value to reference additive effects in non-combination wells
  # (single drug or DMSO)
  no_comb_rows <- apply(
    hsa[, concs], 1,
    function(x) {
      (length(x) - sum(x == 0)) <= 1
    }
  )
  
  hsa$HSA_ref[which(no_comb_rows)] <- hsa$response[which(no_comb_rows)]
  
  # Calculate HSA synergy score
  hsa <- hsa %>%
    dplyr::mutate(HSA_synergy = response - HSA_ref) %>% 
    dplyr::select(-response)
  return(hsa)
}

#' Calculate Loewe Synergy Score
#'
#' \code{Loewe} calculates the synergy score matrix from a dose-response matrix
#' by using a druginteraction reference model introduced by Loewe in 1953.
#'
#' @details Loewe model is a reference model for evaluating the combination
#' effect of two drugs. The basic assumption of this model is "The referece
#' effect of drug combination is the expected effect of a drug combined with
#' itself". \cr
#' \cr
#' The optional arguments \code{drug.col.model}, \code{drug.row.model} are
#' designed for reuse the single drug dose response model fitting results,
#' if it has been down before. Functions \code{\link{FitDoseResponse}} and
#' \code{\link{ExtractSingleDrug}} could be used to calculate these arguments.
#'
#' @param response A data frame. It must contain the columns: "conc1", "conc2",
#'   ..., for the concentration of the combined drugs and "response" for the
#'   observed \%inhibition at certain combination.
#' @param Emax The expected maximum response value in the 4 parameter 
#'   log-logistic model.
#' @param Emin The expected minimum response value in the 4 parameter 
#'   log-logistic model.
#' @param quiet A logical value. If it is \code{TRUE} then the warning message
#'   will not show during calculation.
#'
#' @return A data frame containing the concentrations for drugs, reference
#'   effect and synergy score estimated by Loewe model.
#'
#' @author
#'   \itemize{
#'     \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'     \item Jing Tang \email{jing.tang@helsinki.fi}
#'  }
#'
#' @references \itemize{
#'    \item{Yadav B, Wennerberg K, Aittokallio T, Tang J.(2015).
#'    \href{https://doi.org/10.1016/j.csbj.2015.09.001}{Searching for Drug
#'    Synergy in Complex Dose-Response Landscape Using an Interaction Potency
#'    Model.}Comput Struct Biotechnol J, 13:504– 513.}
#'    \item{[Loewe, 1953] Loewe, S. (1953).
#'    \href{https://www.ncbi.nlm.nih.gov/pubmed/13081480}{The problem of
#'    synergism and antagonism of combined drugs.} Arzneimittelforschung,
#'    3(6):285–290.
#'    }
#' }
#'
#' @export
#'
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' response <- data$response[data$response$block_id == 1, 
#'                           c("conc1", "conc2", "response")]
#' Loewe.score <- Loewe(response)
#'
Loewe <- function(response,
                  Emin = NA,
                  Emax = NA,
                  quiet = TRUE) {
  if (quiet) {
    options(warn = -1)
  }
  
  ndrugs <- ncol(response) - 1

  single_drug_data <- ExtractSingleDrug(response)
  single_drug_model <- lapply(
    single_drug_data,
    function(x) FitDoseResponse(x, Emin = Emin, Emax = Emax)
  )

  single_par <- lapply(single_drug_model, FindModelPar)
  single_type <- lapply(single_drug_model, FindModelType)

  y_loewe <- c()
  dist_loewe <- c()
  for (i in 1:nrow(response)) {
    x <- response[i, c(1:ndrugs)] # concentrations of drugs
    y <- response$response[i] # the observed combination response

    if (length(which(x > 0)) < 2) { # single drugs
      y_loewe[i] <- y
      dist_loewe[i] <- NA
    } else {
      # find the dose of single drugs that achieve the observed combination response
      x_cap <- mapply(function(par, type) .SolveExpDose(y, par, type),
        par = single_par, type = single_type
      )

      if (all(!is.finite(x_cap))) { # if none of drug achieve the combination response
        # max of the single drug response
        y_loewe[i] <- max(mapply(function(model) {
          PredictModelSpecify(model, sum(x))
        },
        model = single_drug_model
        ))
        dist_loewe[i] <- NA
      } else {
        # determine the minimal distance
        tmp <- .SolveLoewe(x, single_par, single_type, nsteps = 100)
        y_loewe[i] <- tmp$y_loewe
        dist_loewe[i] <- tmp$distance
      }
    }
  }
  response$Loewe_ref <- y_loewe

  response$Loewe_synergy <- response$response - y_loewe
  response <- dplyr::select(response, -response)
  # Output results as a long table format
  return(response)
  options(warn = 0)
}


#' Base Line Correction for Drug Combination Matrix
#'
#' \code{CorrectBaseLine} adjusts the base line of drug combination
#' dose-response matrix to make it closer to 0.
#'
#' @param response A drug cobination dose-response matrix. It's column name
#'   and row name are representing the concerntrations of drug added to column
#'   and row, respectively. The values in matrix indicate the inhibition rate to
#'   cell growth.
#' @param method A character value. It indicates the method used for
#'   baseline correction. Available values are:
#'   \itemize{
#'     \item \strong{non} No baseline correction.
#'     \item \strong{part} Adjust only the negative values in the matrix.
#'     \item \strong{all} Adjust all values in the matrix.
#'   }
#'
#' @return A matrix which base line have been adjusted.
#'
#' @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)
#' response <- data$response[data$response$block_id == 1, ]
#' adjusted <- CorrectBaseLine(response, method = "part")
CorrectBaseLine <- function(response, method = c("non", "part", "all")) {
  method <- match.arg(method)
  if (method != "non") {
    single_drug_data <- ExtractSingleDrug(response)
    single_drug_model <- lapply(
      single_drug_data,
      function(x) FitDoseResponse(x)
    )
    
    baseline <- Reduce(min, lapply(single_drug_model, stats::fitted))
    
    if (method == "part") {
      response$response[response$response < 0] <- 
        response$response[response$response < 0] - 
        ((100 - response$response[response$response < 0]) / 100 * baseline)
    } else {
      response$response <- response$response - 
        ((100 - response$response) / 100 * baseline)
    }
  }
  return(response)
}


# Auxiliary Functions -----------------------------------------------------

#' Calculate Distance from a Point to a Plane
#' 
#' This function is used to calculate the distance from a point to a plane. It
#' could also be used in high dimension spaces. The formula comes from
#' https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_plane
#' For two dimension point, the distance to the line w1\*x+w2\*y+b = 0
#' For three dimension point, the distance to the plan w1\*x+w2\*y+w3\*z+b = 0
#'
#' @param w A numeric vector. It contains the parameters for all the coordinates
#'   in the spaces to define the "plan".
#' @param b A numeric value. It is the constant values in the formula which
#'   defines the "plan".
#' @param point A numeric vector. It contains the coordinates in
#'   the spaces to define the "point".
#'  
#' @return A numeric value. It is the distance from point defined by \code{x0}
#'   to the "plane" defined by \code{w} and \code{b}
#'
#' @author
#'   \itemize{
#'     \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'     \item Jing Tang \email{jing.tang@helsinki.fi}
#'  }
#' 
.Distance <- function(w, b, point){
  d <- abs(point %*% w + b) / sqrt(sum(w^2))
  return(d)
}

#' Solve the Expected Dose of Drug to Achieve Given Effect from LL.4 Model
#'
#' This function will solve the fitted four-parameter log-logistic dose-response 
#' model and output the dose of drug at which it could achieve the \% inhibition
#' to cell growth.
#' 
#' @param y The expected effect (\% inhibition) of the drug to cell line
#' @param drug_par The parameters for fitted dose-response model.
#' 
#' @return A numeric value. It indicates the expected dose of drug.
#' 
#' @author
#'   \itemize{
#'     \item Jing Tang \email{jing.tang@helsinki.fi}
#'     \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   }
#' 
.SolveExpDoesLL4 <- function(y, drug_par){
  res <- drug_par[4] * (((y - drug_par[3]) / (drug_par[2] - y)) ^ (1/drug_par[1]))
  # if NAN it means the response cannot be achieved by the drug, the response 
  # can be either too high or too low for the drug to achieve
  if (is.nan(res) == TRUE){
    res <- ifelse(y > max(drug_par[3], drug_par[2]), 1, -1)*Inf
  }
  return(res)
}

#' Solve the Expected Dose of Drug to Achieve Given Effect from L.4 Model
#'
#' This function will solve the fitted four-parameter logistic dose-response 
#' model and output the dose of drug at which it could achieve the \% inhibition
#' to cell growth.
#' 
#' @param y The expected effect (\% inhibition) of the drug to cell line
#' @param drug_par The parameters for fitted dose-response model.
#'
#' @return A numeric value. It indicates the expected dose of drug.
#' 
#' @author
#'   \itemize{
#'     \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'     \item Jing Tang \email{jing.tang@helsinki.fi}
#'  }
#' 
.SolveExpDoesL4 <- function(y, drug_par) {
  res <- exp((drug_par[4] + log((drug_par[3] - y) /
                                 (y - drug_par[2])) / drug_par[1]))
  if (is.nan(res) == TRUE) {
    res <- ifelse(y > max(drug_par[3], drug_par[2]), 1, -1)*Inf
  }
  return(res)
}

#' Solve the Expected Dose of Drug to Achieve Given Effect (\% inhibition)
#'
#' This function will solve the fitted dose-response model and output the dose
#' of drug at which it could achieve the \% inhibition to cell growth.
#' 
#' @param y The expected effect (\% inhibition) of the drug to cell line.
#' @param drug_par The parameters for fitted dose-response model.
#' @param drug_type The type of model was used to fit the dose-response curve.
#'   Available values are "L.4" - four-parameter logistic model; "LL.4" - 
#'   four-parameter log-logistic model.
#' 
#' @return A numeric value. It indicates the expected dose of drug.
#' 
#' @author
#'   \itemize{
#'     \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'     \item Jing Tang \email{jing.tang@helsinki.fi}
#'  }
#' 
.SolveExpDose <- function(y, drug_par, drug_type){
  switch(as.character(drug_type), 
         "LL.4" = .SolveExpDoesLL4(y, drug_par),
         "L.4"  = .SolveExpDoesL4(y, drug_par))
}

#' Solve the Loewe Additive Effect for Concentration Combinations Isobologram
#'
#' @param concs A numeric vector. It contains the concentrations of tested 
#'   drugs.
#' @param drug_par A numeric vector. The parameters in fitted dose response 
#'   curve.
#' @param drug_type The type of model used to fit dose response curve.
#' @param nsteps The total steps to calculate concentration combinations 
#'   approaching to the true Loewe effect.
#'   
#' @return A list contains 3 items:
#'   \itemize{
#'     \item y_loewe the predicted Loewe additive effect which closes to .
#'     \item x_select the expected concentrations for each drug to achieve
#'       y_loewe.
#'     \item distance the smallest distance 
#'   }
#'
#' @author
#'   \itemize{
#'     \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'     \item Jing Tang \email{jing.tang@helsinki.fi}
#'  }
#' 
.SolveLoewe <- function(concs, drug_par, drug_type, nsteps = 100){
  # x is the concentration vector
  x <- as.matrix(concs)
  ndrugs <- length(x)
  dist <- rep(0, nsteps)
  x_test <- mat.or.vec(ndrugs, nsteps)
  # Find the minimal response of all the drugs
  min_y <- min(unlist(lapply(drug_par, function(x) min(x[2], x[3])))) 
  # Find the maximal response of all the drugs
  max_y <- max(unlist(lapply(drug_par, function(x) max(x[2], x[3]))))
  y_test <- seq(min_y, max_y, length.out = nsteps) # test nsteps responses
  # Calculate expected dose at each drug
  for(i in 1:ndrugs){
    tmp <- lapply(y_test, 
                  function(y) {
                    .SolveExpDose(y, 
                                  drug_par = drug_par[[i]], 
                                  drug_type = drug_type[[i]]
                                 )
                  })
    x_test[i,] <- unlist(tmp)
  }
  # Calculate distance between point x to the expected dose plane
  for(j in 1:nsteps){ 
    # note the sign of -1
    dist[j] = .Distance(w = c(1 / x_test[, j]), b = -1, point = x)
  }
  # output the y.loewe corresponding to the minimal distance
  res <- list(y_loewe = y_test[which(dist == min(dist, na.rm = T))], 
             x_select = x_test[,which(dist == min(dist, na.rm = T))], 
             distance = min(dist, na.rm = T)) 
}

#' Bootstraping Sample from Replicates in Response Data
#'
#' @param response A data frame. It contains the dose response information
#'   about one drug combination block with replicates. It must contain the
#'   columns "conc1", "conc2", ... for concentrations of drugs tested and the
#'   "response" column for observed \% inhibition if cell growth.
#'
#' @return A data frame. It contains a full drug combination matrix whose data
#'   points are randomly selected from replicates.
#' 
#' @author
#'   \itemize{
#'     \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'     \item Jing Tang \email{jing.tang@helsinki.fi}
#'  }
#'  
.Bootstrapping <- function(response){
  # unique dose conditions
  concs <- grep("conc\\d+", colnames(response), value = TRUE)
  response <- response %>%
    dplyr::group_by_at(concs) %>%
    dplyr::sample_n(1)
  return(response)
}