R/calculate_sensitivity_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:
#
# CalculateSensitivity: Calculate the synergy scores for drug combinations
# CalculateCSS: Calculate the synergy scores for drug combinations
# CalculateSens: Calculate sensitivity score (relative inhibition)
# ImputeIC50: Impute missing value at IC50 concentration of drug
# PredictResponse: Predict response value at certain drug dose
# CalculateIC50: Transform IC50 from coefficients from fitted dose-response
#                model
#
# Auxiliary functions:
# .ScoreCurve/.ScoreCurve_L4: facility functions for CalculateSens
# .own_log/.own_log2: facility functions for CalculateSens

#' 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) fron one dose-response
#' \strong{matrix}.
#'
#' @param data A list object generated by function \code{\link{ReshapeData}}.
#' @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. 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 columns into \code{data$drug_pairs} table.
#'   The columns are:
#'   \itemize{
#'     \item \strong{ic50_1/2/...} Relative IC50 for drug 1, 2, ...
#'     \item \strong{ri_1/2/...} Relative Inhibition (RI) for drug 1, 2, ...
#'     \item \strong{css1_ic502/...} CSS score of drug 1 while fixing drug 2
#'       at its IC50.
#'     \item \strong{css} Over all CSS score for the whole block. It's the mean
#'       value of the CSS for all drug pairs in the combination. 
#'  }
#'  If there are replicates in the block, this function will add one table named
#'  as "sensitivity_scores_statistics" for the statistics of the values
#'  mentioned about into the input \code{data} list.
#'
#' @export
#' 
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#' 
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' data <- CalculateSensitivity(data)
CalculateSensitivity <- function(data,
                                 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 RI and CSS
  blocks <- unique(response$block_id)
  scores <- NULL
  scores_statistics <- NULL
  for (b in blocks) {
    message("Calculating sensitive scores for block ", b, " ...")
    response_one_block <- response %>% 
      dplyr::filter(block_id == b) %>% 
      dplyr::select(-block_id) %>% 
      dplyr::ungroup()
    concs <- grep("conc\\d", colnames(response_one_block), value = TRUE)
    
    if (data$drug_pairs$replicate[data$drug_pairs$block_id == b]) {
      tmp_iter <- NULL
      set.seed(seed)
      pb <- utils::txtProgressBar(min = 1, max = iteration, style = 3)
      for(i in 1:iteration){
        response_boot <- .Bootstrapping(response_one_block)
        response_boot <- CorrectBaseLine(
          response_boot,
          method = correct_baseline
        )
        # Calculate RI
        single_drug_data <- ExtractSingleDrug(response_boot)
        ri <- as.data.frame(lapply(single_drug_data, CalculateSens))
        colnames(ri) <- sub("conc", "ri_", colnames(ri))
        
        # Calculate IC50 for all drugs
        ic50 <- lapply(
          single_drug_data,
          function(x) {
            model <- FitDoseResponse(x)
            coe <- FindModelPar(model)
            type <- FindModelType(model)
            CalculateIC50(coe, type, max(x$dose))
          }) %>% 
          as.data.frame()
        colnames(ic50) <- sub("ri", "ic50", colnames(ri))
        
        # Calculate CSS
        css <- CalculateCSS(response_boot, ic50 = ic50)
        
        # Assemble data frame
        tmp <- cbind.data.frame(ic50, ri, css)
        tmp_iter <- rbind.data.frame(tmp_iter, tmp)
        utils::setTxtProgressBar(pb, i)
      }
      message("\n")
      SensMean <- colMeans(tmp_iter)
      tmp <- as.data.frame(as.list(SensMean)) %>%
        dplyr::mutate(block_id = b)
      SensSd <- apply(tmp_iter, 2, stats::sd)
      SensSem <- SensSd / iteration
      SensCI95_left <- apply(tmp_iter, 2, 
                             function(x) stats::quantile(x, probs = 0.025))
      SensCI95_right <- apply(tmp_iter, 2, 
                             function(x) stats::quantile(x, probs = 0.975))
      p_value <- apply(
        tmp_iter, 2, 
        function(x) {
          z <- abs(mean(x)) / stats::sd(x)
          p <- exp(-0.717 * z - 0.416 * z ^2)
          p <- formatC(p, format = "e", digits = 2, zero.print = "< 2e-324")
          return(p)
        })
      
      names(SensMean) <- paste0(names(SensMean), "_mean")
      names(SensSem) <- paste0(names(SensSem), "_sem")
      names(SensCI95_left) <- paste0(names(SensCI95_left), "_ci_left")
      names(SensCI95_right) <- paste0(names(SensCI95_right), "_ci_right")
      names(p_value) <-  paste0(names(p_value), "_p_value")
      tmp_sensitivity_statistic <- as.data.frame(as.list(c(SensMean,
                                                           SensSem,
                                                           SensCI95_left,
                                                           SensCI95_right,
                                                           p_value
                                                           ))) %>% 
        dplyr::mutate(block_id = b)
      scores <- rbind.data.frame(scores, tmp)
      scores_statistics <- rbind.data.frame(
        scores_statistics,
        tmp_sensitivity_statistic
      )
    } else{
      response_one_block <- CorrectBaseLine(
        response_one_block,
        method = correct_baseline
      )
      # Calculate RI
      single_drug_data <- ExtractSingleDrug(response_one_block)
      ri <- as.data.frame(lapply(single_drug_data, CalculateSens))
      colnames(ri) <- sub("conc", "ri_", colnames(ri))
      
      # Calculate IC50 for all drugs
      ic50 <-lapply(
        single_drug_data,
        function(x) {
          model <- FitDoseResponse(x)
          coe <- FindModelPar(model)
          type <- FindModelType(model)
          CalculateIC50(coe, type, max(x$dose))
        }) %>% 
        as.data.frame()
      colnames(ic50) <- sub("ri", "ic50", colnames(ri))
      
      # Calculate CSS
      css <- CalculateCSS(response_one_block, ic50 = ic50)
      
      # Assemble data frame
      tmp <- cbind.data.frame(ic50, ri, css) %>% 
        dplyr::mutate(block_id = b)
      scores <- rbind.data.frame(scores, tmp)
    }
  }
  ## 4. Save data into the list
  data$drug_pairs <- data$drug_pairs %>% 
    dplyr::select(
      block_id, 
      setdiff(colnames(data$drug_pairs), colnames(scores))
    ) %>% 
    dplyr::left_join(scores, by = "block_id")
  if (length(scores_statistics) != 0) {
    data$sensitivity_scores_statistics <- dplyr::select(
      scores_statistics,
      block_id,
      dplyr::everything()
    )
  }
  return(data)
}

