Dynamic computation of chunkSize introduces gaps between bins at chunk ends

[dmalzl]

I frequently use deeptools and I am very fond of the versatility that the suite offers, but just recently I encountered a more or less nasty bug that might not be so dramatic for most use cases but is deleterious when an application relies on accurate bin positions.

Using the multiBamSummary tool I counted reads in a 5kb fashion with 16 cores for multiple BAMfiles. These counts per bin were later used to call regions that exhibit a given distribution of counts across these bins. During the development of this caller I discovered something that I already previously noticed. The expected behaviour of tiling a genome into bins would be successive bins of a given binsize e.g. if a chromosome is 100563 bp long I would assume a 5kb binning of it to look like:

chr 0 5000
chr 5000 10000
.
.
chr 95000 100000
chr 100000 100563

However, the deeptools output always exhibited some irregularities such that somewhere along the tiling process some bins seemed to be a little bit larger than the specified binsize, for whatever reason. Since accuracy is crucial for me I dug a little bit into the data and discovered that not the binsize seems to vary from time to time but there are gaps between contiguous stretches of consecutive bins. I dug a little bit into the code and discovered that this most certainly is due to the way the genome is split into chunks for parallel computation.

The chunkSize used for splitting is dynamically computed from some BAM statistic and matches exactly the length of the contiguous bin stretches (in my case 16650712 bp)

Thus, I suspect the following:
The computation of the chunksize results in a number that is not a multiple of the binsize. This scenario seems to be handled by simply omitting bin that are not of a minimum length. Consequently, the last 712 (or whatever binsize dependent residual) are is not included into the output of a given chunk. However subsequent chunks are started at end position of the previous chunk where a little portion was omitted, resulting in a gap of 712 bp (in my case) that is not part of the output.

I suggest fixing this in the next release

[dpryan79]

Indeed, for our own uses we've never all that much cared exactly where the bins are, but since you're probably not the only one for which getting them exactly evenly spaced is important then we should look into addressing this. Thanks for the detailed bug report!

[dpryan79]

Which version of deepTools is giving you non-adjacent bins in this case? With 3.3.1 I get adjacent bins, though they definitely are different sizes around chunk length boundaries.

[dmalzl]

I currently use version 3.3.0 for all my analysis and also for the one where I found the inconsistencies.

[dpryan79]

And you're seeing this with multiBamSummary bins, yes?

[dmalzl]

Exactly

[dpryan79]

Can you provide the exact command you're using? I still don't get the gaps you mentioned with version 3.3.0. In the next release I'll try to add a --genomeChunkLength argument anyway so you can directly modify this for consistency.

[dmalzl]

The used command was:

multiBamSummary bins -b $somebams -l $somelabels -bs 5000 --outRawCounts $rawcountsfile -o $npzfile -p 8 --ignoreDuplicates

[dpryan79]

This should now be fixed in the develop branch. The shorter regions at the end of genomic chunks are no longer skipped so there aren't gaps between bins any more. multiBamSummary also accepts a --genomicChunkSize option now, so you'll be able to tweak for consistency.

[dmalzl]

Thanks