Hello,

Thanks for your effort in producing such a good tool for analysis.

Recently, I am using Diffbind to characterize differential open chromosome regions with ATAC-seq. I have more than 50 samples, equipping with a bam file ranging from 10GB-30GB. However, I found it really hard to handle such a huge amount of data for DiffBind. While I set 'bParallel=TRUE' for dba.count function, it was still too slow to process.

Even when I selected four samples for testing, the tool still broke with two different kinds of error.

library("DiffBind")
samples=read.csv("res.info.csv")
tamoxifen <- dba(sampleSheet=samples)
#sample1 bulk B1 treat1 lib 1 bed
#sample2 bulk B1 treat2 lib 1 bed
#sample3 bulk B1 treat3 lib 1 bed
#sample4 bulk B1 treat4 lib 1 bed
tamoxifen
#4 Samples, 135778 sites in matrix (188638 total):
#         ID Tissue Factor    Condition Treatment Replicate Intervals
#1   sample1   bulk     B1   treat1     lib         1    188126
#2 ...
tamoxifen <- dba.count(tamoxifen,bParallel=TRUE)
#Computing summits...
#Re-centering peaks...
#Reads will be counted as Paired-end.
#Reads will be counted as Paired-end.
#Reads will be counted as Paired-end.
#Reads will be counted as Paired-end.
Type1 error:
#Error: Error processing one or more read files. Check warnings().
#In addition: Warning messages:
#1: In mclapply(arglist, fn, ..., mc.preschedule = TRUE, mc.allow.recursive = TRUE) :
#  scheduled core 2 encountered error in user code, all values of the job will be affected
#2:   error in evaluating the argument 'x' in selecting a method for function 'assay': wrong args for environment subassignment
Type2 error:
#Re-centering peaks...
#Sample: ./input.bam125
#Reads will be counted as Paired-end.
#'BiocParallel' did not register default BiocParallelParam:
#  missing value where TRUE/FALSE needed
#Error in h(simpleError(msg, call)) :
#  error in evaluating the argument 'x' in selecting a method for function 'assay': error in evaluating the argument 'BPPARAM' in selecting a method for function 'bplapply': attempt to select less than one element in get1index

I am wondering if there is any suggestion for my case. Many thanks in advance!

Sincerely yours Xiaoyu

In a Reply, littlefishes20 asks:

In terms of summits setting, it works when I set summits=0! Could you please explain what difference the summit setting could make?

I'm breaking this out into an Answer as it really is a distinct question.

The summits parameter in dba.count() alters the the consensus peakset in multiple, potentially impactful, ways.

In the default case where a consensus peakset has been formed by merging peak calls from many sample replicates, the resulting intervals may get wider and wider. Overlapping peaks are merged by increasing their width to cover all of the overlapping regions, so the more samples in the experiment, the wider the merged peaks can get. This can be problematic in that wider peaks are more likely to include "background" regions. In addition, the distribution of peak widths across the consensus peakset may have some outliers (especially very wide intervals) that can have an impact on the underlying assumptions of the statistical methods doing the modelling.

In general, DiffBind emphasizes confidently identifying regions with different read count distributions, rather than identifying the exact boundaries of those regions. You should be able to rely on the intervals identified as being differentially bound with low FDR, not necessarily that those regions perfectly encompass the precise boundaries of differential enrichment. To do this, we attempt to minimize the inclusion of background regions and focus on intervals where there is most likely to be true enrichment in at least one sample group.

The idea behind summits is then to identify the maximal point of enrichment across the samples in which the peaks were identified, and specify a new peak interval centered on this point. The widths of these intervals are drawn consistently using the value of summits (including summits base-pairs up- and down- stream of the summit, as well as the summit itself). This address both the issue of excess background (so long as the value of summits is appropriate to the assay, specifically not too large) as well as regularizing the distribution of peak widths.

Also worth noting is that if any of the re-centered peaks overlap, they will be merged and again re-centered, resulting in only one peak where there may have been more in the consensus peakset. So there may be cases where the original consensus peakset contains narrow peaks, but the summits parameter is larger, and after re-centering there are overlapping peaks that get merged. There are also edge cases where a summit is calculated to be near the edge of a peak interval and gets merged with an adjacent, but not previously overlapping, peak region after being re-centered.

When summits=FALSE, this re-centering is not performed at all and the (merged) consensus peakset is used as is. Re-centering is also not performed when summits=0, however in this case the summits are computed such that they may be examined, and subsequent calls to dba.count() can be made to see how different values of summits work without having to recompute all of the summits each time.

All of this complexity comes from the fact that in most of these enrichment assays, we do not have a pre-defined set of genomic intervals on which to focus. Windowing approaches that skip peak calling completely, like csaw, avoid some of these issues