R/forest_plot.R

In volcano3D: 3D Volcano Plots and Polar Plots for Three-Class Data

#' Forest plot individual gene from 2x3 factor analysis
#' 
#' Forest plot individual gene from 2x3 factor analysis using either base
#' graphics, plotly or ggplot2.
#' 
#' @param object A 'volc3d' class object from a 2x3 analysis generated by
#'   [deseq_2x3_polar()]
#' @param genes Vector of genes to plot
#' @param scheme Vector of 3 colours for plotting
#' @param labs Optional character vector of labels for the groups
#' @param error_type Either "ci" or "se" to specify whether error bars use 95%
#'   confidence intervals or standard error
#' @param error_width Width of error bars
#' @param facet Logical whether to use facets for individual genes (ggplot2 
#' only)
#' @param gap Size of gap between groupings for each gene
#' @param transpose Logical whether to transpose the plot
#' @param mar Vector of margins on four sides. See [par()]
#' @param ... Optional arguments
#' @return Returns a plot using either base graphics (`forest_plot`), plotly
#'   (`forest_plotly`) or ggplot2 (`forest_ggplot`). `forest_plot` also
#'   invisibly returns the dataframe used for plotting.
#' @seealso [deseq_2x3_polar()]
#' @importFrom graphics abline arrows axis legend par points text
#' @export

forest_plot <- function(object, genes,
                        scheme = c('red', 'green3', 'blue'),
                        labs = NULL,
                        error_type = c("ci", "se"),
                        error_width = 0.05,
                        gap = 1,
                        transpose = FALSE,
                        mar = if (transpose) c(5, 7, 5, 4) else c(5, 5, 5, 3),
                        ...) {
  error_type <- match.arg(error_type)
  data <- forest_df(object, genes, labs, error_type, gap)
  xrange <- range(c(data$val-data$CI, data$val+data$CI, 0), na.rm=TRUE)
  prange <- range(data$pos)
  op <- par(mar=mar)
  on.exit(par(op))
  new.args <- list(...)
  if (transpose) {
    plot.args <- list(x = NA, xlim = xrange, ylim = prange, yaxt="n", bty = "n",
                      xlab = bquote("log"[2]~"FC"), ylab = "")
    if (length(new.args)) plot.args[names(new.args)] <- new.args
    do.call("plot", plot.args)
    arrows(data$val-data$CI, data$pos, data$val+data$CI, code=3, angle=90,
           length=error_width, col=scheme)
    points(data$val, y = data$pos, pch=19, col=scheme)
    abline(v = 0, lty = 2)
    axis(2, data$pos, data$labs, tick = FALSE, las = 1)
    axis(2, data$pos[seq_along(genes)*3-1], genes, tick = FALSE, line=2, las=1)
    text(xrange[2]+ diff(xrange)*0.06, data$pos, data$stars, xpd = NA)
  } else {
    plot.args <- list(x = NA, ylim = xrange, xlim = prange, xaxt="n", bty = "n",
                      ylab = bquote("log"[2]~"FC"), xlab = "", las = 1)
    if (length(new.args)) plot.args[names(new.args)] <- new.args
    do.call("plot", plot.args)
    arrows(data$pos, data$val-data$CI, y1 = data$val+data$CI, code=3, angle=90,
           length=error_width, col=scheme)
    points(data$pos, data$val, pch=19, col=scheme)
    abline(h = 0, lty = 2)
    axis(1, data$pos, data$labs, tick = FALSE, las = 1)
    axis(1, data$pos[seq_along(genes)*3-1], genes, tick = FALSE, line=2)
    text(data$pos, xrange[2]+ diff(xrange)*0.06, data$stars, xpd = NA)
  }
  legend("topright", legend=colnames(object@df[[2]])[1:3], pch=19, col=scheme,
         bty="n", cex=0.8, inset = c(0,-0.25), xpd = NA)
  invisible(data)
}


