Hi,

I'm trying to extract the consensus peakset that DiffBind constructs, with the raw counts in these regions. I have done this two different ways:

1) After running dba.count(), dba.contrast() and dba.analyze() with default settings:

counts_1 <- dba.report(DBA, th=1, bCounts = TRUE, bNormalized = FALSE, DataType = DBA_DATA_FRAME)

2) (I had to run dba.count() normally first, otherwise the second line of code did not run)

DBA <- dba.count(DBA)
DBA <- dba.count(DBA, peaks=NULL, score=DBA_SCORE_READS)
counts_2 <- dba.peakset(DBA, bRetrieve = TRUE, DataType = DBA_DATA_FRAME)

These two methods both yield an object with 104976 consensus peak regions. However, some of these contain 1s in every column (so there are no reads in any of my samples). My first question is, why is this? If these regions have no reads in any sample, why do they make it into the consensus peak set?

My second question is: The number of these regions without reads in different with the two methods. With method 1, I find 11476 such regions; with method 2, there are 9686, even though the total number of regions is the same. What causes this discrepancy, and what is the right way to extract raw read counts in consensus peaks?

Thank you in advance for your answer!

I'm not sure how you called peaks initially? The defaults for determining a consensus peakset are based strictly on peak overlaps, not on counts within peaks.

The default counting method may undercount overlapping regions if you don't supply a fragmentSize and/or have paired-end reads. One thing to try is to call dba.count() with bUseSummarizeOverlaps=TRUE, and compare the results.

The main difference in count scores between the default scoring method DBA_SCORE_TMM_MINUS_FULL and DBA_SCORE_READS, besides the normalization, is that reads in the control track are subtracted (that's the MINUS part). So you have more peaks with no counts because there are more control counts that ChIP counts overlapping that interval.

Also, you should be able to specify score=DBA_SCORE_READS directly to the initial call to dba.count() (but not peaks=NULL, as that tells the function to use the current counts to recompute the score).

If you want to do this outside of DiffBind, there's another way of doing it.

1. Call peaks(e.g. macs2). Convert the peak files to bed files.
2. Merge all the bed files from your samples with bedtools to make a consensus peakset.
3. Covert the merged peak file(*.bed file) into an *.saf file, to input into featureCounts(from SubRead package). Check this link how to make a saf file from bed file using simple awk command https://www.biostars.org/p/228636/.
4. Get the raw read count matrix for your consensus peakset by using that one saf file and all the bam files in featureCounts to get raw count matrix.

You could probably write a bash script to do all these steps in one go.