Hi,

I have been using DiffBind to perform differential binding analysis on ChIP-seq data, comparing two groups with 3 replicates for each. I am now facing an error when using its dba.count function.

I created a DBA object, Expt, using a sample sheet.

Expt <- dba(sampleSheet=”Expt.csv”)

After creating consensus peaksets by

Expt <- dba.peakset(Expt, consensus = -DBA_REPLICATE)

I used the dba.count function as follows.

Expt <- dba.count(Expt, peaks=Expt$masks$Consensus)

I then got an error message as follows.

Error in pv.counts(DBA, peaks = peaks, minOverlap = minOverlap,defaultScore = score,

Can’t count: some peaksets are not associated with a .bam file.

What is weird is that, when I used dba.count function on all the samples, without creating consensus peaksets, it works without showing any error message. In addition, this error message did not appear, when I used an older version of DiffBind. This error happens after I have updated R to ver. 3.3.1 and downloaded the latest version of DiffBind.

Thank you so much.

It looks like you are trying to count reads against a consensus peaksets made up of peaks that are in at least two replicates of at least two sample groups, is that correct? So if you have two sample groups (say WT and AT), and three replicates of each, the call to dba.peakset() with consensus = -DBA_REPLICATE will create two consensus peaksets, one for each sample group, each including peaks that overlap in at least two replicates. When you call dba.count(), you want to get peaks that overlap in both of these consensus sets, is that right?

The way to do this requires a few more steps. First, generate the sample-group consensus sets as you have done, but store it in a new object:

> ExptConsensus <-  dba.peakset(Expt, consensus = -DBA_REPLICATE)

Next limit this to only the consensus peaks:

> ExptConsensus <-  dba(ExptConsensus, mask=ExptConsensus$masks$Consensus)

Now retrieve the consensus of these:

> ConsensusPeaks <- dba.peakset(ExptConsensus, bRetrieve=TRUE)

And pass that in to dba.count():

> Expt <- dba.count(Expt, peaks=ConsensusPeaks)

The key is that when you use dba.peakset() with consensus = -DBA_REPLICATE, it adds new peaksets to the DBA object. These peaksets are not associated with an individual sample, but rather formed from a consensus of multiple samples. Therefore it doesn't make sense to use dba.count() to count reads, as reads are associated with individual samples. The bottom line is that each sample in the DBA object being counted needs to be associated with a set of reads (a bam file); a consensus peakset by itself is no longer associated with any one bam file.

Hope this helps-

Rory