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 > bcl6peaks.txt.centered2000.H3K4me3H3K9Ac

Options are:

-format STR      [bed sam eland] 
-fraglen INT     [0 = no read extension otherwise extend to specified value] 
-nummtnorm INT   [0/1, if 1, perform normalization by number of mapped nucleotides ] 
-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]

• 4. 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")