#' Calculate Combination Sensitivity Score
#' 
#' This function will calculate the Combination Sensitivity Score (CSS) for a
#' drug combination block.
#'
#' @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 ic50 A list. It contains the relative IC50 for all the drugs in the
#'   combination.
#'
#' @return A data frame. It contains the CSS for each drug will one of the other
#'   drugs is at its IC50 and summarized CSS for the whole block.
#'   
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#' 
#' @export
CalculateCSS <- function(response, ic50) {
  concs <- grep("conc\\d", colnames(response), value = TRUE)
  # Calculate CSS for all drug pairs
  conc_zero <- apply(
    dplyr::select(response, - response), 
    1,
    function(x) {
      sum(x != 0) <= 2
    })
  response_two_drugs <- response[conc_zero, ]
  
  conc_pairs <- utils::combn(concs, 2)
  
  css <- vector("list", 2 * ncol(conc_pairs))
  names(css) <- apply(
    conc_pairs, 2, 
    function(x){
      c(paste0(sub("conc", "css", x[1]), sub("conc", "_ic50", x[2])),
        paste0(sub("conc", "css", x[2]), sub("conc", "_ic50", x[1])))
    }
  )
  for (i in 1:ncol(conc_pairs)) {
    pair <- conc_pairs[, i]
    # css_c <- pair[1]
    # ic50_c <- pair[2]
    other_concs <- setdiff(concs, pair)
    if (length(other_concs) > 0) {
      other_concs <- response_two_drugs %>% 
        dplyr::ungroup() %>% 
        dplyr::select(dplyr::all_of(other_concs)) %>% 
        rowSums()
      response_pair <- response_two_drugs[other_concs == 0, ] %>% 
        dplyr::select(dplyr::all_of(pair), response)
    } else {
      response_pair <- response_two_drugs
    }
    for (css_c in pair) {
      ic50_c <- setdiff(pair, css_c)
      tmp_css <- response_pair %>% 
        dplyr::rename(dose = !!ic50_c) %>% 
        dplyr::group_by(!!as.name(css_c)) %>% 
        tidyr::nest(data = dplyr::all_of(c("dose", "response"))) %>% 
        dplyr::mutate(response = furrr::future_map(data, function(x){
          if (nrow(x) == 2) {
            return(x$response[2])
          } else {
            return(PredictResponse(x, dose = ic50[[sub("conc", "ic50_", ic50_c)]]))
          }
        })) %>% 
        dplyr::select(dose = !!as.name(css_c), response) %>% 
        tidyr::unnest(cols = c(response)) %>% 
        CalculateSens()
      css_name <- paste0(
        sub("conc", "css", css_c),
        sub("conc", "_ic50", ic50_c)
      )
      css[[css_name]] <- tmp_css
    }
  }
  css <- as.data.frame(css)
  css$css <- rowMeans(css)
  return(css)
}

