#####
# Code from:
#	Schwartzman A, Jaffe AE, Gavrilov Y, Meyer CA. Multiple testing of local maxima
#	for detection of peaks in ChIP-seq data. Annals of Applied Statistics, Accepted.
# Author: Andrew Jaffe
# Date: 2/7/2013
######

Please find the functions here: http://biostat.jhsph.edu/~ajaffe/code/ChIP-Seq_functions.R
And an example analysis script here: http://biostat.jhsph.edu/~ajaffe/code/ChIP-Seq_analysis_FoxA1.R

The input function needs to be better standardized, but with your own ChIP-seq data, you can jump past the readData() function, assuming you do/have the following:

treat = your reads from the treatment channel
control = your reads from the control channel

You will probably have to "trick" our functions to use your data. the expected columns of each of these (treat, control) are: c("chr","start","end","x","y","strand"), like so (this is the head of 'control', the head of 'treat' looks the same):

   chr  start    end x y strand
1 chr1 233604 233639 0 2      -
2 chr1 559767 559802 0 3      +
3 chr1 742600 742635 0 2      +
4 chr1 742600 742635 0 0      +
5 chr1 744231 744266 0 0      +

where the first 3 columns and and last column come from the read. we do not use columns x or y, so you can just make columns of 0s, or NA. Then you would just make this list:

rawData = list(Treat = treat, Control = control)

then you can use the rest of our functions. I can also send the data used in the paper if you want to play around with our functions.

Let me know if you have any questions or issues.

Andrew Jaffe
andrewejaffe@gmail.com