R/coverage_plot.R

In riccardo-trozzo/BlissR: Package to analyse Bliss-seq data after BLISS-seq preprocessing pipeline

#' A coverage bar plot generator
#'
#' A function that allows you to plot the reads that falls in top n genes
#' @param samples a list of data frames of peaks annotations
#' @param n_top the number of top peaks that you want to plot, default is 100
#' @param plot_dir the directory in which you want the plots, default = "plots"
#' @keywords frequencies
#' @import tidyverse gridExtra ggsci
#' @export
#' @examples
#' coverage_plot(samples, n_top, plot_dir)
coverage_plot <- function(samples, n_top = 100, plot_dir = "plots"){
  count <- 0
  num <- n_top
  dir.create(file.path(plot_dir), showWarnings = FALSE)
  samples_filtered_ordered <- lapply(samples, function(x){
    x %>%
      filter(str_detect(GenomicAnnotation, "Promoter.*|.*Exon.*|Downstream.*")) %>%
      arrange(desc(PeakPvalue))
  })
  pl <- list()
  pl <- lapply(samples_filtered_ordered, function(x){
    count <<- count + 1
    tot_reads <- sum(x$PeakCoverage)
    plot <- x %>%
      select(c("chrom", "PeakCoverage")) %>%
      head(num) %>% mutate(percentage = PeakCoverage / tot_reads * 100)
    plot$order <- rownames(plot)
    plot$order <- factor(plot$order, levels = plot$order)
    ggplot(data = plot, aes( y=percentage, x=order)) +
      geom_bar(stat = "identity", col = pal_npg()(length(samples))[count], fill = pal_npg("nrc", alpha = 0.8)(length(samples))[count]) +
      ggtitle(names(samples_filtered_ordered)[count]) +
      theme(axis.text.x = element_blank()) +
      xlab(paste("top", n_top, "peaks"))
  })
  pdf(paste("coverage_plot.pdf", sep = ""))
  pdf(file.path(plot_dir, "coverage_plot.pdf"))
  do.call(grid.arrange, c(pl, nrow = round(length(samples_filtered_ordered)/2)))
  dev.off()
}