R/plot_expression.R

In leandroroser/RNASeqFunctions: Several functions for rna-seq analysis

#' generate RLE plot
#' @param x matrix with counts (genes in rows, treatments in columns)
#' @param col_factor factor with classes for treatments
#' (with size equal to ncol(x))
#' @export


plot_expression <- function(x, col_factor,
                            method = c("boxplot", "violin"),
                            y_limit = NULL,
                            filter_zero = TRUE) {
  require(ggplot2)
  method <- match.arg(method)

  x_plot <- t(x)

   if(filter_zero) {
    to_filter <- colSums(x_plot) == 0
    cat("Removing ", sum(to_filter), " features with 0 counts for this analysis\n")
    x_plot <- x_plot[ ,-to_filter , drop =FALSE]
  }

  x_plot <- data.frame(sample = rownames(x_plot), Condition = col_factor, x_plot)
  x_plot <- reshape2::melt(x_plot, ids = c("sample", "Condition"))

  # minmax <- quantile(x_plot[, 4],  c(0.25, 0.75))
  # IQR <- minmax[2] - minmax[1]
  # minmax <- c(minmax[1] - 1.5, minmax[2] + 1.5) * IQR
  # minmax <- c(minmax[1]-minmax[1] / 5, minmax[2] + minmax[2] / 5)

  fun_bp <- function(x) {
    out <- as.numeric(boxplot(x, plot = FALSE)$stats)
    names(out) <- c("ymin", "lower", "middle", "upper", "ymax")
    out
  }

  out <- ggplot2::ggplot(data = x_plot , aes(x=sample, y=value)) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    xlab("") +
    ylab("Expression (counts)")


  if(method == "boxplot") {
    cat("Creating box plot...")
    out <- out  +   stat_boxplot(geom="errorbar") +
      stat_summary(fun.data = fun_bp, geom="boxplot", aes(fill = Condition))


    # geom_boxplot(show.legend = TRUE, outlier.colour = NA, aes(fill = Condition)) +
    # geom_errorbar(aes(ymin = min_plot_stats, ymax = max_plot_stats)) +
    # ylim(minmax[1], minmax[2])
  } else {
    cat("Creating violin plot...")
    out <- out  +    geom_violin(show.legend = TRUE, draw_quantiles = c(0.25, 0.5, 0.75),
                                 scale = "width",
                                 aes(fill = Condition))
  }
  if(!is.null(y_limit)) {
    out <- out + ggplot2::ylim(y_limit[1], y_limit[2])
  }
  out
}