#' Calculate Relative Inhibition (RI) for Dose-Response Curve
#'
#' Function \code{CalculateSens} calculates cell line sensitivity to a drug or a
#' combination of drugs from dose response curve.
#'
#' This function measures the sensitivity by calculating the Area Under Curve
#' (AUC) according to the dose response curve. The lower border is chosen as
#' lowest non-zero concentration in the dose response data.
#'
#' @param df A data frame. It contains two variables:
#' \itemize{
#'   \item \strong{dose} the concentrations of drugs.
#'   \item \strong{response} the response of cell lines at corresponding doses.
#'   We use inhibition rate of cell line growth to measure the response.
#' }
#' 
#' @return A numeric value. It is the RI score for input dose-response curve.
#'
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#'
#' @export
#'
#' @examples
#' # LL.4
#' df <- data.frame(dose = c(0, 0.1954, 0.7812, 3.125, 12.5, 50),
#'                  response = c(2.95, 3.76, 18.13, 28.69, 46.66, 58.82))
#' RI <- CalculateSens(df)

CalculateSens <- function(df) {
  df <- df[order(df$dose), ]
  df <- df[which(df$dose != 0),]
  if (nrow(df) == 1) {
    score <- df$response[1]
    res <- score
  } else {
    tryCatch({
      model <- drc::drm(response ~ dose, data = df, fct = drc::LL.4(),
                        control = drc::drmc(errorm = FALSE, noMessage = TRUE,
                                            otrace = TRUE))
      fitcoefs <- model$coefficients
      names(fitcoefs) <- NULL
      # Calculate DSS
      score <- round(.ScoreCurve(d = fitcoefs[3] / 100,
                                 c = fitcoefs[2] / 100,
                                 b = fitcoefs[1],
                                 m = log10(fitcoefs[4]),
                                 c1 = log10(min(df$dose)),
                                 c2 = log10(max(df$dose)),
                                 t = 0), 3)
    }, error = function(e) {
      # Skip zero conc, log, drc::L.4()
      # message(e)
      model <<- drc::drm(response ~ log10(dose), data = df, fct = drc::L.4(),
                         control = drc::drmc(errorm = FALSE, noMessage = TRUE,
                                             otrace = FALSE))
      fitcoefs <- model$coefficients
      names(fitcoefs) <- NULL
      score <<- round(.ScoreCurve_L4(d = fitcoefs[3] / 100,
                                     c = fitcoefs[2] / 100,
                                     b = fitcoefs[1],
                                     e = fitcoefs[4],
                                     c1 = log10(min(df$dose)),
                                     c2 = log10(max(df$dose)),
                                     t = 0), 3)
    })
    res <- score
  }
  return(res)
}

