In amitjavilaventura/seqViewR: Visualization of Omics And Sequencing Data In R

knitr::opts_chunk$set(echo = TRUE, fig.align = "center", fig.width = 4, fig.height = 4)
devtools::load_all("..")
library(dplyr)
library(purrr)

peak_list <- list.files("../testdata", "peak", full.names = T, recursive = T) %>%
  purrr::set_names(c("PeakX", "PeakY")) %>%
  purrr::map(~read.delim(.x))  
upsetPeaks(peak_list)

Load plotmics

# load plotmics
library(plotmics)

Run upsetPeaks()

upsetPeaks() draws an UpSet plot with the intersections between different sets of peaks using the function getVennCounts() (which calls ChIPpeakAnno::makeVennDiagram()) the and the package UpSetR.

Look at the UpSetR package documentation.

Look at the makeVennDiagram() function documentation

Required input

As input, upsetPeaks() takes a named list of data frames with the columns seqnames, start and end.

# read the peak annotation into a list
peak_list <- list.files("../testdata", "peak", full.names = T, recursive = T) %>%
  purrr::set_names(c("PeakX", "PeakY")) %>%
  purrr::map(~read.delim(.x))

peak_list[[1]][1:5, 1:7]

getVennCounts()

getVennCounts() calls ChIPpeakAnno::makeVennDiagram(), retrieves the Venn counts (number of overlaps between different sets of peaks) and builds a matrix of the peaks present in each set.

venn_counts <- getVennCounts(peak_list)

venn_counts$vennCounts
#      PeakX PeakY Counts
# [1,]     0     0      0
# [2,]     0     1     70
# [3,]     1     0    977
# [4,]     1     1     23
# attr(,"class")
# [1] "VennCounts"

venn_counts$matrix[1:5,]
# peak    PeakX   PeakY
# peak1   0       1
# peak2   0       1
# peak3   0       1
# peak4   0       1
# peak5   0       1

upsetPeaks() calls the function getVennCounts() and builds the UpSet plot of the peaks using the Venn counts.

Unexpected intersections

As mentioned, ChIPpeakAnno::makeVennDiagram() is called inside getVennCounts(). This function may have a unexpected outputs when considering the number of overlaps to build the intersection between different sets of regions. Considering the following example:

• We have two sets of regions named A and B: A has 4 regions, 3 of which overlap with regions in B; and B has 6 regions, 5 of which overlap with regions in A.
• In the case we consider the 5 regions from B that overlap with regions in A, the size of the resulting intersection (5) would be greater than the total size of A (4).
• To prevent that, the number of overlaps considered to build the intersection is the lowest possible. In our case, this is 2.
• Hence, with ggVennCounts() we would get 1 unique region from A, one unique region from B and 2 intersections. The sum of any of them is not the total size of A or B.

Default run

upsetPeaks(peak_list)

Customize plot

Change condition names

upsetPeaks(peak_list, conds = c("Condition 1", "Condition 2"))

upsetPeaks(peak_list, conds = c("Condition 1", "Condition 2"), conds_order = c("Condition 2", "Condition 1"))

Different order

upsetPeaks(peak_list, order.by = "freq") # default
upsetPeaks(peak_list, order.by = "degree")
upsetPeaks(peak_list, order.by = c("freq", "degree"))

Change labels

upsetPeaks(peak_list, mainbar.y.label = "This is an Y label for the main barplot", 
           sets.x.label = "This is the X label for the set size")

amitjavilaventura/seqViewR documentation built on Nov. 21, 2023, 10:12 a.m.