Upset plot for overlapping ChIP-seq peaks: an alternative to Venn diagram

ChIP-seq-upset

Ming Tang

January 25, 2016

I am going to demonstrate how to make Upset plot using UpSetR package for ChIP-seq peaks. After running a peak-calling method for 3 samples, you want to compare how many peaks are overlapped in two samples and in three samples. People may draw venn diagram for this kind of problem, but when you have more than 3 samples, the venn diagram is hard to read.

# install it if you have not 
# install.packages("UpSetR")
library(UpSetR)
dummy<- data.frame(peak=c("peak_1", "peak_1", "peak_1","peak_2","peak_2","peak_3","peak_3","peak_4"),
                   sample=c("sample_1","sample_2","sample_3","sample_2","sample_3","sample_2","sample_3","sample_1"))

dummy

##     peak   sample
## 1 peak_1 sample_1
## 2 peak_1 sample_2
## 3 peak_1 sample_3
## 4 peak_2 sample_2
## 5 peak_2 sample_3
## 6 peak_3 sample_2
## 7 peak_3 sample_3
## 8 peak_4 sample_1

To get the format that UpsetR accepts, we have to “spread” the sample column to multiple columns. It is like to change a long format data to wide format. See here as well.

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)

peaks<- dummy %>% mutate(value =1) %>% spread(sample, value, fill =0 )
peaks

##     peak sample_1 sample_2 sample_3
## 1 peak_1        1        1        1
## 2 peak_2        0        1        1
## 3 peak_3        0        1        1
## 4 peak_4        1        0        0

I noticed that one has to convert peaks(local data frame) back to data.frame for UpsetR to work. In this example, I did not convert the peaks to local data frame by tbl_df(peaks).

upset(peaks, order.by = "freq")

find the nearest upstream genes using GRanges

I am using nearest in the GRanges package to annotate my genomic intervals (breakpoints from whole genome sequencing). Those are not ChIP-seq peaks and the strandness for each interval is critical to me. I want to find the closest gene that is upstream of the breakpoint.

Make sure you read my note for the “+” and “-” strandness of bedpe for structural variants. In a word, + means the region is at the 5’ of the breakpoint and - means the region is at the 3’ of the breakpoint. In other words, I need to annotate the breakpoints with the closest genes that are UPSTREAM of the breakpoints. I want the nearest gene that is upstream of the breakpoint no matter what the strand of the gene is. However, when consider which gene is more closer to the breakpoint, the strandness of the gene needs to be considered.

