Back to Elementolab/ChIPseeqer_Tutorial
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
-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:
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")