I have ChIP-seq with two conditions. In condition A, there should a lot of peaks (normal binding). In condition B, there should be none (no binding). Calling peaks with MACS confirms this. When you perform differential binding analysis, you would expect most of the significant peaks to be up in condition A. However, I get half the peaks up and half down. That does not make sense.

I am using DiffBind. My first instinct was that the peaks were being normalized, so in condition B (where there are few reads in the peak regions), those few reads get normalized to a much higher level. However, that should not be the case according to the manual:

When dba.analyze is invoked using the default method=DBA_EDGER, a standardized differential analysis is performed using the edgeR package. ... First, a matrix of counts is constructed for the contrast, with columns for all the samples in the first group, followed by columns for all the samples in the second group. The raw read count is used for this matrix; if the bSubControl parameter is set to TRUE (as it is by default), the raw number of reads in the control sample (if available) will be subtracted (with a minimum final read count of 1). Next the library size is computed for each sample for use in subsequent normalization. By default, this is the total number of reads in the library (calculated from the source BAM//BED file). Alternatively, if the bFullLibrarySize parameter is set to FALSE,the total number of reads in peaks (the sum of each column) is used.

That disproves my initial theory. Is there another explanation?

Hello Igor-

Actually the problem is, indeed, normalization. When using the edgeR package, the default TMM normalization is used. Like most normalization methods, this relies on an assumption that there is a core set of peaks that are similarly enriched in both groups.

You can see the effects of the normalization after performing an analysis by comparing the MA plots with and without normalization:

> par(mfrow=c(2,1)
> dba.plotMA(myDBA, contrast=1)
> dba.plotMA(myDBA,contrast=1, bNormalized=FALSE)

We have noticed that in the most recent releases, the issue you are encountering has gotten more pronounced in experiments where one group has little or no binding. I'm looking into this (I suspect it relates to something that changed in edgeR). We have found that in these cases the normalization used by DESeq2 is less severe and may be more appropriate in this case. You can run both analyses and compare them:

> myDBA <- dba.analyze(myDBA, method=c(DBA_EDGER, DBA_DESEQ2)
> par(mfrow=c(2,2)
> dba.plotMA(myDBA, contrast=1, method=DBA_EDGER, bNormalized=FALSE)
> dba.plotMA(myDBA, contrast=1, method=DBA_EDGER)
> dba.plotMA(myDBA, contrast=1, method=DBA_DESEQ2,bNormalized=FALSE)
> dba.plotMA(myDBA, contrast=1, method=DBA_DESEQ2)

This issue is one of great interest to us and we are looking at alternatives, such as performing only the most basic normalization (for library depth) in these cases.

Cheers-
Rory

The bNormalized parameter to dba.plotMA() and dba.report() has no impact in how the analysis is carried out, it just uses non-normalized counts for the plot and/or report to show some underlying data. There is currently no way with DiffBind to perform an analysis without normalizing (using a normalization method intrinsic to the underlying differential analysis package)..

The TMM normalization used in edgeR is more complicated than just taking into account sequencing depth. For example, if one sample group has a set of enriched peaks A, and the other sample group has the same enriched peaks A plus gains an additional set of enriched peaks B, and they are all sequenced to the same depth, the same number of reads have to be divided between more peaks in the second group. Without normalizing, it will appear that the A peaks will be weaker in the second group when the truth is that the binding affinity hasn't changed. So we need to normalize for more than just library depth. (Another example is differing ChIP efficiencies.) 

I encourage you to read the well-written TMM paper by Robinson and Oshlack.

Have you tried the DESeq2 method? I'm curious to see how you results differ between the two methods.

-R