Elementolab/ChIPseeqerGetReadDensityProfiles

From Icbwiki

Back to Elementolab/ChIPseeqer_Tutorial

ChIPseeqerGetReadDensityProfiles

To run the tools directly from any folder, you need to add the $CHIPSEEQERDIR and $CHIPSEEQERDIR/SCRIPTS to your $PATH variable. Read How to set the CHIPSEEQERDIR variable.

1. If you haven't done it yet, run CS on ChIP-seq data to find peaks. Ignore this step if you have already located peaks.

ChIPseeqer.bin -chipdir BCL6CHIP -inputdir BCL6INPUT -t 15 -fold_t 2 -outfile bcl6peaks.txt

2. Run script to extract regions centered on peak summits (2kb windows used here, ie 1kb on each side)

extract_regions_around_peak_summits.pl --peakfile=bcl6peaks.txt --w=2000 

This will create a file called bcl6peaks.txt.centered2000

3. Extract read densities from other ChIP-seq dataset (eg histone modification; input is ignored here only)

ChIPseeqerGetReadDensityProfiles.bin -intervals bcl6peaks.txt.centered2000 -chipdir H3K4ME3H3K9ACCHIP \
  -format eland -fraglen 0 -outfile bcl6peaks.txt.centered2000.H3K4me3H3K9Ac_density -outepsmap bcl6peaks.txt.centered2000.H3K4me3H3K9Ac_density.eps 

Options are:

-format STR      [bed sam eland] 
-fraglen INT     [0 = no read extension otherwise extend to specified value] 
-rpkmnorm INT    [0/1, if 1, perform normalization by number of reads ]
-uniquereads INT [0/1 if 1 collapose clonal reads; 1 = recommended for TF and histone modification, not for nucleosome positioning ] 
-ws INT          [the window size, can be 10, 100 etc. Default is 10]
-outfile FILE    [the name of the output density file]
-outepsmap FILE  [the name of the output eps file with the 2D density plot]
-xlabel STR      [the label for the x axis of the 2D plot]
-ylabel STR      [the label for the y axis of the 2D plot]

4. See the results. This script will create a .eps file that contains a 2D plot for the densities (averaged per column/bin). The plot will look like this:

ALTERNATIVE:

Load read densities in R, average them, and plot them

m <- read.table("bcl6peaks.txt.centered2000.H3K4me3H3K9Ac", row.names=1)
par(cex=1.3) 
plot(10*(1:200)-1000, colMeans(m), type="l", lwd=5, col="red", ylab="Average H3K4me3+H3K9Ac read density", 
  xlab="Distance to maximum BCL6 peak height (bp)")
dev.copy2eps(file="H3K9Ac-H3K4me3-around-BCL6-peaks.eps")