#' Impute Missing Value at IC50 Concentration of Drug
#'
#' \code{ImputeIC50} uses the particular experiment's values to predict the
#' missing values at the desired IC50 concentration of the drug.
#
#' This function is only called when trying to fix a drug at its selected IC50
#' concentration where the response values have not been tested in experiment.
#'
#' \code{ImputeIC50} fits dose-response models (with \code{\link[drc]{drm}}
#' function) by fixing the concentrations of the
#' \strong{other} drug successively, and uses each fit to predict the missing
#' value at the combination (missing IC50, fixed conc).
#'
#' @param response.mat A matrix. It contains response value of a block of drug
#'   combination.
#' @param row.ic50 A numeric value. The IC50 value of drug added to rows.
#' @param col.ic50 A numeric value. The IC50 value of drug added to columns.
#'
#' @return A data frame contains all response value at the IC50 concentration
#'   of certein drug. It could be directly passed to function
#'   \code{CalculateSens} for scoring.
#'
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#'
#' @export
ImputeIC50 <- function(response.mat, col.ic50, row.ic50) {
  
  colconc <- as.numeric(colnames(response.mat))
  rowconc <- as.numeric(rownames(response.mat))
  n_col <- length(colconc)
  n_row <- length(rowconc)
  
  if (n_row == 2) {
    tempcf_c <- data.frame(dose = colconc, response = response.mat[2, ])
  } else {
    response <- apply(response.mat, 2, function(x){
      df <- data.frame(dose = rowconc, response = x)
      pred <- PredictResponse(df, row.ic50)
      return(pred)
    }
    )
    tempcf_c <- data.frame(dose = colconc, response = response)
  }
  
  if (n_col == 2) {
    tempcf_r <- data.frame(dose = rowconc, response = response.mat[, 2])
  } else {
    response <- apply(response.mat, 1, function(x){
      df <- data.frame(dose = colconc, response = x)
      pred <- PredictResponse(df, col.ic50)
      return(pred)
    }
    )
    tempcf_r <- data.frame(dose = rowconc, response = response)
  }
  
  tempres <- list(tempcf_c = tempcf_c, tempcf_r = tempcf_r)
  return(tempres)
}

#' Calculate Relative IC50 from Fitted Model
#' 
#' This function will calculate the relative IC50 from fitted 4-parameter 
#' log-logistic dose response model.
#'
#' @param coef A numeric vector. It contains the fitted coefficients for
#'   4-parameter log-logistic dose response model.
#' @param type A character value. It indicates the type of model was used for
#'   fitting the dose-response curve. Available values are "L.4" and "LL.4".
#' @param max.conc A numeric value. It indicates the maximum concentration in
#'   the dose-response data
#'
#' @return A numeric value. It is the relative IC50.
#' 
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#' 
#' @export
#'
CalculateIC50 <- function(coef, type, max.conc){
  if (type == "LL.4") {
    ic50 <- coef[["e_EC50"]]
  } else if (type == "L.4") {
    ic50 <- exp(coef[["e_log(EC50)"]])
  } else {
    stop("The input 'type = ", type,
         "' is not available. The available values are 'L.4' and 'LL.4'")
  }
  
  if (ic50 > max.conc) {
    ic50 = max.conc
  }
  
  return (ic50)
}

#' Predict Response Value at Certain Drug Dose
#'
#' \code{PredictResponse} uses \code{\link[drc]{drm}} function to fit the dose
#' response model and generate the predict response value at the given dose.
#'
#' \strong{Note}: Random number generator used in \code{AddNoise} with
#' \code{method = "random"}. If the analysis requires for reproductiblity,
#' plesase set the random seed before calling this function.
#'
#' @param df A data frame. It contains two variable:
#' \itemize{
#'   \item \strong{dose} a serial of concentration of drug;
#'   \item \strong{response} the cell line response to each concentration of
#'   drug. It should be the inhibition rate according to negative control.
#' }
#'
#' @param dose A numeric value. It specifies the dose at which user want to
#'   predict the response of cell line to the drug.
#'
#' @return A numeric value. It is the response value of cell line to the drug at
#'  inputted dose.
#'
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#'
#' @export
PredictResponse <- function(df, dose) {
  if (stats::var(df$response, na.rm = TRUE) == 0) {
    pred <- df$response[1]
  } else {
    model <- FitDoseResponse(df)
    if (model$call$fct[[1]][[3]] == "LL.4") {
      pred <- stats::predict(model, data.frame(dose = dose))
    } else if(model$call$fct[[1]][[3]] == "L.4") {
      pred <- stats::predict(model, data.frame(dose = log(dose)))# NB! use log
    } else {
      stop("Fitted model should be LL.4 or L.4.")
    }
    
    if (pred > 100) {
      pred <- 100 + stats::runif(1, -0.01, 0)
    }
  }
  return(pred)
}

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

