Elementolab/ChIPseeqerGetReadDensityProfiles

From Icbwiki

Revision as of 16:27, 8 November 2010 Eug2002 (Talk | contribs) ← Previous diff		Revision as of 21:28, 21 December 2010 Eug2002 (Talk | contribs) Next diff →
Line 19:		Line 19:
	ChIPseeqerGetReadDensityProfiles.bin -intervals bcl6peaks.txt.centered2000 -chipdir H3K4ME3H3K9ACCHIP \		ChIPseeqerGetReadDensityProfiles.bin -intervals bcl6peaks.txt.centered2000 -chipdir H3K4ME3H3K9ACCHIP \
-	-format eland -fraglen 0 > bcl6peaks.txt.centered2000.H3K4me3H3K9Ac	+	-format eland -fraglen 0 -outfile bcl6peaks.txt.centered2000.H3K4me3H3K9Ac_density -outepsmap bcl6peaks.txt.centered2000.H3K4me3H3K9Ac_density.eps
	Options are:		Options are:
Line 28:		Line 28:
	-uniquereads INT [0/1 if 1 collapose clonal reads; 1 = recommended for TF and histone modification, not for nucleosome positioning ]		-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]		-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. Load read densities in R, average them, and plot them'''	+	'''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:
		+
		+	[[Image:CS_density_plot.png|500px|upright|Average Conservation Profile example]]
		+
		+
		+	ALTERNATIVE:
		+
		+	'''Load read densities in R, average them, and plot them'''
	m <- read.table("bcl6peaks.txt.centered2000.H3K4me3H3K9Ac", row.names=1)		m <- read.table("bcl6peaks.txt.centered2000.H3K4me3H3K9Ac", row.names=1)

---

Revision as of 21:28, 21 December 2010

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] 
-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]
-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")