R/GRbox.R

Defines functions GRbox

Documented in GRbox

#' Boxplots of a given GR metric
#'
#' Given a SummarizedExperiment object created by \code{\link{GRfit}},
#' this function creates boxplots according to the parameters below.
#'
#' @param fitData a SummarizedExperiment object, generated by the GRfit 
#' function.
#' @param metric the GR metric (GR50, GRinf, h_GR, GRmax, GEC50, or GR_AOC) 
#' or traditional metric (IC50, Einf, h, Emax, EC50, or AUC) that will be used
#' for the boxplot.
#' @param groupVariable the name of the variable from data (e.g. drug,
#' cell-line, etc.) to select factors from.
#' @param pointColor a variable that defines the coloring of the points
#' overlayed on the boxplot.
#' @param factors a vector of values of "groupVariable" of data that define
#' which variables to make boxplots for. By default, a separate boxplot is made
#' for each unique value of groupVariable.
#' @param wilA one value or a vector of values from "factors", i.e. 
#' a subset of the boxplots. If specified, a one-sided Wilcoxon rank sum test 
#' (wilcox.test) will be performed between "wilA" and "wilB" and 
#' the results will be displayed on the figure. The null hypothesis 
#' that the values from "wilA" and "wilB" have the same mean will 
#' be tested against the alternative hypothesis that the mean of the 
#' "wilB" values is greater than that of the "wilA" values.
#' @param wilB one value or a vector of values from "factors", i.e. 
#' a subset of the boxplots (not overlapping "wilA").
#' @param plotly a logical value indicating whether to output a ggplot2 graph
#' or an interactive ggplotly graph
#'
#' @return ggplot2 or ggplotly boxplots of the factors along the x-axis, with
#' points colored by the given variable.
#' @author Nicholas Clark
#' @details
#' Given a SummarizedExperiment object created by \code{\link{GRfit}},
#' this function creates boxplots of a given GR metric (GR50, GRmax, etc.) or 
#' traditional metric (IC50, Emax, etc.) 
#' for values of the grouping variable. The results can be viewed in a static
#' ggplot image or an interactive plotly graph.
#'
#' By default, a boxplot is created for all unique
#' values of the grouping variable. The "factors" parameter can be used to
#' specify a smaller subset of values for which to create boxplots.
#' Points are overlayed on the boxplots and
#' they can be colored by the variable specified in the pointColor parameter.
#' If pointColor is set to NULL, the points will all be black. The results can
#' be viewed in a static ggplot image or an interactive plotly graph.
#' @seealso To create the object needed for this function, see
#' \code{\link{GRfit}}. For other visualizations, see \code{\link{GRdrawDRC}}
#' and \code{\link{GRscatter}}. For online GR calculator and browser, see
#' \url{http://www.grcalculator.org}.
#' @examples
#' # Load Case A (example 1) input
#' data("inputCaseA")
#' head(inputCaseA)
#' # Run GRfit function with case = "A"
#' output1 = GRfit(inputData = inputCaseA,
#' groupingVariables = c('cell_line','agent', 'perturbation','replicate',
#' 'time'))
#' GRbox(output1, metric ='GRinf',
#' groupVariable = 'cell_line', pointColor = 'agent' , factors = c('BT20',
#' 'MCF10A'))
#' GRbox(output1, metric ='GRinf',
#' groupVariable = 'cell_line', pointColor = 'cell_line' ,
#' factors = c('BT20', 'MCF10A'), plotly = FALSE)
#' GRbox(output1, metric = 'GR50', groupVariable = 'cell_line',
#' pointColor = 'cell_line', wilA = "BT20", wilB = c("MCF7","MCF10A"))
#' @export

