Another option I use occasionally is the ChIPseeker package in R

The way I would do it is by mapping alignments present in a sample for a given interval.

Split your current bed file by sample name, so that you have one file for sample A, B, and so on:

$ cat sampleA.bed
chr1 1000 1005
chr1 1010 1015
...

$ cat sampleB.bed
chr1 2000 2005
chr1 2010 2015
...

Remember no header in the bed files.

If they are the same intervals for each sample then just generate one bed file for all of the unique intervals, in your case do something like:

$ tail -n+2 data.txt | cut -f 2-4 | sort -u > intervals.bed

Then use samtools to pull out the alignments corresponding to the bed intervals:

samtools view -bLh {sampleA.bed or intervals.bed} sampleA.bam > sampleA.subsampled.sam

Then use seqMINER to plot the sam alignments against all unique intervals. So in seqMINER

Load reference coordinates (i.e., peaks) <- intervals.bed

Load aligned reads <- your subsampled sam files.

Then cluster the peaks however you want. Typically I'll use enrichment linear at 10 clusters.

For the statistically analyses you'll have to do either differential expression, or differential binding analysis. It's not as straight forward as using raw counts for a region. For example, if you are looking at differentially sized regions then you need to normalize raw counts by feature size and sequencing depth. Typically the best way to do this is by using transcripts per million (TPM), then log2 transforming the TPM values. Alternatively, you can use statistically approaches utilizing negative binomial GLM. A couple good packages for this include DeSEQ2, edgeR, and for ChIP sample I'd suggest MACS2 or Homer.