I am trying to use deepTools 3.5 ComputeMatrix + PlotHeatmap to visualise my ATAC-seq data. To prep my data, I merged BAM files of my 4 ATAC-seq replicates, then used code bellow to transform them into read-depth normalised bedgraph/bigwig:

bedtools genomecov -bga -ibam ${bam} -g ${flens} -scale ${scale_factor} > ${bg}

bedGraphToBigWig ${bg} ${flens} ${bw}

I then found my regions of interest by calling peaks from the same bedgraph file using SEACR. I am trying to plot the heatmap over these regions by:

computeMatrix scale-regions -S ATAC-input.bw -R SEACR-peaks.bed  -o matrix.gz -p 4 --missingDataAsZero -a 1000 -b 1000

plotHeatmap -m matrix.gz -o matrix.png --colorMap Purples

I looked at the white regions and there is definitely signal in there. ComputeMatrix seems to be skipping whole chromosomes 10, 11, 12, 13, 14, 15 and a big chunk of chromosome 16. In the screenshots bellow, bottom row is the regions I am trying to plot, middle row is the empty regions in the plot (all values 0 in the matrix), top row is ATAC signal.

Zoomed in on chr16 where the shift happens:

If I zoom into chromosome 16 I don't see any difference between regions that come up as empty or not. I also inspected my input files and I don't see any difference between areas in chr16 that get plotted vs don't. E.g. in my bed regions file, the top 3 regions (bold) come up as empty and the next 3 show signal:

**chr16 73896687    73897570    88032.1 365.705 chr16:73897113-73897117 883**

**chr16 73908008    73908643    37439.1 125.711 chr16:73908262-73908281 635**

**chr16 73997439    73998545    176944  399.99  chr16:73997864-73997865 1106**

chr16   74053025    74053852    94431.9 308.564 chr16:74053304-74053311 827

chr16   74066198    74066840    53233   205.709 chr16:74066718-74066733 642

chr16   74092928    74093735    56078.6 194.281 chr16:74093505-74093513 807

The same area cut from bedgraph (predecesor of bw used for plotting) - the top 3 regions (bold) come up as empty:

**chr16 73998538    73998542    68.5697**

**chr16 73998542    73998544    57.1414**

**chr16 73998544    73998545    11.4283**

chr16   73998545    73998617    0

chr16   73998617    73998627    11.4283

chr16   73998627    73998658    22.8566

I am not getting any error messages. For troubleshooting I have also attempted: using single thread, increasing allocated memory significantly, enabling/disabling sorting, using different deeptools versions (2.5, 3.0, 3.5). Since I am plotting 80k+ regions, I suspected that the issue can be caused by lack of memory, but I tried plotting those regions separately and I still get a completely empty plot:

Any help would be appreciated, I am running out of ideas.