#' CSS Facilitate Function - .ScoreCurve for Curves Fitted by LL.4 Model
#'
#' New function used to score sensitivities given either a single-agent or a
#' fixed conc (combination) columns. The function calculates the AUC of the
#' log10-scaled dose-response curve. \strong{IMPORTANT:} note that with
#' \code{\link[drc]{LL.4}} calls, this value is already logged since the
#' input concentrations are logged.
#'
#' @param b A numeric value, fitted parameter b from \code{\link[drc]{LL.4}}
#'   model.
#' @param c A numeric value, fitted parameter c from \code{\link[drc]{LL.4}}
#'   model.
#' @param d A numeric value, fitted parameter d from \code{\link[drc]{LL.4}}
#'   model.
#' @param m A numeric value, relative IC50 for the curve. log10(e), where e is
#'   the fitted parameter e from \code{\link[drc]{LL.4}} model.
#' @param c1 A numeric value, log10(min conc) (this is the minimal nonzero
#'   concentration).
#' @param c2 A numeric value, log10(max conc) (this is the maximal
#'   concentration).
#' @param t A numeric value, threshold (usually set to zero).
#'
#' @return A numeric value, RI or CSS scores.
#'
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#'
.ScoreCurve <- function(b, c, d, m, c1, c2, t) {
  int_y <- (((((d - c) * .own_log(-b, c2, m)) / ((-b) * log(10))) + c * c2) -
              ((((d - c) * .own_log(-b, c1, m)) / ((-b) * log(10))) + c * c1))
  
  ratio <- int_y / ((1 - t) * (c2 - c1))
  sens <- ratio * 100 # scale by 100
  return(sens)
}

#' CSS Facilitate Function - Log Calculation (nature based) LL.4 Model
#'
#' #' This function calculates ln(1+10^(b*(c-x))) to be used in \code{.ScoreCurve}
#' function
#'
#' @param b A numeric value. It is the fitted parameter b from 
#'   \code{\link[drc]{L.4}} model.
#' @param c A numeric value. It is the fitted parameter c from 
#'   \code{\link[drc]{L.4}} model.
#' @param x A numeric value. It is the relative IC50 for the curve. log10(e), 
#'   where e is the fitted parameter e from \code{\link[drc]{L.4}} model.
#'
#' @return ln(1+10^(b*(c-x)))
#'
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#'
.own_log = function(b, c, x)
{
  arg = 1 + 10^(b*(c-x))
  if(is.infinite(arg)==T) res = b*(c-x)*log(10) else res = log(arg)
  return(res)
}

#' CSS Facilitate Function - .ScoreCurve for Curves Fitted by L.4 Model
#'
#' This function is used to score sensitivities given either a single-agent or a
#' fixed conc (combination) columns. The function calculates the AUC of the
#' log10-scaled dose-response curve.
#'
#' @param b A numeric value, fitted parameter b from \code{\link[drc]{L.4}}
#'   model.
#' @param c A numeric value, fitted parameter c from \code{\link[drc]{L.4}}
#'   model.
#' @param d A numeric value, fitted parameter d from \code{\link[drc]{L.4}}
#'   model.
#' @param e A numeric value, fitted parameter e from \code{\link[drc]{L.4}}
#'   model.
#' @param c1 A numeric value, log10(min conc) (this is the minimal nonzero
#'   concentration).
#' @param c2 A numeric value, log10(max conc) (this is the maximal
#'   concentration).
#' @param t A numeric value, threshold (usually set to zero).
#'
#' @return A numeric value, RI or CSS scores.
#'
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#'
.ScoreCurve_L4 <- function(b, c, d, e, c1, c2, t) {
  int_y <- d * (c2 - c1) + ((c - d) / b) *
    (.own_log2(b * (c2 - e)) - .own_log2(b * (c1 - e)))
  ratio <- int_y / ((1 - t) * (c2 - c1))
  sens <- ratio * 100 # scale by 100
  return(sens)
}

#' CSS Facilitate Function - Log (nature based) Calculation L.4 Model
#'
#' This function calculates ln(1+exp(x)) to be used in \link{.ScoreCurve_L4}
#' function
#'
#' @param x A numeric value. It is relative IC50 for the curve. The fitted
#'   parameter e from \code{\link[drc]{L.4}} model.
#'
#' @return A numeric value. It is ln(1+exp(x))
#'
#' @author
#' \itemize{
#'   \item Shuyu Zheng \email{shuyu.zheng@helsinki.fi}
#'   \item Jing Tang \email{jing.tang@helsinki.fi}
#' }
#'
.own_log2 <- function(x){
  arg = 1 + exp(x)
  if(is.infinite(arg)==T) res = x else res = log(arg)
  return(res)
}