Note that, I am not saying the strandness of the genes (as mentioned here http://ygc.name/2014/01/14/bug-of-r-package-chippeakanno/), but rather the strandness of the intervals.

It looks like distanceToNearest() always returns an obsolute value for distance no matter it is upstream or downstream of the subject hits. We can get around with it.

The bumphunter bioconductor package by Rafael A. Irizarry et.al has a function called annoateNearest can better annotate the nearest distances considering the strandness, but it DID NOT restrict nearest to find only upstream or downstream relative to the query.

Do not make me wrong, bedtools is a great tool and I use it everyday. I was trying to avoid using bedtools and put everything in R for the sake of reproducibility, although bedtools closest documentation has a much better explannation of all the cases I want to use.

In a word, for GRanges, is it possible to find the upstream nearest or downstream nearest as bedtools closest -D -id andbedtools closest -D -iu?

Let me make some toy examples and demonstrate the usage of follow and precede. Read the post and answers there as well.

             -----> start 

————–|————————|——————— geneA in plus strand

——-|———————–|—————————– geneB in minus strand

                 <-----start

                                         |-----|         breakpoint in *plus* strand

suppressMessages(library(GenomicRanges)) 
breakpoint<- GRanges(seqnames = "chr1", ranges=IRanges(start=8, width=2), strand = "+")

genes<- GRanges(seqnames = "chr1", ranges=IRanges(start=c(3,1,12,10), end=c(6,4,15,13)), strand = c("+", "-", "+","+"))
mcols(genes)<- c("geneA", "geneB", "geneC", "geneD")

breakpoint

## GRanges object with 1 range and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1    [8, 9]      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

genes

## GRanges object with 4 ranges and 1 metadata column:
##       seqnames    ranges strand |       value
##          <Rle> <IRanges>  <Rle> | <character>
##   [1]     chr1  [ 3,  6]      + |       geneA
##   [2]     chr1  [ 1,  4]      - |       geneB
##   [3]     chr1  [12, 15]      + |       geneC
##   [4]     chr1  [10, 13]      + |       geneD
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

I want to return geneB which is upstream of the breakpoint. if one uses follow:

follow(breakpoint, genes)

## [1] 1

genes[follow(breakpoint, genes)]

## GRanges object with 1 range and 1 metadata column:
##       seqnames    ranges strand |       value
##          <Rle> <IRanges>  <Rle> | <character>
##   [1]     chr1    [3, 6]      + |       geneA
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

It returns geneA. However, I want to return geneB. geneB is on the minus strand and the transcription start site is closer to the breakpoint.

From the help page:

follow: The opposite of precede, follow returns the index of the range in subject that is directly followed by the range in x. Overlapping ranges are excluded. NA is returned when there are no qualifying ranges in subject.

Orientation and Strand: The relevant orientation for precede and follow is 5’ to 3’, consistent with the direction of translation. Because positional numbering along a chromosome is from left to right and transcription takes place from 5’ to 3’, precede and follow can appear to have ‘opposite’ behavior on the + and - strand. Using positions 5 and 6 as an example, 5 precedes 6 on the + strand but follows 6 on the - strand.

In my case, follow means that the gene is directly followed by the breakpoint considering the strandness of the gene. Because geneA is on plus strand, and the breakpoint is on plus strand, follow(breakpoint, genes) will return geneA, not geneB.

To get geneB, I have to get around by first resize the gene to size 1 (the transcrtiption start site) and then unstrand it.

on the help page:

A range with strand * can be compared to ranges on either the + or - strand. Below we outline the priority when ranges on multiple strands are compared. When ignore.strand=TRUE all ranges are treated as if on the + strand.

x on + strand can match to ranges on both + and * strands. In the case of a tie the first range by order is chosen.

x on - strand can match to ranges on both - and * strands. In the case of a tie the first range by order is chosen.

x on * strand can match to ranges on any of +, - or * strands. In the case of a tie the first range by order is chosen.

follow(breakpoint, unstrand(resize(genes,width=1)))

## [1] 2

genes[follow(breakpoint, unstrand(resize(genes,width=1)))]

## GRanges object with 1 range and 1 metadata column:
##       seqnames    ranges strand |       value
##          <Rle> <IRanges>  <Rle> | <character>
##   [1]     chr1    [1, 4]      - |       geneB
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Now, let’s change the breakpoint to minus strand:

           |-----|                              breakpoint in *minus* strand

———————————-|———|—– geneC in plus strand
——————————-|———–|—— geneD in plus strand

In this case, I want to return gene D. If one uses follow, it will return NA:

breakpoint1<- GRanges(seqnames = "chr1", ranges=IRanges(start=8, width=2), strand = "-")
follow(breakpoint1, genes)

## [1] NA

Because, follow requires that the genes are on the same strand with the breakpoint and the breakpoint follows genes or the genes are followed by the breakpoint. The breakpoint is on the minus strand while geneC and geneD are both on plus strand.

Instead, try

follow(breakpoint1, unstrand(resize(genes,width=1)))

## [1] 4

genes[follow(breakpoint1, unstrand(resize(genes,width=1)))]

## GRanges object with 1 range and 1 metadata column:
##       seqnames    ranges strand |       value
##          <Rle> <IRanges>  <Rle> | <character>
##   [1]     chr1  [10, 13]      + |       geneD
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Let’s change geneD to minus strand:

           |-----|                              breakpoint in *minus* strand

———————————-|———|—– geneC in plus strand
——————————-|———–|—— geneD in minus strand

In this case, I want to return gene C, because geneC is on the plus strand, and its transcription start site is closer to the breakpoint.

genes1<- genes<- GRanges(seqnames = "chr1", ranges=IRanges(start=c(3,1,12,10), end=c(6,4,15,13)), strand = c("+", "-", "+","-"))

mcols(genes1)<- c("geneA", "geneB", "geneC", "geneD")
genes1[follow(breakpoint1, genes1)]

## GRanges object with 1 range and 1 metadata column:
##       seqnames    ranges strand |       value
##          <Rle> <IRanges>  <Rle> | <character>
##   [1]     chr1  [10, 13]      - |       geneD
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

genes1[follow(breakpoint1, unstrand(genes1))]

## GRanges object with 1 range and 1 metadata column:
##       seqnames    ranges strand |       value
##          <Rle> <IRanges>  <Rle> | <character>
##   [1]     chr1  [10, 13]      - |       geneD
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

It returns geneD again, but I want to return gene C.

genes1[follow(breakpoint1, unstrand(resize(genes1, width=1)))]

## GRanges object with 1 range and 1 metadata column:
##       seqnames    ranges strand |       value
##          <Rle> <IRanges>  <Rle> | <character>
##   [1]     chr1  [12, 15]      + |       geneC
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

upstreamGenes<- genes1[follow(breakpoint1, unstrand(resize(genes1, width=1)))]

## distance requires that query and subject on the same strand. 
## Note distance calculation changed in BioC 2.12 to accommodate zero-width ranges in a consistent and intuitive manner

distance(breakpoint1, unstrand(upstreamGenes))

## [1] 2

Conclusions

It seems that using follow(breakpoing, unstrand(resize(genes, width=1))) serves my purpose very well. However, be careful that some breakpoint may not have a upstream gene and it will return NA using follow. GRanges object can not be subscripted by NAs.

Using wget to download specific files from ftp but avoiding the directory structure

I want to download some files from a ftp site, and I only want to download some files with names matching a pattern. How can I do it?
Use wget ! It is a very versatile command and I just got to know several tricks.

When there are many levels of folder, you want to search down to all the folders:

-r --recursive Turn on recursive retrieving.

   -l depth
   --level=depth
       Specify recursion maximum depth level depth.  The default maximum depth is 5.

You can specify what files you want to download or reject using wild cards:

Recursive Accept/Reject Options
-A acclist --accept acclist
-R rejlist --reject rejlist
Specify comma-separated lists of file name suffixes or patterns to accept or reject. Note that if any of the wildcard characters, *, ?, [ or ], appear in an element of acclist or rejlist, it will be treated as a pattern, rather than a suffix.

If you want to save the file to a different name:

-O file --output-document=file
The documents will not be written to the appropriate files, but all will be concatenated together and written to file. If - is used as file, documents will be printed to standard output, disabling link conversion. (Use ./- to print to a file literally named -.)

       Use of -O is not intended to mean simply "use the name file instead of the
       one in the URL;" rather, it is analogous to shell redirection: wget -O file
       http://foo is intended to work like wget -O - http://foo > file; file will be
       truncated immediately, and all downloaded content will be written there.

       For this reason, -N (for timestamp-checking) is not supported in combination
       with -O: since file is always newly created, it will always have a very new
       timestamp. A warning will be issued if this combination is used.

       Similarly, using -r or -p with -O may not work as you expect: Wget won’t just
       download the first file to file and then download the rest to their normal
       names: all downloaded content will be placed in file. This was disabled in
       version 1.11, but has been reinstated (with a warning) in 1.11.2, as there
       are some cases where this behavior can actually have some use.

       Note that a combination with -k is only permitted when downloading a single
       document, as in that case it will just convert all relative URIs to external
       ones; -k makes no sense for multiple URIs when they’re all being downloaded
       to a single file.

If you do not need the folder structure:

-nd --no-directories
Do not create a hierarchy of directories when retrieving recursively. With this option turned on, all files will get saved to the current directory, without clobbering (if a name shows up more than once, the filenames will get extensions .n).

one alternative way is to specify -nH and --cut-dirs=10 together

-nH --no-host-directories
Disable generation of host-prefixed directories. By default, invoking Wget with -r http://fly.srk.fer.hr/ will create a structure of directories beginning with fly.srk.fer.hr/. This option disables such behavior.

--cut-dirs=number
Ignore number directory components. This is useful for getting a fine-grained control over the directory where recursive retrieval will be saved.

       Take, for example, the directory at
       ftp://ftp.xemacs.org/pub/xemacs/.  If you retrieve it with -r, it
       will be saved locally under ftp.xemacs.org/pub/xemacs/.  While the
       -nH option can remove the ftp.xemacs.org/ part, you are still stuck
       with pub/xemacs.  This is where --cut-dirs comes in handy; it makes
       Wget not "see" number remote directory components.  Here are
       several examples of how --cut-dirs option works.

               No options        -> ftp.xemacs.org/pub/xemacs/
               -nH               -> pub/xemacs/
               -nH --cut-dirs=1  -> xemacs/
               -nH --cut-dirs=2  -> .

               --cut-dirs=1      -> ftp.xemacs.org/xemacs/
               ...

       If you just want to get rid of the directory structure, this option
       is similar to a combination of -nd and -P.  However, unlike -nd,
       --cut-dirs does not lose with subdirectories---for instance, with
       -nH --cut-dirs=1, a beta/ subdirectory will be placed to
       xemacs/beta, as one would expect.

If you want to save files to a different folder name:

-P prefix --directory-prefix=prefix
Set directory prefix to prefix. The directory prefix is the directory where all other files and subdirectories will be saved to, i.e. the top of the retrieval tree. The default is . (the current directory).

Continue to download a file:

-c --continue
Continue getting a partially-downloaded file. This is useful when you want to finish up a download started by a previous instance of Wget, or by another program

There are so many different options, just man wget to see all of them! I am impressed on how versatile this command is!

How to make a box-plot with jittered points for multiple groups

Make a box-plot with jittered points for multiple groups using ggplot2

Let me demonstrate it with an example. Use the ToothGrowth data set in the ggplot2 library

library(ggplot2)

ggplot(ToothGrowth, aes(x=as.factor(dose), y=len, color=supp)) + 
        geom_boxplot(position=position_dodge(0.9))+
        geom_jitter(position=position_dodge(0.9)) +
        xlab("dose")

I want to make the points separate from each other rather than on the same vertical line.

ggplot(ToothGrowth, aes(x= as.factor(dose), y=len, color= supp,fill= supp)) + 
        geom_point(position=position_jitterdodge(dodge.width=0.9)) +
        geom_boxplot(fill="white", alpha=0.1, outlier.colour = NA, 
                     position = position_dodge(width=0.9)) +
        xlab("dose")

convert html to pdf by pandoc

Pandoc is a very useful tool to convert common formats.

First install pandoc on mac by:

brew install pandoc

pandoc requires pdflatex to convert to pdfs.

install mactex:
download it and just double click it should install it.

next, put executables including latexpdf to the path.

echo export PATH=$PATH:/usr/texbin/ >> .bashrc
source ~/.bashrc

You can specify margins of the pdf by -V geometry:margin=1in:

pandoc http://quinlanlab.org/tutorials/cshl2013/gemini.html -V geometry:margin=1in -o gemini.pdf`  

when codes in html are too long, they get cut-off

Very thankful, I found the answer in this post:

Save the following as listings-setup.tex

% Contents of listings-setup.tex
\usepackage{xcolor}

\lstset{
    basicstyle=\ttfamily,
    numbers=left,
    keywordstyle=\color[rgb]{0.13,0.29,0.53}\bfseries,
    stringstyle=\color[rgb]{0.31,0.60,0.02},
    commentstyle=\color[rgb]{0.56,0.35,0.01}\itshape,
    numberstyle=\footnotesize,
    stepnumber=1,
    numbersep=5pt,
    backgroundcolor=\color[RGB]{248,248,248},
    showspaces=false,
    showstringspaces=false,
    showtabs=false,
    tabsize=2,
    captionpos=b,
    breaklines=true,
    breakatwhitespace=true,
    breakautoindent=true,
    escapeinside={\%*}{*)},
    linewidth=\textwidth,
    basewidth=0.5em,
}

Then, invoke pandoc:

pandoc https://cran.r-project.org/web/packages/gapmap/vignettes/tcga_example.html --listings -H listings-setup.tex -o gapmap_TCGA.pdf