R/plot_RLE.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_RLE <- function(x, col_factor, method = c("boxplot", "violin"), y_limit = NULL) {
  require(ggplot2)
  method <- match.arg(method)
  x_plot <- log(x + 1)
  x_med <- apply(x_plot, 1, median, na.rm = TRUE)
  x_plot <- t(apply(x_plot, 2, function(y) y - x_med))

  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("Relative Log Expression")+
    geom_hline(yintercept = 0, lty = 2)


  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
}