#' @rdname forest_plot
#' @export
forest_plotly <- function(object, genes,
                          scheme = c('red', 'green3', 'blue'),
                          labs = NULL,
                          error_type = c("ci", "se"),
                          error_width = 4,
                          gap = 1,
                          transpose = FALSE, ...) {
  error_type <- match.arg(error_type)
  data <- forest_df(object, genes, labs, error_type, gap) 
  if (transpose) {
    annot <- list(y = data$pos, x = max(data$val+data$CI, na.rm = TRUE),
                  text = data$stars, showarrow = FALSE, xanchor='left')
    plot_ly(data = data, x = ~val, y = ~pos,
            colors = scheme, color = ~labs,
            type = 'scatter', mode = 'markers',
            error_x = list(array = ~CI, color = ~labs, width = error_width),
            ...) %>%
      layout(xaxis = list(title = "log<sub>2</sub> FC"),
             yaxis = list(title = "",
                          ticktext = genes,
                          tickvals = data$pos[seq_along(genes)*3-1]),
             annotations = annot)
  } else {
    annot <- list(x = data$pos, y = max(data$val+data$CI, na.rm = TRUE),
                  text = data$stars, showarrow = FALSE, yanchor='bottom')
    plot_ly(data = data, y = ~val, x = ~pos,
            colors = scheme, color = ~labs,
            type = 'scatter', mode = 'markers',
            error_y = list(array = ~CI, color = ~labs, width = error_width),
            ...) %>%
      layout(yaxis = list(title = "log<sub>2</sub> FC"),
             xaxis = list(title = "",
                          ticktext = genes,
                          tickvals = data$pos[seq_along(genes)*3-1]),
             annotations = annot)
  }
}


