R/bliss_venn.R

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

#' A Venn diagram generator function
#'
#' This function founds in a list of data frames of peaks regions the one that has more peaks
#' @param peaks a list of data frames of peaks regions
#' @keywords maximum
#' @export
#' @examples
#' find_maximum(list)
find_maximum <- function(peaks){
  max <- 0
  for (condition in peaks) {
    condition <- as.data.frame(condition)
    if (nrow(condition) >= max) {
      max <- nrow(condition)
    }
  }
  return(max)
}

#' A sub function needed to plot venn diagrams of peaks
#'
#' This function allows you to plot venn diagrams of genes
#' @param peaks_list a list of data frames of peaks
#' @keywords peaks
#' @export
#' @import ChIPpeakAnno tidyverse
#' @examples
#' plot_venn(peaks_list, name)
plot_venn <- function(peaks_list, name, plot_dir){
  library(ChIPpeakAnno)
  library(tidyverse)
  pdf(file.path(plot_dir, paste("venn_diagram_", name, ".pdf", sep = "")))
  makeVennDiagram(peaks_list, NameOfPeaks=names(peaks_list),
                  totalTest=find_maximum(peaks_list),scaled=FALSE, euler.d=FALSE,
                  by = "region",
                  category.names = names(peaks_list),
                  height = 3000 ,
                  width = 3000 ,
                  resolution = 600,
                  compression = "lzw",
                  lwd = 1,
                  col = pal_npg()(length(peaks_list)),
                  fill = pal_npg("nrc", alpha = 0.3)(length(peaks_list)),
                  cex = 0.6,
                  fontfamily = "sans",
                  cat.cex = 0.6,
                  cat.default.pos = "outer",
                  cat.fontfamily = "sans",
                  cat.col = pal_npg()(length(peaks_list)),
                  rotation = 1,
  )
  dev.off()
}

#' A sub function needed to plot venn diagrams of genes
#'
#' This function allows you to plot venn diagrams of genes
#' @param genes_lists a list of data frames of gene_names
#' @param_name name of the plot
#' @param plot_dir the directory in which you want the plots, default = "plots"
#' @keywords genes
#' @export
#' @import VennDiagram tidyverse
#' @examples
#' genes_venn(genes_lists, name)
genes_venn <- function(genes_lists, name, plot_dir){
  venn_plot <- venn.diagram(
    x = genes_lists,
    category.names = names(genes_lists),
    filename = file.path(plot_dir, paste('venn_diagram_genes_', name, '.tiff', sep = "")),
    height = 3000 ,
    width = 3000 ,
    resolution = 600,
    compression = "lzw",
    lwd = 1,
    col = pal_npg()(length(genes_lists)),
    fill = pal_npg("nrc", alpha = 0.3)(length(genes_lists)),
    cex = 0.6,
    fontfamily = "sans",
    cat.cex = 0.6,
    cat.default.pos = "outer",
    cat.fontfamily = "sans",
    cat.col = pal_npg()(length(genes_lists)),
    rotation = 1
  )
}

#' A Venn diagram generator function
#'
#' This function allows you to plot venn diagrams choosing between plotting intersection of peaks of of genes
#' @param plot_dir the directory in which you want the plots, default = "plots"
#' @param peaks a list of data frames of genomic coordinates of peaks
#' @param by choose between "peaks" or "genes", peaks if you want to plot the venn of the overlapping peaks and genes if you want to plot the venn of the gene hits
#' @param ... select one ore more parameters from the annotation  table, the venn plot will be restricted to the corresponding regions
#' @keywords venn
#' @export
#' @import GenomicRanges ChIPpeakAnno tidyverse data.table
#' @examples
#' venn_generator(samples, by = "peaks", selection)
venn_generator <- function(plot_dir = "plots", samples, by, ...){
  dir.create(file.path(plot_dir), showWarnings = FALSE)
  args <- rlang::enexprs(...)
  if (by == "peaks") {
    #convert every sample in the list to GRanges
    samples <- lapply(samples, function(condition){
      condition %>%
        distinct() %>%
        GRanges()
    })

    if (args[[1]] == "all") {
      plot_venn(samples, "all", plot_dir)
    }

    else if (args[[1]] == "intergenic") {
      samples <- lapply(samples, function(condition){
        condition %>%
          as.data.frame() %>%
          filter(str_detect(GenomicAnnotation, "Intron.*|.*intergenic.*")) %>%
          GRanges()
      })
      plot_venn(samples, "intergenic", plot_dir)

    }
    else {
      samples <- lapply(samples, function(condition){
        condition %>%
          as.data.frame() %>%
          filter(!!! args) %>%
          GRanges()
      })
      plot_venn(samples, "selection", plot_dir)
    }
  }

  else if (by == "genes") {
    if (args[[1]] == "all") {
      gene_sets <- lapply(samples, function(condition){
        t <- condition %>%
          select(GeneName) %>%
          distinct()
        t <- t[, GeneName] #with this notation "[ ,GeneName]" the dataframe becomes a matrix, that is the format accepted by the venn package
      })
      genes_venn(gene_sets, "all", plot_dir) #plot the venn!
    }
    else {
      gene_sets <- lapply(samples, function(condition){
        t <- condition %>%
          filter(!!! args) %>%
          select(GeneName) %>%
          distinct()
        t <- t[, GeneName] #with this notation "[ ,GeneName]" the dataframe becomes a matrix, that is the format accepted by the venn package
      })
      genes_venn(gene_sets, "selection", plot_dir) #plot the venn!
    }
  }

  else {
    print("option 'by' is mandatory and it can be only 'genes' or 'peaks'")
  }
}