I have a chip-seq data set with three conditions (wt , KO1, KO2) and three time points (0h, 2h, 4h). For each condition and TP I have duplicated samples for IP and Input samples, which looks like that:

I would like not only to do a pairwise comparison of IP vs. Input, but also to be able to do a more complex comparisons. I would like for example to wee if there are any differential binding events based on the differences over time (something like a time-series analysis) within each condition (e.g. comparing wt.0h - wt.2h - wt.4h). 

But to make it even more complicated, I would also like to make a comparison of the different time-series against each other. What I mean by that is for example to compare differential binding events from the wt time-series to KO1 or to KO2. In this case I would like to see if there are peaks (= overlapping genes) in the wt time-series which are significantly different than the peaks (or regions). in the KO sample time-series.

I know this is not so easy, as peaks in one conditions/time-point don't necessarily mean also peaks in a different samples. But these kind of behavior is exactly what i would like to catch. If there is a significant peak in the wt sample, but no peak in the KO I would like to see it. This would a normal peak calling, but when I have multiple crossing samples, I'm not sure how to do the analysis.

I would appreciate any advice.

thanks, Assa

Much of this seems to be doable with judicious use of makeContrasts for the contrast= argument in glmQLFTest (assuming you're using edgeR). If you want to determine whether time has any effect:

con <- makeContrasts(wt.2h - wt.0h, wt.4h - wt.0h, levels=design)

... which asks whether there's any differences between the time points. This is the closest you will get to a time series analysis within each genotype. There's not much point fitting a spline or doing anything more complicated, given that you only have three time points to work with.

If you want to find differential time responses between genotypes:

con <- makeContrasts((wt.2h - wt.0h) - (ko.2h - ko.0h),
                     (wt.4h - wt.0h) - (ko.4h - ko.0h), levels=design)

... which asks whether the WT time effect is that same as the KO time effect. Here, we are looking for genes where the log-fold change between time points is different between genotypes. This test doesn't care about differences in the baseline level of binding (i.e., at time 0) between genotypes, it will only reject if the effect between time points is different between genotypes.

Alternatively, you could do something like:

con <- makeContrasts(wt.0h - ko.0h, wt.2h - ko.2h, wt.4h - ko.4h, levels=design)

... which asks whether there are any differences between WT and KO binding at any time point. Unlike the previous contrast, this will reject the null if there is a difference in the baseline and/or in the magnitude of the time effects between genotypes.

As for your other question: if you're using csaw, you'll be using sliding windows to define genomic regions. Here, there is no concept of a window being "significant" in one sample and not in the others. The most that we do is to filter windows based on their abundance, and this is done using the average count across all samples. This strategy is deliberately blind to the exact differences between conditions, which ensures that we never make choices about which windows to keep or discard based on the presence/absence of binding in different conditions. It means that we don't have to worry about the differences between conditions when defining the genomic regions of interest; this decision is left (correctly) to the hypothesis testing machinery.

Finally, I should add that, depending on the contrast, you may need to normalize differently. If you are comparing ChIP samples to input, you should be normalizing for composition biases. If you are comparing the ChIP samples to each other, you will have to choose between normalization for composition and efficiency biases, depending on the biological system. Check out the user's guide for more details.