get the peaks that shared in multiple samples

Three gotchas when using R for Genomic data analysis

read in all the files in a folder to R

Ming Tang

February 7, 2016

read in all the files in a folder with a loop

see a post
I downloaded some broadpeak data from UCSC and put them in a folder /Users/mtang1/projects/ENCODE_ChIPseq/data/

I cut only first four columns of those files and How can I read in all the bed files as GRanges?

# the base of the path
cancer.base<- "/Users/mtang1/projects/ENCODE_ChIPseq/data"

# list.files will print out all the files in the folder ending with broadPeak
cancer.lines.peak.files<- list.files(cancer.base, pattern="*.broadPeak$", full.names=TRUE)
cancer.lines.peak.files

## [1] "/Users/mtang1/projects/ENCODE_ChIPseq/data/wgEncodeHaibTfbsEcc1Pol2V0416102Dm002p1hPkRep2.broadPeak"
## [2] "/Users/mtang1/projects/ENCODE_ChIPseq/data/wgEncodeHaibTfbsHct116Pol24h8V0416101PkRep2.broadPeak"   
## [3] "/Users/mtang1/projects/ENCODE_ChIPseq/data/wgEncodeHaibTfbsHelas3Pol2Pcr1xPkRep2.broadPeak"         
## [4] "/Users/mtang1/projects/ENCODE_ChIPseq/data/wgEncodeHaibTfbsHepg2Pol2Pcr2xPkRep2.broadPeak"          
## [5] "/Users/mtang1/projects/ENCODE_ChIPseq/data/wgEncodeHaibTfbsK562Pol2V0416101PkRep2.broadPeak"        
## [6] "/Users/mtang1/projects/ENCODE_ChIPseq/data/wgEncodeHaibTfbsPanc1Pol24h8V0416101PkRep1.broadPeak"    
## [7] "/Users/mtang1/projects/ENCODE_ChIPseq/data/wgEncodeHaibTfbsU87Pol24h8V0416101PkRep2.broadPeak"

# using a for loop to read in all the files and assign the cell-line name as the varaible name of the GRanges
# import function in rtracklayer will read in bed file as GRanges

library(rtracklayer)
for (file in cancer.lines.peak.files){
        # extract the cell line name
        cell.line.name <- gsub(".+wgEncodeHaibTfbs(.+)Pol2.+", "\\1", file)
        print (sprintf("reading in %s peak file", cell.line.name))
        file.handler<- file
        print (cell.line.name)
        # assign the cell-line name as the varaible name of the GRanges
        assign(cell.line.name, import(file.handler, format = "BED"))
        
}

## [1] "reading in Ecc1 peak file"
## [1] "Ecc1"
## [1] "reading in Hct116 peak file"
## [1] "Hct116"
## [1] "reading in Helas3 peak file"
## [1] "Helas3"
## [1] "reading in Hepg2 peak file"
## [1] "Hepg2"
## [1] "reading in K562 peak file"
## [1] "K562"
## [1] "reading in Panc1 peak file"
## [1] "Panc1"
## [1] "reading in U87 peak file"
## [1] "U87"

Ecc1

## GRanges object with 35977 ranges and 1 metadata column:
##           seqnames                 ranges strand   |        name
##              <Rle>              <IRanges>  <Rle>   | <character>
##       [1]     chr1       [713816, 714479]      *   |       peak1
##       [2]     chr1       [762384, 763181]      *   |       peak2
##       [3]     chr1       [859112, 859462]      *   |       peak3
##       [4]     chr1       [878092, 878373]      *   |       peak4
##       [5]     chr1       [878472, 878807]      *   |       peak5
##       ...      ...                    ...    ... ...         ...
##   [35973]     chrX [154299651, 154300068]      *   |   peak35973
##   [35974]     chrX [154444181, 154444945]      *   |   peak35974
##   [35975]     chrX [154493584, 154493918]      *   |   peak35975
##   [35976]     chrX [154841987, 154842828]      *   |   peak35976
##   [35977]     chrX [155110716, 155111326]      *   |   peak35977
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths

U87

## GRanges object with 42355 ranges and 1 metadata column:
##           seqnames                 ranges strand   |        name
##              <Rle>              <IRanges>  <Rle>   | <character>
##       [1]     chr1       [ 10081,  10449]      *   |       peak1
##       [2]     chr1       [713834, 714370]      *   |       peak2
##       [3]     chr1       [752419, 752810]      *   |       peak3
##       [4]     chr1       [762601, 763171]      *   |       peak4
##       [5]     chr1       [842211, 842561]      *   |       peak5
##       ...      ...                    ...    ... ...         ...
##   [42351]     chrX [154996980, 154997243]      *   |   peak42351
##   [42352]     chrX [154997284, 154997491]      *   |   peak42352
##   [42353]     chrX [155110793, 155111216]      *   |   peak42353
##   [42354]     chrX [155117874, 155118167]      *   |   peak42354
##   [42355]     chrX [155259636, 155260008]      *   |   peak42355
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths

Use lapply

It is awkward to use for loop in R, let’s just use lapply function. see a post

# full path to every file

## lapply will use import to read in every bed file and keep them in a list.
data<- lapply(cancer.lines.peak.files, import, format = "BED")
names(data)<- gsub(".+wgEncodeHaibTfbs(.+)Pol2.+", "\\1", cancer.lines.peak.files)

data$U87

