R/plot_Pdist.R

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

#' Plot distribution of P-values
#' @param x Output of edgeR::topTags
#' @param multiplot Draw multiple plots?
#' @param plot_title Titles for the plots
#' @param F_abs Use absolute frequency?
#' @export

plot_Pdist <- function(x, multiplot = FALSE, plot_title,
                       F_abs = TRUE, bins = c("scott", "sturges", "FD"),
                       nbins = NULL, nrow, ncol, ...) {

  bins <- match.arg(bins)
  single_plot <- function(y) {
   # y <- data$table
   #  y <- y[!is.na(y$PValue), ]

    if(is.null(nbins)) {
    if(bins == "sturges") {
      classes <- nclass.Sturges(y$PValue)
    } else if(bins == "scott") {
      classes <- nclass.scott(y$PValue)
    } else if(bins == "FD") {
      classes <- nclass.FD(y$PValue)
    }
    } else {
      classes <- nbins
    }

    if(F_abs) {
      ggplot2::ggplot(data=y, ggplot2::aes(x = y$PValue)) +
        ggplot2::geom_histogram(position="identity",  bins = classes, ...) +
        ggplot2::ylab("Frequency") +
        ggplot2::xlab('P-value') +
        ggplot2::xlim(c(-0.03, 1.03))

    } else {
      ggplot2::ggplot(data=y, ggplot2::aes(x = y$PValue)) +
        ggplot2::geom_histogram(aes(y = (..count..)/sum(..count..)),
                                position="identity",  bins = classes, ...) +
        ggplot2::ylab("Relative frequency") +
        ggplot2::xlab('P-value') +
        ggplot2::xlim(c(-0.03, 1.03))
    }
  }

  if(multiplot) {
    merged_table <- list()
    for(i in seq_along(x)) {
      merged_table[[i]] <- x[[i]]$table$PValue
    }
    merged_table <- do.call("cbind", merged_table)
    colnames(merged_table) <- plot_title
    merged_table <- reshape2::melt(merged_table)
    colnames(merged_table)[2:3] <- c("contrast", "PValue")
    out <- single_plot(merged_table) + ggplot2::facet_wrap(~contrast,
                                                           nrow = nrow,
                                                           ncol = ncol,
                                                           scales = "free_y")

  } else {
    out <- single_plot(x$table)
    if(!is.null(title)) {
      out <- out + ggplot2::ggtitle(plot_title)
    }
  }
  out
}