#' @rdname forest_plot
#' @importFrom ggplot2 geom_errorbar geom_vline geom_hline scale_x_continuous
#'   scale_y_continuous theme_minimal xlab ylab facet_grid expand_limits
#'   theme_bw element_line
#' @importFrom rlang .data
#' @export
forest_ggplot <- function(object, genes,
                          scheme = c('red', 'green3', 'blue'),
                          labs = NULL,
                          error_type = c("ci", "se"),
                          error_width = 0.3,
                          facet = TRUE,
                          gap = 1,
                          transpose = FALSE, ...) {
  error_type <- match.arg(error_type)
  data <- forest_df(object, genes, labs, error_type, gap)
  vdiff <- diff(range(c(data$val-data$CI, data$val+data$CI, 0), na.rm=TRUE))
  
  if (facet) {
    if (transpose) {
      xmax <- max(data$val + data$CI, na.rm=TRUE)
      ggplot(data, aes(y=.data$labs, x=.data$val, color=.data$labs)) +
        geom_vline(xintercept = 0, linetype = "dashed", colour = "grey40") +
        geom_errorbar(aes(xmin=.data$val-.data$CI, xmax=.data$val+.data$CI),
                      width=error_width) +
        geom_point() +
        scale_color_manual(values = scheme) +
        facet_grid(gene ~ .) +
        xlab(bquote("log"[2]~"FC")) + ylab("") + labs(color = "") +
        geom_text(y = data$labs,
                  x = xmax +vdiff*0.05,
                  colour = "black",
                  label = data$stars, show.legend = FALSE) +
        expand_limits(x= xmax + vdiff * 0.04) +
        theme_bw() +
        theme(axis.text = element_text(colour = "black"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.spacing = unit(0, "cm"),
              strip.background = element_rect(fill="grey90"))
    } else {
      # not transposed
      ymax <- max(data$val + data$CI, na.rm=TRUE)
      ggplot(data, aes(x=.data$labs, y=.data$val, color=.data$labs)) +
        geom_hline(yintercept = 0, linetype = "dashed", colour = "grey40") +
        geom_errorbar(aes(ymin=.data$val-.data$CI, ymax=.data$val+.data$CI),
                      width=error_width) +
        geom_point() +
        scale_color_manual(values = scheme) +
        facet_grid(. ~ gene) +
        ylab(bquote("log"[2]~"FC")) + xlab("") + labs(color = "") +
        geom_text(x = data$labs,
                  y = ymax +vdiff*0.04,
                  colour = "black",
                  label = data$stars, show.legend = FALSE) +
        expand_limits(y= ymax + vdiff * 0.04) +
        theme_bw() +
        theme(axis.text = element_text(colour = "black"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.spacing = unit(0, "cm"),
              strip.background = element_rect(colour=NA, fill=NA))
    }
  } else {
    # single plot, no facets
    if (transpose) {
      ggplot(data, aes(y=.data$pos, x=.data$val, color=.data$labs)) +
        geom_errorbar(aes(xmin=.data$val-.data$CI, xmax=.data$val+.data$CI),
                      width=error_width) +
        geom_point() +
        scale_color_manual(values = scheme) +
        xlab(bquote("log"[2]~"FC")) + ylab("") + labs(color = "") +
        scale_y_continuous(breaks=data$pos[seq_along(genes)*3-1], 
                           labels=genes) +
        geom_vline(xintercept = 0, linetype = "dashed") +
        annotate(geom = "text",
                 y = data$pos,
                 x = max(data$val + data$CI, na.rm=TRUE) +vdiff*0.05,
                 label = data$stars) +
        theme_minimal() +
        theme(axis.text = element_text(colour = "black"),
              axis.line.x = element_line(color="black"))
    } else {
      ggplot(data, aes(x=.data$pos, y=.data$val, color=.data$labs)) +
        geom_errorbar(aes(ymin=.data$val-.data$CI, ymax=.data$val+.data$CI),
                      width=error_width) +
        geom_point() +
        scale_color_manual(values = scheme) +
        ylab(bquote("log"[2]~"FC")) + xlab("") + labs(color = "") +
        scale_x_continuous(breaks=data$pos[seq_along(genes)*3-1], 
                           labels=genes) +
        geom_hline(yintercept = 0, linetype = "dashed") +
        annotate(geom = "text",
                 x = data$pos,
                 y = max(data$val + data$CI, na.rm=TRUE) +vdiff*0.03,
                 label = data$stars) +
        theme_minimal() +
        theme(axis.text = element_text(colour = "black"),
              axis.line.y = element_line(color="black"))
    }
  }
}


forest_df <- function(object, genes,
                      labs = NULL, error_type = "ci", gap = 1) {
  if(! is(object, "volc3d")) stop("Not a 'volc3d' class object")
  if(!object@df$type %in% c("polar_coords_2x3", "deseq_2x3_polar")) {
    stop("Not a 2x3-way analysis")
  }
  df <- object@df[[2]]
  val <- as.vector(t(df[genes, 1:3]))
  CI <- as.vector(t(df[genes, 4:6]))
  if (error_type == "ci") CI <- CI * 1.96
  ngenes <- length(genes)
  pos <- rep.int(1:3, ngenes) + 
    rep(seq(from=0, by=3+gap, length.out=ngenes), each = 3)
  if (is.null(labs)) labs <- abbreviate(colnames(df)[1:3], 3)
  pval <- as.vector(t(object@padj[genes,]))
  data <- data.frame(gene = factor(rep(genes, each = 3), levels = genes),
                     labs = factor(rep_len(1:3, 3*ngenes),
                                   labels = labs),
                     pos = pos,
                     val = val,
                     CI = CI,
                     pval = pval,
                     stars = stars.pval(pval))
  data
}


# code modified from orphaned package gtools

#' @importFrom stats symnum
stars.pval <- function(p.value) {
  unclass(
    symnum(p.value,
           corr = FALSE, na = FALSE,
           cutpoints = c(0, 0.001, 0.01, 0.05, 1),
           symbols = c("***", "**", "*", "")
    )
  )
}