DiffBind analysis - minCount parameter in dba.count()

Question

[vk]

Why does changing minCount from 0 to 1 in dba.count() (with summits set to false and control reads supplied) change the called status of a peak? With minCount set to 1, a peak is being called in all samples, but with minCount set to 0, the same peak is being called in a few of them?

So, if minCount is set to 1, it will show the peak is called anyway? But that is not the case, cause there are peaks that are not called in half of the samples but called in the other half with minCount set to 1?

Am I missing something basic? What might be happening in my case?

[Rory Stark]

If you can supply a bit more info, I can see if I can reproduce this.

How are you determining how many samples a peak is called in? Can you show a bit of code demonstrating:

• your first call to dba.count()
• show a peak and which samples it is called in
• your second call to dba.count() changing minCount

• show how the peak is appearing to have been called in more/fewer samples

  Also, can you confirm which version of DiffBind you are using?

[vk]

Hi Rory! Thank you for your response. I was referring to the columns "Called1" and "Called2" in the report, they change when I change the minCount.

Diffbind version 3.0.13

Line from the Code:

minCount = 1

exp <- dba.count(exp, bSubControl = TRUE, minCount = 1, summits = FALSE)

minCount = 0

exp <- dba.count(exp, bSubControl = TRUE, minCount = 0, summits = FALSE)

I made no other changes for the two, the following were set:

exp$config$design <- FALSE

exp$config$doBlacklist <- FALSE

exp$config$doGreylist <- FALSE

After count:

exp <- dba.normalize(exp,normalize = DBA_NORM_LIB, library =DBA_LIBSIZE_FULL )

exp <- dba.contrast(exp,categories = DBA_CONDITION)

exp <-dba.analyze(exp)

exp.DB <- dba.report(exp, contrast =1, bCalled = TRUE, bCalledDetail = TRUE, bCounts = TRUE, bNormalized = TRUE)

Please find the link for a google sheet below which has three examples from the report.

Three examples from the report generated

[Rory Stark]

I've tried this on a couple of my own data sets but can not reproduce it (the Called values are always the same when I compare sum(rep0$Called1 != rep1$Called1) and sum(rep0$Called2 != rep1$Called2)).

I'd like to get to the bottom of this. If you could send me a link to your two exp objects immediately after calling dba.count() I can try to see what is going on.

[vk]

Okay, I'm not sure where I have gone wrong. Here you go, exp objects after dba.count

Answer 1

[Rory Stark]

I am able to reproduce this and have identified, and fixed, a bug in how the Called status is maintained after filtering (when minCounts is lowered to 0, some counts are lessened, and additional sites are filtered based on the requirement that at least one sample have an RPKM value of 1 or more).

I have checked in a fix that will appear as DiffBind_3.0.14, available via BiocManager::install() in the next couple of days.

[vk]

oh! Thank you for answering to my question and fixing the bug so quickly! Also, being a newbie to bioinformatics, it is very inspiring to see the time and efforts put into developing such packages for all of us to use. Great work! Thanks, again.

[Rory Stark]

This was actually a non-trivial bug, so I appreciate you giving me the information I needed to reproduce and resolve it!

[vk]

Hi Rory, I analysed new datasets with DiffBind_3.0.15, I am still confused with the results when minCount is set to 0 and 1.

When minCount is set to 0, the results show only those peaks that are called in all the samples (3 WT and 3 KO) without any change to the filter.

In the results obtained when minCount was set to 1, a sample that had almost the same number of reads (shown below) as the other WT's, was not called.

WT_1=792.5
WT_2= 377.89
WT_3=377.26
KO_1=177.93
KO_2= 40.82
KO_3= 69.38
Called1=2
Called2=0

Thank you for your time.