GRbox <- function(fitData, metric, groupVariable, pointColor,
                  factors = "all", wilA = NULL, wilB = NULL, plotly = TRUE) {
  data = cbind(as.data.frame(SummarizedExperiment::colData(fitData)),
               t(SummarizedExperiment::assay(fitData)))
  #bottom_margin = max(nchar(data[[groupVariable]]), na.rm = TRUE)
  data[[groupVariable]] = factor(data[[groupVariable]])
  if(!identical(factors, "all")) {
    if(length(intersect(factors, data[[groupVariable]])) != length(factors)) {
      stop('Factors must be values of the grouping variable')
    }
    data = data[data[[groupVariable]] %in% factors, ]
    #bottom_margin = max(nchar(factors), na.rm = TRUE)
  }
  if(length(intersect(wilA, wilB)) > 0) {
    stop('wilA and wilB must not overlap.')
  }
  if(metric == "GR50") {
    data$`log10(GR50)` = log10(data$GR50)
    metric = "log10(GR50)"
  }
  if(metric == "IC50") {
    data$`log10(IC50)` = log10(data$IC50)
    metric = "log10(IC50)"
  }
  if(metric == "h_GR") {
    data$`log2(h_GR)` = log2(data$h_GR)
    metric = "log2(h_GR)"
  }
  if(metric == "h") {
    data$`log2(h)` = log2(data$h)
    metric = "log2(h)"
  }
  # Get rid of infinite values
  fin = is.finite(data[[metric]])
  data = data[fin,]
  if(!is.null(wilA) & !is.null(wilB)) {
    for(i in 1:length(wilB)) {
      data[[groupVariable]] = stats::relevel(data[[groupVariable]], wilB[i])
    }
    for(i in 1:length(wilA)) {
      data[[groupVariable]] = stats::relevel(data[[groupVariable]], wilA[i])
    }
    # Perform Wilcoxon rank sum test
    rowsA = data[[groupVariable]] %in% wilA
    rowsB = data[[groupVariable]] %in% wilB
    wil_dataA = data[rowsA,metric]
    wil_dataB = data[rowsB,metric]
    wil = stats::wilcox.test(x = wil_dataA, y = wil_dataB,
                             alternative = "less")
    wil_pval = prettyNum(wil$p.value, digits = 2)
  }
  if(plotly == TRUE) {
    p <- ggplot2::ggplot(data, ggplot2::aes_string(x = groupVariable,
                          y = metric, text = 'experiment'))
    p = p + ggplot2::geom_boxplot(ggplot2::aes_string(fill = groupVariable,
                  alpha = 0.3), outlier.color = NA, show.legend = FALSE) +
      ggplot2::geom_jitter(width = 0.5, show.legend = FALSE,
      ggplot2::aes_string(colour = pointColor)) +
      ggplot2::xlab('') + ggplot2::ylab(metric)
    q = plotly::plotly_build(p)
    # Last CRAN version of plotly (3.6.0) uses "q$"
    # Latest github version of plotly (4.3.5) uses "q$x"
    if(is.null(q$data)) {
      # replace q with q$x so code works with new plotly version
      q = q$x
    }
    if(!is.null(wilA) & !is.null(wilB)) {
      top_y = q[[2]]$yaxis$range[2]
      bottom_y = q[[2]]$yaxis$range[1]
      total_y_range = top_y - bottom_y
      # Get top of boxplot whiskers
      whiskers = NULL
      len = length(wilA) + length(wilB)
      for(i in 1:len) {
        whiskers[i] = stats::fivenum(q[[1]][[i]]$y)[5]
      }
      top_whisker = max(whiskers, na.rm = TRUE)
      y_range = (top_y - top_whisker)/total_y_range
      if(y_range < .25) {
        top_y = top_whisker + .25*total_y_range
      }
      lh = top_whisker + total_y_range*(.1)
      bump = total_y_range*(.05)
      ll = lh - bump
      lenA = length(wilA)
      lenB = length(wilB)
      pval = paste("p =", wil_pval)
      if(lenA == 1 & lenB == 1) {
        p = p + ggplot2::annotate("text", x = 1.5, y = lh + bump/2, 
                                  label = pval) + 
          ggplot2::geom_segment(x = 1, y = lh, xend = 2, yend = lh) + 
          ggplot2::geom_segment(x = 1, y = ll, xend = 1, yend = lh) + 
          ggplot2::geom_segment(x = 2, y = ll, xend = 2, yend = lh)
        rm(q)
        q = plotly::plotly_build(p)
      } else if(lenA > 1 & lenB == 1) {
        p = p + ggplot2::annotate("text", x = ((lenA + 1) + ((lenA+1)/2))/2,
                         y = lh + 2*bump, label = pval) +
          ggplot2::geom_segment(x = 1, y = lh, xend = lenA, yend = lh) + 
          ggplot2::geom_segment(x = 1, y = ll, xend = 1, yend = lh) + 
          ggplot2::geom_segment(x = lenA, y = ll, xend = lenA, yend = lh) +
          ggplot2::geom_segment(x = (lenA+1)/2, y = lh + bump, xend = lenA + 1,
                       yend = lh + bump) + 
          ggplot2::geom_segment(x = (lenA+1)/2, y = lh, xend = (lenA+1)/2,
                       yend = lh + bump) + 
          ggplot2::geom_segment(x = lenA+1, y = ll, xend = lenA+1,
                       yend = lh + bump)
        rm(q)
        q = plotly::plotly_build(p)
      } else if(lenA == 1 & lenB > 1) {
        p = p + ggplot2::annotate("text", x = 1.25 + .25*lenB, 
                y = lh + 2*bump, label = pval) +
          ggplot2::geom_segment(x = 1, y = lh+bump, xend = .5*lenB + 1.5,
                                yend = lh+bump) + 
          ggplot2::geom_segment(x = 1, y = ll, xend = 1, yend = lh+bump) + 
          ggplot2::geom_segment(x = 1.5+.5*lenB, y = lh, xend = 1.5+.5*lenB, 
                                yend = lh+bump) +
          ggplot2::geom_segment(x = 2, y = lh, xend = lenB + 1, yend = lh) + 
          ggplot2::geom_segment(x = 2, y = ll, xend = 2, yend = lh) + 
          ggplot2::geom_segment(x = lenB+1, y = ll, xend = lenB+1, yend = lh)
        rm(q)
        q = plotly::plotly_build(p)
      } else if(lenA > 1 & lenB > 1) {
        p = p + ggplot2::annotate("text", x = .25*(lenB-1)+.75*(lenA+1), 
                y = lh + 2*bump, label = pval) +
          ggplot2::geom_segment(x = 1, y = lh, xend = lenA, yend = lh) + 
          ggplot2::geom_segment(x = 1, y = ll, xend = 1, yend = lh) + 
          ggplot2::geom_segment(x = lenA, y = ll, xend = lenA, yend = lh) +
          ggplot2::geom_segment(x = lenA+1, y = lh, xend = lenA+lenB,
                                yend = lh) + 
          ggplot2::geom_segment(x = lenA+1, y = ll, xend = lenA+1, yend = lh) +
          ggplot2::geom_segment(x = lenA+lenB, y = ll, xend = lenA+lenB,
                                yend = lh) +
          ggplot2::geom_segment(x = (lenA+1)/2, y = lh+bump, 
                                xend = (lenA+1)+((lenB-1)/2), yend = lh+bump) +
          ggplot2::geom_segment(x = (lenA+1)/2, y = lh, xend = (lenA+1)/2, 
                       yend = lh+bump) +
          ggplot2::geom_segment(x = (lenA+1)+((lenB-1)/2), y = lh, 
                                xend = (lenA+1)+((lenB-1)/2), yend = lh+bump)
        rm(q)
        q = plotly::plotly_build(p)
      }
    }
    q = plotly::plotly_build(p)
    # Last CRAN version of plotly (3.6.0) uses "q$"
    # Latest github version of plotly (4.3.5) uses "q$x"
    if(!is.null(q$data)) { # old plotly
      bottom_margin = max(nchar(q$layout$xaxis$ticktext), na.rm = TRUE)
      left = nchar(q$layout$xaxis$ticktext[1])
      q$layout$xaxis$tickangle = -45
      q$layout$margin$b = 15 + 6*bottom_margin
      if(left > 10) {
        left_margin = q$layout$margin$l + (left-10)*6
        q$layout$margin$l = left_margin
      }
      return(q)
    } else { # new plotly
      bottom_margin = max(nchar(q$x$layout$xaxis$ticktext), na.rm = TRUE)
      left = nchar(q$x$layout$xaxis$ticktext[1])
      q$x$layout$xaxis$tickangle = -45
      q$x$layout$margin$b = 15 + 6*bottom_margin
      if(left > 10) {
        left_margin = q$x$layout$margin$l + (left-10)*6
        q$x$layout$margin$l = left_margin
      }
      return(q)
    }
    
  } else {
    p <- ggplot2::ggplot(data, ggplot2::aes_string(x = groupVariable,
                                                   y = metric))
    p = p + ggplot2::geom_boxplot(ggplot2::aes_string(
      fill = groupVariable, alpha = 0.3), outlier.color = NA,
      show.legend = FALSE) + ggplot2::geom_jitter(
        width = 0.5, ggplot2::aes_string(colour = pointColor)) +
      ggplot2::xlab('') + ggplot2::ylab(metric) +
      ggplot2::theme_grey(base_size = 14) +
      ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, 
                                                         vjust = 1, hjust=1))
    if(!is.null(wilA) & !is.null(wilB)) {
      q = plotly::plotly_build(p)
      # Last CRAN version of plotly (3.6.0) uses "q$"
      # Latest github version of plotly (4.3.5) uses "q$x"
      if(is.null(q$data)) {
        q = q$x
      }
      # Get y range:
      top_y = q[[2]]$yaxis$range[2]
      bottom_y = q[[2]]$yaxis$range[1]
      total_y_range = top_y - bottom_y
      # Get top of boxplot whiskers
      whiskers = NULL
      len = length(wilA) + length(wilB)
      for(i in 1:len) {
        whiskers[i] = stats::fivenum(q[[1]][[i]]$y)[5]
      }
      top_whisker = max(whiskers, na.rm = TRUE)
      y_range = (top_y - top_whisker)/total_y_range
      if(y_range < .25) {
        top_y = top_whisker + .25*total_y_range
      }
      lh = top_whisker + total_y_range*(.1)
      bump = total_y_range*(.05)
      ll = lh - bump
      lenA = length(wilA)
      lenB = length(wilB)
      pval = paste("p =", wil_pval)
      if(lenA == 1 & lenB == 1) {
        p = p + ggplot2::annotate("text", x = 1.5, y = lh + bump/2, 
                                  label = pval) + 
          ggplot2::geom_segment(x = 1, y = lh, xend = 2, yend = lh) + 
          ggplot2::geom_segment(x = 1, y = ll, xend = 1, yend = lh) + 
          ggplot2::geom_segment(x = 2, y = ll, xend = 2, yend = lh)
      } else if(lenA > 1 & lenB == 1) {
        p = p + ggplot2::annotate("text", x = ((lenA + 1) + ((lenA+1)/2))/2,
                                  y = lh + 2*bump, label = pval) +
          ggplot2::geom_segment(x = 1, y = lh, xend = lenA, yend = lh) + 
          ggplot2::geom_segment(x = 1, y = ll, xend = 1, yend = lh) + 
          ggplot2::geom_segment(x = lenA, y = ll, xend = lenA, yend = lh) +
          ggplot2::geom_segment(x = (lenA+1)/2, y = lh + bump, xend = lenA + 1,
                                yend = lh + bump) + 
          ggplot2::geom_segment(x = (lenA+1)/2, y = lh, xend = (lenA+1)/2,
                                yend = lh + bump) + 
          ggplot2::geom_segment(x = lenA+1, y = ll, xend = lenA+1,
                                yend = lh + bump)
      } else if(lenA == 1 & lenB > 1) {
        p = p + ggplot2::annotate("text", x = 1.25 + .25*lenB, 
                                  y = lh + 2*bump, label = pval) +
          ggplot2::geom_segment(x = 1, y = lh+bump, xend = .5*lenB + 1.5,
                                yend = lh+bump) + 
          ggplot2::geom_segment(x = 1, y = ll, xend = 1, yend = lh+bump) + 
          ggplot2::geom_segment(x = 1.5+.5*lenB, y = lh, xend = 1.5+.5*lenB, 
                                yend = lh+bump) +
          ggplot2::geom_segment(x = 2, y = lh, xend = lenB + 1, yend = lh) + 
          ggplot2::geom_segment(x = 2, y = ll, xend = 2, yend = lh) + 
          ggplot2::geom_segment(x = lenB+1, y = ll, xend = lenB+1, yend = lh)
      } else if(lenA > 1 & lenB > 1) {
        p = p + ggplot2::annotate("text", x = .25*(lenB-1)+.75*(lenA+1), 
                                  y = lh + 2*bump, label = pval) +
          ggplot2::geom_segment(x = 1, y = lh, xend = lenA, yend = lh) + 
          ggplot2::geom_segment(x = 1, y = ll, xend = 1, yend = lh) + 
          ggplot2::geom_segment(x = lenA, y = ll, xend = lenA, yend = lh) +
          ggplot2::geom_segment(x = lenA+1, y = lh, xend = lenA+lenB,
                                yend = lh) + 
          ggplot2::geom_segment(x = lenA+1, y = ll, xend = lenA+1, yend = lh) +
          ggplot2::geom_segment(x = lenA+lenB, y = ll, xend = lenA+lenB,
                                yend = lh) +
          ggplot2::geom_segment(x = (lenA+1)/2, y = lh+bump, 
                                xend = (lenA+1)+((lenB-1)/2), yend = lh+bump) +
          ggplot2::geom_segment(x = (lenA+1)/2, y = lh, xend = (lenA+1)/2, 
                                yend = lh+bump) +
          ggplot2::geom_segment(x = (lenA+1)+((lenB-1)/2), y = lh, 
                                xend = (lenA+1)+((lenB-1)/2), yend = lh+bump)
      }
    }
    return(p)
  }
}
uc-bd2k/GRmetrics_old documentation built on May 3, 2019, 2:13 p.m.