R/summariseScreen.R

In lujunyan1118/DrugScreenExplorer_dev: A pipeline for processing drug sensitivity screen data

#' Summarise the drug effect across concentrations
#'
#' This function summarises the normalized drug effect (viability) across multiple concentrations into a single value.
#' The summarisation is based on the normalized viability values, and therefore plate normalisation needs to be performed first.
#' If the incubation/edge effect corrected viability values are present, the normalisation will be performed on both uncorrected and corrected viability values and
#' a suffix, '.cor' will be added to the columns that contain the summarised values.
#'
#' @param screenData the data frame containing screening data generated by \code{readScreen()} function
#' @param method a character string or a vector of character strings specifying the methods for summarising the effect. Currently three methods are supported:
#' average, which simply calculates the averaged drug effect across all concentrations;
#' AUC, which calculates the normalized area under dose-response curve using the trapezoidal method;
#' IC50, which will first fits a robust four parameter log-logistic model using dr4pl package and then summarise the drug effect
#' using the four parameters (upper limit, lower limit, slope and IC50/EC50) and area under the sigmoid curve.
#' Multiple methods can be used at the same time by specifying a vector method names.
#' @export
#' @import dr4pl tidyr
#' @return Depends on the user-specified summarisation method, the following columns will be added to the input data frame:
#'
#' for \code{method = 'average'}
#' \item{meanViab}{the averaged viability across all concentrations for a certain drug}
#'
#' for \code{method = 'AUC'}
#' \item{AUC}{the area under the dose-response curve calculated by using the trapezoidal rule.}
#'
#' for \code{method = 'IC50'}
#' \item{UpperLimit}{the upper limit of the fitted logistic model}
#' \item{IC50}{the half maximal inhibitory concentration (IC50), or the half maximal effective concentration (EC50), depends on the type of input data}
#' \item{Slope}{the slope of the fitted logistic model}
#' \item{lowerLimit}{the upper limit of the fitted logistic model}
#' \item{AUC}{the area under the fitted sigmoid doese-response curve}
#' If the logistic model can not be fitted for a certain drug, the above four values will be set to NA.
#'
#' If a column that contains the edge/incubation effect corrected viability values is present in the input data frame,
#' the above columns with a suffix \code{.cor}, which represent the summarised values based
#' on edge/incubation corrected viabilities, will also be added to the input data frame.
#'
#' @examples
#' data(screenData_normalized)
#' screenData <- summariseScreen(screenData_normalized)
#' # Please see the vignette for more information.

summariseScreen <- function(screenData, method = "average") {

  # firstly check if necessary columns are presented
  if (!all(c("name", "concentration", "sampleID") %in% colnames(screenData))) {
    stop("No information of name, concentration or sampleID found")
  }

  if (!"normVal" %in% colnames(screenData)) {
    stop("Column 'normVal' not found. Plate normalization should be performed first")
  }

  ifCor <- "normVal.cor" %in% colnames(screenData)

  # only summarise for non-control wells
  subData <- screenData[screenData$wellType %in% "sample", ]

  if (ifCor) {
    subData <- group_by(subData, sampleID, name, concentration) %>%
      summarise(viab = mean(normVal, na.rm = TRUE), viab.cor = mean(normVal.cor, na.rm = TRUE)) %>%
      ungroup()

  } else {
    subData <- group_by(subData, sampleID, name, concentration) %>%
      summarise(viab = mean(normVal, na.rm = TRUE)) %>%
      ungroup()

  }

  if ("average" %in% method) {
    if (ifCor) {
      sumTab <- group_by(subData, sampleID, name) %>%
        summarise(meanViab = mean(viab, na.rm = TRUE), meanViab.cor = mean(viab.cor)) %>%
        ungroup()

    } else {
      sumTab <- group_by(subData, sampleID, name) %>%
        summarise(meanViab = mean(viab, na.rm = TRUE)) %>%
        ungroup()

    }

    if ("meanViab" %in% colnames(screenData)) screenData[["meanViab"]] <- NULL

    screenData <- left_join(screenData, sumTab, by = c("sampleID", "name"))
  }

  if ("AUC" %in% method) {
    if (ifCor) {
      sumTab <- group_by(subData, sampleID, name) %>%
        summarize(AUC = calcAUC(viab, concentration), AUC.cor = calcAUC(viab.cor, concentration)) %>%
        ungroup()

    } else {
      sumTab <- group_by(subData, sampleID, name) %>%
        summarize(AUC = calcAUC(viab, concentration)) %>%
        ungroup()

    }
    if ("AUC" %in% colnames(screenData)) screenData[["AUC"]] <- NULL

    screenData <- left_join(screenData, sumTab, by = c("sampleID", "name"))
  }

  if ("IC50" %in% method) {

    sumTab <- group_by(subData, sampleID, name) %>% nest() %>%
      dplyr::mutate(model = purrr::map(data, ~sumIC50(viab ~ concentration, .))) %>%
      unnest("model") %>% ungroup()

    sumTab$data <- NULL

    if (ifCor) {
      sumTab.cor <- group_by(subData, sampleID, name) %>% nest() %>%
        dplyr::mutate(model = purrr::map(data, ~sumIC50(viab.cor ~ concentration, .))) %>%
        unnest("model") %>% ungroup()

      sumTab.cor$data <- NULL
      newColName <- colnames(sumTab.cor)
      newColName[!newColName %in% c("sampleID", "name")] <-
        paste0(newColName[!newColName %in% c("sampleID", "name")], ".cor")

      colnames(sumTab.cor) <- newColName
      sumTab <- left_join(sumTab, sumTab.cor, by = c("sampleID", "name"))

    }

    for (colN in c("UpperLimit", "IC50", "Slope", "LowerLimit", "AUC")) {
      if (colN %in% colnames(screenData)) screenData[[colN]] <- NULL
    }

    screenData <- left_join(screenData, sumTab, by = c("sampleID", "name"))
  }

  return(screenData)
}