Tutorial: Visualization of ChIP-Seq peak overlaps using HOMER mergePeaks and UpSetR
11
gravatar for steve
2.8 years ago by
steve1.9k
United States
steve1.9k wrote:

EDIT: I've aggregated some of the scripts here, with example output here

I had previously made a post here on how to create venn diagrams in R using the venn output from HOMER's mergePeaks tool. However this is limited in scope to comparisons of 2 to 5 sets of peaks. For more flexibility and scalability, I created a script that uses the UpSetR package in R to create UpSet plots of the overlaps. The script itself can be found here with a sample basic workflow & output. I've also copied the script below:

#!/usr/bin/env Rscript

## USAGE: multi_peaks_UpSet_plot.R <sampleID> /path/to/venn.txt
## This script will process venn output from HOMER mergePeaks and make an UpSet plot
## this is currently only for R version 3.3.0; module unload r; module load r/3.3.0; R

# resources:
# http://www.caleydo.org/tools/upset/
# https://cran.r-project.org/web/packages/UpSetR/vignettes/basic.usage.html#example-5-alternative-input-formats 
# https://github.com/hms-dbmi/UpSetR

# ~~~~~ LOAD PACKAGES ~~~~~~~ #
library("UpSetR"); library("ggplot2"); library("grid"); library("plyr")
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # 


# ~~~~~ GET SCRIPT ARGS ~~~~~~~ #
args <- commandArgs(TRUE); cat("Script args are:\n"); args
SampleID<-args[1]
venn_table_file<-args[2]

plot_outdir<-dirname(venn_table_file)
plot_filepath<-paste0(plot_outdir,"/",SampleID,"_UpSetR_plot.pdf") 
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # 


# ~~~~~ PARSE THE VENN TABLE ~~~~~~~ #
# read in the venn text
venn_table_df<-read.table(venn_table_file,header = TRUE,sep = "\t",stringsAsFactors = FALSE,check.names = FALSE)
# > venn_table_df
# Sample2-H3K27AC.bed Sample1-H3K27AC.bed Total                                Name
# 1                                   X 29005                   Sample1-H3K27AC.bed
# 2                 X                   19523                   Sample2-H3K27AC.bed
# 3                 X                 X 35344 Sample2-H3K27AC.bed|Sample1-H3K27AC.bed


# get the venn categories
venn_categories<-colnames(venn_table_df)[!colnames(venn_table_df) %in% c("Total","Name")] 
cat("Venn categories are:\n"); venn_categories
# > cat("Venn categories are:\n"); venn_categories
# Venn categories are:
#   [1] "Sample2-H3K27AC.bed" "Sample1-H3K27AC.bed"


# venn_categories
num_categories<-length(venn_categories)
cat("Num categories are:\n"); num_categories
# > num_categories<-length(venn_categories)
# > cat("Num categories are:\n"); num_categories
# Num categories are:
#   [1] 2


# make a summary table
venn_summary<-venn_table_df[!colnames(venn_table_df) %in% venn_categories]
cat("Venn summary table is categories are:\n"); venn_summary
# > venn_summary<-venn_table_df[!colnames(venn_table_df) %in% venn_categories]
# > cat("Venn summary table is categories are:\n"); venn_summary
# Venn summary table is categories are:
#   Total                                Name
# 1 29005                   Sample1-H3K27AC.bed
# 2 19523                   Sample2-H3K27AC.bed
# 3 35344 Sample2-H3K27AC.bed|Sample1-H3K27AC.bed


# write summary table
plot_filepath<-paste0(plot_outdir,"/",SampleID,"_UpSetR_plot.pdf") 
write.table(venn_summary,file = paste0(plot_outdir,"/","venn_summary.tsv"),quote = FALSE,row.names = FALSE)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # 


# ~~~~~ SET UP THE PLOT ~~~~~~~ #
# convert the summary table to a numeric/int vector, with element names as the combination names
# # swap the | character with &; for passing to UpSet fromExpression
upset_expression<-setNames(venn_summary[['Total']],gsub("|","&",venn_summary[['Name']],fixed = TRUE))
# > upset_expression
# Sample1-H3K27AC.bed                   Sample2-H3K27AC.bed
# 29005                               19523
# Sample2-H3K27AC.bed&Sample1-H3K27AC.bed
# 35344

# save the plot into a PDF
cat("Plot output file is:\n"); plot_filepath

pdf(plot_filepath,width = 8,height = 8)
upset(fromExpression(upset_expression),order.by = "degree", decreasing = F,mainbar.y.label = "Overlapping Peaks", sets.x.label = "Peaks per Category") # , group.by = "sets"
dev.off()

Example output: screen shot 2017-10-30 at 4 31 15 pm

homer chip-seq peaks tutorial R • 2.6k views
ADD COMMENTlink modified 16 months ago • written 2.8 years ago by steve1.9k

Looks great steve. Can you post an example figure?

ADD REPLYlink written 16 months ago by Kevin Blighe39k
1

yes I've updated the original post and included example output here

ADD REPLYlink written 16 months ago by steve1.9k

Pretty creative plot - thanks!

ADD REPLYlink written 16 months ago by Kevin Blighe39k

Really cool way to visualise the overlapping peaks! Thanks Steve!

ADD REPLYlink written 16 months ago by Sej Modha4.1k
0
gravatar for morovatunc
2.5 years ago by
morovatunc400
Turkey
morovatunc400 wrote:

Hi steve, I would like to ask you something related the interpretation of the HOMER Venn content.

I think venn matrix should be diagonal but when I check the values it seems they are not. If you have any idea about this, could you enlighten me ?

I think homer is a great tool of identifying overlapping peaks but this is the only disadvantage that which makes me confused a lot.

Please check this link for a example. I have asked this a while ago but yet to be solved.

HOMER mergepeaks matrix interpreatation

ADD COMMENTlink modified 2.5 years ago • written 2.5 years ago by morovatunc400
1

That is a great question, but unfortunately I do not have an answer. I have never understood the purpose or usage of the matrix output, so I usually just ignore it / don't output it.

ADD REPLYlink written 2.5 years ago by steve1.9k

To be honest, I also kept it unused. The resulting file of the overlaps is very good. Especially where you can get the information about how many sub peaks are in a big merged peak.

Anyways, thank you for the anwser!

ADD REPLYlink written 2.5 years ago by morovatunc400
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2231 users visited in the last hour