## GRanges object with 42355 ranges and 1 metadata column:
##           seqnames                 ranges strand   |        name
##              <Rle>              <IRanges>  <Rle>   | <character>
##       [1]     chr1       [ 10081,  10449]      *   |       peak1
##       [2]     chr1       [713834, 714370]      *   |       peak2
##       [3]     chr1       [752419, 752810]      *   |       peak3
##       [4]     chr1       [762601, 763171]      *   |       peak4
##       [5]     chr1       [842211, 842561]      *   |       peak5
##       ...      ...                    ...    ... ...         ...
##   [42351]     chrX [154996980, 154997243]      *   |   peak42351
##   [42352]     chrX [154997284, 154997491]      *   |   peak42352
##   [42353]     chrX [155110793, 155111216]      *   |   peak42353
##   [42354]     chrX [155117874, 155118167]      *   |   peak42354
##   [42355]     chrX [155259636, 155260008]      *   |   peak42355
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths

Using Awk to join two files based on several columns

I was reading a thread on stackoverflow and found that this post was very interesting. I will go through the problem and the awk solution. Again, Awk is awesome!

Note:

It looks a bit messy when I copied the ipython notebook directly into Blogger. Ideally, one can write the blog using ipython notebook directly, but it does not work well with Blogger. see here http://blog.fperez.org/2012/09/blogging-with-ipython-notebook.html

I do not want to use other blogging platforms for now. So, just bear it and you can go to here to see the clean ipython notebook. Github now randers ipython notebook.

I created some dummy files.
file_a is a tab-delimited bed file with 6 colums:

In [1]:

cat file_a.bed

chr1 123 aa b c d
chr1 234 a b c d
chr1 345 aa b c d
chr1 456 a b c d

file_b is the file that contain additional infomation, which we want to add to file_a:

In [2]:

cat file_b.bed

xxxx abcd chr1 123 aa c d e
yyyy defg chr1 345 aa e f g

we want to annotate file_a based on the fact that columns 3,4,5 in file_b are the same as columns 1,2,3 in file_a.
To do this, we are going to use Awk associated array. see a link

Let me execute the awk one-liner first and then explain what's going on here:

In [7]:

awk 'NR==FNR{a[$3,$4,$5]=$1OFS$2;next}{$6=a[$1,$2,$3];print}' OFS='\t' \
file_b.bed file_a.bed

chr1 123 aa b c xxxx abcd
chr1 234 a b c 
chr1 345 aa b c yyyy defg
chr1 456 a b c 

we annotated file_a using file_b. Aka, we added first two columns from file_b to file_a.
There are several things happening here:
we see built-in variables in awk: NR and FNR. NR is the line number of the
current processing line.
when awk read in multiple files, awk NR variable will give the total number of records
relative to all the input file. Awk FNR will give you number of records for each inpu
file. see a link
here
for all the built-in variables in awk.
Let's deomonstrate the difference between NR and FNR:

In [5]:

awk '{print FILENAME, NR}' file_a.bed file_b.bed

file_a.bed 1
file_a.bed 2
file_a.bed 3
file_a.bed 4
file_b.bed 5
file_b.bed 6

FILENAME is another built-in variable for the input file name of awk.
There are 4 lines in file_a and 2 lines in file_b, and NR increments for the total lines.
compare with FNR:

In [6]:

awk '{print FILENAME, FNR}' file_a.bed file_b.bed

file_a.bed 1
file_a.bed 2
file_a.bed 3
file_a.bed 4
file_b.bed 1
file_b.bed 2

Now, awk prints out the line numbers in respect to each file.

From the awk code, we are reading file_b first. NR==FNR means when NR equals to FNR
(this is true only for file_b) do the following: {a[$3,$4,$5]=$1OFS$2;next}.
We created an associated array named a using columns 3,4,5 in file_b as keys and
the columns 1 and 2: $1"\t"$2 as values. we set OFS="\t" in the end of the command.
next means to proceed for the next line, rather than execute the following { } code block.

when awk reads in the second file (file_a.bed), NR==FNR is not true, awk
exectues the second { } code block: {$6=a[$1,$2,$3];print}
we look up the associated array a we created from file_b.bed using
the first three columns in file_a.bed as keys, and assign column 6 to
the looked-up values and print it out the whole line.

Conclusion: Awk is very powerful in text wrangling. Once get used to the
syntax, you can do fairly complicated formatting in an awk one-liner.
I strongly recommand you to learn it.

how to get a genome-wide motif bed file

Someone was asking this question on Seqanswers http://seqanswers.com/forums/showthread.php?t=40762&highlight=genome+motif+bed

For motif analysis, the most popular program is MEME http://meme.nbcr.net/meme/
There are a bunch of tools in the suites including some useful ones for ChIP-seq
I also saw RAST http://rsat01.biologie.ens.fr/rsa-tools/index.html and oPOSSUM were mentioned http://opossum.cisreg.ca/oPOSSUM3/

I analyze ChIP-seq data a lot. Usually, one gets a bed file containing the positions of the peaks. However, many motif analysis programs require fasta file as input.

There are many ways to get fasta file based on coordinates.
See:
http://www.biostars.org/p/7481/
and my previous post http://crazyhottommy.blogspot.com/2013/04/batch-converting-coordinates-to.html

Going back to the question, I've found several ways to get a bed file containing the coordinates for the motifs.

1. Homer homepage
http://homer.salk.edu/homer/
at the very bottom, there are several links for human and mouse .The file contains all the known motif coordinates for the whole genome. Files are big (several Gb). one can get the specific motifs occurrences by grep.

2.  from UCSC software http://genome.ucsc.edu/ENCODE/analysisTools.html
ENCODE-motifs  http://compbio.mit.edu/encode-motifs/
at the bottom of the page:

only the human data are available.

3. Motif-map http://motifmap.ics.uci.edu/
click motif search

Several more organisms are supported

search the motif you want ( I use CTCF as an example)

Click save on the bottom right.

Click my motifs on the upper right.

export to bed or other format files.